An integrative analysis of DNA methylation in osteosarcoma.

BACKGROUND: The study aimed to analyze aberrantly methylated genes, relevant pathways and transcription factors (TFs) in osteosarcoma (OS) development.
METHODS: Based on the DNA methylation microarray data GSE36002 that were downloaded from GEO database, the differentially methylated genes in promoter regions were identified between OS and normal samples. Pathway and function enrichment analyses of differentially methylated genes was performed. Subsequently, protein-protein interaction (PPI) network was constructed, followed by identification of cancer-associated differentially methylated genes and significant differentially methylated TFs.
RESULTS: A total of 1379 hyper-methylation regions and 169 hypo-methylation regions in promoter regions were identified in OS samples compared to normal samples. The differentially hyper-methylated genes were significantly enriched in Neuroactive ligand-receptor interaction pathway, and Peroxisome proliferator activated receptor (PPAR) signaling pathway. The differentially hypo-methylated genes were significantly enriched in Toll-like receptor signaling pathway. In PPI network, signal transducers and activators of transcription (STAT3) had high degree (degree=21). MAX interactor 1, dimerization protein (MXI1), STAT3 and T-cell acute lymphocytic leukemia 1 (TAL1) were significant TFs enriched with target genes in OS samples. They were found to be cancer-associated and hyper-methylated in OS samples.
CONCLUSION: Neuroactive ligand-receptor interaction, PPAR signaling, Toll-like receptor signaling pathways are implicated in OS. MXI1, STAT3, and TAL1 may be important TFs involved in OS development.

Introduction

Osteosarcoma (OS) is the most common primary malignant bone tumor in children and adolescents. More than 15% of patients with OS develop metastases, frequently to lung [1]. For patients with metastasis or recurrence, the long-term survival rate is less than 20% [2], [3]. Understanding of the molecular mechanisms of OS would provide a basis for developing new therapeutic strategies.

It has been demonstrated that genetic alternations in the status of DNA methylation count among the most common molecular alterations in human neoplasia [4]. DNA methylation usually results in obstruction of the promoter region, hampering gene transcription and causing gene silencing [5]. Quite a few studies have reported findings related to DNA methylation in OS. It has been reported that methylation of frizzled-related proteins (SFRPs) may promote Wnt signaling pathway, thereby enhancing OS cell invasion [6]. Hyper-methylation of p14ADP-ribosylation factor (ARF) and estrogen receptor 1 (ESR1) have been found in OS as well, and may have implications in the prognosis of OS patients [7]. Moreover, Lu et al. [8] have reported that iroquois homeobox 1 (IRX1) enhances OS metastasis and may be a potential molecular marker. Recent study [9] suggests that promoter hyper-methylation of reversion-inducing cysteine-rich protein with Kazal motifs (RECK) is a causative factor in metastasis of OS. Despite some advances have been achieved in this field, certain mechanisms underlying OS remain largely unknown. Microarray analysis has been applied to identify the gene alterations and screen potential targets in human OS cell lines [10]. Kresse et al. [11] integrated genomewide genetic and epigenetic profiles from the EuroBoNeT panel, a European Network of Excellence on bone tumors (http://www.eurobonet.eu), of 19 human osteosarcoma cell lines based on microarray technologies, and deposited the DNA methylation dataset in the Gene Expression Omnibus (GEO) data repository (accession number GSE36002). In their study, they have comprehensively analyzed the relationships of DNA copy number, DNA methylation and mRNA expression in OS. Additionally, they screened out the differentially methylated genes and performed functional enrichment analysis. However, the interaction between differentially methylated genes, and the classification of these genes have not been analyzed.

Since methylation of CpG islands in promoter regions is a mechanism for inactivating genes, we estimated the differentially methylated genes in promoter regions between OS and normal samples from GSE36002 in this study. Pathway and functional enrichment analyses of differentially methylated genes was then performed. Subsequently, protein-protein interaction (PPI) network was constructed to analyze the interactions between differentially methylated genes. Furthermore, the classifications of differentially methylated genes, including tumor suppressor (TS) genes, oncogenes and TFs were investigated.

Materials and methods

DNA methylation microarray data

The GSE36002 DNA methylation microarray data [11] were downloaded from GEO database in National Center of Biotechnology Information (NCBI) (http://www.ncbi.nlm.nih.gov/gds/), based on the platform of Illumina Human Methylation 27 BeadChip (Illumina Inc., California, USA). GSE36002 dataset is comprised of 25 samples: 19 osteosarcoma cell lines samples and 6 normal control samples (four normal bone samples and two osteoblast cultures).

Identification of differentially methylated genes

DNA methylation microarray data were processed with simple scaling normalization using Lumi package (http://bioconductor.org/packages/release/bioc/html/lumi.html) [12], [13] in Bioconductor [14]. Subsequently, the differentially methylated regions (DMR) between osteosarcoma and normal samples were identified using methyAnalysis with M-value>1 and false discovery rate (FDR)<0.01.

M-values were calculated by the formula:

FDR method is also called Benjamini and Hochberg (BH) method [15], used to adjust p value. In detail, the original p values of all genes were ranked in descending order. The maximum p value was assigned as n, and the minimum was assigned as 1. The adjusted p value (FDR) was calculated as followed:where n represents the number of all genes; i represents the ith p value (from minimum to maximum).

The identified DMR were performed gene annotation to screen the differentially methylated sites located in promoter region. Dkhil et al. [16] have reported that most of the promoters display the changes of DNA methylation in their Ups-regions, which are between +500 and +2000 bp upstream from the transcription start site of genes. Therefore, in consideration of numerous genes in the dataset, the promoter region was defined as 2000 bp upstream of the transcription start site in this study.

Functional enrichment analysis of differentially methylated genes

Functional enrichment analyses included Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, Reactome pathways and Gene Ontology (GO) terms analyses. The differentially methylated genes were performed functional enrichment analysis using GOFunction (http://www.bioconductor.org/packages/release/bioc/html/GOFunction.html) in R package. P value<0.05 and count (number of significantly enriched genes)≥2 were used as thresholds.

PPI network construction

PPI analysis can provide new insights into protein function, and uncover the generic organization principles of functional cellular networks [17]. Therefore, we constructed PPI network to further analyze the functions of differentially methylated gene. The Search Tool for the Retrieval of Interacting Genes (STRING) database (http://www.mybiosoftware.com/pathway-analysis/7789) provides the information of both experimental and predicted interactions [18]. With this tool, the differentially methylated gene pairs with significant interactions (combine score>0.9) were identified and used for construction of PPI network which was visualized using Cytoscape (http://cytoscape.org/) [19]. In the network, the nodes with higher degrees were defined as hub nodes.

Gene classification analysis

Based on tumor suppressor (TS) genes database [20] and tumor associated genes (TAG) database [21], the known TS genes and oncogenes were selected from the obtained differentially methylated genes.

TF analysis

On the basis of the TF-target gene pairs in Encode [22], [23], the differentially methylated genes corresponding to the differentially methylated TFs were identified. Then according to the obtained TF-target gene interaction pairs, significant TFs (p value<0.01) were screened using hypergeometric analysis. The formula of hypergeometric distribution was shown as follows:where N represents the total number of genes; n represents the number of target genes regulated by TF; M represents the number of differentially methylated genes; k represents the number of differentially methylated genes regulated by TF.

Results

In the study, a total of 2845 differentially methylated loci were identified in OS samples relative to normal samples, including 1379 hyper-methylation regions and 169 hypo-methylation regions in the promoter region. Additionally, the numbers of hyper- and hypo-methylated loci in enhancers are shown in Table 1, and in exons and other regions are shown in Table 2.

Table 1

The numbers of hyper- and hypo-methylated loci in enhancer region.

Enhancer	Hyper	Hypo	Total
TRUE	43	163	206
FALSE	1490	1071	2561

Table 2

The numbers of hyper- and hypo-methylated loci in exon, promoter and other regions.

Regions	Hyper	Hypo	Total
1stExon	57	489	546
3'UTR	3	8	11
Body	84	619	703
IGR	12	25	37
Promoter	1379	169	1548

UTR: untranslated regions; IGR: intergenic region.

Functional enrichment analyses of differentially methylated genes

Pathway enrichment analysis revealed that the differentially hyper-methylated genes were significantly enriched in KEGG pathways including Pathways in cancer, Bladder cancer, Neuroactive ligand-receptor interaction, and Peroxisome proliferator-activated receptor (PPAR) signaling pathways (Table 3). The differentially hypo-methylated genes were significantly related to Cytokine-cytokine receptor interaction and Toll-like receptor signaling pathways (Table 4). With regard to Reactome pathways, the differentially hyper-methylated genes were primarily enriched in Gastrin-CREB signaling pathway via PKC and MAPK, and G alpha (q) signaling events pathway, while the differentially hypo-methylated genes were mainly enriched in Defensins, and Alpha-defensins pathways (Table 5). For GO function, the differentially hyper-methylated genes were predominately enriched in developmental process and anatomical structure development (Table 6), while the differentially hypo-methylated genes were primarily enriched in extracellular region (Table 7).

Table 3

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (top 10) of the differentially hyper-methylated genes.

KEGG ID	Name	Gene count	P value
5200	Pathways in cancer	40	0.002712
4080	Neuroactive ligand-receptor interaction	34	0.004135
3320	PPAR signaling pathway	12	0.007642
4974	Protein digestion and absorption	13	0.009805
4512	ECM-receptor interaction	13	0.014496
5219	Bladder cancer	8	0.014911
4920	Adipocytokine signaling pathway	11	0.015917
5414	Dilated cardiomyopathy	13	0.022575
4916	Melanogenesis	14	0.025452
4340	Hedgehog signaling pathway	9	0.028947

Gene count represents the number of gene enriched in the pathway.

Table 4

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (top 10) of the differentially hypo-methylated genes.

KEGG ID	Name	Gene count	P value
4740	Olfactory transduction	10	0.004797
5320	Autoimmune thyroid disease	3	0.015042
4060	Cytokine-cytokine receptor interaction	7	0.016163
4620	Toll-like receptor signaling pathway	4	0.018845
4950	Maturity onset diabetes of the young	2	0.025692
830	Retinol metabolism	3	0.026061
4622	RIG-I-like receptor signaling pathway	3	0.034039
982	Drug metabolism - cytochrome P450	3	0.036526
4140	Regulation of autophagy	2	0.045367

Gene count represents the number of gene enriched in the pathway.

Table 5

The reactome pathways (top 5) of the differentially hyper- and hypo-methylated genes.

Reactome ID	Name	Gene count	P value
Differentially hyper-methylated genes
881907	Gastrin-CREB signaling pathway via PKC and MAPK	34	6.26E−06
416476	G alpha (q) signaling events	29	8.10E−05
163685	Integration of energy metabolism	19	0.000224877
500792	GPCR ligand binding	51	0.000330922
1474244	Extracellular matrix organization	35	0.000489987
Differentially hypo-methylated genes
1461973	Defensins	7	4.84E−07
1462054	Alpha-defensins	3	0.000113
1592389	Activation of Matrix Metalloproteinases	4	0.000245
381753	Olfactory Signaling Pathway	11	0.000577
1461957	Beta defensins	4	0.000602

Gene count represents the number of gene enriched in the pathway.

Table 6

The gene ontology (GO) of the differentially hyper-methylated genes.

GO categories	Name	Gene count	P value
BP	Developmental process	483	1.11E−16
	Anatomical structure development	436	2.22E−16
	System development	387	4.44E−16
	Cell differentiation	323	2.00E−15
CC	Extracellular region part	141	1.79E−10
	Plasma membrane part	211	1.14E−09
	Proteinaceous extracellular matrix	60	1.29E−09
	Extracellular matrix	64	1.22E−08
MF	Sequence-specific DNA binding	92	3.03E−08
	Sequence-specific DNA binding RNA Polymerase II transcription factor activity	46	8.37E−07
	Voltage-gated cation channel activity	28	9.45E−07
	Sequence-specific DNA binding transcription factor activity	118	1.16E−06

BP: biological process; CC: cellular component; MF: molecular function; Gene count represents the number of gene enriched in GO term.

Table 7

The gene ontology (GO) of the differentially hypomethylated genes.

GO categories	Name	Gene count	P value
BP	Keratinocyte differentiation	9	3.54E−07
	Epidermal cell differentiation	10	7.98E−07
	Keratinization	6	1.88E−06
	Defense response	29	3.62E−06
CC	Extracellular region	50	3.73E−11
	Extracellular space	23	1.18E−06
	Extracellular region part	27	1.81E−06
	Extracellular region	50	3.73E−11
MF	Serine-type endopeptidase activity	8	2.01E−05
	Serine-type peptidase activity	8	5.12E−05
	Serine hydrolase activity	8	5.56E−05
	Monooxygenase activity	6	8.15E−05

BP: biological process; CC: cellular component; MF: molecular function; Gene count represents the number of gene enriched in GO term.

In this study, a PPI network with 442 nodes and 693 PPI pairs was constructed (Fig. 1). In this network, adenylate cyclase 2 (ADCY2) (degree=34), proopiomelanocortin (POMC) (degree=29), signal transducers and activators of transcription (STAT3) (degree=21), neuropeptide Y (NPY) (degree=20), and somatostatin (SST) (degree=18) with degrees more than 17 were regarded as hub nodes. The sub-network modules that contained the five hub nodes are shown in Fig. 2A–E.

Fig. 1

A protein-protein interaction network of differentially methylated genes. Red nodes stand for differentially hyper-methylated genes, while green nodes stand for differentially hypo-methylated genes.

Fig. 2

The networks that contained five hub genes of adenylate cyclase 2 (ADCY2), proopiomelanocortin (POMC), signal transducers and activators of transcription (STAT3), neuropeptide Y (NPY), and somatostatin (SST). Red nodes stand for differentially hyper-methylated genes, while green nodes stand for differentially hypo-methylated genes.

Among the identified differentially hyper-methylated genes, 27 were found to be oncogenes, 78 were TS genes and 22 were tumor-associated genes (oncogenes or TS genes). The identified differentially hypo-methylated genes contained one oncogene, five TS genes and one tumor-associated gene.

A total of 76 TFs that were significantly enriched in target genes were identified (Supplementary material Fig. 1). Among these TFs, MAX interactor 1, dimerization protein (MXI1), STAT3, and T-cell acute lymphocytic leukemia 1 (TAL1) were also cancer-associated TFs. Interestingly, STAT3 was a hub node in the PPI network. Furthermore, the methylation levels of the three TFs were up-regulated in OS samples relative to control samples (Fig. 3).

Fig. 3

The methylation levels of MXI1, STAT3 and TAL1 in all samples. Horizontal axis represents 25 samples, and vertical axis represents the M value of methylation level of MXI1, STAT3 and TAL1.

Discussion

OS is a rare malignant tumor with a high tendency of metastasis. Aberrant methylation of tumor-related genes in the promoter region is associated with human tumors [24]. In the current study, a total of 2845 differentially methylated loci were identified in OS samples compared to normal samples, including 1379 hyper-methylation regions and 169 hypo-methylation regions in the promoter region. Although so many differentially methylated loci were identified, not all of them may involve in OS progression. Further analyses revealed that the differentially hypo-methylated genes were significantly enriched in Toll-like receptor signaling pathway, while differentially hyper-methylated genes were significantly enriched in Neuroactive ligand-receptor interaction pathway, and PPAR signaling pathway. Additionally, ADCY2 and POMC had the top two highest degrees in the PPI network. Moreover, MXI1, STAT3, and TAL1 were significant TFs enriched with target genes, and they were also tumor-associated TFs.

Toll-like receptors play important roles in the innate immune system, especially in inflammatory response, which is considered to be an important epigenetic factor contributing to neoplasia and tumor progression [25], [26]. Triggering of toll-like receptor 4 (TLR4) on metastatic breast cancer cells has been found to promote adhesion and invasive migration of tumor cells [27]. There is evidence that inflammatory cytokines, such as TNF-α and IL-1 are required for promoting the tumorigenesis of osteosarcoma [28]. The present study found that the differentially hypo-methylated genes in OS samples were related to Toll-like receptor signaling pathway. These findings led to a speculation that the differentially hypo-methylated genes and the Toll-like receptor signaling pathway might be involved in OS-related inflammation. A recent microarray analysis reported that Neuroactive ligand-receptor interaction pathway and PPAR signaling pathway were associated with OS metastasis [29]. Moreover, it has been established that PPAR agonists can suppress OS proliferation and growth, and induce apoptosis [30]. Results of this study suggest that the roles of Neuroactive ligand-receptor interaction pathway and PPAR signaling pathway in OS may be partly due to hyper-methylation of genes.

It has been reported that PPI network can organize all protein-coding genes into a large network which provides a better understanding of the functional organization of the proteome [17]. In the present study, PPI analysis found that ADCY2 and POMC had the top two highest degrees in the constructed PPI network. Genome-wide studies have shown that deletion of a hub protein is more likely to be lethal than deletion of a non-hub protein. As a result, hub nodes may play more important roles in OS than the other nodes. Presently, we focused on the top two nodes with higher degrees for discussion. Interestingly, they have been reported to be implicated in the progression of OS. Recently, Sun et al. [31] performed gene expression profiling analysis of OS cell lines, and found that ADCY2 was involved in purine metabolism pathway in OS. POMC has been found in osteoclasts [32], which is a precursor of β-endorphin. Baamonde et al. [33] reported that endogenous β-endorphin involved in the initial stages of murine OS. Taken together, our study may further suggested the role of ADCY2 and POMC in OS.

STAT3 is a member of STAT protein family and plays a critical role in cell growth and apoptosis [34]. Increasing studies have established that STAT3 activation promotes the development of OS [35], [36]. STAT3 was found to be hyper-methylated in OS samples in the present study, suggesting that STAT3 might play a role in OS. It has been reported that STAT3 had an oncogenic or a tumor suppressor role in human brain tumor depending on the mutational profile of the tumor [37]. It indicates that STAT3 may also play dual roles in promoting or discouraging the development of OS. Moreover, STAT3 was a hub node (degree=21) in the PPI network of this study, suggesting that the role of STAT3 in OS might be associated with its interactions with other proteins.

MXI1 is a member of the MAD gene family, which plays a key role in regulating cell proliferation and growth [38], [39]. It is located at 10q24-25, where loss of heterozygosity occurs in several human cancers, including endometrial cancers, prostate tumors, gliomas, renal cell carcinomas, and small-cell lung cancers [40]. Importantly, numerous experiments have verified that MXI1 functions as a tumor suppressor gene via negatively modulating the promoter of oncogene c-Myc, which is involved in the control of cell proliferation and apoptosis [41], [42], [43], [44]. To our knowledge, the role of MXI1 in OS remains unclear. The study revealed hyper-methylation of MXI1 in OS. In light of these results, we speculate that MXI1 may inhibit progression of OS by compromising cell proliferation and inducing apoptosis. TAL1 is a basic helix-loop-helix TF required for blood cell development [45]. In the study of Kresse et al. [46], TAL1 was found to be significantly differentially methylated between OS and normal control samples. In agreement with the study above, the present study also revealed that it was a differentially methylated gene and a tumor-associated TF. Therefore, TAL1 might be an important TF associated with OS.

In this study, some limitations have to be acknowledged. First, the sample size was a little small. Second, no experiment was conducted to validate the expression level of these methylated genes, including MXI1, STAT3 and TAL1. Therefore, further studies with more samples and experimental validations are needed to validate our results.

Taken together, microarray analysis of epigenetic alterations revealed that Neuroactive ligand-receptor interaction pathway, PPAR signaling pathway and toll-like receptor signaling pathway may be implicated in tumorigenesis of OS. MXI1, STAT3, and TAL1 might be critical TFs in development of OS. Results of the study provide more enlightening insights concerning the association of DNA methylation with tumorigenesis of OS. More studies are warranted to verify these results.

Competing interests

The authors declare that they have no competing interests.

Authors contribution

Jie Xu, Deng Li and Ruofan Ma Wrote Manuscript, Designed Research,

Zhiqing Cai, Yingbin, Zhang Yulin Huang and Baohua Su help Analyzed Data.