Bioinformatic identification of candidate genes induced by trichostatin A in BGC-823 gastric cancer cells.

The aim of the present study was to identify the candidate genes induced by trichostatin A (TSA) in BGC-823 gastric cancer (GC) cells and to explore the possible inhibition mechanism of TSA in GC. Gene expression data were obtained through chip detection, and differentially expressed genes (DEGs) between GC cells treated with TSA and untreated GC cells (control group) were identified. Gene ontology analysis of the DEGs was performed using the database for annotation, visualization and integrated discovery. Then sub-pathway enrichment analysis was performed and a microRNA (miRNA) regulatory network was constructed. We selected 76 DEGs, among which 43 were downregulated genes and 33 were upregulated genes. By sub-pathway enrichment analysis of the DEGs, the propanoate metabolism pathway was selected as the sub-pathway. By constructing a miRNA regulatory network, we identified that DKK1 and KLF13 were the top hub nodes. The propanoate metabolism pathway and the genes DKK1 and KLF13 may play significant roles in the inhibition of GC induced by TSA. These genes may be potential therapeutic targets for GC. However, further experiments are still required to confirm our results.

Introduction

Gastric cancer (GC) is one of the most common malignancies worldwide and a leading cause of cancer-related mortality (1,2). Little is known about GC-associated genes or its diverse clinical properties, including metastatic status, invasiveness, histological type and responsiveness to chemotherapy (3). The carcinogenic mechanism of GC is still not fully understood, although previous studies have demonstrated a number of genetic alterations in this disease (3).

Studies have revealed that trichostatin A (TSA) has an inhibitory role in GC (4–7). In eukaryotic transcription, chromatin modifications by histone acetyltransferase or histone deacetyltransferase (HDAC) represent a fundamental mechanism of transcriptional regulation (8). Evidence from further studies suggests that histone acetylation alters nucleosomal structures and facilitates the accessibility of transcription-associated factors to chromatin DNA through the disruption of interactions between histones and DNA (9–11). As an HDAC inhibitor, TSA has been known to cause a variety of phenotypic changes, including cell-cycle arrest in the G1/G2 phase, apoptosis and differentiation in cultured transformed cells (12–14). Although significant efforts have been made, the candidate genes and inhibition mechanism of TSA in GC remain unclear.

In the present study, microarray data were obtained through chip detection, and the differentially expressed genes (DEGs) between TSA-treated GC cells and untreated GC cells were identified. In addition, gene ontology (GO) analysis of the DEGs was performed. Then sub-pathway enrichment analysis was performed, and a microRNA (miRNA) regulatory network was constructed. Through the identification of GC-associated genes and biological changes, possible molecular mechanisms and potential therapeutic targets for GC were explored.

Materials and methods

Affymetrix microarray data

The GC cell line, BGC-823, was grown in RPMI-1640 supplemented with 10% fetal bovine serum and antibiotics in 5% CO2. The cells were grown to ~70–80% confluence and treated with 330 nM TSA (Sigma-Aldrich, St. Louis, MO, USA) for 24 h. Then total RNA was isolated from TSA-treated cells using TRIzol (Invitrogen Life Technologies, Gaithersburg, MD, USA). The construction of a fluorescence probe and hybridization were performed using a 3DNA array detection kit (Genisphere, Hatfield, PA, USA). cDNA microarrays were performed on a cDNA chip containing 21,522 cDNA clones selected from the Human Genome Oligo set version 2.1 (Operon, Ebersberg, Germany). Each slide contained 12 control genes to normalize the signal intensities of the different fluorescent dyes. Each hybridization array was scanned on a LuxScan 10KA (CapitalBio, Beijing, China). The intensity of each hybridization signal was calculated by the GenePix Pro 4.0 program (Axon, Schaumburg, IL, USA). Three repeats were included in this study, and each repeat contained TSA treated GC cells and untreated cells.

Data processing and DEG analysis

The signal intensities of TSA treated GC cells and untreated GC cells in each repeat were obtained through chip detection and the gene expression ratio (TSA-treated GC cells / untreated GC cells) in each repeat was calculated. Then, the DEGs between the TSA-treated group and the control group were analyzed, and the score was calculated by the significant analysis of microarray (SAM) algorithm (15). Each gene was assigned a difference score (defined as “d”) based on the significance of its gene expression changes between TSA treated GC cells and untreated GC cells. The multiple testing correction was performed using a false discovery rate (FDR) (16). The fold change of the expression of individual genes was also observed for the differential expression test. DEGs with an FDR<0.05 were considered to be significant. DEGs with a score of (d) ≥1.25 and a ratio (TSA-treated GC cells / untreated GC cells) >1.5 in at least two groups were defined as upregulated DEGs, and DEGs with a score of (d) ≤-1.25 and a ratio (TSA-treated GC cells / untreated GC cells) <0.66 in at least two groups were defined as downregulated DEGs.

GO analysis

GO analysis has become a commonly used approach for functional studies of large-scale genomic or transcriptomic data (17). The database for annotation, visualization and integrated discovery (DAVID) (18) consists of an integrated biological knowledge base and analytic tools aimed at systematically extracting biological meaning from large gene or protein lists. The GO function of DEGs was analyzed using DAVID. P<0.05 was considered to indicate a statistically significant difference.

Sub-pathway enrichment analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (19) contains information on how molecules or genes are networked, and is complementary to most of the existing molecular biological databases that contain information on individual molecules or genes. The KEGG pathways of the DEGs were analyzed using DAVID. P<0.05 was considered to indicate a statistically significant difference. The closer the components within the metabolic pathways, the greater the similarity in their biological functions. Therefore, the identification of the sub-pathway of disease is crucial. K-clique was used to divide the metabolic pathway into sub-pathways through the iSubpathwayMiner package in R (https://www.r-project.org/). The sub-pathway with P<0.05 was considered to be significant.

miRNA regulatory network analysis

miRNAs are a type of endogenous non-coding RNA with regulatory functions, and their size is ~20–25 nucleotides. Mature miRNAs are assembled into RNA-induced silencing complexes, then target genes are identified and regulated by them. Thus, the identification of the miRNA regulatory network is essential. We predicted the interactions between target genes and miRNAs based on the TarBase (http://microrna.gr/tarbase) (20), TargetScan (www.targetscan.org) (21) and miRecord (http://mire-cords.biolead.org) (22) databases. The miRNA regulatory network was constructed using Cytoscape (http://cytoscape.org/) (23).

Literature mining analysis of key genes

The online tool GenCLip2.0 was used to mine DEGs in literature (24). A list of DEG names was input into GenCLip2.0, then the Gene Cluster with Literature Profiles module was used to cluster based on the documented frequency of the input genes. GenCLip2.0 also provided GO and pathway enrichment analysis of genes. The GC research-related DEGs were mined from the literature and these DEGs were submitted to GO and pathway enrichment analyses.

Results

Screening and GO analysis of DEGs

Three groups of gene expression ratio (TSA-treated GC cells / untreated GC cells) about TSA-treated group and control group were obtained through chip detection. A total of 76 DEGs (43 downregulated genes and 33 upregulated genes) were obtained. The results of GO analysis revealed that the upregulated DEGs were mainly enriched in biological processes including skeletal muscle contraction, response to organic substance and regulation of cell adhesion, and the downregulated DEGs were mainly enriched in processes including the DNA catabolic process, DNA metabolic process and response to the tumor necrosis factor (Table I).

Table I.

Gene ontology analysis results of differentially expressed genes.

Gene type	Category	Term	Count	P-value
Upregulated	GOTERM_BP_FAT	GO:0003010~voluntary skeletal muscle contraction	2	0.004135
	GOTERM_BP_FAT	GO:0014721~twitch skeletal muscle contraction	2	0.004135
	GOTERM_BP_FAT	GO:0031444~slow-twitch skeletal muscle fiber contraction	2	0.004135
	GOTERM_BP_FAT	GO:0010033~response to organic substance	6	0.015009
	GOTERM_BP_FAT	GO:0003009~skeletal muscle contraction	2	0.024566
	GOTERM_BP_FAT	GO:0030155~regulation of cell adhesion	3	0.032410
	GOTERM_BP_FAT	GO:0050881~musculoskeletal movement	2	0.038627
	GOTERM_BP_FAT	GO:0050879~multicellular organismal movement	2	0.038627
	GOTERM_BP_FAT	GO:0009725~response to hormone stimulus	4	0.039305
	GOTERM_BP_FAT	GO:0050873~brown fat cell differentiation	2	0.048550
	GOTERM_CC_FAT	GO:0005576~extracellular region	10	0.010773
	GOTERM_CC_FAT	GO:0005615~extracellular space	5	0.042191
	GOTERM_MF_FAT	GO:0003735~structural constituent of ribosome	3	0.041104
Downregulated	GOTERM_BP_FAT	GO:0006308~DNA catabolic process	3	0.004610
	GOTERM_BP_FAT	GO:0006259~DNA metabolic process	5	0.009728
	GOTERM_BP_FAT	GO:0006301~post-replication repair	2	0.013524
	GOTERM_BP_FAT	GO:0034612~response to tumor necrosis factor	2	0.021888
	GOTERM_BP_FAT	GO:0051129~negative regulation of cellular component organization	3	0.023964
	GOTERM_BP_FAT	GO:0065003~macromolecular complex assembly	5	0.024316
	GOTERM_BP_FAT	GO:0006974~response to DNA damage stimulus	4	0.024466
	GOTERM_BP_FAT	GO:0070647~protein modification by small protein conjugation or removal	3	0.029893
	GOTERM_BP_FAT	GO:0043933~macromolecular complex subunit organization	5	0.030046
	GOTERM_BP_FAT	GO:0006309~DNA fragmentation involved in apoptosis	2	0.031835
	GOTERM_BP_FAT	GO:0044265~cellular macromolecule catabolic process	5	0.032124
	GOTERM_BP_FAT	GO:0006921~cell structure disassembly during apoptosis	2	0.038412
	GOTERM_BP_FAT	GO:0009057~macromolecule catabolic process	5	0.040630
	GOTERM_BP_FAT	GO:0000737~DNA catabolic process, endonucleolytic	2	0.043317
	GOTERM_BP_FAT	GO:0030262~apoptotic nuclear changes	2	0.043317
	GOTERM_BP_FAT	GO:0009314~response to radiation	3	0.044895
	GOTERM_CC_FAT	GO:0043232~intracellular non-membrane-bounded organelle	9	0.025314
	GOTERM_CC_FAT	GO:0043228~non-membrane-bounded organelle	9	0.025314
	GOTERM_CC_FAT	GO:0070013~intracellular organelle lumen	7	0.038954
	GOTERM_CC_FAT	GO:0043233~organelle lumen	7	0.042966
	GOTERM_CC_FAT	GO:0031974~membrane-enclosed lumen	7	0.046708

BP, biological process; CC, biological process; MF, molecular function.

The sub-pathway of disease was screened. The propanoate metabolism pathway was the sub-pathway selected in our study. Propionyl-coenzyme A carboxylase, beta polypeptide (PCCB) was the DEG enriched in this sub-pathway (Fig. 1).

Figure 1.

Sub-pathway enrichment analysis of differentially expressed genes (DEGs). Digital nodes refer to enzymes; letter nodes refer to genes; the red-bordered node refers to DEGs enriched in the pathway; the lines represent the interactions.

According to the miRNA regulatory network, the upregulated genes dickkopf wnt signaling pathway inhibitor 1 (DKK1), follistatin (FST), clusterin (CLU), inhibitor of DNA binding 3, dominant negative helix-loop-helix protein (ID3), cysteine-rich, angiogenic inducer 61 (CYR61), x-box binding protein 1 (XBP1) and small nuclear ribonucleoprotein polypeptide N (SNRPN) were the hub nodes, and among them DKK1 was the top hub node (Fig. 2A). From the miRNA regulatory network, the downregulated genes kruppel-like factor 13 (KLF13), breast cancer 1, early onset (BRCA1), caspase 3, apoptosis-related cysteine peptidase (CASP3), 2′-5′-oligoadenylate synthetase-like (OASL) and histone cluster 1, H2ac (HIST1H2AC) were the hub nodes, and among them KLF13 was the top hub node (Fig. 2B).

Figure 2.

miRNA regulatory network. (A) miRNA regulatory network constructed with miRNA and upregulated differentially expressed genes (DEGs). Red nodes refer to upregulated DEGs; blue nodes refer to miRNAs; the lines represent the interactions. (B) miRNA regulatory network constructed with miRNA and downregulated DEGs. Green nodes refer to downregulated DEGs; blue nodes refer to miRNAs; the lines represent the interactions.

A heatmap of DEGs was constructed (Fig. 3). From the heatmap, the upregulated genes were mainly enriched in functions including epidermal growth factor, inflammation and tumor necrosis factor; the downregulated genes were mainly enriched in biological processes including squamous cell carcinoma entries, apoptosis and the DNA replication process. The results of GO analysis revealed that the upregulated genes were mainly enriched in functions including extracellular space and response to organic substance (Table II); the downregulated genes were mainly enriched in the DNA metabolic process (Table II).

Figure 3.

Literature mining results of the key genes. (A) Results for upregulated differentially expressed genes (DEGs) in the heatmap; (B) results for downregulated DEGs in the heatmap.

Table II.

Gene ontology analysis results of the key genes.

Gene type	Pathway	Hit enrichment score	P-value	Q-value	GO ID
Upregulated	Cluster1	4.97
	Extracellular space	8.00	4.72E-06	0.002702511	GO:0005615
	Extracellular region	13.00	2.42E-05	0.004620421	GO:0005576
	Single1	4.90
	Response to organic substance	13.00	1.27E-05	0.003630953	GO:0010033
	Single2	4.13
	BMP signaling pathway	4.00	7.42E-05	0.010635069	GO:0030509
Downregulated	Single1	6.14
	DNA metabolic process	8.00	7.25E-07	0.000332722	GO:0006259

GO, gene ontology.

Discussion

GC is one of the most common malignant diseases and a leading cause of cancer-related mortality worldwide (1,2). The mechanism of TSA inhibition in GC is unclear. In the present study, we used bioinformatics approaches to explore the candidate genes induced by TSA in GC, and identified 76 DEGs between TSA-treated GC cells and untreated GC cells. Among these DEGs, 43 genes were upregulated and 33 were downregulated. In addition, sub-pathway enrichment analysis was performed, and the propanoate metabolism sub-pathway was identified. Then the miRNA regulatory network was constructed; from this, DKK1 was identified as the top hub node among the upregulated DEGs, while KLF13 was the top hub node among the downregulated DEGs.

The propanoate metabolism pathway was the sub-pathway enriched in our study. To date, there have been no studies on the correlation between the propanoate metabolism pathway and GC. However, certain studies have reported that this sub-pathway is associated with other cancer types. Zhu et al observed that bladder cancer-related long noncoding RNAs (lncRNAs) demonstrated a significant association with the propanoate metabolism pathway (25). Furthermore, certain studies have reported that dysregulation of lncRNAs is a primary feature of several human cancers, including prostate cancer, breast cancer, GC, bladder cancer and hepatocellular carcinoma (26–32). The kidney cancer-related proteins were demonstrated to be involved with a high degree of confidence in propanoate metabolism, pyruvate metabolism, the urea cycle and arginine pathways (33). Gmeiner et al observed that the downregulated DEGs in colon cancer were involved with propanoate metabolism and the isoleucine degradation pathway (34). Yang et al revealed that valine relating to lung cancer was involved in a number of metabolic processes including the propanoate metabolism process (35). Gonzalez-Angulo et al demonstrated that the metabolic pathways (glutamate metabolism, chondroitin sulfate biosynthesis and propanoate metabolism) were the most consistently upregulated in basal-like residual cancers (36). Considering the studies above, the propanoate metabolism pathway may play a significant role in TSA inhibition in GC, and may be a therapeutic target for GC.

DKK1 and KLF13 were observed to be hub nodes in our miRNA regulatory networks. We revealed that DKK1 was an upregulated DEG. Hirata et al suggested that DKK1 induced cell apoptosis and inhibited cell growth in renal cell carcinoma (RCC) (37). The levels of DKK1 were decreased in RCC cell lines, but they increased following treatment with TSA. These authors also revealed that TSA induced histone acetylation at the DKK1 promoter, which resulted in the reactivation of DKK1 expression in human RCC. This result demonstrated that the silencing of DKK1 is caused by histone modification (37). Lee et al demonstrated that DKK1 was repressed by histone deacetylation in cervical cancer cells (38). The DKK family was demonstrated to suppress cell growth and induce cell apoptosis in gastrointestinal and colon cancer cells (39,40). Studies have revealed that DKK1 is preferentially expressed in lung cancer, and it was identified as a novel prognostic biomarker as well as a therapeutic target for lung cancers (41,42). Certain studies demonstrated that DKK1 significantly reduced tumor growth (39,43). DKK1 levels may aid in diagnosing GC, and the gene may be a novel prognostic marker for GC (35). Our results revealed that DKK1 was an upregulated DEG in GC cells treated with TSA. Furthermore, the literature mining results from our study revealed that DKK1 was mainly enriched in functions including cell growth, cell differentiation and tumor necrosis factor; these functions have been demonstrated to be correlated with certain cancers (44,45). In our study, DKK1 was overexpressed in GC cells induced by TSA. Considering the studies above, one of the possible inhibition mechanisms of TSA in GC may be that TSA induces histone acetylation at the DKK1 promoter, which results in the reactivation of DKK1 expression (37,38). In other words, DKK1 may be a significant therapeutic target for GC.

KLF13 is a zinc finger transcription factor known to play a role in proliferation, differentiation, cell cycle progression and apoptosis (46–48). KLF13 is required for the expression of several oncogenes including cyclin D1, which is a known oncogene in oral squamous cell carcinoma (47,49). Furthermore, KLF13 is believed to play a role in cancer, and altered expression of it contributes to tumorigenesis (50–53). Henson et al revealed that KLF13 was overexpressed in oral cancer cells, and artificially reducing its cellular levels decreases cell proliferation and malignancy; therefore, KLF13 may be a useful biomarker for early detection and a possible target for therapy (49). KLF13 is also involved in the proliferation and differentiation of the heart (47,54). KLF13 was downregulated in our study. All of the studies above indicated that decreased expression of KLF13 may be a possible inhibition mechanism of GC induced by TSA. In summary, KLF13 may be utilized as a biomarker for detection and a new therapeutic target for GC.

In conclusion, the propanoate metabolism pathway and the associated genes DKK1 and KLF13 may play significant roles in the inhibition of GC induced by TSA. They may be possible therapeutic targets for GC. However, further studies are necessary to verify the clinical applications of these genes as biological targets for GC treatment.