Integrated bioinformatics identifies the dysregulation induced by aberrant gene methylation in colorectal carcinoma.

Colorectal carcinoma (CRC) is one of the most common cancers, and is associated with a poor clinical outcome. The key genes and potential prognostic markers in colorectal carcinoma remain to be identified and explored for clinical application. DNA expression/methylation profiles were downloaded from the Gene Expression Omnibus (GEO) database to identify differentially expressed/methylated genes (DEGs and DEMs). A total of 255 genes and 372 genes were identified as being up-regulated and down-regulated, respectively, in GSE113513, GSE81558, and GSE89076. There were a total of 3350 hypermethylated genes and 443 hypomethylated genes identified in GSE48684. Twenty genes were found to be hypermethylated as well as down-regulated, and a functional enrichment analysis revealed that these genes were mainly involved in cancer-related pathways. Among these 20 genes, GPM6A, HAND2 and C2orf40 were related to poor outcomes in cancer patients based on a survival analysis. Concurrent decreases of GPM6A, HAND2 and C2orf40 protein expression were observed in highly-differentiated colorectal carcinoma tissues, and higher expression levels were found in undifferentiated or minimally-differentiated colorectal carcinoma tissues. In conclusion, 20 genes were found to be downregulated and hypermethylated in CRC, among which GPM6A, HAND2 and C2orf40 were explored for their potential prognostic value.

Introduction

Colorectal carcinoma is one of the most common cancer types worldwide, accounting for 9% and 7% all new cancer diagnoses in men and women, respectively, and ranking third in estimated new cancer cases and deaths in the United States. Many risk factors may contribute to the development of CRC, such as age, excessive alcohol use, a sedentary lifestyle, and obesity. In addition, specific genes containing sporadic, inherited and familial mutations can lead to the development of CRC., The choice off treatment for CRC patients varies according to the clinical characteristics, including the tumor size, presence of metastases, stage and personal status. Surgical resection and chemotherapy are often combined to treat patients. Monoclonal antibodies and agents targeting vascular endothelial growth factor (VEGF) and the epidermal growth factor receptor (EGFR) also show great prospects in CRC treatment. Although therapies for CRC have improved over the last few decades, the overall survival (OS) and quality of life for patients remains relatively poor, with a median OS of approximately 30 months.

Aberrant methylation of genes has been widely reported in cancers, and not only affects carcinogenesis, but also contributes to metastasis., In this study, we performed analyses that integrated differentially-methylated genes and differentially-expressed genes to better understand the impact of aberrant methylation of gene expression and the importance in CRC. A functional enrichment analysis and PPI network were applied to retrieve core pathways and genes. Finally, a survival analysis was performed on the downregulated and hypermethylated genes, and immunohistochemistry was used to validate the expression levels of gene products with prognostic value. A flow chart summarizing this study is provided as Fig. 1.

Figure 1

A flow chart of the present study.

Methods

Data acquisition and processing

DNA methylation data for GSE48684 were obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), which includes 41 normal colon samples, 42 colon adenomas, and 64 colorectal cancers. The methylation data were based on the Illumina Human Methylation 450 Bead Chip platform. Gene expression data were obtained from the GEO database under accession numbers GSE89076, GSE81558 and GSE113513., GSE89076 includes 41 colorectal carcinoma tissues and 30 adjacent normal tissues, with data based on the Agilent-039494 SurePrint G3 Human GE v2 8 × 60 K Microarray 039381. GSE81558 includes 23 sporadic colorectal adenocarcinomas and 19 liver metastases from 23 patients with liver metabolism and 9 non-tumoral colorectal tissues. The 23 sporadic colorectal adenocarcinomas and 9 non-tumoral tissues were selected for further analysis. GSE113513 contains 14 matched colorectal cancer tissue and normal tissue samples from 14 patients. Both GSE81558 and GSE113513 were based on the Affymetrix Human Gene Expression Array.

Data processing

After RMA background correction, log2 logarithmic transformation and quantile normalization, the probe ID of the gene expression microarrays was matched to the gene symbol, and the average values of genes with several probes were calculated as the gene expression levels, as fulfilled by applying the R software Affy package. The R software limma package was used to screen for aberrantly expressed genes. For the DNA methylation data, after calculating normalized betas from data downloaded from the GEO database using the R software wateRmelon package, limma packages were applied to compare the differences between colorectal carcinoma tissues and normal tissues.

Functional enrichment analysis

Functional enrichment analysis is an effective method to obtain general knowledge about a large number of genes, including gene ontology (GO) term enrichment and KEGG pathway enrichment. The GO enrichment analysis consists of three parts: the biological process, cellular component and molecular function. The Gene Ontology (http://geneontology.org/) and KEGG Pathway (https://www.genome.jp/kegg/pathway.html) databases were used to detect the enriched terms and pathways related to DEGs. A P value was calculated based on the cumulative hypergeometric distribution. Terms and pathways with values of P < 0.05 were considered to be statistically significant.

PPI network construction

To detect the interactions among DEGs, the STRING v11 online database was applied to predict the protein-protein interactions, which include not only experimentally-validated interactions, but also interactions among functionally-grouped proteins maintained in KEGG pathways. Protein-protein interactions with a combined score higher than 0.7 were selected for PPI network construction, and the network was visualized using online software (REF). The degree of connectivity was calculated to identify the most important nodes in the network, namely hub genes.

Survival analysis

The TCGA online database has accrued more than 10,000 primary samples from patients 52 types of tumors (https://portal.gdc.cancer.gov/). Different types of data are included that were generated using several different techniques, such as gene expression profiling, microRNA profiling, DNA methylation profiling, copy number variation profiling, SNP genotyping and so on. Gene expression profiling of colorectal carcinoma and matched clinical traits data were downloaded, and samples with both of these datasets were selected for further survival analysis. A total of 616 samples were obtained, the patients were classified into two groups: those with expression values greater than the median were assigned to the high expression group, while those with values lower than the median were assigned to the low expression group.

Immunohistochemistry

Formalin-fixed paraffin-embedded blocks of 50 colorectal carcinoma tissue samples and 50 adjacent normal tissue samples were obtained from the Department of General Surgery of the Second Affiliated Hospital of Soochow University. Written informed consent was obtained for the collection and analysis of these samples. The study was approved by the Human Ethics Committee of the Second Affiliated Hospital of Soochow University. A Novus antibody (NBP1-81799) for GPM6A, an R&D antibody (AF3876-SP) for HAND2 and an ImmunoWay antibody for C2orf40 (YT1457) were used at a 1:50 dilution for immunostaining of the paraffin-embedded colorectal carcinoma tissues.

The protein staining was evaluated by board-certified pathologists who had no knowledge of the patients’ clinical status. The expression scores for the proteins were given separately for the staining intensity and the area of staining. The intensity of staining was scored as follows: weak, 1; strong, 2. The area of staining was quantified as follows: <50%, 1; >50%, 2. Each sample received a final score that was the product of the intensity and area scores.

Results

Identification of differentially expressed genes and methylated genes

Comparisons of the gene expression of tumor and normal tissues from the 3 GEO gene expression profiles were performed using the R software limma package, and mRNAs with a |log2 (fold change)|>1, p-value < 0.05 and q-value < 0.05 were identified as differentially expressed genes. Heatmaps and volcano plots are presented in Figure 2, Figure 3. A total of 3420 DEGs were identified in GSE89076, consisting of 1733 up-regulated genes and 1687 down-regulated genes. A total of 1316 aberrantly expressed genes were identified in GSE113513, including 524 up-regulated genes and 792 down-regulated genes. There were 428 up-regulated genes and 604 down-regulated genes in GSE81558, with 1032 genes found to be aberrantly expressed in total. Further, 255 and 372 genes were detected to be up- and down-regulated in all three profiles, as shown in Fig. 4, and were used for further analysis. After applying exclusion criteria as |log2 (fold change) |> 0.263, p-value < 0.05 and q-value < 0.05, 3350 genes were considered to be hypermethylated and 443 genes were hypomethylated in GSE48684, as shown in Fig. 2D and Fig. 3D.

Figure 2

The differentially expressed/methylated mRNAs (A) volcano plot of GSE89076; (B) volcano plot of GSE81558; (C) volcano plot of GSE13513; (D) volcano plot of GSE48684.

Figure 3

Heatmaps for the DEGs and DEMs in the cancer versus normal tissue. (A) GSE89076; (B) GSE81558; (C) GSE13513; (D) GSE48684.

Figure 4

A Venn plot of the differentially expressed genes. (A) A Venn plot of the upregulated genes in the three expression profiles; (B) a Venn plot of the downregulated genes in the three expression profiles.

Functional annotation of DEGs

We utilized the Gene Ontology database and KEGG pathway database to find terms and pathways associated with the development of colorectal carcinoma. Notably, the upregulated DEGs were mainly involved in homeostasis-related terms, while the enriched pathways involved the cell cycle, p53 signaling, TGF-beta signaling and Wnt signaling, which were closely related to the development of cancer, as shown in Fig. 5A and B. The downregulated DEGs mainly involved metabolism-related pathways at the KEGG level, including metabolic pathways, pancreatic secretion, mineral absorption, nitrogen metabolism, caffeine metabolism and tyrosine metabolism, as shown in Fig. 5C and D. The KEGG pathway DEMs mainly participated were also predicted (Fig. 5E and F). Notably, hypermethylated genes mainly participated in neuroactive ligand receptor interactions, cell adhesion, CAMP signaling, calcium signaling and cholinergic synapses at the KEGG pathway level.

Figure 5

Results of the GO and KEGG pathway enrichment analyses for DEGs and DEMs. (A) GO enrichment for upregulated DEGs; (B) KEGG pathway enrichment for upregulated DEGs; (C) GO enrichment for downregulated DEGs; (D) KEGG pathway enrichment for downregulated DEGs; (E) KEGG pathway enrichment for upregulated DEMs; (F) KEGG pathway enrichment for downregulated DEMs.

Constructing a PPI network of DEMs

Using a PPI network to visualize protein-protein interactions involved in the development of colorectal carcinoma is a reliable way to screen for potential hub genes. When we set the required confidence (combined score) > 0.7, there were a total of 960 interactions among 136 upregulated DEGs that were identified. The genes with the highest degree of interactions within the network were CDK1, CCNB1, BUB1, TOP2A, and MAD2L1. Similar processing was performed on the downregulated DEGs, and 267 interactions among 136 nodes were identified, with GNG7, LPAR1, SST, GNAI1, NPY, CXCL12, and CCL5 having the highest degree of connectivity (Fig. 6A and B).

Figure 6

The PPI networks. (A) PPI network of the upregulated genes; (B) PPI network of the downregulated genes. The nodes represent dysregulated genes, the edges represent interactions. PPI: protein-protein interactions.

Identifying hypermethylated genes with low expression, and defining the related terms and pathways

By comparing 1106 promoter-associated hypermethylated genes with the 372 down-regulated genes, we identified 20 genes that were both hypermethylated and had low expression (Hyper-LEGs), including ACADS, ANO5, BEND5, BMP2, C2orf40, CHL1, EIF4E3, EPB41L3, FOXF2, GPM6A, HAND2, LAMA1, LIFR, NPY, NR3C1, NR5A2, PRKCB, PTPRR, STOX2 and TCF21. These were used for further GO and KEGG enrichment analyses, as shown in Fig. 7A. These 20 Hyper-LEGs had functions related to amoebiasis, cancer development, focal adhesion, cytokine-cytokine receptor interactions, MAPK signaling and neuroactive ligand-receptor interactions, as presented in Fig. 7B and C.

Figure 7

The downregulated and hypermethylated genes (Hyper-LEGs). (A) A Venn plot of the downregulated DEGs and hypermethylated DEMs; (B) The GO enrichment analysis of downregulated and hypermethylated genes; (C) A KEGG pathway enrichment analysis of downregulated and hypermethylated genes.

Identification of genes related to the prognosis of CRC and validation by immunohistochemistry

A survival analysis was performed to evaluate the effects of these Hyper-LEGs on the prognosis of CRC. The upregulation of GPM6A, HAND2 and C2orf40 were closely related to a poorer OS, as shown in Fig. 8. An immunohistochemical analysis was performed on 50 colorectal carcinoma tissues and the 50 corresponding adjacent normal tissues. Concurrent decreases in GPM6A, HAND2 and C2orf40 protein expression were observed in highly-differentiated colorectal carcinoma tissues, and higher expression was observed in the poorly-differentiated colorectal carcinoma tissues, as shown in Fig. 9.

Figure 8

The results of the survival analysis. (A) GPM6A; (B) HAND2; (C) C2orf40.

Figure 9

Concurrent decreases in the protein expression of GPM6A, HAND2 and C2orf40 were observed in colorectal carcinoma tissue samples by immunohistochemistry. (A) GPM6A; (B) HAND2; (C) C2orf40.

Discussion

In this study, we integrated multiple expression profiles and methylation profiles from CRC patients to identify genes that play a crucial role in the onset and development of CRC. Twenty genes that were both downregulated and hypermethylated were identified, and were found to be mainly enriched in cancer formation- and progression-related pathways, focal adhesion and MAPK signaling. To evaluate the prognostic impact of these Hyper-LEGs, a survival analysis was conducted. Among these, GPM6A, HAND2 and C2orf40 were found to be significantly related to the survival of CRC patients. The expression levels of the three gene products was confirmed by immunohistochemistry.

By comparing three databases of CRC expression profiles (GSE113513, GSE81558, GSE89076), we finally identified 255 and 372 genes that were consistently up- or down-regulated in all three profiles. A functional enrichment analysis demonstrated that upregulated DEGs were mainly related to pathways involving the cell cycle, p53 signaling, TGF-beta signaling and Wnt signaling at the KEGG pathway level. The downregulated DEGs were enriched in metabolism-related pathways. We detected a total of 3350 hypermethylated and 443 hypomethylated genes. The hypermethylated genes mainly participated in neuroactive ligand receptor interactions, cell adhesion, CAMP signaling, calcium signaling and cholinergic synapses at the KEGG pathway level. Of note, impaired cell adhesion may promote cancer migration and metastasis (REF). The pathways identified have all been reported to play a crucial role in the onset of CRC or in driving its progression.,, Dysregulated of the WNT pathway leads to aberrant cell growth and stem cell differentiation (REF). It also impairs cellular adhesion, and contributes to cell migration and metastasis. MAPK signaling plays an important role in cell proliferation and differentiation. Abnormal MAPK signaling pathway promotes cell proliferation and survival, thus enhancing CRC development (REF). Activation of the MAPK signaling pathway also weakens cell-cell adhesions via induction of the epithelial-mesenchymal transition (EMT) and TGF-beta signaling, thus promoting metastasis. The MAPK signaling pathway is also involved in inflammation, which may contribute to the initiation and progression of CRC.,

A PPI network was constructed to visualize the protein-protein interactions among DEGs and to pinpoint hub genes in the network. GNG7, LPAR1, SST, GNAI1, NPY, CXCL12, CCL5, CCL19, CXCL13, CCL21, CCL28, INSL5 and P2PY14 were the top 13 genes with the highest connectivity in the network.

We further detected genes which were both downregulated and hypermethylated. Twenty genes (ACADS, ANO5, BEND5, BMP2, C2orf40, CHL1, EIF4E3, EPB41L3, FOXF2, GPM6A, HAND2, LAMA1, LIFR, NPY, NR3C1, NR5A2, PRKCB, PTPRR, STOX2 and TCF21) were identified. A KEGG analysis indicated that these genes were related to amoebiasis and cancer-related pathways, such as focal adhesion and MAPK signaling. The upregulation of GPM6A, HAND2 and C2orf40 were found to contribute to poorer prognosis in CRC patients.

GPM6A is a transmembrane protein which is abundant on the surface of neuronal cells. GPM6A has also been reported to be an oncogene in human lymphoid leukemia, inducing lymphocyte proliferation. Heart and neural crest derivatives express 2 (HAND2), which encodes transcription factors that play a crucial role in the development various organs, especially the heart. We found that HAND2 was both downregulated and hypermethylated. Hypermethylation of HAND2 has previously been reported in cervical cancer and lung cancer, as well as colon cancer and rectal cancer.23, 24, 25 HAND2 methylation has been confirmed to affect the development of endometrial cancer, but how it exerts its effects on colorectal cancer remains largely unknown. It has been reported that c2orf40 (chromosome 2 open reading frame 40, also called esophageal cancer-related gene 4 (ECRG4) or augurin) is a tumor suppressor, which is hypermethylated in various types of cancer, including esophageal cancer, colorectal carcinoma, glioma and breast cancer.,, Overexpression of C2orf40 impairs colon cancer cell proliferation, and re-expression of the silenced C2orf40 gene can restore its inhibition of colon cell growth. Intriguingly, although GPM6A, HAND2 and C2orf40 were all downregulated and hypermethylated in tumor tissues, low expression levels of these three genes were all associated with significantly longer survival, indicating that the functions of these genes may differ during cancer initiation and development. Immunohistochemical staining showed that concurrent decreases of GPM6A, HAND2 and C2orf40 protein expression were observed in highly-differentiated colorectal carcinoma tissues, with higher expression in poorly-differentiated colorectal carcinoma tissues, which may explain the phenomenon. However, the exact mechanisms underlying the effects on the prognosis need further exploration.

Conclusion

In this study, we integrated expression and methylation profiles to identify key genes involved in CRC development. We identified 20 genes that were both downregulated and hypermethylated in CRC. Among these, GPM6A, HAND2 and C2orf40 were of potential prognostic value. However, how GPM6A, HAND2 and C2orf40 function in the onset and development of CRC requires further study.

Availability of data and materials

The dataset analyzed in the current study is available on GEO database.

https://www.ncbi.nlm.nih.gov/search/all/?term=GSE113513

https://www.ncbi.nlm.nih.gov/search/all/?term=GSE81558.

https://www.ncbi.nlm.nih.gov/search/all/?term=GSE89076.

https://www.ncbi.nlm.nih.gov/search/all/?term=GSE48684.

Authors contribution

Wei Chen and Chungen Xing designed the study; Zhenyu Ye and Yecheng Li wrote the paper; Jiaming Xie, Zhenyu Feng, Xiaodong Yang and Yong Wu collected and analyzed data; Yuwei Pu, Jiawei Gao, Xiangrong Xu, Zhaobi Zhu and Wei Li created the color figures in the paper; all authors contributed to the intellectual contents and approved the final version of the manuscript.

Ethics approval and consent to participate

The experimental protocol was established according to the ethical guidelines of the Helsinki Declaration and was approved by the Human Ethics Committee of the Second Affiliated Hospital of Soochow University. Written informed consent was obtained from all individuals or their guardians.

Conflict of Interests

The authors declare that they have no competing interests related to this work.

Funding

The present study was partially supported by the [grant number 81672970], the [grant number BK20160338], the projects of [grant numbers SYS201552, SS201753 and SYS2018054], the [grant number SZYJTD201803], the [grant number KJXW2017013], Project [grant number KYCX19_1986] and the Project Funding.