Quantification of Transcriptome Responses of the Rumen Epithelium to Butyrate Infusion using RNA-seq Technology.

Short-chain fatty acids (SCFAs), such as butyrate, produced by gut microorganisms, play a critical role in energy metabolism and physiology of ruminants as well as in human health. In this study, the temporal effect of elevated butyrate concentrations on the transcriptome of the rumen epithelium was quantified via serial biopsy sampling using RNA-seq technology. The mean number of genes transcribed in the rumen epithelial transcriptome was 17,323.63 ± 277.20 (±SD; N = 24) while the core transcriptome consisted of 15,025 genes. Collectively, 80 genes were identified as being significantly impacted by butyrate infusion across all time points sampled. Maximal transcriptional effect of butyrate on the rumen epithelium was observed at the 72-h infusion when the abundance of 58 genes was altered. The initial reaction of the rumen epithelium to elevated exogenous butyrate may represent a stress response as Gene Ontology (GO) terms identified were predominantly related to responses to bacteria and biotic stimuli. An algorithm for the reconstruction of accurate cellular networks (ARACNE) inferred regulatory gene networks with 113,738 direct interactions in the butyrate-epithelium interactome using a combined cutoff of an error tolerance (ɛ = 0.10) and a stringent P-value threshold of mutual information (5.0 × 10(-11)). Several regulatory networks were controlled by transcription factors, such as CREBBP and TTF2, which were regulated by butyrate. Our findings provide insight into the regulation of butyrate transport and metabolism in the rumen epithelium, which will guide our future efforts in exploiting potential beneficial effect of butyrate in animal well-being and human health.

Introduction

Short-chain fatty acids (SCFAs, also known as volatile fatty acids or VFAs) are produced in the gastrointestinal tract by microbial fermentation of carbohydrate substrates of both exogenous (dietary fiber) and endogenous origin (mucus). SCFAs contribute up to 75% of the total energy requirement in ruminants and ~10% for humans.1 A broad range of biological functions have been ascribed to SCFAs, specifically butyrate, in addition to their value in energy and nutrition metabolism.1,2 While as a preferred energy source for ruminal epithelial cells, butyrate plays a critical role in the development of the rumen epithelium.3 As an indirect factor, butyrate is able to affect blood flow in the gastrointestinal (GI) tract as well as insulin and glucagon secretion.1 Furthermore, butyrate reinforces intestinal barriers and modulates motility and visceral sensitivity of the intestine.4 Butyrate modulates cell differentiation and motility, induces apoptosis, and inhibits cell proliferation in vitro5,6 in contrast to the observed effects on cell proliferation documented in vivo7 Butyrate is known to be anti-inflammatory and anti-tumorigenic.4 Recently, butyrate has been shown to reduce bacterial translocation and internalization resulting from metabolic stress in human colon-derived cells.8 Beneficial effects of butyrate on human health and animal well-being have been extensively reviewed.4,9

Steady-state butyrate concentration in the rumen represents a balance between butyrate biosynthesis by ruminal microorganisms and absorption and utilization by the rumen epithelium. Among factors that affect butyrate output of the gut microbial ecosystem, such as redox potential and pH, gut turnover and motility, and mucosal transport,10 substrates in the gut play an essential role in butyrate production, possibly due to their influence on microbial population composition and dynamics.10–12 The community-level response of rumen microorganisms to exogenous butyrate has recently been characterized.12 Dietary supplements, such as butyrogenic resistant starch, and probiotic administration of butyrate-producing bacteria have been exploited in human medicine to boost gut butyrate concentration for its potential health benefit.13–15 In ruminants, butyrate produced by gut microbes is rapidly removed by the rumen or gut epithelium, and subsequently metabolized to ketone bodies (acetoacetate and β-hydroxybutyrate) or oxidized to CO2 prior to being released into circulation such that only 25% of the butyrate enters the systemic circulation.16 While substantial amounts of butyrate are absorbed across the rumen epithelium via passive diffusion, either trans- or para-cellularly,17 other forms of butyrate transport may exist. Indeed, a butyrate transporter, SLC5A8, has been identified.18 Much of our knowledge relating to the transcriptional effects of butyrate is derived from in vitro studies.2 However, butyrate has long been known to have an opposite effect in vivo, especially on cell proliferation. In the present study, we quantified the changes in the rumen epithelial transcriptome of dairy cows in response to butyrate infusion using RNA-seq technology and bioinformatic tools.

Materials and Methods

Animal experiment

Four ruminally-cannulated Holstein cows used in this study were previously described.12 Briefly, the cows in their mid-lactation were fed ad libitum standard lactation rations with 50% corn silage and 50% concentrate on a dry matter basis. All animal care and handling were conducted according to the guidelines approved by the USDA Beltsville Area Institutional Animal Care (Protocol # 09-008). Cows were moved to a tie stall barn for adaptation and acclimation 7 days prior to the experiment. Ruminal infusion of butyrate was initiated immediately following 0h sampling (baseline controls) and thereafter continued for 168 h at a rate of 5.0 L/day of a 2.5 M solution (representing >10% of daily anticipated metabolizable energy intake to support lactation) in a buffered saliva solution (pH 7.0; 3.8% KHCO3, 7.3% NaHCO3) as a continuous infusion. Following the 168 h infusion, infusion was stopped and cows were maintained on the basal lactation ration for an additional 168 h for sampling. Rumen epithelial samples were serially collected via biopsy through rumen fistulae at 0, 24, 72, and 168 h of infusion, and 24 and 168 h post infusion (post 24 h and post 168 h, respectively). Ruminal pH was monitored using a standard pH meter and recorded at every sampling. Rumen epithelial samples were saved in RNAlater and stored in −80C until RNA extraction. Additionally, a bovine epithelial cell line treated with 10 mM of sodium butyrate (N = 4) and PBS, respectively, for 24 h was used as an in vitro model for comparison.2

RNA extraction and sequencing using RNA-seq

Total RNA from 24 rumen epithelial samples and 8 cell line samples was extracted using Trizol (Invitrogen, Carlsbad, CA, USA) followed by DNase digestion and Qiagen RNeasy column purification (Qiagen, Valencia, CA, USA) as previously described.19,20 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 pooled based on their respective sample-specific 6-bp adaptors and sequenced at 50 bp/sequence read using an Illumina HiSeq 2000sequencer as described previously.21 Approximately 25.7 million sequences per sample (Mean ± SD = 25,746,641 ± 7,637,237.4) were generated.

Data analysis and bioinformatics

Raw sequence reads were first checked using our quality control pipeline. Nucleotides of each raw sequence read were scanned for low quality and trimmed using SolexaQA.22 Trimming removed ~24.9% of the raw reads; and mean length of resultant input reads were 48.96 bp (±0.46, SD). Input reads were then aligned to the bovine reference genome (UMD3.1) using TopHat.23 Approximately 96.68% of input reads were aligned to the bovine genome whereas 71.58% of input reads were uniquely mapped to the gene space in the bovine genome. Each SAM output file 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 files24 to test for differential gene expression. Mapped reads were normalized based on upper-quartile normalization method. Cuffdiff models the variance in fragment counts across replicates using the negative binomial distribution as described.25

Differentially expressed genes in the transcriptome identified were further analyzed using Gene Ontology (GO) analysis (https://github.com/tanghaibao/goatools). Enrichment of certain GO terms was determined based on Fisher’s exact test. A multiple correction control (permutation to control false discovery rate, FDR) was implemented to set up the threshold to obtain the lists of significantly over-represented GO terms.

Co-regulatory gene networks were analyzed using Algorithm for the Reconstruction of Accurate Cellular Networks (ARACNE).26 Raw mapped sequence reads from 32 samples, including 24 rumen epithelial samples and 8 samples from a bovine epithelial cell line were first normalized (per million mapped reads) and log2 transformed. Only 15,469 genes with mean mapped reads greater than 10.0 across all 32 samples were used for subsequent analysis. Relative expression values of these genes were used as input files to infer global regulatory networks.

Results

A snapshot of the rumen epithelial transcriptome

Collectively, a total of 20,371 genes were detected at least once in at least one of the 24 rumen epithelial samples. The number of genes expressed per sample was 17,323.63 ± 277.20 (Mean ± SD). The number of genes with mean hits ≥ 5.0 were 15,015 whereas the core transcriptome of the rumen epithelium consisted of 15,025 genes, ie, genes represented by at least one sequence hit in each of the 24 samples tested.

According to the transcript abundance classification described previously,21 approximately 70.78% of genes transcribed in the rumen epithelium can be classified into “very rare” with a relative abundance ≤15 molecules per cell (Table 1), followed by “not expressed” (0 molecule per cell) at 17.24% and “rare” (16–99 molecules per cell) at 10.41%. “ Moderately abundant” genes accounted for 1.32% while the abundant transcripts with ≥500 molecules per cell accounted for only 0.25% of genes detected in the transcriptome. Unlike what was observed in the bovine muscle tissue where the transcriptome was dominated by a handful of genes (Baldwin et al 2012, personal communication), the most abundant gene in the rumen epithelial transcriptome was a mitochondrial gene, cytochrome c oxidase subunit I, represented by ~6.23% of sequence reads. The 20 most abundant genes in the rumen epithelial transcriptome, which accounted for ~23.8% of sequences, included 6 genes of mitochondrial origin, such as cytochrome b, and 3 keratin genes, such as keratin, type I cytoskeletal 17 (Ensembl ID: ENSBTAG00000006806) and keratin, type II cytoskeletal 59 kDa, component IV ( ENSBTAG00000039425). This is in contrast to the bovine small intestine transcriptome,21 in which the 20 most abundant genes accounted for ~13.3% of sequences. Among the 20 most abundant transcripts, only 3 genes were shared in the transcriptome of both the small intestine and the rumen epithelium, including cytochrome c oxidase subunit I (ENSBTAG00000043561), elongation factor 1-alpha 1 (ENSBTAG00000014534), and a novel protein-coding gene (ENSBTAG00000034185).

Table 1

Transcript abundance in the bovine rumen epithelium.

Transcript category	Rumen epithelium	Small intestine
Not expressed	17.24%	15.23%
Very rare	70.78%	68.27%
Rare	10.41%	14.91%
Moderately abundant	1.32%	1.41%
Abundant	0.25%	0.18%

Genes significantly impacted by butyrate

A total of 80 genes were found to be affected in the rumen epithelium by butyrate infusion at ≥ one time point at a stringent cutoff, FDR < 5% (Table 2). Compared to pre-infusion at 0h, 24 h infusion affected 15 genes, including down-regulation of solute carrier family 5 (iodide transporter), member 8 (SLC5A8), leukocyte surface antigen (CD47), lipopolysaccharide binding protein (LBP), immunoglobulin lambda-like polypeptide 1 (IGLL1), and lactoperoxidase (LPO). The maximal effect of butyrate was observed at 72 h infusion when 59 genes were influenced, including lymphocyte antigen 6D (LY6D), keratin 75 (KRT75), and tubulin beta-3 chain (TBB3). While the impact on some genes, such as a novel protein (ENSBTAG00000039035), which was strongly induced by butyrate at only the 24 h infusion sampling, was seemingly transient, transcriptional responses of butyrate on many genes were long-lasting. The ruminal butyrate concentration returned to its pre- infusion basal level 168 h after the infusion withdrawal (Fig. 1). However, as many as 13 genes remained significantly affected at 168 h post-infusion, compared to 0h, and these included cadherin 26 (CDH26), desmoglein 1 (DSG1), and kallikrein-related peptidase 10 (KLK10). 53 of the 80 differentially expressed genes were also significantly impacted by butyrate in the in vitro model. However, the remaining 27 genes were only significantly regulated by butyrate in the rumen epithelium, such as dual oxidase maturation factors 1 and 2 (DUOXA1 and DUOXA2) and keratin 36 and keratin 75 (KRT 36 and KRT75). Additionally, of the 80 genes affected by butyrate, a majority were down-regulated by butyrate treatment, which is consistent with the effect of butyrate observed in bovine epithelial cells in vitro (Wu et al 2012, personal communication). Eight genes up-regulated by butyrate at ≥one time point were heat shock 70 kDa protein 6 (HSPA6, 24 h), transglutaminase 2 (TGM2, 72 h and 168 h), alcohol dehydrogenase 1C (class I), gamma polypeptide (ADH1C, 72 h and post 168 h), small nucleolar RNA SNORA31 (SNORA31, 168 h), farnesyl pyrophosphate synthase (FPPS, post 124 h), keratin 36 (KRT36, post 168 h), ornithine decarboxylase (ODC, post 168 h), and peptidoglycan recognition protein 2 (PGLYRP2, post 168 h).

Table 2

Genes significantly impacted by butyrate in the rumen epithelium detected by RNA-seq technology.

Gene_id	Symbol	Gene name	Locus	24 h	72 h	168 h	Post 24 h	Post 168 h
ENSBTAG00000001100	IL22RA1	Interleukin 22 receptor, alpha 1	2:129351697–129378160	0.57	0.37**	0.41	0.83	0.92
ENSBTAG00000001497	MRAS	Muscle RAS oncogene homolog	1:131765942–131834001	0.96	0.76	0.86	0.81	0.64***
ENSBTAG00000002088		Uncharacterized protein	11:46575177–46584984	0.51***	0.93	0.97	0.96	1.25
ENSBTAG00000002404	TMX2	Thioredoxin-related transmembrane protein 2	15:82363843–82372527	0.87	0.73**	0.82**	0.98	0.95
ENSBTAG00000003192	TBC1D16	TBC1 domain family, member 16	19:53181378–53245548	0.70	0.56**	0.79	0.85	0.71
ENSBTAG00000003585	CD47	CD47 molecule	1:53103995–53169038	0.77***	0.73***	0.89	1.06	1.21
ENSBTAG00000003948 FPPS		Farnesyl pyrophosphate synthase	3:15356046–15374937	1.01	0.85	0.88	1.27**	1.46***
ENSBTAG00000004208	PLEKHG6	Pleckstrin homology domain containing, family G member 6	5:104415804–104431128	0.73	0.36**	0.41*	0.72	0.82
ENSBTAG00000004256	ODC	Ornithine decarboxylase	11:87178097–87182503	0.74	0.56	0.65	1.29	2.96**
ENSBTAG00000004322	FOS	FBJ murine osteosarcoma viral oncogene homolog	10:86883738–86887169	0.80	0.43**	0.55	0.89	0.75
ENSBTAG00000004884	SLC4A5	Solute carrier family 4, member 5	11:10283281–10403977	0.57	0.39**	0.48	0.58	0.67
ENSBTAG00000004971		GRAM domain containing 1C	1:58963949–59063256	0.68***	0.68***	0.82**	1.13	0.72***
ENSBTAG00000005397	CD46	Membrane cofactor protein	16:77483016–77522665	0.82	0.74**	0.85	0.88	0.81
ENSBTAG00000005466	C8orf58	Chromosome 8 open reading frame 58	8:70406943–70412090	0.74	0.70**	0.84	0.75	0.75
ENSBTAG00000005580	SZT2	Seizure threshold 2 homolog (mouse)	3:103160271–103201060	0.91	0.79*	0.95	0.92	0.92
ENSBTAG00000006464		ADP-ribosylation factor-binding protein GGA1	5:109960752–109979896	0.91	0.70***	0.82***	0.86	0.87
ENSBTAG00000006984	CD55	CD55 molecule, decay accelerating factor for complement	16:5105576–5124229	0.84	0.68**	0.84	0.95	1.08
ENSBTAG00000007399	LAMP2	Lysosomal-associated membrane protein 2	X:4846777–4879916	0.88	0.78*	0.82	0.83	0.80
ENSBTAG00000008442	NSF	N-ethylmaleimide-sensitive factor	19:46081985–46228565	0.74	0.57**	0.73	0.78	0.82
ENSBTAG00000008700	CRCT1	Cysteine-rich C-terminal 1	3:18222997–18224799	0.55	0.40*	0.39*	0.91	0.68
ENSBTAG00000009575	C6orf170	Chromosome 6 open reading frame 170	9:30195128–30400435	0.86	0.76	0.75	0.68*	0.69
ENSBTAG00000009995	CPEB4	Cytoplasmic polyadenylation element binding protein 4	20:5512160–5584540	0.75	0.67*	0.79	0.84	0.75
ENSBTAG00000010670	VSIG8	V-set and immunoglobulin domain containing 8	3:9942946–9950485	0.58	0.27***	0.34***	0.62	0.72
ENSBTAG00000011482	SLC43A2	Solute carrier family 43, member 2	19:23261579–23301008	0.80	0.47*	0.49**	0.62	0.82
ENSBTAG00000011525	SLC5A8	Solute carrier family 5 (iodide transporter), member 8	5:65387536–65453818	0.32*	0.33*	0.38	0.90	1.49
ENSBTAG00000012280	PGLYRP2	Peptidoglycan recognition protein 2	7:8702089–8707156	1.09	0.64	0.82	2.27	7.00***
ENSBTAG00000012780	LPO	Lactoperoxidase	19:9441748–9471498	0.37***	0.72	0.76	1.14	0.72
ENSBTAG00000013205		Interleukin 1 receptor accessory protein	1:77198181–77343645	0.85	0.64***	0.75	0.87	1.00
ENSBTAG00000013507	SPRR3	Small proline-rich protein 3	3:17796521–17798340	0.61	0.46**	0.66	0.89	1.41
ENSBTAG00000013578		Uncharacterized protein	3:32269407–32284496	0.23***	0.61	0.62	0.88	0.73
ENSBTAG00000013831	DSG1	Desmoglein 1	24:26089488–26133405	0.64	0.26***	0.21***	0.35**	0.34**
ENSBTAG00000013999	FAM20B	Family with sequence similarity 20, member B	16:61810680–61852109	0.82	0.68***	0.83	0.90	0.80
ENSBTAG00000014289	ZFYVE9	Zinc finger, FYVE domain containing 9	3:94600423–94713143	0.79	0.67*	0.75	0.80	0.79
ENSBTAG00000014367	PRKX	Protein kinase, X-linked	X:138715323–138747777	0.87	0.70***	0.86	0.88	0.87
ENSBTAG00000014945	IQSEC2	IQ motif and Sec7 domain 2	X:96082196–96158317	0.82	0.60***	0.76	0.80	0.83
ENSBTAG00000014984	RPRD2	Regulation of nuclear pre-mRNA domain containing 2	3:20254978–20345588	0.93	0.80	0.93	0.84	0.72*
ENSBTAG00000015129	KLK10	Kallikrein-related peptidase 10	18:57448224–57454197	0.58	0.45*	0.36**	0.68	0.35**
ENSBTAG00000015796	DNAJC5	DnaJ (Hsp40) homolog, subfamily C, member 5	13:54382946–54408528	0.85	0.66***	0.78	0.96	1.05
ENSBTAG00000015915	SBSN	Suprabasin	18:46358552–46363200	0.66	0.24***	0.16***	0.32*	0.40
ENSBTAG00000016121	KRT75	Keratin 75	5:27630614–27640849	0.81	0.45*	0.47	0.59	0.77
ENSBTAG00000016208	TGM2	Transglutaminase 2	13:67663137–67697607	1.40	2.05*	2.23*	0.94	1.62
ENSBTAG00000016239	DUOXA2	Dual oxidase maturation factor 2	10:65525529–65540482	0.50	0.48*	0.55	1.11	1.22
ENSBTAG00000016242	DUOXA1	Dual oxidase maturation factor 1	10:65525529–65540482	0.50	0.48*	0.55	1.11	1.22
ENSBTAG00000016864	LBP	Lipopolysaccharide binding protein	13:67874473–67910095	0.47*	0.54	0.63	0.86	0.92
ENSBTAG00000016968	CHST4	Carbohydrate (N-acetylglucosamine 6-O) sulfotransferase 4	18:39665072–39671884	0.38*	0.39*	0.55	0.94	0.71
ENSBTAG00000018252	ARRDC4	Arrestin domain containing 4	21:8924762–8939028	0.49*	0.46**	0.41**	0.96	0.87
ENSBTAG00000018937		Heme oxygenase (decycling) 2	25:3621446–3649596	0.88	0.72***	0.81*	0.90	0.94
ENSBTAG00000019107	GAS7	Growth arrest-specific 7	19:29670743–29876650	0.96	0.69*	0.87	0.86	0.81
ENSBTAG00000019822	TPPP3	Tubulin polymerization-promoting protein family member 3	18:35121242–35124820	1.73	0.29	0.20*	0.69	0.71
ENSBTAG00000019876	DSC1	Desmocollin 1	24:26283414–26321590	1.07	0.78	0.66*	0.63	0.76
ENSBTAG00000020261	CDH26	Cadherin 26	13:57056796–57091107	0.21	0.20	0.16	0.20	0.13*
ENSBTAG00000020384	WDR47	WD repeat domain 47	3:34446313–34506800	0.82	0.57**	0.62**	0.75	0.75
ENSBTAG00000023026	SERPINB2	Serpin peptidase inhibitor, clade B (ovalbumin), member 2	24:12747968–12765314	1.10	0.44**	0.43**	0.64	0.99
ENSBTAG00000023730	TBB3	Tubulin beta-3 chain	18:14761259–14769914	0.79	0.49*	0.51	0.69	0.81
ENSBTAG00000025531	PSORS1C2	Psoriasis susceptibility 1 candidate 2	23:27793272–27794712	0.50	0.49*	0.63	1.19	0.97
ENSBTAG00000026585		Elongin A	2:129738382–129752889	0.87	0.70*	0.81	1.00	1.01
ENSBTAG00000027225	LAP	Lingual antimicrobial peptide	27:5124122–5125987	0.35***	0.36***	0.38**	0.55	0.76
ENSBTAG00000030301	MAP7D1	MAP7 domain containing 1	3:110233981–110243621	0.70	0.48*	0.62	0.83	1.00
ENSBTAG00000030484	KLK12	Kallikrein-related peptidase 12	18:57461616–57465567	0.53	0.47	0.29**	0.60	0.24**
ENSBTAG00000031160	IGLL1	Immunoglobulin lambda-like polypeptide 1	17:72826679–73152862	0.61**	0.80	1.10	1.12	0.77
ENSBTAG00000031551	PRSS53	Protease, serine, 53	25:27450355–27455369	0.46*	0.40***	0.38**	0.67	0.76
ENSBTAG00000033107	OSMR	Oncostatin M receptor	20:35521409–35588186	0.71	0.56**	0.69	0.93	1.08
ENSBTAG00000033398	HKDC1	Hexokinase domain containing 1	28:25752420–25791456	0.63	0.31**	0.30**	0.48	0.77
ENSBTAG00000033669	IL17RC	Interleukin 17 receptor C	22:16873473–16885025	0.91	0.70*	0.80	0.96	1.16
ENSBTAG00000034189	ADH1C	Alcohol dehydrogenase 1C (class I), gamma polypeptide	6:26656287–26672428	0.86	4.57*	3.46	0.94	5.55***
ENSBTAG00000034498	LY6D	Lymphocyte antigen 6D	14:2801382–2803020	0.56	0.19*	0.37	0.34	0.52
ENSBTAG00000035171	SERPINB5	Serpin peptidase inhibitor, clade B (ovalbumin), member 5	24:62319471–62329373	0.89	0.69	0.52	0.30**	0.43
ENSBTAG00000038540		Uncharacterized protein	26:25060218–25073503	0.97	0.66*	0.72	0.80	0.81
ENSBTAG00000038878	Aox4	Aldehyde oxidase 4	2:89699123–89771871	0.68	0.46	0.32**	0.16***	0.44
ENSBTAG00000039035		Uncharacterized protein	3:8027844–8029776	13.40***	0.98	1.10	2.05	0.85
ENSBTAG00000040356	SPIRE2	Spire homolog 2 (Drosophila)	18:14694451–14723123	0.78	0.65**	0.80	0.90	0.90
ENSBTAG00000040394	PRSS27	Protease, serine 27	25:2126943–2133828	0.46*	0.50	0.50	0.62	0.63
ENSBTAG00000040507	ARHGEF18	Rho/Rac guanine nucleotide exchange factor (GEF) 18	7:17513880–17565624	0.74	0.58*	0.75	0.94	0.90
ENSBTAG00000043553	GPX3	Glutathione peroxidase 3 (plasma)	7:64286947–64295116	0.51	0.49*	0.52	0.71	0.67
ENSBTAG00000044673	SNORA31	Small nucleolar RNA SNORA31	4:15035753–15183208	0.99	2.17	2.70**	0.98	1.10
ENSBTAG00000046124	KRT36	Keratin 36	19:42325913–42329903	0.91	2.96	1.34	2.37	4.41*
ENSBTAG00000046158	CFB	Complement factor B	23:27225269–27231292	0.61	0.50*	0.62	1.00	1.83
ENSBTAG00000046482		Uncharacterized protein	10:26113953–26114424	0.35*	0.34*	0.52	0.69	0.70
ENSBTAG00000046917	GPR22	G protein-coupled receptor 22	4:48576431–48851596	0.83	0.74**	0.76***	0.92	0.82***
ENSBTAG00000047345		Uncharacterized protein	18:52984271–53016520	0.86	0.67*	0.81	0.98	1.05

Notes: Numbers denote fold changes (treated/control at 0 h, N = 4).

FDR < 5%;

FDR < 1%;

FDR < 0.1%.

Figure 1

Ruminal concentrations (mM) of short-chain fatty acids in response to butyrate infusion.

Notes: 0 h = baseline, immediately prior to initiation of continuous infusion with butyrate. 24 h, 72 h, and 168 h = 24, 72, and 168 h after continuous infusion, respectively. Post 24 h and post 168 h = 24 h and 168 h after infusion withdrawal, respectively. Red: acetate; Blue: butyrate; Green: propionate. Error bars = standard deviation, N = 4.

Gene Ontology (GO) analysis

GO terms significantly enriched in the rumen epithelium in response to butyrate infusion were identified using GOATOOLs. At 24 h infusion, five biological processes were significantly enriched at a Bonferroni- corrected P value < 0.05 (Table 3A), including responses to bacteria and biotic stimuli. When all 80 genes significantly impacted by butyrate at one or more time points were analyzed, biological processes such as responses to bacteria and other organisms were enriched, consistent with the 24 h observation. In contrast, cellular components related to “organelle” were significantly depleted by butyrate infusion when assessed relative to all 15,015 genes expressed in the rumen epithelium with mean mapped reads ≥5.0 (Table 3B).

Table 3

Gene Ontology (GO) terms significantly impacted in the rumen epithelium by butyrate. (A) At 24 h infusion. (B) Collectively in all time points.

GO_id	Status	GO description	Ratio	P_unadjusted	P_bonferroni	FDR
(A)
GO:0009617	E	Response to bacterium	4/67	4.85E-07	1.72E-04	0.0000
GO:0051707	E	Response to other organism	4/137	8.52E-06	3.02E-03	0.0100
GO:0042742	E	Defense response to bacterium	3/56	2.19E-05	7.77E-03	0.0270
GO:0051704	E	Multi-organism process	4/208	4.41E-05	1.56E-02	0.0390
GO:0009607	E	Response to biotic stimulus	4/219	5.39E-05	1.91E-02	0.0430
GO:0050729	E	Positive regulation of inflammatory response	2/25	0.000278	9.85E-02	0.1220
(B)
GO:0043231	P	Intracellular membrane-bounded organelle	8/4872	3.47E-06	0.0029	0.0030
GO:0043227	P	Membrane-bounded organelle	9/4877	1.89E-05	0.0158	0.0120
GO:0043229	P	Intracellular organelle	12/5577	2.07E-05	0.0174	0.0130
GO:0009617	E	Response to bacterium	5/68	3.03E-05	0.0254	0.0550
GO:0043226	P	Organelle	13/5583	6.13E-05	0.0514	0.0740
GO:0004252	E	Serine-type endopeptidase activity	5/81	7.06E-05	0.0593	0.0780
GO:0051707	E	Response to other organism	6/140	1.01E-04	0.0845	0.1540

Notes:

Ratio = the number of significantly regulated genes that can be assigned to this GO term/all genes assigned to this GO term [among 15,015 genes expressed in the rumen epithelium with a mean hit count ≥ 5.0]. Under the status, “E” (enriched) = GO term is significantly higher in the study group than those in the population while “P” (depleted) = GO terms significantly lower in the study group than in the population. The number of permutation = 1000.

Regulatory gene networks

Using a combined stringent cutoff of an error tolerance ɛ = 0.10 and a P-value threshold of mutual information (MI) at 5.0 × 10−11, ARACNE inferred global gene networks with 113,738 direct interactions in the butyrateepithelium interactome. Top 67 genes controlling >5% of all direct interactions and can be considered as hubs or master regulators (Supplementary File 1). Among these genes, 46 were impacted by butyrate in in vivo and/or in vitro models. Of the 46, 7 were known transcription factors or translation regulators, including CREB binding protein (CREBBP), eukaryotic translation termination factor 1 (ETF1), hypoxia inducible factor 1, alpha (HIF1A), lysine (K)-specific demethylase 5A (KDM5A), regulatory factor X, 2 (RFX2), transcription elongation factor A (3TCEA3), and transcription termination factor, RNA polymerase II (TTF2). The master regulator with the greatest number of direct interactions (or the first neighbors) was ETF1 (ENSBTAG00000011415). The product of this gene acts as an omnipotent protein release factor (RF) and directs the termination of nascent peptide synthesis (translation) in response to termination codons. ETF1 was down-regulated by butyrate in vitro (P value = 3.97 × 10−8, FDR = 1.00 × 10−5). ARACNE inferred 123 direct interactions (1st neighbors) and 6,195 2nd neighbors for ETF1. Among the indirect interactions (2nd neighbors), 1052 of the 1472 genes involved were significantly impacted by butyrate, including a total of 111 transcription factors in vitro. The predominant GO terms identified as enriched among its 1st neighbors relate to the cellular component organelle (GO:0043226), including membrane-bounded organelle (GO:0043227) and intracellular organelle (GO:0043229). CREBBP (ENSBTAG00000026403) encodes a protein with intrinsic histone acetyltransferase activity, which acetylates both histone and non-histone proteins. CREBBP acts as a transcriptional coactivator of RNA polymerase II-mediated transcription and is known to play a critical role in cell growth and human diseases. This gene was significantly down-regulated by butyrate infusion and had 87 1st neighbors (Fig. 2). An additional 189 genes were involved in 3,327 of indirect interactions of CREBBP (2nd neighbors). Of these CREBBP 2nd neighbors, 140 were regulated by butyrate including 19 transcription factors. Not coincidently, the GO terms predominantly enriched in the direct interaction of CREBBP were related to biological processes such as histone modification (GO:0016570) and covalent chromatin modification (GO0016569), cellular components such as histone acetyltransferase complex (GO:0000123) and NuA4 histone acetyltransferase complex (GO:0035267), and molecular functions such as histone acetyltransferase activity (GO:0004402) and nucleic acid binding (GO:0003676).

Figure 2

A regulatory gene network controlled by CREB binding protein (CREBBP).

Notes: The network was inferred using ARACNE at a combined stringent cutoff of an error tolerance ɛ = 0.10 and a P-value threshold of mutual information (MI) at 5.0 × 10−11. CREBBP had 87 direct interactions (the first neighbors) and 4139 indirect interactions (the second neighbors, not shown). The expression of the CREBBP gene, a transcription factor, at the mRNA level was significantly regulated by butyrate. The color cycle with gene symbols represented genes significantly regulated by butyrate.

Discussion

Multiple effects of butyrate on local (intestinal) and systemic physiological function are presumably due to its actions as a potent regulator of gene expression.4–6 Butyrate metabolism, both biosynthesis by gut microorganisms and utilization in the gut, is tightly regulated and thus, imbalance could play a critical role in the pathology of humans and animals. Typically, increased production of SCFA, especially butyrate, results in lower luminal pH, which, in turn, creates favorable niche for and stimulates proliferation of butyrate-producing bacteria.12 In modern intensive production systems used for ruminants, especially during physiologically critical stages, such as weaning and the transition or periparturient period in dairy cattle, high concentrate rations are fed which can result in abrupt increases in SCFA production. This circumstance often results in production in excess of the utilization capacity, subsequently leading to the development of ruminal acidosis, a prominent digestive disorder with a significant economical impact.27,28 Individual variation in the capacity for uptake of acetate and butyrate among animals fed the same rations can explain the degree of acidosis observed.29 In fact, considerable effort has focused on development of practical nutritional intervention strategies for use during transition periods to enhance rumen epithelial surface area prior to the introduction of high VFA producing rations to facilitate better absorptive capacity and mitigate negative effects of acidosis.30 To date, these strategies have yielded mixed results with respect to increasing papillae surface area and total absorptive capacity (Andersen, personal communication). The development of a better understanding of the regulation of both rumen epithelial growth promotion and VFA transport capacity remains a concern to the industry.29 In this study, we assessed the rumen epithelial transcriptome dynamics when butyrate concentration was elevated due to direct ruminal infusion with buffered butyric acid in order to assess transcriptional effects of butyrate and its potential role in regulating butyrate transport in cattle.

The majority of our knowledge on regulatory impact of butyrate on global gene expression is derived from in vitro studies and observations.6,31,32 However, caution in interpretation of in vitro data and application of knowledge gained in vitro to in vivo models is warranted. SCFAs, such as butyrate, are known to promote rumen development and stimulate the proliferation of rumen epithelial cells in vivo.33 However, butyrate inhibits the proliferation of epithelial cells of the large intestine, rumen and kidney by down- regulating genes controlling cell proliferation in vitro.2,32 Butyrate also induces apoptosis and differentiation of tumor cells.34,35 Moreover, an in vitro study of epithelial cells of different origin (rat small intestine vs. human colon) has demonstrated that the cell type affects butyrate uptake characteristics.36 In fact, opposite effects of butyrate are observed for many cellular processes, such as cell proliferation and division, between in vitro and in vivo models and are clearly reflected in transcriptome characteristics.33,37 A number of genes related to cell proliferation and cell cycle progress were significantly down-regulated by a 24-h 10 mM butyrate incubation of established bovine rumen epithelial cells in long-term culturing (Wu et al 2012, personal communication). In contrast, expression of these same genes in the rumen epithelium in the present data set was not altered despite a 2-fold increase from 19.5 mM to 38.5 mM in intraruminal butyrate concentration (Fig. 1). Interestingly, the abundance of a butyrate transporter, SLC5A8, in the rumen epithelial transcriptome was significantly reduced concomitant with intra-ruminal butyrate concentration increases. However, SLC5A8 expression at the mRNA level was significantly increased ~21 fold by butyrate exposure in vitro (Fig. 3). This apparent opposite effect of butyrate on the expression of its transporter is suggestive of alternative regulatory mechanisms relating to butyrate uptake control and transport by the intact rumen epithelium and cells in culture. Moreover, cellular butyrate metabolism may be different between in vivo and in vitro models due to changes in the rate of removal of end product as well as changes between cell-cell interactions and micro-environments present in vivo, but disrupted in vitro. To this point, epithelial metabolism of butyrate, especially the pathways leading to ketogenesis, helps to maintain a butyrate concentration gradient in vivo, which in turn facilitates butyrate uptake and affects butyrate intracellular concentrations.30 Rate-limiting enzymes in the ruminal ketogenic process, such as acetyl-CoA acetyl transferase (ACAT) and 3-hydroxy-3-methylglutaryl CoA synthases (HMGCS) 1 (cytoplasmic) and 2 (mitochondrial), play an important role in regulating butyrate metabolism at the substrate level. As depicted in Table 4, key enzymes in butyrate metabolic pathways exhibited a different expression pattern between the cell line and the rumen epithelium. As expected, HMGCS2 of mitochondrial origin was significantly up-regulated by butyrate in vitro, in response to increased butyrate concentration. However, expression of HMGCS2 remained unchanged in vivo.

Figure 3

Relative expression of a butyrate transporter, solute carrier family 5 (iodide transporter), member 8 (SLC5A8).

Notes: The number denotes the relative abundance of the transcript in both the bovine epithelial cell and in the bovine rumen epithelium. CT: control, cells treated with PBS; BT: cells treated with 10 mM butyrate for 24 h in vitro. *False discovery rate (FDR) < 0.05; ***FDR < 0.001.

Table 4

Expression levels of key enzymes involved in butyrate metabolism in the rumen epithelium detected by RNA-seq technology.

Gene_id	Annotation	Symbol	In vitro	24 h	72 h	168 h	Post 24 h	Post 168 h
ENSBTAG00000002827	Acetyl-CoA acetyltransferase 2	ACAT2	0.32***	0.95	0.76	0.76	1.11	1.24
ENSBTAG00000011839	3-hydroxy-3-methylglutaryl-CoA synthase 1	HMGCS1	0.63***	0.91	0.67	0.67	1.49	1.69
ENSBTAG00000003898	3-hydroxy-3-methylglutaryl-CoA synthase 2	HMGCS2	1.30*	0.95	0.86	0.94	0.84	0.91
ENSBTAG00000015107	Monocarboxylate transporter 1	MCT1	0.52***	0.98	0.84	0.95	0.90	0.84

Notes: The number denotes fold changes (treated/control).

FDR < 0.05;

FDR < 0.001.

The biological interpretation of high-throughput expression data generated using microarrays or RNA-seq technology requires both differential expression and differential network analyses.38 Many transcriptional regulators exert their impact on biological functions via post-transcriptional mechanisms with subtle or no apparent changes at mRNA expression level detectable by tools for assessing differential expression alone. However, differential network analysis relies on powerful computational tools to extract accurate regulatory gene networks reflecting causal interactions underlying biological processes or phenotypes. Of several algorithms available, those based on information theory, including estimating mutual information values, such as ARACNE26,39 perform well in inferring global gene networks, especially for smaller sample sizes.40,41 ARACNE is based on the assumptions that the expression level of a given gene is a random variable and the mutual relationships between them can be derived by statistical dependences.42 Our results provided further support for the utility of this approach in constructing regulatory gene networks that depict phenotypes and regulation of biological processes. An example from this current study is the regulatory network controlled by FBJ murine osteosarcoma viral oncogene homolog, or c-fos (FOS), which was significantly down-regulated by butyrate in both in vitro and in vivo models. As a transcription factor, FOS dimerizes with another oncogene JUN to form the AP-1 complex, which regulates transcription of a diverse range of genes and is implicated in many biological processes including cell proliferation and differentiation as well as tumor transformation and progression. In the current data set, ARACNE inferred a network of four direct interactions (1st neighbors) for FOS and 32 indirect interactions (Fig. 4). All four direct interactions, cingulin (CGN), heparin-binding epidermal growth factor-like growth factor (HBEGF), intermediate filament family orphan 2 (IFFO2), and jun proto-oncogene (JUN), were up-regulated by butyrate (at both P < 10−6 and FDR < 10−5). Moreover, of 24 genes in the 2nd neighbors category, all were also regulated by butyrate, including three transcription factors, JUN, upstream transcription factor 2, c-fos interacting (USF2), and REST corepressor 1 (RCOR1). GO analysis identified GO terms significantly enriched in this network, including SMAD binding (GO:0046332 and GO:0070412), SMAD protein signal transduction (GO:0060395), and transforming growth factor (TGF)-β receptor signaling (GO:0007179). FOS binding at the TGF-β1 promoter proximal AP-1 site is required for TGF-β1 production by colon carcinoma cells.43 Indeed, previous studies have shown that SMAD proteins cooperate with FOS/JUN complex to mediate TGF-β-induced transcription.44 Thus, ARACNE correctly inferred a direct interaction between FOS and JUN (Fig. 4) as well as interaction between JOS and HBEGF. HBEGF plays a pivotal role in mediating the early cellular response to intestinal injury by serving as a potent cytoprotective factor.45 Other experimental evidence also provides a strong support of a direct interaction between FOS and HBEGF.45,46

Figure 4

A regulatory gene network controlled by FOS.

Notes: The network was inferred using ARACNE at a combined stringent cutoff of an error tolerance ɛ = 0.10 and a P-value threshold of mutual information (MI) at 5.0 × 10−11. FOS had 4 direct interactions (the first neighbors) and 32 indirect interactions. The expression of all genes in this network at the mRNA level was significantly regulated by butyrate. The green/yellow color with gene symbols represents transcription factors.

Tight junctions between epithelial cells regulate the permeability of molecules via para-cellular pathways as well as bacterial translocation across the gut epithelial layer, thus having the potential to strengthen intestinal barrier function and being involved in intestinal pathology.47 In addition to cingulin that was induced by butyrate, at least 16 other tight junction related genes were regulated by butyrate in vitro. A strong up-regulation of the major macromolecular components, such as claudins (CLDN1, CLDN3, CLDN4, CLDN7, CLDN12, and CLDN23), tight junction protein 3 (TJP3), and junctional adhesion molecules 2 and 3 (JAM2 and JAM3), by butyrate supports the contention that these genes may play an important role in maintaining and/or restoring intestinal barrier function. Indeed, evidence indicates that probiotics, such as Lactobacillus plantarum MB452, enhance intestinal barrier function via increasing expression of genes encoding proteins involved in tight function formation.48 It is known that elevated ruminal butyrate results in a profound change in the ruminal microbial composition, including a potential stimulating effect on butyrate-producing bacteria.12 It is foreseeable that an elevated concentration of butyrate in the lumen of the gut could play a regulatory role in the maintenance of intestinal barrier function via the expression of genes involved in tight junctions. Our future work will include experimental verification of the global regulatory gene networks inferred by computational tools.