Transcriptome Characterization by RNA-seq Unravels the Mechanisms of Butyrate-Induced Epigenomic Regulation in Bovine Cells

Short-chain fatty acids (SCFAs), especially butyrate, affect cell differentiation, proliferation, and motility. Butyrate also induces cell cycle arrest and apoptosis through its inhibition of histone deacetylases (HDACs). In addition, butyrate is a potent inducer of histone hyper-acetylation in cells. Therefore, this SCFA provides an excellent in vitro model for studying the epigenomic regulation of gene expression induced by histone acetylation. In this study, we analyzed the differential in vitro expression of genes induced by butyrate in bovine epithelial cells by using deep RNA-sequencing technology (RNA-seq). The number of sequences read, ranging from 57,303,693 to 78,933,744, were generated per sample. Approximately 11,408 genes were significantly impacted by butyrate, with a false discovery rate (FDR) <0.05. The predominant cellular processes affected by butyrate included cell morphological changes, cell cycle arrest, and apoptosis. Our results provided insight into the transcriptome alterations induced by butyrate, which will undoubtedly facilitate our understanding of the molecular mechanisms underlying butyrate-induced epigenomic regulation in bovine cells.


Introduction
Short-chain fatty acids (SCFAs), such as acetate, propionate, and butyrate, are important nutrients in ruminants. SCFAs are produced during the microbial fermentation of dietary fiber in the gastrointestinal tract and are directly absorbed at the site of production and oxidized for cell energy production and use [1]. In humans, colonic microbiota convert dietary fiber into prodigious amounts of SCFAs that benefit the human host through numerous metabolic, trophic, and chemopreventative effects [2]. The SCFA butyrate, in particular, also serves as an inhibitor of histone deacetylases (HDACs), which are critical epigenetic regulators [3,4,5]. Therefore, butyrate could act to reactivate epigenetically silenced genes by increasing global histone acetylation [6]. Epigenetic modifications play a key role in the regulation of gene expression, and HDAC activity contributes significantly to epigenetic modification. The HDACs are part of a transcriptional co-repressor complex that influences various tumor suppressor genes. HDACs also play significant roles in several human cancers, making HDAC inhibitors an important emerging class of chemotherapeutic agents.
Chromatin modification has evidently evolved to be a very important mechanism for the epigenetic regulation of the transcriptional status of a genome [4]. Butyrate is not only important for its nutritional impact. It also has profound impacts at the gene level, altering cell differentiation, proliferation, and motility and inducing cell cycle arrest and apoptosis [3]. The foremost biochemical change induced by butyrate and other HDAC inhibitors is the global hyper-acetylation of histones [3,7]. Clear evidence has linked modifications in chromatin structure to cell cycle progression, DNA replication, and overall chromosome stability [8,9]. Cultured bovine cells respond to the hyperacetylation of histones induced by butyrate at physiological concentrations by arrest in the early G1 phase and the cessation of DNA synthesis. Butyrate at a relatively high concentration also induces apoptosis in an established bovine cell line, the Madin-Darby bovine kidney epithelial cell line (MDBK) [3]. The modulation of genome expression through chromatin structural changes by processes such as histone acetylation is considered a major genetic control mechanism.
Histone lysine acetylation has emerged as an essential regulator of genome organization and function. As a HDAC inhibitor (HDACi), butyrate is a strong inducer of the hyper-acetylation of histone in cells and provides an excellent in vitro model for the study of the epigenomic regulation of gene expression induced by histone acetylation. An investigation of the global gene expression profiles of MDBK cells and their regulation by sodium butyrate has recently been conducted using a high-density oligonucleotide microarray [10]. The profound changes observed in gene expression in bovine cells following butyrate treatment demonstrate the pleiotropic effects of histone acetylation [5]. As nutrition research shifts from epidemiology and physiology to the study of molecular interactions with the genome and the elucidation of these less-obvious nutritional effects, a detailed knowledge of changes in gene expression becomes necessary as a basis for understanding these molecular mechanisms.
In the present study, we report our findings on the function and pathways induced by butyrate in MDBK cells. We used deep RNA sequencing to provide a significant amount of novel gene information for bovine cell transcription, which can then be used for further transcriptomic studies or to gain a deeper understanding of the bovine genome and transcriptome. This study also provides a significant amount of information for the epigenetic regulation induced by butyrate. Our data show that butyrateinduced histone acetylation results in subsequent changes in the accessibility of the DNA to transcription activities. Transcriptomic characterization using deep RNA sequencing facilitates the identification of the potential mechanisms underlying gene expression and the epigenomic regulation of cellular functions induced by butyrate.

Butyrate treatment induces changes in cell morphology and cell cycle arrest
We previously reported that butyrate induces cell cycle arrest in MDBK cells. In preparation for deep RNA sequencing, we first endeavored to confirm that the butyrate induced cell cycle arrest. When cells were treated with 10 mM butyrate for 24 hours, cell morphology became distorted. Cells with large vacuoles, with ragged membranes, lacking distinct intracellular organelles, and having increased spaces between cells were readily visible and recurrent. Flow cytometry analysis of the cell population profiles for DNA content and BrdU labeling also confirmed that the cells were arrested at the G1 and G1/S boundary. The incorporation of the BrdU label suggested that DNA synthesis was blocked by butyrate treatment. Western blotting also confirmed that butyrate induced the hyper-acetylation of H3 ( Figure 1).

RNA-seq provided a comprehensive view of the bovine cell transcriptome
In total, 57,303,693 to 78,933,744 sequence reads were generated per sample (Table 1), and 24,526 genes had at last one sequence hit in at least one sample. Of these, 16,212 genes were shared by all samples and can be considered to be the core transcriptome of the bovine epithelial cell. A mean value of 19,4776155 (mean6SD) genes was detected in the butyratetreated group, while 17,6266125 (mean6SD) genes were detected in the control group. Table 2 summarizes the alignment results. Among these genes, 11,408 genes showed a significant differential expressed at a strict false discovery rate (FDR) ,0.05 (Table S1).
Previous gene expression profiling in MDBK cells and the induction of histone acetylation by butyrate has been analyzed by using bovine oligonucleotide microarray. In this previous study, 30 genes representing different expression levels and functional classes were selected and validated by real-time RT-PCR [5]. We were able to confirm over 70% of the differentially regulated genes that were identified by the microarray experiment using RNA-seq. However, RNA-seq allowed us to identify a significantly greater numbers of genes that were induced by butyrate, but which had not been previously associated with the biological effects of butyrate. Transcriptome characterization by RNA-seq also identified 587 genes that were uniquely expressed in butyratetreated cells, but had not been previously detected by a microarray experiment in cells given similar treatments [5].

Functional annotation of differentially expressed genes induced by butyrate
The biological relevance of butyrate-induced gene expression in bovine cells was explored by the Gene Ontology (GO) classification. Table 3 lists 65 GO terms that were significantly perturbed by butyrate treatment. The most-represented biological processes and molecular functions, sorted by statistical significance in both terms of p-value and FDR, included nucleic acid metabolic process, DNA metabolic processes, the regulation of the cell cycle, and DNA replication.
Global function and pathway analyses identified the mechanism of butyrate-induced cell cycle arrest The functional category and pathway analysis of differentially expressed genes in cells treated with butyrate were explored using the IPA (Ingenuity Pathways Analysis, IngenuityH Systems, www. ingenuity.com) Knowledge Base. Of the 24,525 genes in the data set, 13,885 genes were mapped, and 10,637 genes were not mapped in the database. These genes were uploaded for IPA. Among the 13,885 mapped genes, 8,862 genes were identified with matched gene symbols and were used in pathway analysis. Of these, 5,542 genes were significantly up-regulated, while 3,320 genes were significantly down-regulated by butyrate. In comparison, the earlier microarray reports [10] identified only 371 genes (285 genes down-and 86 up-regulated genes) for the IPA analysis.
Functional analysis identified the biological functions and/or diseases that were most significantly enriched in the dataset. When the functional category analysis of the genes was performed, genes from the datasets that were associated with biological functions and/or diseases in the Ingenuity Pathways Knowledge Base were considered for analysis. Fischer's exact test was used to calculate the P values. The top five molecular and cellular functions, as determined by P-value, are listed in Table 4. These five functional categories may represent the mechanisms underlying the essential biological effects of butyrate treatment, including cell morphology changes, cell cycle arrest, and apoptosis. The number of genes defined in each function category was greatly extended by RNAseq to include 2,257 genes involved in cell death and 2,322 genes involved in cellular growth and proliferation (Table 4).
We illustrated the functional changes induced by butyrate treatment by separately comparing the functional categories that were up-or down-regulated. Figure 2 shows the top fifteen functional categories that were significantly enriched in either upor down-regulated genes. Cell cycle; DNA replication, recombination, and repair; and RNA post-transcriptional modification were among the functional categories that were significantly impaired by butyrate. In contrast, cell death, cellular growth and proliferation, molecular transport, and cellular signaling categories were enhanced by butyrate.
Four canonical pathways (Cell cycle G2/M DNA damage checkpoints, purine metabolism, pyrimidine metabolism, and G1/ S checkpoint regulation) previously identified by the microarray experiment were also confirmed by the RNA-seq analysis. In addition, many other pathways were significantly impacted by butyrate treatment, including those directly related to cell cycle regulation, DNA replication, and cell cycle control of chromosomal replication; these finding were consistent with the observed phenotypic changes in cell cycle arrest and the blockage of DNA synthesis induced by butyrate. Signaling pathways, including NF-kB, IGF-1, p53, TGF-b, and apoptosis signaling, were also significantly induced by butyrate ( Figure 3 and Table S2).
Butyrate induced extensive deregulation of genes related to cell cycle progression IPA analysis identified 1,117 genes associated with cell cycle progression that were differentially regulated by butyrate (p-value: up to 5.86E 228 ) ( Table 4). These genes were involved in various checkpoint pathways and selected examples of these pathways were analyzed in further detail. A complete list of pathways is presented in the supplementary material (Table S3).
Regulation of the cell cycle: The G1/S checkpoint control is vital for normal cell division. Deregulation of the expression of checkpoint proteins can lead to apoptosis or tumorigenesis. This pathway highlights the key components of G1/S checkpoint regulation. Our data indicated that the G1/S checkpoint regulation pathway is one of the significantly down-regulated   Figure 4). It is very interesting to note that HDACs were among these deregulated genes; for example, HCACs 1, 4, 7, 9, and 10 were downregulated, while HDACs 2, 3, 5, 6, and 11 were significantly upregulated by butyrate. Regulation of DNA replication: cell cycle control of chromosomal replication is another canonical pathway closely related to cell cycle progression. The top functions of these pathways included DNA replication, recombination, and repair; cell cycle regulation; cellular assembly; and cellular organization. The stable propagation of genetic information requires that the entire genome of an organism be faithfully replicated only once in each cell cycle. Therefore, chromosomal DNA replication in eukaryotic cells entails a series of complex events that includes the recognition of origins, the firing of replication origins, the loading of DNA polymerases onto the origins, and the elongation of newly synthesized DNA. The initiation of DNA replication takes place only at specific loci on the chromosomal DNA, which are termed replication origins. The Origin Recognition Complex (ORC) includes six components (ORC1 to ORC6), which are specifically associated with replication origin throughout the cell cycle. ORC serves as a hallmark of the origins and is highly conserved. ORC1 is the largest subunit of the origin recognition complex and the association of ORC1 with chromatin appears to be the ratelimiting step in the assembly of a functional pre-replication complex [11]. Our data revealed that 23 genes from a total of 30 genes involved in this pathway were regulated by butyrate. These genes (such as CDC45, CDC6, CDC7, CDK, CDT1, CHEK2, DBF4, DNA Polymerase, MCM, ORC, ORC-CDC45-CDT1-MCM-RPA, ORC1, ORC2, ORC3, ORC4, ORC5, ORC6, which are the important components for the formation of the prereplication complex, as well as RC and RPA) were all significantly down-regulated ( Figure 5).
The canonical pathway of cell cycle regulation by BTG proteins may also play an important role in butyrate-induced cell cycle arrest. As shown in Figure 6, both BTG1 and BTG2 were upregulated by butyrate treatment. BTG1 expression reaches a maximum in the G0/G1 phases of the cell cycle and then begins to undergo down-regulation as cells progress through G1. BTG1 negatively regulates cell proliferation [12]. BTG2 proteins are anti-proliferation proteins involved in cell cycle regulation, growth arrest, and differentiation. The activation of BTGs may lead to the down-regulation of the cyclin E/CDK2 complex and other members of the cyclin family that are essential for the progression of the cell cycle from G1 to the S phase and that are responsible for the regulation of cyclin-dependent kinases. All of these differentially expressed genes and their functions are in agreement with our results and with the observed biological effects of butyrate.
Our data also demonstrated that cytokinesis was significantly down-regulated, with a p-value of 5.22E 209 and a Z-score of 22.781 (Table 5). A total of 48 genes in this pathway were downregulated, including Aurora kinases A, B and C. The KIF (kinesin superfamily of microtubule-associated motors) members, such as KIF 4A, C1, 20A, 23, were also significantly down-regulated (from 24.0 to 27.6 fold). These findings confirm the earlier microarray results that showed that butyrate induced changes in the expression of genes related to cytokinesis [5].
Transcription factors: Transcription factors are a group of proteins that bind to specific DNA sequences and control the transcription of genetic information from DNA to mRNA [13]. Transcription factors either promote (as activators) or block (as repressors) the recruitment of RNA polymerase to specific genes. Table 6 lists the major transcription factors identified by RNA-seq and IPA analysis that were involved in the regulatory effect of butyrate. Only the genes with a predicted activation state, either activated or inhibited, are listed.
TP53, one of the most important transcription factors, was found in the center of a down-regulated network in the microarray profiling of butyrate-induced regulation. However, microarrays failed to detect changes in TP53 gene expression. In the present experiment with RNA-seq, TP53 was clearly down-regulated by butyrate (,4 fold). TP53 targets 518 genes in the entire dataset of differentially expressed genes induced by butyrate (Table S4). In addition to TP53, butyrate also induced the expression of TP53BP1, TP53BP2, TP53I13, TP53INP1 (tumor protein p53 inducible nuclear protein 1), and TP53I11. TP53INP1 was upregulated (2.5 fold), while functionally-associated gene TP73 was up-regulated almost 24-fold. All of these changes in gene expression suggest a cell-cycle regulation network that may enhance cell cycle arrest.
The expression of non-coding RNAs (ncRNA) was disrupted by butyrate treatment: Deep RNA-seq also reveals a significant amount of information regarding ncRNA. There are 24 ncRNAs that are differentially expressed due to the butyrate treatment. Those ncRNAs belong to different types of ncRNAs, including snoRNA (small nucleolar RNA), snRNA (splicesomal RNA), and some miscRNAs (Table 7). Particularly, the expression of 10 snoRNAs (5 down-regulated and 5 up-regulated) was found to be disrupted by the butyrate treatment. snoRNAs are intermediatesized ncRNAs (60-300 bp). They are components of small nucleolar ribonucleoproteins (snoRNPs), which are complexes that are responsible for the modification and processing of ribosomal RNA [14]. More importantly, a large proportion of snoRNAs have been found to be further processed into smaller molecules, such as microRNAs (miRNAs) [15]. Surprisingly, only one miRNA with differential expression was detected. We suspect that the RNA purification protocols may exclude small RNAs. It is certainly interesting to follow-up this finding to look into the functionality of the disruption of the expression of ncRNA induced by butyrate.

Discussion
During the last few years, several publications have reported the use of HDAC inhibitors to study histone acetylation and gene regulation. An important question to be addressed by the study of histone modification is how modifications affect not only chromatin dynamics but also various processes (e.g., DNA replication, RNA transcription) along the DNA-template. These processes can be influenced by a number of post-translational modifications of histones, including acetylation, methylation, phosphorylation, and ubiquitination. These modifications may not act alone, but in concert and in a context-dependent manner to facilitate or repress chromatin-mediated processes [6].
Our previous studies [3,5,10,16] revealed that VFAs, especially butyrate, participate in metabolism, both as nutrients and as regulators of histone modification, thereby regulating the 'epigenomic code.' These findings implicate histone modifications induced by butyrate as determinants of bovine phenotype and in bovine ruminal development.
Epigenomics is an emerging area of scientific investigation that is confirming the complexity of the mechanisms used to determine the how, when, and where of gene expression in order to ensure the normal development, health, and homeostasis of the animal. The recently completed profiling of global gene expression used a high-density oligonucleotide microarray [5,10] to identify 450 genes in bovine kidney epithelial cells that were significantly regulated by sodium butyrate at a very stringent false discovery rate (FDR) of 0%. The functional category and pathway analyses of the microarray data revealed that four canonical pathways (cell cycles: G2/M DNA damage checkpoint, pyrimidine metabolism, G1/S checkpoint regulation, and purine metabolism) were significantly perturbed. The biologically relevant networks and pathways of these genes were also identified, including genes such as IGF2, TGFB1, TP53, E2F4, and CDC2, which were established as central to these networks. However, because they are restricted to probes designed to target the genes in a given species' genome, hybridization-based microarray technologies offer a limited ability to fully catalogue and quantify the diverse RNA molecules that are expressed from genomes over a wide range of levels [17], and they often fail to capture the full catalogue of transcripts and their variations.
The development of the next-generation sequencing (NGS) has provided novel tools for expression profiling and genome analysis [17,18,19]. As a vital step towards a comprehensive understanding of the molecular mechanism of butyrate-induced acetylation, as well as its biological effects, the present study was designed to utilize next-generation sequencing technology in order to provide a more complete characterization of the RNA transcripts of MDBK cells. This study also focused on the comparison between the control group (without butyrate treatment) and the cells    The significance value associated with a function in Global Analysis is a measure for how likely it is that genes from the dataset file under investigation participate in that function. The significance is expressed as a p-value, which is calculated using the right-tailed Fisher's Exact Test. doi:10.1371/journal.pone.0036940.g002 treated with 10 mM butyrate for 24 hours. With technical replicates (four lanes for controls and four lanes for butyratetreated samples), the samples were deep-sequenced, with an average of more than 67 million reads per sample, and the results were used to estimate the differences induced by butyrate treatment. Therefore, our results show a very reliable and detailed profiling of the changes in gene expression induced by butyrate. This study has generated comprehensive information on an experimental system that can be used in many functional genomics studies of bovine cells. To the best of our knowledge, this is the first study that has used NGS and IPA to identify the influences of butyrate on transcriptomic characterization in a normal bovine cell line. IPA analysis revealed that butyrate exerts a very broad range of effects on many biological pathways through its inhibitory action on HDACs in the MDBK cell line. Our NGS results, with a comparative transcriptomic profiling approach, extended far beyond the findings reported using microarray technologies [10,20]. The phenomenal number of genes we identified that fall within a broad range of functional categories appear to provide a very detailed molecular basis for the butyrate-induced biological effects.
The stable propagation of genetic information requires that the entire genome of an organism be faithfully replicated only once in each cell cycle. In eukaryotes, this replication is initiated at hundreds to thousands of replication origins distributed over the genome, each of which must be prohibited from re-initiating DNA replication within a single cell cycle [21]. Initiation of DNA replication is a two-step process: First, initiation proteins are assembled onto the replication origin in a stepwise fashion to develop a pre-replication complex. Second, the initiation complex is activated by protein kinases, resulting in the establishment of replication forks. This process is tightly regulated, such that initiation at a given replication origin occurs only once per cell cycle. In addition, initiation is down-regulated in response to agents that damage DNA or block DNA replication.
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 the cell-division cycle is accurately regulated to ensure that normal cell division is synchronous with cell growth and that the initiation of DNA synthesis (the S phase) is timed precisely to avoid inappropriate DNA amplification. The G1/S checkpoint control is vital for normal cell division and involves the key components that include cell cycle kinases, CDK4/6-cyclin D and CDK2-cyclin E, and the transcription complex composed of the retinoblastoma protein (Rb) and transcription factor E2F. The activation of E2F is necessary for the G1-S transition. In the present report, CDK4/6 and cyclins E and E2F were significantly down-regulated by butyrate-induced histone acetylation. In contrast, p21, a cell cycle inhibitor protein, was significantly up-regulated. All of these perturbations of gene expression in the G1/S cell cycle checkpoint pathways are consistent with the observed biological effects of butyrate, which induces cell cycle arrest at the G1/S boundary [3].
Butyrate is able to inhibit all class I HDACs. It also seems to affect many other epigenetic-related enzymes by regulating the expression of genes. The missing link is why this inhibition of enzymatic activities, in turn, regulates their own expression at the mRNA level. In this report, we found a vastly complicated depiction of the expression of HDACs induced by butyrate treatment. Whereas the expression of HCACs 7, 8, and 9 are down-regulated, HDACs 5 and 11 are up-regulated, and HDACs 1, 2, 4, and 6 are unchanged (Table S1). HDAC inhibitors that affect the expression of the HDACs themselves have been observed in mouse neural cells [22]. In that report, both TSA and SB indeed elevated the expression of HADC1, HDAC3, HDAC5, and HDAC6, whereas the mRNA levels for HDAC 2 and HDAC7 did not change. The mRNA levels of HDAC8 and HDAC10 were not detectable in these cells. The mechanism and biological relevance of HDAC inhibitors in the regulation of the 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.
Butyrate, as a histone deacetylase inhibitor, can also decrease histone methylation [23], suggesting an interplay between histone acetylation and histone methylation. An emerging possibility is that histone modifications can influence one another. In other words, there may be ''crosstalk among histone modification'' [24]. Consistently, KDM5B, a specific histone demethylase (H3trimethyl-K4), was significantly up-regulated by butyrate treatment (Table S1). However, JSRID2, which is directly related to histone methylation and responsible for maintaining the methylation level on histone H3 lysine 27 trimethylation (H3K27me3) [25], was also significantly up-regulated. JARID2 possesses an in vitro methyl-protective activity, stabilizing Polycom Repressive Complex 2 (PRC2)-catalyzed H3K27me3 by protecting it from the activity of H3K27 demathylases [26]. These data may indicate that different histone marks (modifications) are differentially regulated and that in turn, differentially regulated histone marks regulate different biological functions [27]. On the other side, a reversal of DNA methylation by butyrate has also recently been reported to occur by the regulation of DNA (cytosine-5-)methyltransferase 1 (DNMT1) through ERK signaling [28]. In this report, we found that three DNA methyltransferases (DNMTs), DNMT1, DNMT3A, and DNMT3B, were significantly down-regulated by the butyrate treatment (Table S1). While DNMT1 functions in the establishment and regulation of tissuespecific patterns of methylated cytosine residues, DNMT3A and DNMT3B function in the de novo methylation of DNA [29,30]. These DNMTs are regulated by several mechanisms in terms of their expression and catalytic activity. However, for the first time, our data directly indicated that histone modification has a role in the regulation of the expression of DNMTs, thereby affecting the level of DNA methylation.
The first clear evidence that a six-subunit ''origin recognition complex's'' (ORC) activity in mammalian cells is regulated by cell cycle-dependent changes in the affinity of the largest subunit (Orc1) for chromatin has been reported [11,21]. Evidence has since confirmed these findings and extended them to show that mammalian Orc1 is selectively ubiquitinated and phosphorylated during the S-to-M-phase transition, while ORC subunits 2 to 5, which constitute a stable core complex, remain tightly bound to chromatin throughout cell division [31]. In addition, a second mechanism prevents the assembly of a functional ORC until the completion of mitosis: the selective association of Orc1 with Cdk1 (Cdc2)/cyclin A during the G2/M phase of cell division. This association accounted for the appearance in M-phase cells with hyperphosphorylated Orc1 that was subsequently dephosphorylated during the M-to-G1 transition [32]. The rebinding of Orc1 to chromatin follows the same time course as the degradation of  cyclin B, suggesting that the exit from mitosis triggers Orc1 binding to chromatin. In fact, the inhibition of Cdk activity in metaphase cells resulted in the rapid binding of Orc1 to chromatin, and NGS profiling shows that all six subunits of ORC are down-regulated by butyrate-induced histone acetylation, adding yet another layer of regulation of ORC activities via the modified expression of those genes. In our previous microarray profiling [10], some of the components of this pathway were found to be perturbed by butyrate-induced gene regulation; however, ORC1 was the only one of the six ORC complex genes that was detected to be a down-regulated gene. In the present report, ORC1 is still the most significantly down-regulated gene, but the other ORC components (ORC2 to ORC6) are all also identified as down-regulated. This result certainly indicates the superb sensitivity of deep RNA sequencing. We also found significant up-regulation of both BTG1 and BTG2. The BTG family member-2 (BTG2) has antiproliferative activity, and the expression of BTG2 in cycling cells induces the accumulation of hypophosphorylated, growth-inhibitory forms of retinoblastoma protein (Rb) and leads to G1 arrest through the impairment of DNA synthesis. These up-regulated antiproliferation activities are strengthened by the extensive repression of cyclin-dependent kinase and cell cycle-related genes that are clearly associated with the cell growth arrest induced by butyrate.
Tumor protein p53 (TP53, a nuclear protein), transcription factor E2F4, and many other transcription factors, were deregulated by butyrate treatment in the present study. TP53 plays an essential role in the regulation of the cell cycle, specifically in the transition from G0 to G1. It is found in very low levels in normal cells; however, in a variety of transformed cell lines, it is expressed in high amounts and is believed to contribute to transformation and malignancy. P53 is a DNA-binding protein that contains DNA-binding, oligomerization, and transcription activation domains. P53 is postulated to bind as a tetramer to a p53-binding site and activate the expression of downstream genes that inhibit growth and/or invasion, thereby functioning as a tumor suppressor. P53 has been extensively studied for its function and involvement in butyrate-induced biological effects [33,34,35]. Butyrate efficiently suppresses the growth of WT p53-containing cells. It leads to a major G2/M arrest of cells in the presence of p53, while cells without wild-type p53 accumulate mainly in the G1 phase of the cell cycle. Apoptosis induction by butyrate is also greatly reduced in the absence of p53, suggesting that a p53 pathway mediates, in part, growth suppression by butyrate and that p53 status may be an important determinant of chemosensitivity to butyrate [36]. Our data also indicate that the TP53 genes may have different responses and different roles to play in normal and transformed cells. In our dataset, 518 genes were potential targets for TP53 regulation. Among these 518 genes, 238 genes showed expression directions consistent with the activation of TP53. However, one remaining question is why the expression of TP53 was down-regulated, even as its function was more active. As an extremely regulated gene, two major factors may contribute to this complexity of TP53. First, the expression of TP53 is subject to multiple regulations at transcriptional, post-transcriptional, and translational levels, with very complex expression patterns of alternative splicing, alternative promoter usage, and alternative translation. Secondly, the regulation of p53 function is extremely complex and occurs at many levels. Post-translational modifications of p53 (phosphorylation, methylation, acetylation, etc.) alter the functions of p53 (recognition of DNA sequences, interactions with transcription factors at promoters of target genes, etc.) [37]. Indeed, deep RNA-seq and IPA analysis revealed significant changes in the expression of genes related to the molecular function of protein post-translational modification (Figure 2). There are 333 genes related to the phosphorylation of proteins, 80 genes related to the tyrosine phosphorylation of proteins, and 106 genes related to the activation of protein kinase, which is upregulated by butyrate. The possibility exists that the modification of p53 is affected by butyrate, directly or indirectly. Clearly, more studies are still required to understand the exact roles that TP53 plays in butyrate-induced biological effects.
In conclusion, the acetylation of histone tails is essential for diverse cellular processes, such as DNA replication and cell cycle progression. Butyrate-induced histone hyper-acetylation, however, has some divergent activities, including the induction of cell cycle arrest, gene expression, and apoptosis [3,10]. The transcriptome characterization of bovine cells using RNAseq identified transcriptional control mechanisms via butyrate. Our results extended our knowledge of the regulatory effects of butyrate on gene expression and will undoubtedly provide insight into the molecular mechanisms of in vivo butyrate-induced epigenomic regulation.

Cell culture and butyrate treatment
Madin-Darby bovine kidney epithelial cells (MDBK, American Type Culture Collection, Manassas, VA., and Catalog No. CCL-22) were cultured in Eagle's minimal essential medium and supplemented with 5% fetal bovine serum (Invitrogen, Carlsbad, CA) in 25 cm 2 flasks, as described in our previous report [3]. At approximately 50% confluence, during the exponential phase, the cells were treated for 24 hours with 10 mM sodium butyrate (Calbiochem, San Diego, CA). A butyrate concentration of 10 mM was selected as it represents a physiologically relevant dose and has previously been successfully used to evoke desired changes in cell cycle dynamics [3]. Four replicate flasks of cells for both treatment and control groups (i.e., a total of 8 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 [5]. The RNA integrity was verified using an Agilent Bioanalyzer 2100 (Agilent, Palo Alto, CA, USA). High-quality RNA (RNA Integrity number or 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 [38]. Approximately 67.5 million reads per sample (mean 6 sd = 67,527,11168,330,388.6) were generated.

Data analysis and bioinformatics
Raw sequence reads were first checked using our quality control pipeline. Nucleotides of each raw read were scanned for low quality and trimmed using SolexaQA [39]. Trimmed reads were aligned to the bovine reference genome (Btau 4.0) using TopHat [40]. Each SAM output file per sample from TopHat alignment, along with the GTF file from ENSEMBL bovine genebuild v65.0, were used in the Cuffdiff program in the Cufflink package (v1.3.0) as input files [41] to test for 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 [42].
Differentially-expressed genes in the transcriptome were further analyzed using Gene Ontology (GO) analysis (GOseq) [43]. 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, IngenuityH Systems, and 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 [44,45,46,47,48].

Canonical pathway analysis of data sets
Analysis of canonical pathways identified the pathways from the IPA library of canonical pathways that were most significant to the data set. Genes from the data set that were associated with a canonical pathway in the Ingenuity Pathways Knowledge Base were considered for the analysis. The significance of the association between the data set and the canonical pathway was measured in two ways: 1) a ratio of the number of genes from the data set that map to the pathway divided by the total number of genes that map to the canonical pathway was displayed. 2) Fischer's exact test was used to calculate a p-value determining the probability that the association between the genes in the dataset and the canonical pathway was explained by chance alone.

Functional analysis of data sets
The Functional Analysis identified the biological functions and/ or diseases that were most significant to the data set. Genes from the datasets that were associated with biological functions and/or diseases in the Ingenuity Pathways Knowledge Base were considered for the analysis. Fischer's exact test was used to calculate a p-value determining the probability that each biological function and/or disease assigned to that data set was due to chance alone.

Pathways analysis and network generation
A data set containing gene identifiers and corresponding expression values was uploaded into in the application. Each gene identifier was mapped to its corresponding gene object in the Ingenuity Pathways Knowledge Base. These genes, called Focus Genes, were overlaid onto a global molecular network developed from information contained in the Ingenuity Pathways Knowledge Base. Networks of these Focus Genes were then algorithmically generated based on their connectivity.

Functional analysis of a network
The Functional Analysis of a network identified the biological functions and/or diseases that were most significant to the genes in the network. The network genes associated with biological functions and/or diseases in the Ingenuity Pathways Knowledge Base were considered for the analysis. Fischer's exact test was used to calculate a p-value determining the probability that each biological function and/or disease assigned to that network was due to chance alone.

Network/pathways graphical representation
A network pathway is a graphical representation of the molecular relationships between genes/gene products. Genes or gene products were represented as nodes, and the biological relationship between two nodes were represented as an edge (line). All edges were supported by at least 1 reference from the literature, from a textbook, or from canonical information stored in the Ingenuity Pathways Knowledge Base. The intensity of the node color indicated the degree of up-(red) or down-(green) regulation.