An integrated analysis of DNA promoter methylation, microRNA regulation, and gene expression in gastric adenocarcinoma.

BACKGROUND: Gastric adenocarcinoma (GAC), a common type of gastric cancer, poses a significant public health threat worldwide. This study aimed to determine the transcriptional regulatory mechanisms of GAC.
METHODS: HTSeq-FPKM raw data were obtained from The Cancer Genome Atlas Stomach Adenocarcinoma data collection. Subsequently, the limma package in R was used to identify differentially expressed genes (DEGs). Differentially methylated genes (DMGs), DEGs, and differentially expressed microRNAs (miRNAs) in normal, and tumor tissues of the same patients were screened and compared using R software tools. A functional enrichment analysis was performed using Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) for various DEGs, DMGs, promoter methylation, and miRNAs. DEG-specific methylation and transcription factors were analyzed using ENCODE ChIP-seq.
RESULTS: DEGs were centrally modified by the histone trimethylation of lysine 27 on histone H3 (H3K27me3). Upstream transcription factors of DEGs were enriched in different ChIP-seq clusters, such as Forkhead Box M1, E2F Transcription Factor 4, and suppressor of zest 12. Integrated regulatory networks of DEGs, promoter methylation, and miRNAs were constructed. Two miRNAs (hsa-mir-1 and hsa-mir-133a) and four DEGs (A disintegrin and metalloproteinase domain 12, transcription factor AP-2 alpha, solute carrier family 5 member 7, and cadherin 19) separately played important roles in the integrated regulatory network. Therefore, these DEGs, DMGs, promoter methylation, and miRNAs may play an important role in GAC pathogenesis.
CONCLUSIONS: In summary, the present study results provide insights into the oncogenesis and progression of GAC, thus accelerating the development of novel targeted GAC therapies. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Gastric cancer (GC) is the fifth most common malignancy worldwide, with the third-highest incidence and mortality rates of all types of cancer (1). In 2018, more than one million new GC-related cases were diagnosed with an estimated 783,000 global deaths (1). Moreover, incidence rates are markedly elevated in East Asian countries, particularly in Mongolia, Japan, and the Republic of Korea, whereas those in North America, Northern Europe, and African regions are generally low (1). Surgery is considered a potentially curative therapy for early GC; however, the prognosis for patients with advanced adenocarcinoma remains poor, despite improvements in chemotherapy (2).

Epigenetic alteration is an important event during carcinogenesis; it plays a critical role in the transcriptional silencing of tumor suppressor genes (3). Detecting the regulators of gene expression associated with cancer progression is one of the fundamental challenges in cancer research. DNA methylation occurs almost exclusively in the context of CpG islands,defined as regions of more than 200 bases with a G+C content of at least 50% and a ratio of statistically expected CpG frequencies of at least 0.6 (4). Moreover, the methylation of these CpG island shores referred to as regions of lower CpG density that lie close (~2 kb) to CpG islands, is likely to occur in promoter regions near the transcription start site. Additionally, microRNAs (miRNAs) are a class of endogenous small RNAs containing 20–25 nucleotides that regulate mRNA translation and degradation by perfect, or near-perfect, complementarity to the 3' untranslated regions (UTRs) of the target mRNA (5). miRNAs act as regulators of post-transcriptional levels and transcriptome changes and are also major components of the epigenome. Active regulatory functions of nuclear miRNAs have been described in gene promoter regions and exert further effects via epigenetic pathways (6). Hypermethylation and hypomethylation of miRNAs frequently occur in human cancer, representing a new level of complexity in gene regulation. Furthermore, in this research, we found that miRNAs and DNA methylation were in mutual regulation in GAC. Specifically, miRNAs regulate DNA methylation by targeting DNA methyltransferases or methylation-related proteins (7). Nevertheless, the mechanism of the synergistic interactions between DNA methylation and miRNAs as epigenetic regulators of transcriptomic changes, as well as their association with clinical outcomes, have remained largely unexplored in cancer research (8).

Previous studies have identified various genes related to gastric adenocarcinoma (GAC). For instance, Pilehchian Langroudi et al. (9) analyzed FAT4 expression and promoter methylation in GC and reported that FAT4 is a tumor suppressor gene in cell adhesion. This indicated that a reduced expression of FAT4 and increased methylation of its promoter might be a key mechanism associated with tumor growth (9). Furthermore, targeted therapies, such as trastuzumab, an antibody against HER2 (also known as ERBB2), and the VEGFR-2 antibody ramucirumab, have recently been initiated (10).

GAC is a major subtype of GC, but the pathogenesis of GAC has yet to be characterized. To reduce the significant morbidity and mortality associated with GAC, it is critical to identify GAC-associated genes and mechanisms. Tumorigenesis is regulated by epigenetic phenomena, including nucleosome remodeling by histone modifications, DNA methylation, and miRNA-mediated targeting of genes associated with various biochemical pathways. Epigenetic abnormalities are heritable as part of gene transcription and collaborate with genetic changes to cause cancer evolution. General hypomethylation and focal hypermethylation of the noncoding region, especially the promoter-related CpG island associated with gene silencing, are the characteristics of cancer cells. This process may be related to the acquisition of histone inhibitory markers (11). Furthermore, mutated genes controlling epigenetic factors have further strengthened the importance of epigenetics in cancer (12). It is worth mentioning that mutations in genes caused by epigenetic changes can be frequently observed (13). Most cancers harbor frequent mutations in genes that encode components of the epigenetic machinery, resulting in abnormalities in the epigenome, which can further affect gene expression patterns and genomic stability (14,15). Hence, conducting an integrated analysis of differentially expressed genes (DEGs) and the expression of regulatory factors, such as methylation, mRNA splicing, transcription factors (TFs), and miRNAs, is an effective strategy for investigating GAC pathogenesis. The present study used microarray technology and bioinformatics analyses to investigate the regulatory network of GAC-associated miRNAs and their target genes in an effort to further explore the function of promoter methylation and miRNAs in the pathogenesis of GAC. We hope these findings will serve as a foundation for future studies. We present the following article in accordance with the MDAR reporting checklist (available at https://dx.doi.org/10.21037/atm-21-3211).

Methods

Datasets

We downloaded the GAC dataset from The Cancer Genome Atlas (TCGA) data portal (https://gdc.cancer.-gov/), which included the DNA methylation dataset, the mRNA expression profile dataset, and the miRNA expression profile dataset. In total, we obtained 397 DNA tissue samples (395 GAC, 2 normal), 407 mRNA tissue samples (375 GAC, 32 normal), and 491 miRNA tissue samples (446 GAC, 45 normal). The three datasets were then combined using R software (3.5.2 version, http://cran.r-project.org/src/base/R-3/R-3.5.2.tar.gz). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Differential expression analysis of mRNA and miRNA

We obtained the differential expressions of mRNAs and miRNAs by using the limma package (16) in R [logFC >2, false discovery rate (FDR) <0.01]. The target genes of miRNAs were acquired from mirWalk (http://zmf.umm.uniheidelberg.de/apps/zmf/mirwalk/micrornapredictedtarget.html).

Screening of DEGs and differentially methylated genes (DMGs)

We used the limma package in R (16,17) to compare DEGs in GAC and normal tissue samples after adjusting the FDR <0.01 and |log2fold-change (FC)| ≥1 threshold. Differentially methylated regions (DMRs) were discovered using the minfi R software (17), with a q value <0.05. By applying gplots in the limma package (https://cran.r-project.org/web/packages/gplots/), we produced volcano plots to analyze the DEGs with adjusted P values <0.01 and |log2FC| ≥2.

Functional analysis of DEGs

Gene Ontology (GO), consisting of three categories (biological process, BP; cellular component, CC; and molecular function, MF), was used to predict the potential functions of gene products. The GO functional enrichment and annotation of DEGs were computed using Enrichr (18) (https://amp.pharm.mssm.edu/Enrichr/). We used the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, which connects genomic information with functional coding information to systematically analyze gene functions (19). We also used Cytoscape (http://apps.cytoscape.org/apps/keggscape) to reproduce equivalent detailed hand-drawn pathway diagrams (20). The GO and KEGG enrichment analyses were performed for the DEGs, DMGs, and the target genes of miRNAs, as well as promoter methylation. The P values of the enriched terms were corrected by the Holm Bonferroni method (21). An adjusted P value <0.05 was applied as the cut-off criterion. ChIP-X enrichment analysis (ChEA) (22) and ENCODE ChIP-seq (23) were used to screen the enriched TFs located upstream of the DEGs. Lists of the mammalian gene symbols, for which the program computes over-represented TFs from the ChIP-X database, were then imported using ChEA.

Correlation analysis of DEGs, differential miRNAs, and promoter methylation

We used a network graph to determine the correlations between GAC DEGs, differential miRNAs, and promoter methylation. Funrich (24) (http://www.mybiosoftware.com/funrich-functional-enrichment-analysis-tool.html), a stand-alone software tool, was used to produce Venn diagrams to show the functional enrichment and interaction network analysis of the identified genes (25). The data visualization tool MEXPRESS (https://mexpress.be/) was used to analyze further the precise methylation locus and expression (26).

Construction of the risk score-based prognostic model

The GAC tumor samples from the TCGA-STAD cohort (n=350) were enrolled as the training cohort to construct the prognostic model, and the trimmed mean of M-values (TMM) normalized expression matrix of 17 interested genes were obtained from . Ten-fold cross-validation were applied to the multivariate cox regression model to train the minimum lambda value, with 17 genes as the parameters. The minimum lambda value (λmin≈0.046) was applied to the original multivariate cox model, returning a subset of variables containing 5 variables (ADAM12, CDH19, RIMS1, SHOX2, and IGF2BP1). Furthermore, the risk score model was built on the Lasso-Cox regression result, with 5 genes selected based on the subset with λmin applied. The coefficients of the risk score model were weighted by the coefficients of the lasso multivariate Cox regression model. The formula for calculating the risk score is: risk score = 0.058316123 × EXP + 0.014146433 × EXP + 0.015680240 × EXP + 0.032075667 × EXP + 0.007903263 × EXP.

Figure 1

Interrelated networks in the context of DEGs, differentially expressed miRNAs, and gene expression. Methylation and microRNA regulatory network of 31 DEGs. The downregulated and upregulated genes are shown in green and red circles, respectively. The upregulated and downregulated miRNAs are displayed in red and green rectangles, respectively. DEGs, differentially expressed genes

Validation of the prognostic model in the training and validation cohort

In the training cohort, time-dependent receiver operating characteristic (ROC) analysis was adopted to evaluate the performance of the risk score-based prognostic signature to predict 1-, 3- and 5-year survival. Kaplan-Meier (KM) survival curves combined with a log-rank test were used to test the differences in prognosis between the high- and low-risk groups. Gastric tumor samples from the GSE15459 dataset (n=192) were extracted as the validation cohort. Similarly, time-dependent ROC analysis and KM survival analysis were also performed on the validation cohort.

Construction and assessment of a predictive nomogram

Univariate and multivariate analyses were successively performed on the 5 genes-based risk score prognostic model alone and combined with other clinical characteristics (including TNM stage, age, and sex) of patients with GAC to test their predictive efficacy. In the training cohort, univariate Cox regression analysis was used to identify clinical characteristics that were significantly associated with overall survival (OS). Characteristics with P<0.05 in the univariate analysis were included in the multivariate Cox regression model. Subsequently, nomogram for 1-, 3-, and 5-year OS were constructed using the rms R package (27) as a visualizing aid to obtain predicted values manually from the multivariate analysis. The performance of the nomogram was evaluated using Harrell’s concordance index (c-index) and calibration curves.

Overall survival analysis

To identify which interactions between miRNAs and methylation are significantly associated with GAC prognosis, we used the Kaplan-Meier (KM) overall survival analysis from UALCAN (http://ualcan.path.uab.edu/index.html) (28). This program aids in investigating gene expression patterns related to clinical parameters. A set of KM plots (29) from both the univariate analysis (considering only gene expression) and multivariate analysis (considering clinical parameters along with gene expression) were provided from the output page of the survival analysis after using the “scan by gene” option in UALCAN. A KM survival plot was generated for every gene in each TCGA cancer type, using the “survival” (https://cran.r-project.org/web/packages/survival/index.html) and “survminer” (http://mirrors.dotsrc.org/pub/cran/web/packages/survminer/index.html) packages. The survival curves of samples with high gene expressions and low/medium gene expressions were compared using the log-rank test (30). Based on the KM plots, we then ran survival analyses for the two extreme subgroups; the first subgroup comprised the lowest quantiles of miRNAs and methylation profiles, and the second subgroup comprised the highest quantiles of both profiles. Observation of the KM plots allowed us to assess the relationship between gene expression and overall survival in GAC patients.

Statistical analysis

Statistical analyses were performed using R software (https://mirrors.tuna.tsinghua.edu.cn/CRAN/). Kaplan-Meier survival analysis was conducted together with the log-rank test to compare the prognostic risk of GAC patients with different average expression levels. Unless otherwise noted, a P value <0.05 was considered statistical significance.

Results

Identification of DEGs, target genes of miRNAs, and promotor methylation

The overall study design is outlined in a flow chart in . A total of 397 DNA samples, 407 mRNA samples, and 491 miRNA samples from TCGA were analyzed using the minifi package in R with a q value <0.05 to identify DMRs; in turn, 8,615 DMGs were identified using annover. Finally, 2,059 differential promoter methylation genes were identified. Using the limma package in R, 1,160 downregulated target genes and 3,480 upregulated target genes of miRNAs were identified. Moreover, 507 downregulated genes and 919 upregulated genes were identified. The DEGs were selected by volcano plot filtering (P value <0.05 and |log2FC| ≥2) ().

Figure 2

Flow chart representing the detailed analysis process applied in this work. TCGA-STAD, The Cancer Genome Atlas Stomach Adenocarcinoma; DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology.

Figure 3

Volcano plot depicting the differentially expressed genes between normal gastric and GAC samples. The red dots indicate upregulated genes, the green dots indicate downregulated genes, and the white block represents the remaining genes—|log2FC| <2. The Y-axis represents FDR, and the X-axis represents the value of log2FC. DEGs, differentially expressed genes; GAC, gastric adenocarcinoma.

GO enrichment analysis of methylation-related genes

To further investigate the function of methylation-related genes in GAC, we used the GO enrichment analysis in Enrichr, and the KEGG pathway enrichment analysis in Cytoscape. Four specific variables were used in the analyses: DEGs, DMGs, the target genes of miRNAs, and promoter methylation. The top seven GO pathways are shown in . Specifically, DEGs were enriched in the extracellular matrix organization (ECM) (BP), the spindle of cell components (CC), and endopeptidase activity (MF) (). For DMGs, the predominantly enriched genes included the development of the nervous system (BP), exon (CC), and transcriptional activator and factor activity (MF) (). The target genes of miRNAs were enriched in the regulation of transcription from RNA polymerase II promoter (BP), nuclear body (CC), and protein kinase activity (MF) (). Promoter methylation genes showed enrichment in the regulation of transcription and DNA-template (BP), an integral component of the plasma membrane (CC), and transcription regulatory region DNA binding (MF) ().

Figure 4

GO enrichment analysis of methylation-related genes. GO enrichment results for the (A) DEGs, (B) DMGs, (C) miRNA target genes, (D) and promotor methylation. The data are represented as the mean ± SEM. The Student’s t-test was used to analyze significant differences, P<0.05 vs. shCtrl. shCtrl, negative control cells; GO, Gene Ontology; DEGs, differentially expressed genes; DMGs, differentially methylated genes; BP, biological process; CC, cellular component; MF, molecular function.

KEGG pathway enrichment analysis of methylation-related genes

Based on the KEGG pathway enrichment analysis, DEGs were significantly enriched in the IL-17 signaling pathway, cytokine-cytokine receptor interaction, and cell cycle (Figure S1), while other pathways exhibited little correlation with GAC. DMGs were associated with the neuroactive ligand-receptor interaction, cAMP signaling pathway, RAP1 signaling pathway, exon guidance, cell adhesion molecules, ribosome biogenesis in eukaryotes, and ribosomes (Figure S2). The target genes of miRNAs were correlated with endocytosis, signaling pathways regulating pluripotency of stem cells, proteoglycans in cancer, the RAP1 signaling pathway, olfactory transduction, axon guidance, the mitogen-activated protein kinase (MAPK) signaling pathway, pathways in cancer, and RAS signaling pathways (Figure S3). Promotor methylation genes were strongly related to neuroactive ligand-receptor interactions and cell adhesion molecules (Figure S4).

Transcriptional analysis of DEGs

ChEA2 analyses obtained from the ENCODE database () indicated that the screened DEGs were modified and regulated by multi-cancer cell line histones, including trimethylation of lysine 27 on histone H3 (H3K27me3). The binding patterns of the upstream TFs were not as clustered as those of the histone modification; however, they were enriched in different ChIP-seq TF clusters in different cell lines (), such as forkhead box M1 (FOXM1) in MCF7 cells and endometrial cells (ECC1), E2F transcription factor 4 (E2F4) in HeLaS3 cells and murine erythroleukemia (MEL) cells, and the suppressor of zest 12 (SUZ12) in NTERA2 clone D1 (NT2D1).

Figure 5

Transcriptional analysis of DEGs. Enriched TFs are located upstream of the DEGs, as identified by ChIP-X enrichment analysis (A) and ENCODE ChIP-seq of histone modifications and TFs (B). TFs, transcription factors; DEGs, differentially expressed genes.

Integrated analysis of DEGs, differential miRNAs, and promoter methylation

In the Venn diagram shown in , a total of 31 genes appeared in the overlapped area between differentially expressed miRNA target genes, DEGs, and promoter methylated genes. In addition, based on the hypermethylation and hypomethylation levels, potential regulatory networks associated with methylation, miRNA expression, and gene expression were generated (). Within this network were 17 upregulated and 3 downregulated miRNAs and 13 upregulated and 4 downregulated genes. In addition, 6 hypermethylated genes and 11 hypomethylated genes were included. In total, 20 miRNAs were found in the network from which miR-133a (degree =4) had the highest degree. In addition, A disintegrin and metalloproteinase domain 12 (ADAM12) (31), orthodenticle homeobox 1 (OTX1), and transcription factor AP-2 alpha (TFAP2A) (32) had relatively high degrees. Moreover, ADAM12 was associated with upregulation of hsa-mir-501, hsa-mir-19a, and hsa-mir-130b, and downregulation of hsa-mir-1 and hsa-mir-133a (). TFAP2P correlated with upregulation of hsa-mir-429, hsa-mir-200b, hsa-mir-200c, and hsa-mir-135b, as well as the downregulation of hsa-mir-133a.

Figure 6

Correlation analysis of the DEGs, differentially expressed miRNAs, and promotor methylation. The Venn diagram representing the functional enrichment and interaction network formed in the context of the DEGs, miRNA target genes, and promoter methylation are also shown. DEGs, differentially expressed genes.

Establishment of the risk score and prognostic model of the nomogram

The method section has presented the formula to calculate the risk score based on 5 genes (ADAM12, CDH19, RIMS1, SHOX2, and IGF2BP1), the method section has presented the formula to calculate the risk score. Based on the median f the risk score, the 350 GAC patients in the training cohort were divided into a high-risk group (n=175) and a low-risk group (n=175). Time-dependent ROC curves () were drawn to evaluate the predictive efficacy of the constructed risk score-based prognostic model, with AUCs for the prediction of 1 year, 3 years, and 5 years of OS were 0.607, 0.632, and 0.735, respectively. Kaplan-Meier survival analysis () indicated that patients in the high-risk group had worse prognoses than the low-risk group (P=0.001). The GSE15459 dataset (n=192) was used as the validation cohort to evaluate the performance of the risk score-based model; it comprised 96 patients in the high-risk group and 96 patients in the low-risk group, respectively. The time-dependent ROC curve and Kaplan-Meier survival curve were also plotted. As shown in , the AUC for the prediction of OS at 1 year, 3 years, 5 years in the validation cohort were 0.568, 0.587 and 0.615, respectively; worse prognoses were revealed in patients in the high-risk group than in the low-risk group (P=0.0089).

Figure 7

Exploration of the risk score and prognostic model of the nomogram in the training and validation cohort. (A) ROC curve analysis of risk score-based prognostic model in the training cohort. (B) Kaplan-Meier survival curve of the high- and low-risk subgroups in the training cohort. (C) ROC curve analysis of risk score-based prognostic model in the validation cohort. (D) Kaplan-Meier survival curve of the high- and low-risk subgroups in the validation cohort. (E) Nomogram for predicting the OS rates at 1-, 3- and 5-year in the training cohort. (F) Calibration curve of the nomogram for predicting 1-, 3-, and 5-year OS in the training cohort. AUC, area under the curve; OS, overall survival; ROC, receiver operating characteristic.

To construct a prognostic model of nomogram, the univariate and multivariate analyses of OS regarding canonical clinicopathologic characteristics were performed in the training cohort (Table S1). In the univariate analysis, OS was significantly correlated with risk score (P=0.001), TNM-T stage (P=0.018), TNM-N stage (P=0.005), TNM-M stage (P=0.014), and age at diagnosis (P=0.023). However, in the multivariate analysis, only the risk score (P=0.007), TNM-N stage (P=0.017), TNM-M stage (P=0.017), and age at diagnosis (P=0.012) remained significantly associated with OS. As shown in , a nomogram was further constructed to predict 1-, 3-, and 5-year OS in the training cohort (TCGA-STAD dataset) to provide an accurate prediction of the prognosis. The estimated total points were calculated from clinical factors significantly associated with OS in the univariate analysis. For example, a high-risk score patient (67.5 points) aged older than 60 years old (65 points) with pathologic T3/T4 (27.5 points), N1/2/3 (75 points), M1 stage (100 points) would score 335 points in total for OS, and his/her predicted 1-, 3-, 5-year OS would be approximately 50%, 12.5%, and less than 10%respectively. Additionally, the prognostic nomogram model was internally validated by assessing both discriminations using the c index and plotting calibration curves (). The c-index of the combined nomogram (0.640, 95% CI: 0.615–0.666) was higher than that of other clinicopathological characteristics alone.

Analysis of the relationships between key gene methylation and expression

To further analyze the precise methylation loci and expression, MEXPRESS was used to investigate the relationships between ADAM12 and TFAP2A methylation and expression (Figures S5,S6). The dark blue regions in Figure S5 demonstrate a negative correlation with ADAM12 gene expression (Pearson’s correlation coefficients for each probe are shown on the right), indicating that ADAM12 methylation silences gene expression. The dark blue regions in Figure S6 demonstrate a positive correlation with TFAP2A gene expression, indicating that TFAP2A methylation activates gene expression.

Survival outcomes

Survival outcomes were analyzed using UALCAN to investigate whether the gene expression levels of the hypermethylation key drivers ADAM12, SLC5A7, and CDH19, and the hypomethylation driver TFAP2A, play important roles in GAC patient prognosis. The log-rank P values of ADAM12, SLC5A7, CDH19, and TFAP2A were 0.0044, 0.013, 0.036, and 0.0066, respectively. Therefore, regulation of oncogenes caused by aberrant methylation could result in different overall survival. Specifically, patients with a high TFAP2A expression () had better prognoses, while those with a low/medium expression of ADAM12 (), CDH19 (), and SLC5A7 () had poorer prognoses.

Figure 8

Survival analyzes to investigate the impact of the key methylation driver genes for the prognoses of GAC patients. Kaplan-Meier survival curves demonstrating the prognostic value of (A) TFAP2A, (B) ADAM12, (C) CDH19, (D), and SLC5A7. A two-sided P value of <0.05 was regarded as statistically significant; of note, the P value was FDR-corrected. P<0.05 vs. shCtrl. shCtrl, negative control cells. ADAM12, A disintegrin and metalloproteinase domain 12; TFAP2A, transcription factor AP-2 alpha; SLC5A7, solute carrier family 5 member 7; CDH19, cadherin 19; FDR, false discovery rate.

Discussion

In this study, we performed an integrated analysis of DEGs, DMGs, the target genes of miRNAs, and promoter methylation. Several DEGs, miRNAs, and TFs were screened with respect to metabolism and cell apoptosis found in various diseases, and which are known to play a significant role in the progression of GAC. Zhao et al. (33) identified several novel genetic aberrations, gene network modules, and miRNA target interactions within TCGA. However, their results differed slightly from ours, likely due to the bioinformatics tools used. Additionally, as the pathogenesis of GAC remains unclear, the results presented here may help to clarify the molecular mechanisms of this disease.

Our GO analysis demonstrated enrichment of DEGs in the ECM. Extracellular composition and organization are regulated to control cell behavior and differentiation; however, dysregulation of ECM dynamics leads to the development of various diseases, such as cancer (34). Moreover, the DMGs in our study were enriched in transcriptional activator and TF activity. For example, the promoter region contains putative binding sites for multiple TFs, including the signal transducer and activator of transcription 6 (STAT6), a downstream effector of IL-4 (35). The GO enrichment analysis also revealed that changes in the target genes of miRNAs were primarily enriched in protein kinase activity, which is known to be related to cancer proliferation, cell invasion, and migration (36-40). In terms of promoter methylation, GO was enriched in the transcription regulatory region for DNA binding. Changes in DNA methylation associated with GC have been previously studied and provide additional clues into the pathogenesis of the disease (41).

The KEGG analysis showed that DEGs were primarily enriched in the IL-17 signaling pathway, cytokine-cytokine receptor interaction, and cell cycling, which is consistent with previous reports (42). It has been reported that IL-17 plays an important role in the upregulation of VEGF to promote tumor angiogenesis. We also identified that DMGs were enriched in the cAMP signaling pathway, RAP1 signaling pathway, and cell adhesion molecules. cAMP interacts with the HERG protein by binding to the cAMP-binding domain of HERG protein and subsequently impacts HERG in GAC (43). Meanwhile, reduced RAP1 signaling, which has been reported in prostate cancer, is related to cell adhesion, migration, and survival associated with metastasis. Cell adhesion molecules were in accordance with the studies above (44). Furthermore, the target genes of miRNAs were enriched in proteoglycans in cancer, the RAP1 signaling pathway, MAPK signaling pathway, pathways in cancer, and RAS signaling pathways. Proteoglycans perform multiple functions, including the regulation of tumor cell growth, survival, adhesion, metastasis, and angiogenesis through interactions between their charged glycosaminoglycan chains and effector proteins (45).

The core protein of proteoglycans can also work with other proteins, such as integrins, to accommodate their signaling (46). The MAPK signaling pathway is thought to be involved in a myriad of mechanisms associated with eukaryotic cell regulation by coordinating with diverse receptor families, such as growth factors closely related to the immune system (47). In addition, MAPK regulates the activation of gene transcription, protein synthesis, cell cycle machinery, cell death, and differentiation (47). RAS also plays a key role in malignant cells, including the deregulation of tumor cell growth, programmed cell death, and the ability to invade and induce new blood vessel formation, making it a new therapeutic target for GAC (48,49). Promoter methylation was also significantly enriched in cell adhesion molecules. These results provide the foundation for further GAC research; however, additional studies are required to verify our findings.

TF expression and histone modifications are under strict control. Hence, during dynamic packaging, maintaining stable and ordered chromatin is vital to ensure normal cellular homeostasis. Factors influencing the number of chromatin-associated cellular events, including transcription, warped histones, and DNA, are subject to covalent post-translational modifications (50). In the present study, a comparative analysis of histone modifications in tumor and normal tissues was conducted, which revealed that the DEGs were primarily regulated by H3K27me3 in several cancer cell lines. H3K27me3, a post-translational modification, is highly correlated with genomic silencing (51). Many studies have reported that cytosine methylation of DNA regulatory sequences is associated with the transcriptional inactivation of genes, while hypomethylation supports transcriptional activation (52,53). Zhang et al. explored the relevance of H3K27me3 and DNA methylation in terms of molecular mechanisms to better understand the relationship between DNA methylation and histone modifications in cancer-associated gene silencing (53). Our results suggest that abnormal modification of H3K27me3 may play an important role in GAC with regulation predominantly executed by H3K27me3. Other studies have reported that hypomethylated genes in H3K27me3 activate transcription, whereas hypermethylated genes inactivate transcription (52,54).

The enrichment analysis performed with ChIP-seq revealed that DEGs were enriched in different TF clusters, including E2F4, FOXM1, and SUZ12. E2F4, a member of E2F, regulates various cellular functions related to the cell cycle and apoptosis (55). In fact, GAC frequently exhibits a mutation of the adenosine-guanine-cytosine (AGC) repeat of the E2F4 gene. Furthermore, E2F4 appears to promote tumor progression in GAC (56,57). FOXM1 participates in cell cycling by regulating both the transition from G1 to S and progression to mitosis (58). Moreover, FOXM1 mediates the promotion of human GAC angiogenesis, growth, and metastasis, and enhanced regulation of FOXM1 leads to the acceleration of GAC (58,59). SUZ12 is primarily involved in histone modification (60), whereas HOTAIR expression positively correlates with SUZ12 expression levels and, therefore, may affect the epigenetic state of cancer tissues. A high expression of HOTAIR is associated with higher stages and lymph node metastasis in GAC (61). Consequently, a high expression of SUZ12 may contribute to GAC progression.

The constructed differential miRNAs and promoter methylation networks showed that the hub nodes included SLC5A7 and CDH19. Among these hub genes, ADAM12 and TFAP2A showed the highest node degrees. Significant changes in the expression of key upstream genes may affect a large number of downstream target genes. Indeed, altered DNA methylation in the promoter region of genes has been shown to cause the inactivation of tumor suppressors and other cancer-related genes and is regarded as the most well-defined epigenetic characteristic in GC (62). ADAM12 is significantly upregulated in GAC, inhibiting cancer proliferation, invasion, and metastasis (63). In addition, the upregulation of hypomethylated ADAM12 has also been reported in breast, pancreatic, and ovarian cancers, and may serve as a valuable biomarker for diagnosis, treatment, and prognosis (64,65). In the cellular context, TFAP2A is individually related to cell differentiation and development and cancer progression/regression. Wang et al. suggested that a reduced expression of TFAP2P is associated with GAC prognosis and may be a potential marker for improved prognosis in GAC patients (66). Moreover, ZNF471 recruits KAP1 to the promoter of target genes, thereby inducing H3K9me3 enrichment and repressing oncogenic TFAP2A (67). However, few studies have examined SLC5A7 in the context of GAC. The downregulated expression of SLC5A7 in various cancers is reportedly associated with poor clinical outcomes (68). For instance, SLC5A7 was found to be downregulated and hypermethylated in lung adenocarcinoma (LUAD). However, it was shown to be upregulated and hypomethylated in squamous cell carcinoma (LUSC) (69). The overall survival (OS) of patients with high promotor region methylation of LUSC was better than for those with low methylation; however, no significant correlation was observed between methylation and OS in LUAD (69). CDH19 encodes a tyrosine phosphatase associated with adherent junctions. As GC develops from the same glands as colorectal tumors, CDH19 may also be related to GAC (70). Furthermore, new human antibodies targeting CDH19 have been applied to prevent, treat, and ameliorate melanoma and metastatic melanoma disease.

The generated network of DEGs, differential miRNAs, and promoter methylation constructed in this study identified two significant differentially expressed miRNAs, miR-1 and miR133a, with miR-133a showing the highest node degree (five). miRNAs can distinguish between cancer and normal tissues or various cancers (71,72), and may regulate diverse protein-coding genes and participate in different molecular pathways associated with tumor evolution and progression (73-75). Recent studies have claimed that miR-1 and miR-133a function as tumor suppressors, and their oncogenes have been reported in different cell environments (76) . Indeed, overexpression of miR-1 can inhibit cell growth, differentiation, and migration, while its downregulation promotes cell proliferation and migration in endothelial cells of GC and activates proangiogenesis-associated signaling (77). Moreover, the potential anti-tumorigenic function of miR-1 in primary HCCs is related to the methylation-mediated silencing of miR-1 and upregulation of the potential oncogenic target genes of miR-1 (78). Recently, miR-1 has been suggested as a therapeutic target for prostate cancer and NSCLC (79,80).

Furthermore, reduced miR-133a expression, caused by hypermethylation of its promoter or negative regulation of TFs, was associated with larger and more invasive GAC tumors (81,82), resulting in metastasis of peripheral organs (83). In addition, altered miR-133a expression has been reported in various types of cancer, including esophageal squamous cell carcinoma and bladder cancer, mainly regulating the expression of FSCN1 as a tumor suppressor (84,85). DNA methylation could, directly and indirectly, affect the drug response. Directly, DNA methylation involves gene transcription and thus impacts the methylation level of epigenome targets. Target genes and the expression of proteins affected by epigenome targets is the indirect way that is imposed on the response to agents of cancer cells (86). Therefore, as the key DEGs, miRNAs, and TFs mentioned above and based on the relationship with methylation, the outcome of treatment with GAC could be predicted by detecting hypermethylated or hypomethylated levels of the corresponding targets.

Clinical or/and pathological staging at diagnosis and treatment are practical components that determine the prognosis of GAC patients (87). Methylation-driven gene expression could be predicted by hypomethylation and hypermethylation of corresponding genes exactly based on bioinformatic analysis (88). Nowadays, Gastrointestinal endoscopy (GIS) has served as one of the most widely used diagnostic tools for pathological analysis, is neither suitable for early diagnosis nor free from the risk of morbidity, and detecting miRNAs in body fluids can be a potential method in the diagnosis of GAC. It has been demonstrated that endogenous miRNA in serum or plasma remains stable. Whether in extreme conditions such as boiling, very low/high pH levels, extended storage time, and multiple freeze-thaw cycles (89,90), or stored for more than a decade in archival tissues and human serum specimens, miRNA can be easily detected (91). The stability and resistance to degradation of miRNAs make them more useful biomarkers for cancer diagnosis. Q-RT-PCR is widely utilized to identify previously unknown new miRNAs, which is easier to perform and is a cost-effective technique. Blood-based (92) and gastric juice-based (93) miRNAs are also considered potential biomarkers for early detection of GAC. miRNA-125b was identified as a crucial miRNA in GC development from miRNA-gene network analyses related to gastric oncogenesis (93). Another study found that hsa-miR-421/hsa-miR-29b-1-5p targeted CREBZF and could play an important role in MKN-74 cell migration, suggesting that increased CREBZF by inhibition of hsa-miR-421/hsa-miR-29b-1-5p could be pivotal in preventing gastric cancer in its early stage (94). Moreover, proliferation, as well as migration and invasion abilities of GAC cells, were restricted by the overexpression of HAND2-AS1 and HIF3A and enhanced by miR-184. Furthermore, HAND2-AS1 rescued enhanced GAC cell proliferation, migration, and invasion abilities as well as the glycolytic process caused by hypoxia via miR-184/HIF3A (95). Therefore, once preclinical screening and validation are resolved, it will promote the development of miRNA biomarkers in body fluids for clinical applications.

Conclusions

Our study reports on an integrated analysis of high-throughput sequencing data for DEGs, differential miRNAs, and promoter methylation associated with GAC. Four possible GAC-related DEGs, including ADAM12, TFAP2A, SCL5A7, and CDH19, and two possible GAC-related miRNAs, miR-133a and miR-1, were identified. Moreover, DMRs were screened and integrated. GAC-related promoter methylation may be under comparable transcriptional regulation. Furthermore, several TFs and miRNAs, as well as promoter methylation, may play critical roles in GAC tumorigenesis. These results provide the foundation for further GAC research and will need to be confirmed by additional experiments.

The article’s supplementary files as