Aberrantly methylated-differentially expressed genes and related pathways in cholangiocarcinoma.

ABSTRACT: This study aimed to explore aberrantly methylated-differentially expressed genes and related pathways in cholangiocarcinoma (CCA).The mRNA expression data (GSE26566) and methylation profiling data (GSE44965) were collected from the Gene Expression Omnibus (GEO) Datasets. Differentially expressed genes and differentially methylated genes were identified using GEO2R. Gene ontology analysis and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses were performed using clusterprofiler in R. MCODE clustering tool was used to screen modules of the protein-protein interaction network in Cytoscape. Related pathways of hub gene by using gene set enrichment analysis.Eighty-one hypermethylated, lowly expressed genes (Hyper-LGs) and 76 hypomethylated, highly expressed genes (Hypo-HGs) were identified in this study. Hyper-LGs were enriched in ion channel binding and transcription factor activity, which was associated with Mineral absorption and Cell adhesion molecules. Hypo-HGs were enriched in cysteine-type endopeptidase activity, which was associated with Sphingolipid signaling pathway and T cell receptor signaling pathway. Based on protein-protein interaction networks, MYC and VWF were identified as hub genes for Hyper-LGs, and no hub genes for Hypo-HGs.This study found methylated-differentially expressed genes and signaling pathways that are connected with the CCA by using a series of bioinformatics databases and tools. MYC and VWF act as hub genes of CCA, which can be used as biomarkers based on aberrant methylation for the accurate diagnosis and treatment of CCA.

Introduction

Cholangiocarcinoma (CCA) is the second most common liver malignancy after hepatocellular carcinoma. It is a malignant tumor of the bile duct epithelium that occurs from the left and right hepatic ducts to the lower end of the common bile duct.[ It can be divided into intrahepatic cholangiocarcinoma, perihilar cholangiocarcinoma and common bile duct cancer according to the location.[ Perihilar cholangiocarcinoma is more common, accounting for 60% to 80%. Cholangiocarcinoma predominantly occurs in people over the age of 50, and there is no difference between men and women.[ Because there are no obvious symptoms in the early stage of cholangiocarcinoma, and neither serological examination nor imaging examination can accurately diagnose gallbladder cancer. Most patients are found at an advanced stage and have a poor prognosis. Therefore, we need a kind of biomarker that can quickly and accurately diagnose cholangiocarcinoma.[

Epigenetics studies the variety of heritability on gene expression without changing the nucleotide sequence of a gene. DNA methylation is one of the widest studies among epigenetic modifications, which mainly includes DNA hypomethylation and DNA hypermethylation.[ They induce cancer by activating oncogene expression and inactivating suppressor genes. Abnormal DNA methylation at CpG island is found in most cancers, and it can lead to cancer cells proliferation by downregulating the expression of suppressor.[ Thus it can be seen that aberrant DNA methylation plays an important role in the occurrence and development of tumor. In order to discover methylated-differentially expressed genes (MDEGs) and pathways in CCA, we used the mRNA expression data (GSE26566) and methylation profiling data (GSE44965) in Gene Expression Omnibus (GEO) Datasets. By comparing the normal tissue and cholangiocarcinoma tissue, we identified hypermethylated, lowly expressed genes (Hyper-LGs) and hypomethylated, highly expressed genes (Hypo-HGs). And then we explored the functions of the MDEGs by GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. Finally, we found the hub genes and pathways by protein–protein interaction (PPI) network and Gene Set Enrichment Analysis (GSEA). We aimed to provide new ideas for the diagnosis and treatment of cholangiocarcinoma.

Materials and methods

Microarray data

Our study analyzed one mRNA microarray data sets and one methylation profiling data sets to explore MDEGs between cholangiocarcinoma samples and normal tissue samples. We collected the mRNA expression data (GSE26566) and methylation profiling data (GSE44965) from the Gene Expression Omnibus (GEO) Datasets (https://www.ncbi.nlm.nih.gov/gds/). 4 normal tissue samples and 104 cholangiocarcinoma tissue samples were contained in the GSE26566 series, while 9 cholangiocarcinoma samples and 9 adjacent normal samples were found in the GSE44965 series. The study was approved by the Ethics Committee of the First Hospital of China Medical University.

Data processing

We selected differentially expressed genes and methylation genes between cholangiocarcinoma and normal samples separately in each of the data sets by using GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/). GEO2R is an interactive tool for comparing different samples groups to identify differential expression in a GEO series. FDR <0.05 and |logFC| >1 were considered statistically significant in expression datasets. While in methylation datasets, FDR <0.05 was considered statistically significant. Duplicate gene probes and unspecific probes will be removed.

Gene ontology and pathway enrichment analyses

Gene ontology analysis (GO) and KEGG are important bioinformatics tool. KEGG is equivalent to collection of databases, which contains diseases, information of genomes, biological pathways, and chemical substances. Go analysis can be used to annotate genes and gene products. GO and KEGG pathway enrichment analyses were performed using clusterprofiler in R. P <.05 was considered statistically significant.

Integration of PPI network and module analysis

We used The Retrieval of Interacting Genes (STRING) database tool (string-db.org) to analyse interactive relationships among the MDEGs. It was considered statistically significant when interactions with a combined score >0.4. Hub genes were considered when genes with connection numbers > 5. According to its results, PPI networks were constrcted by using the Cytoscape software. Then we used the plug-in Molecular Complex Detection (MCODE) clustering tool to filter modules of the PPI network in Cytoscape. MCODE scores ≥4 and number of nodes ≥4 were defined by a module. STRING was performed GO and KEGG pathway enrichment analyses of the genes in the selected MCODE modules. P < .05 was indicated significant differences.

Gene set enrichment analysis for hub gene

Co-expression analysis of high-dimensional expression data has proven effective for the study of gene functions. Here, using COEXPEDIA (http://www.coexpedia.org/), we searched for co-expressed genes of hub gene. Gene Set Enrichment Analysis (GSEA) is a computational method that assesses whether an a priori defined set of genes shows statistically significant, concordant differences between two biological states. Using GSEA, we searched for related pathways of hub gene.

Results

Identification of MDEGs in cholangiocarcinoma

1663 genes were extracted from the GSE26566 and 2573 gened were collected from GSE44965 data sets. A total of 807 genes and 856 genes were down-regulated and up-regulated in the expression profiling data sets, respectively. And 1146 genes were hypermethylated and 1427 genes were hypomethylated in the methylation profiling data sets. 81 Hyper-LGs and 76 Hypo-HGs were found by comparing the differentially expressed genes with the differentially methylated genes. Flowchart of bioinformatics analysis was seen in Figure 1. And the Volcano plot of GSE26566 was shown in Figure 2

Figure 1

Flowchart of bioinformatics analysis. DMGs = differentially methylated genes, DEGs = differentially expressed genes, Hypo-HGs = hypomethylated, highly expressed genes, Hyper-LGs = hypermethylated, lowly expressed genes.

Figure 2

The Volcano plot of GSE26566.

GO and pathway functional enrichment analyses

We used clusterprofiler to perform GO and KEGG pathway enrichment analyses for exploring the functions of the MDEGs. Hypo-HGs were enriched mainly in cysteine-type endopeptidase activity involved in apoptotic signaling pathway, and cell adhesion molecule binding (Supplementary Digital Content Figure 1). Hyper-LGs were enriched mainly in ion channel binding and transcription factor activity (Fig. 3). In the KEGG analysis, Hyper-LGs was associated with Mineral absorption and Cell adhesion molecules. While Hypo-HGs was mainly associated with Sphingolipid signaling pathway and T cell receptor signaling pathway (Fig. 4).

Figure 3

GO pathway enrichment analyses. GO = gene ontology analysis.

Figure 4

KEGG pathway enrichment analyses. KEGG = Kyoto Encyclopedia of Genes and Genomes.

PPI network construction and module selection

81 nodes and 43 edges were found in the Hyper-LGs network, while 75 nodes and 27 edges were in the Hypo-HGs network by uploading the MDEGs into the PPI network tool (Fig. 5 and Supplementary Digital content Figure 2). In figures, common genes were indicated by orange circles and hub genes were indicated by blue dots. The size increases with degree. The hub genes among the Hyper-LGs were annotated as MYC and VWF. No genes with connection numbers >5 were found among the Hypo-HGs.

Figure 5

The PPI network for the MDEGs. MDEGs = methylated-differentially expressed genes.

Module analysis

There was one module in the Hyper-LG network and one module in the Hypo-HG network. By the enrichment analysis, the genes in modules of Hyper-LG were associated with positive regulation of peptidase activity and cellular response to stress (Fig. 6A). However, no term was found in the Hypo-HG module (Fig. 6B)

Figure 6

Module analysis.

Gene set enrichment analysis for bub gene

Using COEXPEDIA, we found 572 genes associated with VMF and 408 genes related with MYC. We then used clusterprofiler to do GSEA. After analysis. MYC was associated with TNF signaling pathway and Jak-STAT signaling pathway. While, VMF was associated with Protein digestion and absorption and cGMP-PKG signaling pathway. (Fig. 7)

Figure 7

Gene set enrichment analysis for bub gene.

Discussion

In recent years, epigenetic studies on cancer have gradually increased, and gene methylation changes may play an important role in the occurrence and development of cancer. DNA methylation is an important epigenetic change. DNA methylation is a chemical modification mediated by DNA methyltransferase, which binds methyl groups to the 5’carbon position of CpG dinucleotide cytosine. CpG dinucleotides are mostly present in the CpG island region of the 5’ promoter of the gene. The methylation changes of the promoter CpG island are important to regulate gene expression.[ As a result, more and more studies explore the connection between the aberrant DNA methylation and the development of cancer. Lots of MDEGs are found and can become biomarkers in the diagnosis of cancer, especially in digestive system. Liu[ found 10 hub genes in 411 Hypo-HGs and 239 Hyper-LGs of colorectal cancer. 445 Hyper-LGs and 129 Hypo-HGs are found as MDEGs in gastric cancer, and CASR, CXCL12, SST were identified hub genes.[ In terms of hepatocellular carcinoma,[ Sang discovered 5 hub genes in 266 Hyper-LGs and 5 hub genes in 161 Hypo-HGs. And Cai[ use different data set to find the new hub gene CDCA5, which can change the cell cycle and make cells proliferate. In addition, many MDEGs are discussed in other systems, such as breast cancer,[ osteosarcoma,[ nasopharyngeal carcinoma,[ and ovarian cancer.[ Because the early symptoms are not obvious, most CCA patients have been diagnosed at the terminal stage. In terms of treatment, CCA is not sensitive to radiotherapy and chemotherapy, and surgery is the only way at present. In order to find new diagnosis and treatment method, there are many studies talking about the relationship between aberrant DNA methylation and CCA. Nakaoka[ summaried MLH, DCLK1, CDO1, ZSCAN18, ZNF331, p14 (ARF), p16 (INK4a, CDKN2A), DAPK, CCND2, CDH13, GRIN2B, RUNX3, TWIST1, EGFR, LKB1 that are frequently methylated genes in cholangiocarcinoma. Afterwards, Zhang[ identified 98 Hyper-LGs and 93 Hypo-HG. And then they found 9 hub genes (F2, AHSG, RRM2, AURKB, CCNA2, TOP2A, BIRC5, PLK1, ASPM). In order to find more biomarkers, we use different mRNA expression data (GSE26566) and methylation profiling data (GSE44965) In our study. Finally, we identified 81 Hyper-LGs and 76 Hypo-HGs that may be associated with the molecular regulation of CCA. Two hub genes, MYC and VWF, were both from Hyper-LGs. MYC is the common oncogenes, which can encode nucleoprotein (C-myc, N-myc, L-myc).[ C-myc protein can act on transcriptional regulatory region, and promote the production of cell cycle regulatory protein. It can not only promote excessive cell proliferation, but also inhibit the cell apoptosis. Study has shown that the abnormal expression of C-myc can lead to the activation of the serine synthesis pathway, which can promote the metabolic transformation of cancer cells. It has an important role in maintaining the survival and proliferation of cancer cells and promoting cancer development.[ VWF gene encode plasma protein vWF, whose function is to help platelet adhesion and hemostasis. Some studies pointed out that VWF gene can lead to the development of tumor by inhibiting angiogenesis and apoptosis. It may also promote the spread of cancer cells.[ Abnormal expression of VWF has been found in thyroid adenoma,[ colon carcinoma,[ non-small cell lung cancer,[ and clear cell renal cell carcinoma.[

According to the result of GSEA, MYC was associated with TNF signaling pathway and Jak-STAT signaling pathway in CCA. VMF was associated with protein digestion and absorption and cGMP-PKG signaling pathway. TNF signal[ combine with TNF receptor and make it multimerization, which bring an interaction with TRADD molecular. This process can induce the activation of NF-kB by TNFR-associated factor and receptor interacting protein. NF-kB enter into nucleus, combine with specific DNA sequence, start and enhance the transcription. Continuous and abnormal activation of NF-kB will cause cell proliferation out of control and inhibit the apoptosis. Aberrant JAK/STAT pathway has been found in many kinds of cancers, such as prostate cancer, hematopoietic malignancies and lung cancer.[ JAK (Janus Kinase) is activated after receiving signal and catalyse receptor tyrosine phosphorylation. STAT (signal transducer and activator of transcription) combine with receptor and enter into nucleus, which can promote gene expression. JAK/STAT signaling pathway plays an important role in cholangiocarcinoma. It is connected with Pr1, IL-6, EGF, and HGF that regulate the cholangiocellular function and enhance the survival and proliferation of malignant bile duct cells.[ STAT3 can also increase cell adhesion molecule to cause the tumor metastasis.[ In addition, the support of nutrients are necessary for the progression of tumors. So protein digestion and absorption is closely connected with the development of CCA.

According to our Gene ontology analysis, Hyper-LGs were mainly enriched in ion channel binding and transcription factor activity. Hypo-HGs were enriched in cell adhesion molecule binding and cysteine-type endopeptidase activity involved in apoptotic signaling pathway. These results are reasonable. As we all know, ion channel (K+, Cl-, Ca2+) is a switch for ions to enter and exit the cell. And abnormal ion channels can cause ion imbalance inside and outside the cell, which can act as signal factors to regular cell proliferation and cause different kinds of cancer.[ In addition, cell adhesion molecules are closed related to spread and metastasis of tumor. They can promote the adhesion of tumor cells to the basement membrane and extracellular matrix, and make tumor cells to infiltrate the basement membrane to the interstitial and form metastasis. However, the relationship between cysteine-type endopeptidase activity and the development of tumor should be explore further. According to the KEGG analysis, Hyper-LGs was associated with mineral absorption and cell adhesion molecules. Hypo-HGs was associated Sphingolipid signaling pathway and T cell receptor signaling pathway. Some studies[ found mineral supplements could prevent cancer, but definite mechanism was not clear. So mineral malabsorption may have an effect on the cancer. Sphingolipid signaling pathway is a complex process, which refers to a lot of enzymes. Sphingosine kinase 1 (SPHK1) generate sphingosine-1-phosphate which make the inhibitor of NF-κB kinase phosphorylate and NF-κB inhibitor degrade by binding to the TNF receptor-associated factor 2. SPHK1 can activate NF-kB signaling pathway to participate in cancer cell proliferation, invasion, migration and metastasis.[ It has been found SPHK1 has an important effect on various cancers.[

Module analysis showed KLF2, CAV1, MYC, and SNCA in the modules of Hyper-LG were associated with positive regulation of peptidase activity and cellular response to stress. What puzzles us is KLF2 because it belongs to tumor suppressor gene. Although there are few studies to explore the relationship between KLF2 and CCA, KLF2 plays an important role in orther cancer. KLF2 methylation reduce the expression of KLF2 in NSCLC cells, which is significantly related to lymph node metastasis and advanced TNM staging. This result may be related to the small sample size. So the relationship between KLF2 and CCA should be explored further.[ Overexpression of CAV1 promotes cancer cell growth and increases invasiveness. Most of CAV1 overexpression occur in hepatocellular carcinoma,[ and more studies should be discuss the relationship between CAV1 and CCA. There are many types of peptidases and they are highly expressed in many cancers. For example, kallikrein-related peptidases overexpressed in most digestive system tumor by regulating cell growth, angiogenesis, and metastasis.[ In endometrial carcinoma, dipeptidyl peptidase IV overexpression induce hypoxia-inducible factor 1a-vascular endothelial growth factor A signaling pathway to cause the cancer.[ So we speculated that the positive regulation of peptidase activity might promote the development of CCA. A growing number of evidence points show the contribution of SNCA to the etiology of Parkinson's disease.[ However, there are few studies on the relationship between genes SNCA and tumors. In colorectal cancer, SNCA were significantly hypermethylated in tissue samples and stool samples.[ But its effect on CCA is still unclear, and further research is needed.

Our research results are affected by some factors. On the one hand, the sample size of the study is small, and we only explore one dataset. On the other hand, CCA is related to some risk factors, such as clonorchis sinensis infection and intrahepatic bile duct stones. But these factors can not be ruled out. These factors may cause the difference of the results.

Some aspects are restricted in this study. First of all, the study lacks further experimental verification of the effect of abnormal methylation gene expression and function. Secondly, due to the use of bioinformatics arrays and accessible data tools, we did not investigate clinical parameters and prognosis. Third, because only two microarray patterns were analyzed, the sample size was not large enough. Therefore, a large sample study is needed to verify the research results. The etiology was not analyzed in our study. Therefore, supplementary molecular experiments should be encouraged to further validate our research results

In summary, we found MDEGs and signaling pathways that are connected with the CCA through the combined analyses of gene expression and methylation profile microarray data. New information of MDEGs in CCA was revealed in this study. This information may help to improve the understanding of the epigenetic regulation mechanism of CCA occurrence and development. In addition, MYC and VWF act as hub genes in the network module related to MDEG. These genes may have potential predictive and prognostic value, and can be used as methylation-based biomarkers for accurate diagnosis and treatment of CCA.

Conclusion

In conclusion, we found MDEGs and signaling pathways that are connected with the CCA by using a series of bioinformatics databases and tools. MYC and VWF act as hub genes of CCA in 81 Hyper-LGs and 76 Hypo-HGs. These hub genes may help to accurately diagnose and treatment of cholangiocarcinoma. We need to explore further to confirm the function of identified genes in cholangiocarcinoma.

Author contributions

CReDIT: Conceptualization, Y.L.; Data curation, G.L.; Formal analysis, H.T.; Funding acquisition, Y.L.; Investigation, X.Z.; Methodology, X.Z.; Project administration, G.L.; Resources, G.L.; Software, H.T.; Supervision, Y.L.; Validation, Y.L.; Visualization, Y.L.; Writing - original draft, X.Z.; Writing - review & editing, G.L. and X.Z.

Conceptualization: Yiling Li.

Data curation: Hae Kim, Lin Guan.

Formal analysis: Haoyu Tian.

Funding acquisition: Yiling Li.

Investigation: Xinhe Zhang.

Methodology: Hae Kim, Xinhe Zhang.

Project administration: Lin Guan.

Resources: Lin Guan.

Software: Hae Kim, Haoyu Tian.

Supervision: Yiling Li.

Validation: Yiling Li.

Visualization: Hae Kim, Yiling Li.

Writing – original draft: Hae Kim, Xinhe Zhang.

Writing – review & editing: Lin Guan, Xinhe Zhang.