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%.[1,2] 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.[3–5] 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.[6,7] 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.
2.1 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.
2.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.
2.3 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.
2.4 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.
3.1 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).
3.2 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.
3.3 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.
3.4 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).
3.5 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.
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.
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.
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.[25,26] Repression of MTHFD2 decreased the invasion and migration of breast cancer cell lines.[27,28] 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.[29–31] 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).[36,37] Extensive roles for moesin in tumors have been uncovered, and moesin overexpression is correlated with metastasis and poor prognosis in different cancer types.[38–41] 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.[44–46] 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.
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.
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.
. Chevrollier GS, Giugliano DN, Palazzo F, et al. Patients with non-response to neoadjuvant chemoradiation for esophageal cancer have no survival advantage over patients undergoing primary esophagectomy. J Gastrointest Surg 2020;24:28898.
. National Cancer Institute - Surveillance, Epidemiology, and End Results Program. Based on data from SEER 2009-2015. Available at: https://seer.cancer.gov/statfacts/html/esoph.html
. Accessed April 01, 2020.
. Esteller M. Epigenetics in cancer. N Engl J Med 2010;358:114859.
. Zhang H, Dong S, Feng J. Epigenetic profiling and mRNA expression reveal candidate genes as biomarkers for colorectal cancer. J Cell Biochem 2019;120:1076776.
. Wang XX, Xiao FH, Li QG, et al. Large-scale DNA methylation
expression analysis across 12 solid cancers reveals hypermethylation in the calcium-signaling pathway. Oncotarget 2017;8:1186876.
. Jones PA, Baylin SB. The fundamental role of epigenetic events in cancer. Nat Rev Genet 2002;3:41528.
. Jammula S, Katz-Summercorn AC, Li X, et al. Identification of subtypes of Barrett's Esophagus and esophageal adenocarcinoma based on DNA methylation
profiles and integration of transcriptome and Genome Data. Gastroenterology 2020;158:1682.e197.e1.
. Liu Y, Zhang M, He T, et al. Epigenetic silencing of IGFBPL1 promotes esophageal cancer growth by activating PI3K-AKT signaling. Clin Epigenetics 2020;12:22.
. Ma Y, He S, Gao A, et al. Methylation silencing of TGF-β receptor type II is involved in malignant transformation of esophageal squamous cell carcinoma. Clin Epigenetics 2020;12:25.
. Chen Y, Liao LD, Wu ZY, et al. Identification of key genes by integrating DNA methylation
and next-generation transcriptome sequencing for esophageal squamous cell carcinoma. Aging (Albany NY) 2020;12:133265.
. Huang AM, Kao YT, Toh S, et al. UBE2M-mediated p27(Kip1) degradation in gemcitabine cytotoxicity. Biochem Pharmacol 2011;82:3542.
. Couto-Vieira J, Nicolau-Neto P, Costa EP, et al. Multi-cancer V-ATPase molecular signatures: a distinctive balance of subunit C isoforms in esophageal carcinoma
. EBioMedicine 2020;51:102581.
. Lee JJ, Natsuizaka M, Ohashi S, et al. Hypoxia activates the cyclooxygenase-2-prostaglandin E synthase axis. Carcinogenesis 2010;31:42734.
. Li X, Zhou F, Jiang C, et al. Identification of a DNA methylome profile of esophageal squamous cell carcinoma and potential plasma epigenetic biomarkers for early diagnosis. PLoS One 2014;9:e103162.
. Diboun I, Wernisch L, Orengo CA, Koltzenburg M. Microarray analysis after RNA amplification can detect pronounced differences in gene expression using limma. BMC Genomics 2006;7:252.
. Aryee MJ, Jaffe AE, Hector CB, et al. Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation
microarrays. Bioinformatics 2014;30:13639.
. Bar-Joseph Z, Gifford DK, Jaakkola TS. Fast optimal leaf ordering for hierarchical clustering. Bioinformatics 2001;17:S229.
. Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 2009;4:4457.
. Ogata H, Goto S, Sato K, et al. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res 1999;27:2934.
. Szklarczyk D, Franceschini A, Wyder S, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res 2015;43:D44752.
. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 2003;13:2498504.
. Liang B, Li C, Zhao J. Identification of key pathways and genes in colorectal cancer using bioinformatics analysis. Med Oncol 2016;33:016829.
. Lin DY. Cox regression analysis of multivariate failure time data: the marginal approach. Stat Med 1994;13:223347.
. Xu X, Qiao M, Zhang Y, et al. Quantitative proteomics study of breast cancer cell lines isolated from a single patient: discovery of TIMM17A as a marker for breast cancer. Proteomics 2010;10:137490.
. Nilsson R, Jain M, Madhusudhan N, et al. Metabolic enzyme expression highlights a key role for MTHFD2 and the mitochondrial folate pathway in cancer. Nat Commun 2014;5:3128.
. Liu F, Liu Y, He C, et al. Increased MTHFD2 expression is associated with poor prognosis
in breast cancer. Tumor Biol 2014;35:868590.
. Selcuklu SD, Donoghue MT, Rehmet K, et al. MicroRNA-9 inhibition of cell proliferation and identification of novel miR-9 targets by transcriptome profiling in breast cancer cells. J Biol Chem 2012;287:2951628.
. Lehtinen L, Ketola K, Makela R, et al. High-throughput RNAi screening for novel modulators of vimentin expression identifies MTHFD2 as a regulator of breast cancer cell migration and invasion. Oncotarget 2013;4:4863.
. Mattaini KR, Sullivan MR, Vander Heiden MG. The importance of serine metabolism in cancer. J Cell Biol 2016;214:24957.
. Shin M, Momb J, Appling DR. Human mitochondrial MTHFD2 is a dual redox cofactor-specific methylenetetrahydrofolate dehydrogenase/methenyltetrahydrofolate cyclohydrolase. Cancer Metab 2017;5:0173.
. Wei Y, Liu P, Li Q, et al. The effect of MTHFD2 on the proliferation and migration of colorectal cancer cell lines. Onco Targets Ther 2019;12:636170.
. Yan Y, Zhang D, Lei YY, et al. MicroRNA-33a-5p suppresses colorectal cancer cell growth by inhibiting MTHFD2. Clin Exp Pharmacol Physiol 2019;46:192836.
. Rollin R, Alvarez-Lafuente R, Marco F, et al. The ubiquitin-proteasome pathway and viral infections in articular cartilage of patients with osteoarthritis. Rheumatol Int 2009;29:96972.
. Selcuklu SD, Donoghue MTA, Rehmet K, et al. MicroRNA-9 inhibition of cell proliferation and identification of novel mir-9 targets by transcriptome profiling in breast cancer cells. J Biol Chem 2012;287:2951628.
. Xu T, Zhang K, Shi J, et al. MicroRNA-940 inhibits glioma progression by blocking mitochondrial folate metabolism through targeting of MTHFD2. Am J Cancer Res 2019;9:25069.
. Jiang L, Phang JM, Yu J, et al. CLIC proteins, ezrin, radixin, moesin and the coupling of membranes to the actin cytoskeleton: a smoking gun? Biophys Acta 2014;1838:64357.
. Deng W, Cho S, Li R. FERM domain of moesin desorbs the basic-rich cytoplasmic domain of l-selectin from the anionic membrane surface. J Mol Biol 2013;425:354962.
. Charafe-Jauffret E, Monville F, Bertucci F, et al. Moesin expression is a marker of basal breast carcinomas. Int J Cancer 2007;121:177985.
. Schlecht NF, Brandwein-Gensler M, Smith RV, et al. Cytoplasmic ezrin and moesin correlate with poor survival in head and neck squamous cell carcinoma. Head Neck Pathol 2012;6:23243.
. Wang X, Liu M, Zhao C. Expression of ezrin and moesin related to invasion, metastasis and prognosis
of laryngeal squamous cell carcinoma. Genet Mol Res 2014;13:800213.
. Barros FBA, Assao A, Garcia NG, et al. Moesin expression by tumor cells is an unfavorable prognostic biomarker
for oral cancer. BMC Cancer 2018;18:53.
. Wu Q, Chen D, Luo Q, et al. Extracellular matrix protein 1 recruits moesin to facilitate invadopodia formation and breast cancer metastasis. Cancer Lett 2018;437:4455.
. Lan S, Zheng X, Hu P, et al. Moesin facilitates metastasis of hepatocellular carcinoma cells by improving invadopodia formation and activating β-catenin/MMP9 axis. Biochem Biophys Res Commun 2020;524:8618.
. Jin W, Chang M, Sun S-C. Peli: a family of signal-responsive E3 ubiquitin ligases mediating TLR signaling and T-cell tolerance. Cell Mol Immunol 2012;9:11322.
. Humphries F, Moynagh PN. Molecular and physiological roles of Pellino E3 ubiquitin ligases in immunity. Immunol Rev 2015;266:93108.
. Medvedev AE, Murphy M, Zhou H, et al. E3 ubiquitin ligases Pellinos as regulators of pattern recognition receptor signaling and immune responses. Immunol Rev 2015;266:10922.
. Park H-Y, Go H, Song HR, et al. Pellino 1 promotes lymphomagenesis by deregulating BCL6 polyubiquitination. J Clin Investig 2014;124:497688.
. Choe J-Y, Park M, Yun JY, et al. PELI1 expression is correlated with MYC and BCL6 expression and associated with poor prognosis
in diffuse large B-cell lymphoma. Mod Pathol 2016;29:131323.
. Jeon YK, Kim CK, Hwang KR, et al. Pellino-1 promotes lung carcinogenesis via the stabilization of Slug and Snail through K63-mediated polyubiquitination. Cell Death Differ 2017;24:46980.
. Jeon YK, Kim CK, Koh J, et al. Pellino-1 confers chemoresistance in lung cancer cells by upregulating cIAP2 through Lys63-mediated polyubiquitination. Oncotarget 2016;7:4181124.