Abstract
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 [35]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.
Keywords: Differentially methylated gene, Osteosarcoma, Protein-protein
interaction, Transcription factor
1. 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 [36][1]. For patients with metastasis or
recurrence, the long-term survival rate is less than 20% [37][2],
[38][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 [39][4]. DNA methylation usually results in obstruction of
the promoter region, hampering gene transcription and causing gene
silencing [40][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 [41][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 [42][7]. Moreover, Lu et al. [43][8] have reported that
iroquois homeobox 1 (IRX1) enhances OS metastasis and may be a
potential molecular marker. Recent study [44][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 [45][10]. Kresse et al. [46][11] integrated
genomewide genetic and epigenetic profiles from the EuroBoNeT panel, a
European Network of Excellence on bone tumors
([47]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 [48]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 [49]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.
2. Materials and methods
2.1. DNA methylation microarray data
The [50]GSE36002 DNA methylation microarray data [51][11] were
downloaded from GEO database in National Center of Biotechnology
Information (NCBI) ([52]http://www.ncbi.nlm.nih.gov/gds/), based on the
platform of Illumina Human Methylation 27 BeadChip (Illumina Inc.,
California, USA). [53]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).
2.2. Identification of differentially methylated genes
DNA methylation microarray data were processed with simple scaling
normalization using Lumi package
([54]http://bioconductor.org/packages/release/bioc/html/lumi.html)
[55][12], [56][13] in Bioconductor [57][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:
[MATH: M−value=log2methylated probe
intensityunmethylated probe
intensity :MATH]
(1)
FDR method is also called Benjamini and Hochberg (BH) method [58][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:
[MATH: FDR=originalpvalue*(n/i) :MATH]
(2)
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. [59][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.
2.3. 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
([60]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.
2.4. PPI network construction
PPI analysis can provide new insights into protein function, and
uncover the generic organization principles of functional cellular
networks [61][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
([62]http://www.mybiosoftware.com/pathway-analysis/7789) provides the
information of both experimental and predicted interactions [63][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
([64]http://cytoscape.org/) [65][19]. In the network, the nodes with
higher degrees were defined as hub nodes.
2.5. Gene classification analysis
Based on tumor suppressor (TS) genes database [66][20] and tumor
associated genes (TAG) database [67][21], the known TS genes and
oncogenes were selected from the obtained differentially methylated
genes.
2.6. TF analysis
On the basis of the TF-target gene pairs in Encode [68][22], [69][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:
[MATH: P(X=k)=CMkCN−Mn−kCNnk<
/mi>∈{0,1,2,…m} :MATH]
(3)
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.
3. Results
3.1. Identification of differentially methylated genes
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 [70]Table 1, and in
exons and other regions are shown in [71]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
[72]Open in a new tab
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
[73]Open in a new tab
UTR: untranslated regions; IGR: intergenic region.
3.2. 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 ([74]Table 3). The differentially
hypo-methylated genes were significantly related to Cytokine-cytokine
receptor interaction and Toll-like receptor signaling pathways
([75]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 ([76]Table 5). For
GO function, the differentially hyper-methylated genes were
predominately enriched in developmental process and anatomical
structure development ([77]Table 6), while the differentially
hypo-methylated genes were primarily enriched in extracellular region
([78]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
[79]Open in a new tab
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
[80]Open in a new tab
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
[81]Open in a new tab
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
[82]Open in a new tab
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
[83]Open in a new tab
BP: biological process; CC: cellular component; MF: molecular function;
Gene count represents the number of gene enriched in GO term.
3.3. PPI network construction
In this study, a PPI network with 442 nodes and 693 PPI pairs was
constructed ([84]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 [85]Fig. 2A–E.
Fig. 1.
[86]Fig. 1
[87]Open in a new tab
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.
[88]Fig. 2
[89]Open in a new tab
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.
3.4. Gene classification analysis
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.
3.5. TF analysis
A total of 76 TFs that were significantly enriched in target genes were
identified ([90]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 ([91]Fig. 3).
Fig. 3.
[92]Fig. 3
[93]Open in a new tab
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.
4. 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 [94][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 [95][25], [96][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 [97][27]. There is
evidence that inflammatory cytokines, such as TNF-α and IL-1 are
required for promoting the tumorigenesis of osteosarcoma [98][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 [99][29]. Moreover, it has
been established that PPAR agonists can suppress OS proliferation and
growth, and induce apoptosis [100][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 [101][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. [102][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 [103][32],
which is a precursor of β-endorphin. Baamonde et al. [104][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 [105][34]. Increasing studies have
established that STAT3 activation promotes the development of OS
[106][35], [107][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 [108][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 [109][38], [110][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 [111][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 [112][41], [113][42], [114][43], [115][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 [116][45]. In the study of
Kresse et al. [117][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.
Acknowledgements