Five Hub Genes Can Be The Potential DNA Methylation Biomarkers For Cholangiocarcinoma Using Bioinformatics Analysis

Introduction

Cholangiocarcinoma is featured with biliary epithelial differentiation. As the second most primary hepatic cancer, cholangiocarcinoma is rising globally.¹ Nowadays, we classify CCA into three types based on anatomical location: intrahepatic cholangiocarcinoma (ICC), hilar cholangiocarcinoma and extrahepatic cholangiocarcinoma (ECC).² ICC is a special cholangiocarcinoma kind which is located proximally to the bile ducts, identified to be a highly aggressive malignancy worldwide. Though the overall incidence is low, treatment options for ICC patients are limited.³ For advanced CCA patients, the chemotherapy regimen of gemcitabine and cisplatin is the main treatment. Unfortunately, most CCAs showed resistant to chemotherapy and relapse or metastasize quickly after chemotherapy treatment.^4,5 Thus, it is for us to ascertain novel therapeutic targets and effective treatment options for this disease.

In carcinogenesis, methylation could be a symbolic event through epigenetic studies.^6,7 Several kinds of methylation, such as promoter methylation and CpG island methylation, are strongly associated with CCA.⁸ Therefore, we tried to discover novel DNA methylation biomarkers in CCA patients and provide potential therapeutic targets for this indicating disease. Through analyzing GEO expression profiling data, we detected the downregulated genes between CCA and normal liver tissues. Gene Ontology (GO) functional annotation analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed for the screened DEGs. Next, we used a protein–protein interaction (PPI) network to identify hub suppress genes related to CCA. Then, through MEXPRESS database, we performed methylation levels detection to find out that F2, AHSG, ALDH8A1, SERPIND1 and AGXT were higher methylated in promoter in CCA compared to normal samples. Besides, these results were also further proved in CCA cell lines using RT-PCR. Therefore, we put forward that these five genes may be the DNA methylation biomarkers in CCA.

Methods And Materials

Cell Culture

In this study, we used four kinds of CCA cell lines (CCLP-1, HuCCT-1, HCCC-9810, and RBE) and a normal human biliary epithelial cell line (HIBEC). We purchased these cell lines from Cell Bank of Type Culture Collection of Chinese Academy of Sciences (Shanghai, China). According to their supplier’s protocol, all cell lines were cultured in RPMI-1640 complete medium (Biological Industries, Kibbutz Beit-Haemek, Israel) supplied with 10% fetal bovine serum (FBS; Moregate Biotech, Brisbane, Australia) and were cultured in an incubator of 37°C and 5% CO₂ environment.

Data Source

Datasets we analyzed in this research were obtained mainly from GEO database (http://www.ncbi.nlm.nih.gov/geo/) after searching for keywords related to CCA. Totally, we obtained 101 series about CCA from GEO database. Then, we selected three separate gene expression profiles (GSE107943, GSE26566, and GSE119337) for our study, about 228 human cholangiocarcinoma samples information were retrieved from these profiles. Among them, GSE107943 was based on the GPL18573 platform, GSE26566 was based on GPL6104, and GSE119337 was based on platform GPL11154.

Data Processing Of DEGs

In order to find out the DEGS between CCA and normal liver samples, we used the GEO2R online analysis tool in NCBI (https://www.ncbi.nlm.nih.gov/geo/geo2r/). We also used this tool to calculate the adjusted P-value and |logFC|. We considered the genes which meet cutoff standard requirements (adjusted P<0.05 and |logFC|≥2.0) as DEGs. Statistical analysis was carried out for each dataset. The Venn diagram was performed by using the web tool (bioinformatics.psb.ugent.be/webtools/Venn/).

GO And KEGG Pathway Analysis Of DEGs

Gene Ontology (GO) is used widely in functional annotation and enrichment analysis; biological process (BP), molecular function (MF), and cellular component (CC) are the three major components of gene function. KEGG is a database resource for collecting large number of data about molecular-level information, biological pathways, chemical substances which is generated by high-throughput experimental technologies. Through the Database for Annotation, Visualization and Integrated Discovery (DAVID) tools (https://david.ncifcrf.gov/), we carried out the GO annotation analysis and KEGG pathway enrichment analysis of DEGs. P<0.01 and gene counts≥10 were considered statistically significant.

PPI Network Construction And Hub Gene Identification

In order to obtain a PPI map, we upload the data of DEGs to the database analysis platform named Search Tool for the Retrieval of Interacting Genes (STRING). If the PPI pairs have the combined score >0.4, we extracted these pairs and visualize the PPI network using Cytoscape software (www.cytoscape.org/). Nodes with higher degree of connection, which were more essential for maintaining the stability of the entire network. The top 10 genes in the central index are considered to be core candidate genes.

RNA Isolation And Quantitative Real-Time PCR (qRT-PCR)

We isolated the total RNA using TRIzol reagent (TaKaRa, China) and then it was transcribed into cDNA using the PrimeScript RT Reagent Kit (TaKaRa, China). With an SYBR Green PCR Kit (Takara, China), the quantitative real-time PCR was carried out by an ABI 7500 FAST RT-PCR system (Applied Biosystems, USA). Relative quantification of mRNA expression was based on the 2ΔΔCt method after normalization to the endogenous reference GAPDH. Table S1 presented the primers we used in this research.

The Correlation Between Gene Expression And Methylation Around In The Promoter Region

MEXPRESS is a web tool which could offer the visualized analysis of clinical data, the expression (normalized RNASeqV2 value), and methylation TCGA and detect the relationship between them for one single gene in the specific tumor type. In this webtool, it executed the Pearson correlation to evaluate the difference between expression value and methylation level. When MEXPRESS faced with clinical parameter contains only two levels, the system will use a P value to compare the difference. The false discovery rate was used to correct for multiple comparisons.

Results

Identification Of DEGs

In this study, we selected three groups of gene expression profiles (GSE107943, GSE26566, and GSE119337). As shown in Table 1, GSE107943 contained 31 CCA specimens and 30 normal specimens, GSE26566 contained 104 CCA specimens and 65 normal liver specimens, and GSE119337 included 30 CCA samples and 30 normal liver samples. Matched with the criteria of P<0.05 and |logFC|≥2, 576 upregulated genes and 851 downregulated genes were identified through comparing the CCA specimens to liver normal specimens. In GSE26566 expression profile, 1076 DEGs were identified; 463 genes were upregulated and 613 genes were downregulated. And from GSE65194 dataset, 586 DEGs including 212 upregulated genes and 374 downregulated genes were discovered. The DEGs levels were illustrated by Heatmap in Figure S1. Subsequently, we carried out Venn analysis to get the intersection of the DEG profiles. On the whole, 115 DEGs were downregulated in all three groups (Figure 1A).

Table 1 Statistics Of The Three Microarray Databases Derived From The GEO Database

Figure 1 (A) Venn diagram of downregulated DEGs common to all three GEO datasets. (B) Protein–protein interaction network constructed with the differentially expressed genes.

Functional Enrichment Analyses Of Downregulated DEGs

Then, we evaluate the downregulated DEGs biological functions in CCAs, GO function and KEGG pathway enrichment analysis were carried out by DAVID (Table 2). CC, BP, and MF ontologies are the three components of enriched GO terms. The GO analysis results showed that downregulated DEGs were mainly enriched in CCs, including extracellular exosome, extracellular space, mitochondrion, and blood microparticle. The results of MF analysis indicated that the downregulated DEGs were significantly enriched in serine-type endopeptidase activity, calcium ion binding, and serine-type endopeptidase inhibitor activity. In addition, terms associated with complement and coagulation cascades, metabolic pathways, peroxisome, and biosynthesis of antibiotics were obtained from the results of KEGG pathway enrichment analysis.

Table 2 Significantly Enriched GO Terms And KEGG Pathways Of DEGs

PPI Network Construction And Hub Gene Identification

As shown in Table 2, 28 terms for KEGG pathway enrichment analysis were significantly over-represented among these downregulated DEGs. Then, the probability of relationships between pathways was evaluated by STRING tools. The PPI network was presented in Figure 1B. Besides, the top 10 genes evaluated by connectivity degree in the PPI network were identified (Table 3). From Table 3, we found that the degree of Coagulation factor II⁹ (F2) was 27, which ranked first among these 10 genes. Then, the following genes were plasminogen¹⁰ (PLG; degree= 21), enoyl-CoA hydratase and 3-hydroxyacyl CoA dehydrogenase¹¹ (EHHADH; degree= 20), alanine-glyoxylate and serine–pyruvate aminotransferase¹² (AGXT; degree= 15), alpha 2-HS glycoprotein (AHSG; degree= 15),¹³ carboxypeptidase B2¹⁴ (CPB2; degree= 14), aldehyde dehydrogenase 8 family member A1¹⁵ (ALDH8A1; degree= 14), serpin family D member 1 (SERPIND1; degree= 12), serpin family F member 2¹⁶ (SERPINF2; degree= 11), serine hydroxyl-methyltransferase 1¹⁷ (SHMT1; degree= 11). All of these hub genes were downregulated in CCA.

Table 3 Top 10 Hub Genes With Higher Degree Of Connectivity

F2, AHSG, ALDH8A1, SERPIND1 And AGXT Might Be The Potential DNA Methylation Biomarkers For CCA

To further clarify the expression of these hub genes in CCA and normal tissues, we performed the differential expression analysis by using GEPIA tools. From the box-profiles in Figure 2, these results indicated that all these 10 hub genes were lower expressed in CCA than normal liver tissues (p<0.05), in other words, 10 hub genes could be the suppressor gene in CCA. DNA promoter hyper-methylation is an important pattern which impact tumorigenesis of cancers, and this pattern could be found in various tumor suppressor genes.¹⁸ Therefore, we speculated whether these 10 hub genes were associated with DNA promoter methylation in CCA. Then, the methylation levels of hub genes in CCA were detected by MEXPRESS database. The results revealed that the promoter region of F2, AHSG, ALDH8A1, SERPIND1 and AGXT showed significantly higher methylation levels (p<0.05) in CCA compared to normal tissues (Figure 3). However, the left 5 genes showed no significant methylation level change between tumor and normal tissues (Figure S2). In addition, we detected the relationship between the five hub genes and the survival of CCA patients. Using GEPIA database, as Figure 4 showed, higher expression of these five genes would predict the better prognosis of CCA patients.

Figure 2 The expression of 10 hub genes in CCA samples and normal tissue samples through GAPIA database (*p<0.05).

Figure 3 Visualization of the TCGA data for F2, AHSG, ALDH8A1, SERPIND1, and AGXT in cholangiocarcinoma using MEXPRESS. Samples were ordered by their expression value. This figure showed the negative correlation between hub gene expression and promoter methylation (highlighted part with red box), with the Pearson correlation coefficients on the right (*p<0.05, **p<0.01, ***p<0.001).

Figure 4 The five hub genes (F2, AHSG, ALDH8A1, SERPIND1, and AGXT) impacted the survival rate of CCA patients.

Then, we evaluated F2, AHSG, ALDH8A1, SERPIND1 and AGXT expressions level in a panel of five cell lines: four CCA (CCLP-1, HuCCT-1, HCCC-9810 and RBE) cell lines and one non-tumorigenic or benign (HIBEC) cell line. The mRNA expression measured by RT-PCR revealed that F2, AHSG, ALDH8A1, SERPIND1 and AGXT transcription levels in cancerous cell lines were lower than those in the non-tumorigenic or benign control cells (Figure 5A). Moreover, the expression of five genes above in CCLP-1, HuCCT-1, HCCC-9810 and RBE cells was significantly increased after treatment with 5-Aza, a DNA methyltransferase inhibitor (Figure 5B). And, we also analyzed the methylation sequencing data form TCGA, we found out that compared with the adjacent normal liver tissues the CCA have the higher methylation levels in the promoter (Figure S3). All these results indicated that F2, AHSG, ALDH8A1, SERPIND1 and AGXT may be the potential DNA methylation biomarkers for CCA.

Figure 5 (A) Results of quantitative real-time PCR revealed that F2, AHSG, ALDH8A1, SERPIND1, and AGXT transcription levels in four cancerous cell lines and benign cells. Error bars represent the SD from three independent experiments. (B) Quantitative real-time PCR indicated that 5-Aza treatment for 96 hrs increased the expression of F2, AHSG, ALDH8A1, SERPIND1, and AGXT in CCA cells (*p<0.05, **p<0.01).

Discussion

In order to identify the potential key genes associated with CCA, gene expression and protein–protein interaction analysis were executed in our study. We found 115 downregulated DEGS in CCA compared with normal tissue by comparing gene expression profiling data from the three GEO dataset. Then, we constructed the PPI network to detect the inner-relationship of the downregulated DEGs, and 10 hub genes were identified, including F2, PLG, EHHADH, AGXT, AHSG, CPB2, ALDH8A1, SERPIND1, SERPINF2 and SHMT1.

In various cancers, DNA methylation plays an important role in the gene expression regulation. Several methylated genes are often shared among the gastrointestinal tract.^6,19–21 Then we identified the relationship between DEGs expression and DNA methylation using MEXPRESS. As Figure 4 showed, default MEXPRESS plot for hub genes in CCA with the samples sorted based on their expression value. The MEXPRESS analysis results have indicated that there was a negative correlation between the expression and methylation around the promoter region which is highlighted with red box, which showed that the expressions of F2, AHSG, ALDH8A1, SERPIND1 and AGXT might be regulated through promoter DNA methylation (the left 5 hub genes showed no significant methylation level change between tumor and normal tissues).

Totally, we have shown the panel in this research included F2, AHSG, ALDH8A1, SERPIND1, and AGXT might be the novel biomarkers for CCA. Although the role of these genes in cancer is not fully understood, the high performance of the identified five‐gene biomarker panel suggests that it could be suitable for discriminating the DNA methylation expression changes during the carcinogenetic process in bile tract epithelium. Some research have shown that DNA methylation always occurred in the early stage of cancers.²⁰ This panel may be of early diagnostic value for CCA patients, and we can combine this biomarker panel of this study with bile duct brush biopsy technology, which has the potential to diagnose CCA at an early stage that can be cured by surgery or transplantation.²²

Conclusion

Through our bioinformatics analysis, we identified 112 genes were downregulated in CCA. Based on the GO analysis and PPI network, 10 hub genes were found out from these DEGs, including F2, PLG, EHHADH, AGXT, AHSG, CPB2, ALDH8A1, SERPIND1, SERPINF2, and SHMT1. Epigenetics is the main factor which causes the reduction of tumor suppressor gene expression in tumor, and DNA methylation plays an important role in the epigenetic modification. Using MEXPRESS and RT-PCR, we found that F2, AHSG, ALDH8A1, SERPIND1 and AGXT could be methylated in the promoter. These 5 genes may be the potential DNA methylation biomarkers and therapy targets for CCA.