Identification of MTHFD2 as a novel prognosis biomarker in esophageal carcinoma patients based on transcriptomic data and methylation profiling.

DNA methylation is an important epigenetic regulatory mechanism in esophageal carcinoma (EC) and is associated with genomic instability and carcinogenesis. In the present study, we aimed to identify tumor biomarkers for predicting prognosis of EC patients.We downloaded mRNA expression profiles and DNA methylation profiles associated with EC from the Gene Expression Omnibus database. Differentially expressed and differentially methylated genes between tumor tissues and adjacent normal tissue samples were identified. Functional enrichment analyses were performed, followed by the construction of protein-protein interaction networks. Data were validated based on methylation profiles from The Cancer Genome Atlas. Candidate genes were further verified according to survival analysis and Cox regression analysis.We uncovered multiple genes with differential expression or methylation in tumor samples compared with normal samples. After taking the intersection of 3 differential gene sets, we obtained a total of 232 overlapping genes. Functional enrichment analysis revealed that these genes are related to pathways such as "glutathione metabolism," "p53 signaling pathway," and "focal adhesion." Furthermore, 8 hub genes with inversed expression and methylation correlation were identified as candidate genes. The abnormal expression levels of MSN, PELI1, and MTHFD2 were correlated with overall survival times in EC patients (P < .05). Only MTHFD2 was significantly associated with a pathologic stage according to univariate analysis (P = .037) and multivariate analysis (P = .043).Our study identified several novel EC biomarkers with prognostic value by integrated analysis of transcriptomic data and methylation profiles. MTHFD2 could serve as an independent biomarker for predicting prognosis and pathological stages of EC.

Introduction

Esophageal carcinoma (EC) is one of the most common malignant cancers. EC is prevalent in Chinese populations and is responsible for a large proportion of all diagnosed EC cases worldwide. Despite some recent therapeutic improvements, most EC patients are diagnosed in the advanced stages, and the overall 5-year survival rate is <20%.[ Therefore, the development of novel therapeutic targets and molecular biomarkers for early detection and clinical management of EC is urgently needed.

Genomic alterations and epigenetic regulation play pivotal roles in carcinogenesis and development. DNA methylation is a dominant form of epigenetic modification involved in the regulation of the cancer genome.[ Hypermethylation at CpG islands or promoter regions often leads to the silencing of tumor suppressor genes, while hypomethylation in these regions is associated with upregulation of oncogenes.[ Large-scale and multi-omics profiling can provide a comprehensive analysis of the genomic and epigenomic aberrations in various cancers. In esophageal squamous cell carcinomas (ESCC), many cancer-related genes function in modulating DNA methylation. For example, the complete methylation of IGFBPL1 found in EC cells can suppress EC cell growth by inhibiting PI3K-AKT signaling.[ Promoter hypermethylation of TGFBR2 was identified in ESCC samples, and lentiviral mediated overexpression of TGFBR2 could inhibit ESCC cell proliferation.[ Through integrated analysis of DNA methylation and RNA sequencing data, Chen et al[ identified 16 key ESCC genes such as SIX4, CRABP2, and EHD3, which demonstrate an inverse correlation between DNA methylation and mRNA expression. Thus, integrated multi-omics data can accurately screen potential prognostic biomarkers and provide novel insights into precision medicine for cancer.

In the current study, through integrated analysis of transcriptomic and methylation data from the Gene Expression Omnibus (GEO) database, we screened differentially expressed and differentially methylated genes in EC samples compared with normal tissue samples. Functional enrichment analysis was conducted to identify enriched Gene Ontology (GO) terms and pathway categories of candidate genes, followed by protein–protein interaction (PPI) network construction. After further data validation, we focused on 3 genes with inverse correlations between expression level and DNA methylation status (MSN, MTHFD2, and PELI1). Survival and Cox regression analyses were performed to explore their prognostic values in EC patients. A schematic of the bioinformatics pipeline used to analyze transcriptomic data and DNA methylation profiles is presented in Fig. 1. Our systemic analysis provides new insights into the pathobiology of EC and may aid in the development of novel prognostic biomarkers for EC patients.

Figure 1

A schematic of the bioinformatics analysis pipeline used for transcriptome analysis and generation of DNA methylation profiles.

Methods

Data resources and screening of differentially expressed genes (DEGs)

Transcriptome data (accession numbers GSE20347,[ GSE75241,[ and GSE17351[) and DNA methylation profiles (GSE52826[) associated with EC including corresponding probe annotation information were downloaded from the GEO database. The GSE20347 mRNA microarray datasets included 34 clinical samples and were tested on an Affymetrix Human Genome U133A 2.0 Array platform, while GSE75241 consisted of 30 biopsy specimens and was verified on an Affymetrix Human Exon 1.0 ST Array platform. GSE17351 included 5 paired primary ESCC tumor and normal tissues tested on an Affymetrix Human Genome U133 Plus 2.0 Array platform. EC methylation data from The Cancer Genome Atlas database were used for data validation. Moreover, GSE52826 DNA methylation profiles were generated using Infinium methylation 450K BeadChips, and the datasets consisted of 12 clinical samples from ESCC patients and healthy individuals.

Additionally, the limma package[ was used to screen DEGs between tumor tissues and adjacent normal tissues. The minfi package[ was used to analyze the GSE52826 DNA methylation dataset. The criteria for DEG classification and differentially methylated genes were as follows: |Fold Change| >2 and adjusted P < .05, and |beta| >2 and adjusted P < .05, respectively. Volcano plots were generated to visualize upregulated and downregulated genes using the pheatmap package.[ Differences were deemed statistically significant when P < .05 and |log2fold (FC)| ≥2.

Functional enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/)[ is a bioinformatics resource for performing batch functional annotation according to GO terms relating to biological processes (BP), cellular components (CC), and molecular function (MF). Moreover, Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/)[ enrichment analyses were performed to identify potential pathways associated with these DEGs. Statistical significance was defined as P < .05.

Construction of PPI network

Protein interactions were searched for using the online tool Search Tool for the Retrieval of Interacting Genes (STRING) (http://www.string-db.org/, version 11.0)[ with a threshold correlation coefficient >0.4. Cytoscape software[ was used to construct the PPI network. The functional modules were screened from large protein interaction networks using the MCODE plugin.[ Statistical significance was defined as P < .05.

Survival analysis

Pearson correlation was used to evaluate correlations between mRNA levels and methylation levels of each gene. Coefficient values < 0 were considered statistically significant. In addition, we further performed survival analysis[ and constructed Cox regression models[ to screen candidate genes associated with EC progression. Based on univariate and multivariate analyses, we identified several genes with prognostic value.

Results

Identification of differentially expressed and methylated genes in EC

Transcriptome analysis and DNA methylation profiling were conducted using the limma and minfi packages. After data normalization, hierarchical clustering analysis of GSE20347 and GSE75241 showed that multiple genes were differentially expressed in EC tissue compared with normal samples (Fig. 2A and B). Subsequently, we found numerous differentially methylated genes from the GSE52826 dataset, and randomly selected 100 genes for clustering analysis. Heat map results indicate that these genes exhibited differential methylation patterns in EC samples compared with normal samples. The large number of points located on each side of the volcano plot indicates that many probes exhibited differential methylation status (Fig. 2C). We extracted the intersection of the 3 differential gene sets, and finally obtained a total of 232 overlapping genes with differential expression or methylation (Fig. 2D).

Figure 2

Identification of candidate genes in EC samples according to microarray data. A–C (left). Heat maps representing the hierarchical clustering results of the GSE20347, GSE75241, and GSE52826 datasets. The columns and rows represent clinical samples and differentially expressed genes, respectively. Red indicates upregulated or hypermethylated genes, while green refers to downregulated or hypomethylated genes. A–C (right). Volcano plots representing changes in gene expression at the transcriptional and methylation levels. The genes with | log2fold (FC)| ≥2 and P < .05 were respectively labeled with red or green. The green dots on the left side indicate significantly downregulated or hypomethylated genes, while the red dots on the right side indicate upregulated or hypermethylated genes. D. Venn diagram analysis showing the 232 overlapping genes that were selected as candidate genes for further analysis based on the 3 datasets. EC = esophageal carcinoma.

Functional enrichment analysis of genes related to EC

GO analysis was conducted to illustrate that EC-related proteins are involved in many BP, CC, and MF categories. By evaluating gene count numbers and FDR values, we identified the top 20 BPs such as “cytoskeleton,” “cell cycle,” “biological adhesion,” “cell adhesion,” “extracellular region part,” “M phase,” and “response to wounding” (Fig. 3A). KEGG enrichment analysis indicated that these candidate genes were associated with several signaling pathways (Fig. 3B) such as “small cell lung cancer” (hsa05222, count = 5), “p53 signaling pathway” (hsa04115, count = 7), “glutathione metabolism” (hsa00480, count = 7), “focal adhesion” (hsa04510, count = 9), “fatty acid metabolism” (hsa00071, count = 4), “ECM-receptor interaction” (hsa04512, count = 9), “DNA replication” (hsa03030, count = 4), “cell cycle” (hsa04110, count = 8), and “biosynthesis of unsaturated fatty acids” (hsa01040, count = 3). The signaling pathways are visualized in Fig. 3 using clueGO plugin.

Figure 3

Functional enrichment analysis of differentially expressed genes and differentially methylated genes in ECs. A. The top 20 functional terms (FDR < 0.05) were determined by GO term analysis evaluating the BP, CC, and MF. Red dots indicate a smaller P-value. The node size corresponds to the gene count, and larger dots indicate a higher number of genes with enriched GO terms. The intensity of dot color corresponds to the logarithmic P-value (a stronger red indicates higher statistical significance). B. KEGG pathway categories enriched in the list of differentially expressed genes. The columns represent pathway categories. The red color of the columns corresponds to logarithmic P-values. C. Significantly enriched pathway categories as determined by ClueGO plugin enrichment analysis. BP = biological process; CC = cellular component; EC = esophageal carcinoma; GO = Gene Ontology; KEGG = Kyoto Encyclopedia of Genes and Genomes; MF = molecular function.

PPI network analysis

We searched the STRING database and explored the potential interactions between the EC proteins. PPI networks were visualized using Cytoscape software (Fig. 4). The network consisted of multiple edges and nodes which represent the genes and interactions, respectively. We further identified several hub genes from the PPI network and constructed a regulatory network using the MCODE plugin (Fig. 4B). These hub genes were selected as candidate genes associated with EC development, and included MTHFD2, MCM6, CENPE, NCAPH, RRM2, PELI1, KIF15, CCNG2, and MGST2.

Figure 4

PPI network associated with differentially expressed genes in EC. A. Results of PPI network analysis. B. Regulatory network of hub genes related to EC. EC = esophageal carcinoma; PPI = predicted protein–protein interaction.

Eight genes are differentially methylated in EC

We analyzed EC methylation and transcriptomic data downloaded from the TCGA database. Using coefficient and P thresholds of <–.3 and <.05, respectively, we identified 8 genes with inverse correlations between mRNA expression and DNA methylation (Fig. 5). These 8 genes were hypermethylated and downregulated, or hypomethylated and upregulated in tumor tissue, and are related to EC progression (CCNB1, CENPF, KIF15, MCM7, MSN, MTHFD2, NCAPH, and PELI1).

Figure 5

Hub genes with inverse methylation and expression status (coefficient values < −0.3 and P < .05). The axes represent DNA methylation (X) and gene expression levels (Y). Rpearson is the Pearson correlation coefficient, and n is the number of clinical samples.

MTHFD2 is an independent risk factor for EC development

Survival analysis revealed that patients with high expression of MSN or PELI1 displayed longer overall survival times compared with those with low expression level of these genes (Fig. 6A; MSN, hazard ratio [HR]: 0.672, 95% confidence interval [CI] 0.169–2.681, P = .0468; PELI1, HR: 0.228, 95% CI 0.037–1.405, P = .0452). In contrast, EC patients with lower expression of MTHFD2 exhibited better prognosis than individuals with high MTHFD2 (HR: 2.394, 95% CI 0.312–18.345, P = .0361). These results suggest that MSN and PELI1 may function as tumor suppressor genes, while MTHFD2 may act as an oncogene.

Figure 6

Survival analysis and data validation identifying MTHFD2 as a prognostic biomarker in EC. A. Kaplan–Meier plots depicting survival analysis of three candidate genes: MSN, MTHFD2, and PELI1. B. Data validation to identify MTHFD2 expression in normal and tumor tissues based on the TCGA and GSE17351 datasets. EC = esophageal carcinoma.

In addition, we analyzed the overall survival of EC patients (paying attention to several clinicopathological characteristics such as age, sex, pathologic stage, and neoplasm type; Table 1). Univariate and multivariate analyses showed that only the expression of MTHFD2 expression was significantly associated with pathologic stage (HR: 1.971, 95% CI 1.032–3.764, P = .037; HR: 1.885, 95% CI 0.965–3.681, P = .043), indicating that MTHFD2 is an independent prognostic factor in EC progression.

Table 1

Univariate and multivariate analyses of clinicopathological characteristics and important genes with overall survival in TCGA esophageal cancer cohort.

MTHFD2 expression was further analyzed in tumors and adjacent normal tissues by data validation of the TCGA and GSE17351 datasets (Fig. 6B). MTHFD2 was significantly upregulated in esophageal cancer samples compared with normal tissues (GSE17351 dataset, n = 10, P = .008; TCGA dataset, n = 173, P < .001), which was in line with our prediction results. Thus, MTHFD2 can serve as an independent EC biomarker associated with poor prognosis.

Discussion

In this present study, we identified 232 genes with aberrant expression and methylation status based on the analysis of transcriptome and methylation datasets. Functional enrichment analysis of these genes revealed enrichment of pathways such as “ECM-receptor interaction,” “focal adhesion,” “p53 signaling pathway,” “cell cycle,” and “glutathione metabolism.” Eight genes were identified as hub genes in the PPI network such as MTHFD2, MCM6, and CENPE. According to survival analysis, deregulation of 3 genes (MSN, MTHFD2, and PELI1) was correlated with overall survival time in EC. Cox regression analysis showed that only MTHFD2 expression status was significantly correlated with pathologic stage (P < .05). Finally, data validation using the TCGA database demonstrated that MTHFD2 was significantly upregulated in tumor tissue compared with adjacent normal samples, indicating that MTHFD2 might have prognostic value in EC development.

MTHFD2 encodes a mitochondrial enzyme with methylenetetrahydrofolate dehydrogenase and cyclohydrolase activities.[ This protein is expressed in the developing embryo and is absent in most adult tissues. Interestingly, recent studies have focused on the functions of MTHFD2 in cancer biology. Overexpression of the MTHFD2 protein was detected in 19 cancer types; high MTHFD2 mRNA expression was correlated with poor prognosis in breast cancer patients.[ Repression of MTHFD2 decreased the invasion and migration of breast cancer cell lines.[ However, little is known about the regulatory mechanism(s) of MTHFD2 in EC. Based on metabolic enzyme expression analysis, MTHFD2 has a key role in the mitochondrial one-carbon folate pathway in various cancers.[ Folate metabolism is central to cell proliferation, and MTHFD2 was reported as a folate-coupled enzyme broadly required for cancer cell proliferation. It was shown to be responsible for mitochondrial NADPH production in proliferating cancer cells and involved in biochemical reactions including deoxythymidylate, purine nucleotide production, and amino acid inter conversion.[ Thus, the effects of MTHFD2 on EC progression exerted through the regulation of the 1 carbon metabolism and folate pathways require further validation.

Furthermore, MTHFD2 may be regulated by several microRNAs. In colorectal cancer, MTHFD2 expression promotes cancer cell growth, and MTHFD2 is targeted by miR-33a-5p, a microRNA with a major role in tumorigenesis.[ In addition, miR-92a inhibits cancer cell proliferation and induces apoptosis by regulating MTHFD2 in acute myeloid leukemia.[ Similarly, miR-9 can directly target MTHFD2 to exert anti-proliferative and pro-apoptotic effects in breast cancer cells.[ A recent study demonstrated that miRNA-940 interacts with MTHFD2 to induce mitochondrial folate metabolism dysfunction and suppression of glioma progression.[ Thus, the mitochondrial isozyme MTHFD2 is upregulated and associated with cancer prognosis. Targeting MTHFD2 and the mitochondrial folate pathway might be a potential novel therapeutic strategy for cancer.

As for other genes, MSN encodes the moesin protein belonging to the ERM family (ezrin, radixin, and moesin).[ Extensive roles for moesin in tumors have been uncovered, and moesin overexpression is correlated with metastasis and poor prognosis in different cancer types.[ Wu et al[ revealed that ECM1 interacts with MSN and facilitates invadopodia formation, which is required for stromal invasion and metastasis in breast cancer. Similarly, in hepatocellular carcinoma cells, Lan et al[ identified that moesin promotes metastasis by improving invadopodia formation and activating the β-catenin/MMP9 axis. However, the potential role of MSN in EC remained unknown. Our results now suggest that MSN expression is correlated with poor prognosis in EC patients, indicating MSN plays a vital role in EC progression. In addition, PELI1 (pellino-1) functions as an E3 ubiquitin ligase that is degraded by the proteasome. It is involved in inflammatory and autoimmune diseases through mediating post-translational modifications in animal cells.[ A previous study on B-cell lymphomas found that PELI1 promotes lymphomagenesis by regulating BCL6 polyubiquitination.[ PELI1 is upregulated in some high-grade B-cell lymphomas, and relatively downregulated in low grade B-cell lymphomas or T-cell lymphomas.[ Abnormal expression of PELI1 is correlated with MYC and BCL6, and is associated with poor prognosis in B-cell lymphoma patients. However, the potential function(s) of PELI1 in solid tumors remains largely unknown. Overexpression of PELI1 increases cell proliferation and survival, and contributes to lung tumorigenesis by promoting the epithelial to mesenchymal transition.[ Furthermore, PELI1 promotes cell survival and chemoresistance by upregulating the expression of cIAP2 in lung cancer.[ Taken together, these findings suggest that PELI1 might be a potential therapeutic target in some tumor types.

Conclusion

In conclusion, our study identified 8 hub genes with aberrant expression and methylation levels in EC patients through integrated analysis of transcriptome and methylation data. MSN, MTHFD2, and PELI1 were identified as potential tumor biomarkers with prognostic value in EC. Expression of MTHFD2 was significantly correlated with pathological stage and poor prognosis, highlighting its important role in tumor development, and ability to predict EC prognosis.

Author contributions

Data curation: Jianlin Wang, Judong Luo.

Formal analysis: Zhiqiang Sun, Fei Sun.

Funding acquisition: Jianlin Wang, Ze Kong, Jingping Yu.

Resources: Fei Sun, Ze Kong.

Software: Jianlin Wang, Judong Luo, Zhiqiang Sun.

Writing – original draft: Jianlin Wang.

Writing – review & editing: Jianlin Wang, Jingping Yu.