Bioinformatic Dissecting of TP53 Regulation Pathway Underlying Butyrate-induced Histone Modification in Epigenetic Regulation.

Butyrate affects cell proliferation, differentiation, and motility. Butyrate inhibits histone deacetylase (HDAC) activities and induces cell-cycle arrest and apoptosis. TP53 is one of the most active upstream regulators discovered by ingenuity pathways analysis (IPA) in our RNA-sequencing data set. TP53 signaling pathway plays key role in many cellular processes. TP53 pathway and their involvement in cellular functions modified by butyrate treatment were scrutinized in this report by data mining the RNA-sequencing data using IPA (Ingenuity System(®)). The TP53 mechanistic pathway targets more than 600 genes. Downstream analysis predicted the activation of the TP53 pathway after butyrate treatment. The data mining also revealed that nine transcription factors are downstream regulators in TP53 signaling pathways. The analysis results also indicated that butyrate not only inhibits the HDAC activities, but also regulates genes encoding the HDAC enzymes through modification of histones and epigenomic landscape.

Introduction

Epigenomics is one of the most rapidly expanding areas in biology, and epigenomic mechanisms, such as histone modification and DNA methylation, play important roles in nearly every aspect of molecular and cellular biology.1 Epigenomics is the study of both heritable and non-heritable changes in the regulation of gene activity and expression that occur without an alteration in the DNA sequence. This dynamic and rapidly developing discipline is making its significant impact across the biomedical sciences. Epigenomic studies have been utilized to discover gene regulatory mechanisms, to identify and detect disease biomarkers, as well as for growing implication for strategies to cancer control.2 Histone modifications, especially histone acetylation and methylation, play a dominant role in epigenomic regulation. After the discovery of histone acetylation in the 1960s,3,4 it was suggested that this post- translational modification of histone could provide an enzymatic mechanism for modulating the interactions between histones and DNA in ways that affect the structure and function of chromatin.5 Correlations between gene activation for RNA synthesis and increased acetylation of the histones were first reported during the 1970s.6 The simple postulation of such activity of histone modification was that acetylation of the lysine residues in the basic amino-terminal regions of the histones neutralizes their positive charges and would be expected to weaken their interactions with the phosphate groups of the DNA strand enveloping the nucleosome core.7 Now, after the discovery of histone acetylation and methylation over 50 years ago,4 modifications of histones are strongly implicated in virtually all genomic processes.

Beyond their nutritional impact, short-chain fatty acids (SCFAs) (especially butyrate) modulate cell differentiation, proliferation, motility, and in particular, induce cell-cycle arrest and apoptosis. All three major components of short-chain fatty acids (acetate, propionate, and butyrate) induce apoptosis and inhibit cell proliferation. However, butyrate was the most potent in induction of apoptosis and inhibition of cell proliferation.8,9 It is also proved to be a strong regulator of genomic activities. During the late 1960s and early 1970s, histone deacetylation, the reverse action of histone acetylation; and the enzyme responsible for histone deacetylation, histone deacetylase (HDAC), were discovered.10,11 Thereafter, butyrate became the first chemical compound discovered with inhibitory activity against histone deacetylation during the late 1970s.7,12,13 Butyrate treatment enhances histone acetylation by inhibiting HDAC activities. Butyrate-mediated regulation of apoptosis, gene expression, cell cycle, transport, growth, and proliferation were shown on a cellular level in the Madin-Darby bovine kidney (MDBK) cell system.14,15 Butyrate-induced epigenomic regulation was evaluated with both microarray and RNA-sequence technologies in our laboratory.15,16 Therefore, this SCFA provides an excellent in vitro model for studying the epigenomic regulation of gene expression induced by histone acetylation. We reported earlier that there are about 11,408 genes significantly differentially expressed after butyrate treatment in bovine cells.15 TP53, one of the most important transcription factors, was found in the center of butyrate-induced regulation.15 For a deeper understanding of the transcriptome alterations induced by butyrate and the identification of the potential mechanisms underlying gene expression and the epigenomic regulation of cellular functions induced by butyrate, we conducted a bioinformatic analysis of the functional category, pathway and integrated network, using Ingenuity Pathways Knowledge Base (http://www.ingenuity.com) by mining the RNA-sequencing dataset to dissect TP53 regulation pathway. The major findings are presented in this report.

Method and Materials

Cell culture and butyrate treatment

As described in our earlier publication,15 MDBK epithelial cells were cultured in Eagle’s minimal essential medium and supplemented with 5% fetal bovine serum (Invitrogen, Carlsbad, CA, USA) in 25 cm2 flasks. During the exponential phase, the cells were treated for 24 hours with 10 mM sodium butyrate (Calbiochem, San Diego, CA, USA). Four replicate flasks of cells for both treatment and control groups (ie a total of eight samples) were used for the RNA-seq experiments.

RNA extraction and sequencing using RNA-seq

Total RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA, USA) followed by DNase digestion and Qiagen RNeasy column purification (Qiagen, Valencia, CA, USA), as previously described.17 High-quality RNA (RNA integrity number (RIN) > 9.0) was processed using an Illumina TruSeq RNA sample prep kit following the manufacturer’s instruction (Illumina, San Diego, CA, USA). After quality control procedures, individual RNA-seq libraries were then pooled based on their respective 6-bp adaptors and sequenced at 50 bp/sequence, read using an Illumina HiSeq 2000 sequencer, as described previously.18 Approximately 67.5 million reads per sample (mean ± sd = 67,527,111 ± 8,330,388.6) were generated.

Data analysis and bioinformatics

RNA-seq reads were aligned to the bovine reference genome (Btau 4.0) using TopHat.19 Each SAM (sequence alignment/map) output file per sample from TopHat alignment, along with the GTF (general transfer format) file from ENSEMBL bovine gene-build v65.0, was used in the Cuffdiff program in the Cufflink package (v1.3.0) as input files20 to test differential expression. Mapped reads were normalized based on the upper-quartile normalization method, and Cuffdiff modeled the variance in fragment counts across replicates using the negative binomial distribution as described previously.15

Differentially expressed genes in the transcriptome were further analyzed using gene ontology (GO) analysis (GOseq).21 Enrichment of certain GO terms was determined based on Fisher’s exact test. A multiple correction control (permutation to control false discovery rate) was implemented to set up the threshold to obtain the lists of significantly over-represented GO terms.

The molecular processes, molecular functions, and genetic networks following butyrate treatment were further evaluated by analyzing differentially expressed genes using ingenuity pathways analysis (IPA) (Ingenuity Systems®, www.ingenuity.com). IPA is a software application that enables biologists to identify the biological mechanisms, pathways, and functions most relevant to their experimental datasets or genes of interest.22–26

Results and Discussion

Upstream analysis reveals TP53 as one of the most active upstream regulator

IPA analysis can help to answer the questions such as what biological processes and cellular functions are affected in the experiment. IPA also reveals what functions are expected to increase or decrease upon gene expression changes in the experiment dataset and what genes are involved in the processes. IPA up/down stream analysis also predicts the effect of gene expression changes in the dataset on biological processes and disease or on cellular functions. IPA uses the regulation z-score algorithm to make predictions. The z-score algorithm is designed to reduce the chance that random data will generate significant predictions. Analysis with Ingenuity Pathways Knowledge Base on our RNA-sequencing datasets15 revealed that 1010 genes are targets of TP53. Among these 1010 genes, 658 genes are differentially expressed after butyrate treatment. Functional analysis by IPA indicated that TP53 gets activated on butyrate treatment with the highest z-score of 5.116 and overlap P-value of 6.39E-42. A total of 317 genes from 658 differentially expressed genes have expression direction consistent with activation of TP53. The 658 genes and their expression values are listed in Supplement Table 1. However, the mechanistic network shows the complicated interaction between the TP53 genes and its targets (Fig. 1, see Supplement Table 1 for the detailed gene list). Downstream analysis is warranted to understand the regulational pathways. The data sets supporting the results of this article have been deposited in National Center for Biotechnology Information (NCBI)’s Gene Expression Omnibus in Sequence Read Archive (SRA). This dataset (SRA051007) included eight separate RNA-seq data with a total of 540.2 million sequences.

Figure 1

TP53 mechanistic network and targets of TP53. More than 600 genes are targets of TP53 revealed by RNA sequencing and IPA analysis. Those highlighted with red color are upregulated genes and those with green dots are downregulated genes. The names of these genes are listed in Supplement Table 1.

Nine downstream regulators are the center of the TP53 mechanistic network of gene regulation

Nine downstream regulators (TP53, TP63, TP73, MYC, RBL1, RB1, E2f, SIN3A, and HDAC) are revealed by IPA analysis as the center of the TP53 mechanistic regulation network (Fig. 2A). Several members of this network are also regulating or regulated by a network of their own. For example, the E2 factor (E2f) family (Fig. 2B) of transcription factors is a unique transcriptional network. E2f factors have been known for long time to be the prominent regulators of S-phase entry. Recent advances in the studies indicated that E2f has functions that reached beyond G1/S control and impacted cell proliferation in several different ways.27 The E2f family of transcription factors, as the downstream effectors of the retinoblastoma protein (RB) pathway, functions as both an activator and repressor of transcription. It also acts in cell proliferation and differentiation, and as a regulator of cell death. Clear evidence has linked modifications in chromatin structure to cell-cycle progression, DNA replication, and overall chromosome stability.28,29 In our previous reports, as a direct result of hyper-acetylation of histones induced by butyrate at physiologic concentrations, cultured bovine cells are arrested in the early G1 phase, and DNA synthesis is blocked. Butyrate also induces apoptosis in an established bovine cell line, specifically MDBK epithelial cell line, at a relatively high concentration.14,30

Figure 2

Dataset was analyzed with the IPA software (Ingenuity Systems®, www.ingeuity.com). (A) Nine transcription factors involved in TP53 mechanistic pathways and the prediction of the regulatory functions in the dataset. The number under the gene’s name is the Log2 fold changes of gene expression induced by butyrate. (B) E2f transcription factor family in the TP53 mechanistic pathway.

In eukaryotic cells, cell-cycle checkpoint regulation assures the fidelity of cell division. The G1 (first gap phase)/S cell-cycle checkpoint controls the passage of eukaryotic cells from the G1 into the S phase. Mitogen-dependent progression through the G1 of cell-division cycle is accurately regulated in an order such that normal cell division is synchronous with cell growth and the initiation of DNA synthesis (S phase) is timed precisely to avoid inappropriate DNA amplification. The G1/S checkpoint control is vital in normal cell division. The key components involved in the G1/S checkpoint are the cell-cycle kinases, CDK4/6-cyclin D and CDK2-cyclin E, as well as the transcription complex that includes RB and transcription factor E2F. The activation of E2F is necessary for the G1–S transition. In the present report, E2F is significantly down regulated by the butyrate-induced histone acetylation. Furthermore, the activities of the family of E2f transcription factors are depressed as predicted in IPA downstream analysis. All these perturbations of the gene expression in the G1/S cell-cycle checkpoint pathways are consistent with the biological effects of butyrate, which induces cell-cycle arrest at the G1/S boundary.14

Cross-talk between TP53 and HDACs

Among SCFAs, butyrate induces the highest acetylation. It was believed that butyrate inhibits class I and class IIa but not class IIb HDACs (HDAC6 and HDAC10). However, more recent studies established more evidence that butyrate not only inhibits the HDAC activities but also affects many other epigenetic-related enzymes by regulating the expression of the respective genes encoding HDACs attributable to the potential of butyrate to alter epigenetic landscapes. In addition, cross-talking between upstream and downstream regulators such as TP53 and HDACs sets up another layer of regulation network. HDACs are active components of transcriptional regulation complexes. HDACs downregulate the functions of p53, the TP53 gene product. It was suggested that interactions of p53 and HDACs likely result in p53 deacetylation, thereby reducing its transcriptional activity.27 Inhibitors of HDAC can restore the p53 pathway.31 Our data are consistent with the findings that butyrate, a HDAC inhibitor, activates the TP53 pathway. In return, activation of TP53 may directly or indirectly regulate the expression of HDACs (Fig. 3). Our RNA-sequencing data revealed that the expression of HDACs was modulated by butyrate. Especially, HDACs 7, 8, 9, and 10 were downregulated and HDACs 3, 5, and 11 were upregulated. Figure 3 shows the interaction between TP53 and HDACs. The downstream analysis also indicates inhibition of the activities of HDACs 2, 4, 6, and 9. The mechanism and biological relevance of HDAC inhibitors in the regulation of expression of HDACs is not clear, but may possibly indicate the existence of an auto-regulatory feedback loop for the expression of several HDACs after their activities are inhibited. It will be of great interest to understand the connections between the inhibition of HDAC activities and the regulation of expression of these genes.

Figure 3

HDAC family as downstream regulators of TP53: dataset was analyzed with the IPA software (Ingenuity Systems®, www.ingeuity.com). The number under the gene’s name is the Log2 fold changes of gene expression induced by butyrate.

Canonical pathway of TP53 signaling

The dynamic canonical pathways are well-characterized metabolic and cell signaling pathways. Canonical pathways are distinct from networks. The pathways (metabolic and signaling pathways) are generated in IPA before data input, based on the literature, and do not change upon data input. The p53 (TP53) tumor suppressor protein is a key transcriptional regulator that responds to a variety of cellular stresses such as DNA damage, UV irradiation, and hypoxia. P53 regulates key cellular processes such as DNA repair, cell-cycle progression, angiogenesis, and apoptosis. These p53-dependent pathways shut down damaged cells, either through apoptosis or cell-cycle arrest. Loss of p53 function is thought to be a contributing factor in the majority of cancer cases. Figure 4 represents the TP53 signaling pathway and butyrate-induced differentially expressed genes involved in TP53 signaling pathway. The RNA-sequencing data strongly support that TP53 signaling pathways are crucially involved in the cellular functions that were moderated by the butyrate treatment. These differentially expressed genes underlie the butyrate-induced modifications of the cellular functions, such as cell-cycle arrest, cell survival, and apoptosis.

Figure 4

Canonical pathway of TP53 signaling: dataset was analyzed with the IPA software (Ingenuity Systems®, www.ingeuity.com). The number under the gene’s name is the Log2 fold changes of gene expression induced by butyrate. The color indicates the expression level of the genes (red indicating upregulated and green indicating downregulated genes).

Conclusions

Bioinformatic dissecting and data mining the RNA-sequencing data using IPA reveals that TP53 plays an important role in the butyrate-induced epigenetic regulation of gene expression. The TP53 mechanistic pathway targets more than 600 genes. Downstream analysis predicted the activation of the TP53 pathway after butyrate treatment. Data mining also revealed that a total of nine transcription factors are involved in TP53 signaling pathways. These transcription factors such as E2f affect the gene functions through the network of their own. Butyrate not only inhibits HDAC activities but also regulates the expression of the genes that encode these enzymes, conceivably through histone modification and changes of epigenomic landscape. In addition, butyrate, as an HDAC inhibitor, can also alter histone methylation, implying the interplay between histone acetylation and histone methylation. These findings will facilitate our understanding of the molecular mechanisms underlying butyrate-induced epigenomic regulation in bovine cells.

Supplementary table 1. List of target genes of TP53. The table was generated through the use of IPA (Ingenuity® Systems, www.ingenuity.com).