Skip to content

Advertisement

BMC Genomics

Research article
Open Access

Transcriptome analysis of ruminal epithelia revealed potential regulatory mechanisms involved in host adaptation to gradual high fermentable dietary transition in beef cattle

BMC Genomics201718:976

https://doi.org/10.1186/s12864-017-4317-y

©  The Author(s). 2018

Received: 28 November 2016

Accepted: 17 November 2017

Published: 19 December 2017

Abstract

Background

The transition from a high forage to a highly fermentable diet can induce digestive disorders in the rumen. To date, the host mechanisms that regulate the adaption to such dietary transition are largely unknown. To understand the molecular mechanisms involved in such phenomena, RNA-sequencing was performed to identify the changes in the transcriptome of ruminal epithelia during gradual transition from a diet containing 0% to 89% grain.

Results

In total, the expression of 11,044, 11,322 and 11,282 genes were detected in ruminal epithelia of beef heifers (n = 15) fed 0%, 72% and 89% barley grain diet, respectively. The transcriptome profiles of rumen epithelia differed between low grain diet (LGD) (0% grain) and high grain diet (HGD) (72% and 89%), and HGD tended to reduce the expression of genes involved in epithelial catalytic and binding activities. When diet was changed from 72% to 89% grain, the mean ruminal pH change was significantly different among individual heifers with five of them decreased (down group (DG); from 6.30±0.09 to 5.87±0.15, P < 0.01) and five of them increased (up group (UG); from 5.84±0.42 to 6.35±0.37, P < 0.05). The functional analysis of differentially expressed (DE) genes revealed inhibited “Immune response of leukocytes”, “Attraction of phagocytes”, and “Cell movement of leukocytes” (P < 0.05) functions (Z-score = −2.2, −2.2 and −2.0, respectively) in DG, and inhibited “Concentration of lipid” and “Proliferation of epithelial cells” functions in UG (Z-score = −2.0, and −1.8, respectively). In addition, the expression of genes involved in ketogenesis (HMGCL) and lipid synthesis (SREBF2, FABP4) was increased in DG, while the expression of ketogenesis (ACAT2, HMGCS) and cholesterol synthesis related genes (HMGC and FDPS) were deceased in UG. Furthermore, the upstream regulators were found to be involved in the regulation of immune response and cell cycle progress, and SNP (g.46834311A > G) in FABP4 was identified between two groups of animals (P < 0.1).

Conclusion

The identified genes, upstream regulators, and SNP could be potential genetic markers that may account for the varied individual ruminal pH responses to the dietary transition stress.

Keywords

TranscriptomeRuminal epitheliaAdaptationDietary transitionBeef cattle

Background

To meet the global demand for meat consumption, it has become a common practice in the beef industry to use intensive feeding strategies such as high energy and high concentrate diets to finish cattle [1]. However, feeding high concentrate diets has been reported to be associated with digestive disorders that could lead to ruminal acidosis [2], laminitis [3], liver abscesses [4], and hindgut acidosis [5] in cattle.

When cattle are fed high fermentable diets, the increase and accumulation in short chain fatty acid (SCFA) production together with subsequent decrease in ruminal pH are usually observed due to increased microbial fermentation [6]. In the meantime, ruminal epithelium plays an important role in response to the ruminal pH change and regulation via SCFA absorption [7–9]. Therefore, highly fermentable diets could partly elicit an adaptive response by rumen epithelium to maintain the normal function of rumen. Indeed, when a gradual transition strategy is applied, ruminants have strong ability to adapt to highly fermentable diet [2, 10]. However, the individual variation in the adaptation to high concentrate diet has been widely observed in both beef cattle [2, 11] and dairy cows [12, 13]. It is suggested that the differences in the rate and ways of SCFA absorption by the animal may explain the individual variation in the severity of subacute ruminal acidosis in sheep [8]. Recently, Schlau et al. [11] reported that the expression of sodium hydrogen exchanger isoform 3 (NHE3) gene was different between acidosis-resistant and acidosis-susceptible steers during rapid high grain diet transition, indicating the differences in intracellular proton removal could be attributed to variation in the host response.

To date, the underlying regulatory mechanisms for the host response to ruminal pH change has not been well defined. In this study, we aimed to identify the molecular mechanisms for variation in the response to a gradual transition to a high grain diet within the same animal by characterizing the global gene expression pattern of ruminal epithelia using RNA-seq based transcriptome profiling. The identified mechanisms may help to explain the observed animal variation in maintaining a balanced ruminal pH, which provide new insight into decreasing risk of ruminal acidosis for the beef industry.

Results

Rumen epithelial transcriptomes fed 0%, 72% and 89% grain diets

A total of 1130 million (25.11 ± 2.89 million reads per sample) high-quality 100-bp paired-end reads were obtained from 45 ruminal papillae samples collected from 15 heifers fed diets containing 0%, 72%, and 89% grain, respectively. Of these reads, ~86.1% of them were mapped to the bovine genome (UMD 3.1) and the expression of 11,044, 11,322 and 11,282 genes were detected (with reads per million (RPM) > 1 in 15 heifers fed each diet) in the ruminal epithelial tissue under each dietary condition, respectively. Among them, the expression of 10,880 genes was commonly detected from three dietary conditions (Fig. 1a). The most relevant gene ontology (GO) terms of these commonly expressed genes were “catalytic activity” and “bindings”, followed by “nucleic acid binding transcription factor activity”, “structural molecular activity”, “enzyme regulator activity”, and “receptor activity” (Fig. 1b). When the transcriptome profiles were further compared among three diets, the dietary-dependent expression of genes was detected, with 43, 119 and 102 genes expressed only in 0%, 72% and 89% grain, respectively (Fig. 1a). Further analysis on these diet-dependent genes using DAVID (Database for Annotation, Visualization and Integrated Discovery) showed different functional annotation clusters (Additional file 1 worksheet 2). For example, the enriched gene ontologies for 43 genes under 0% grain diet were “immune response” and “inflammatory response”, and for the 119 genes under 72% grain diet 102 genes under 89% grain, their enriched functions were “extracellular matrix organization”, “protein catabolic process”, respectively. Although the most of these dietary dependent genes were expressed at low level (Fig. 1c) with RPM < 5, the expression of 13 genes showed relatively high abundance (RPM > 5) when heifers were fed diets containing 0% (4 genes), 72% (3 genes), and 89% grain (6 genes), respectively (Table 1).
Figure 1
Fig. 1

Dietary differences of the transcriptome profiles among different diets (0%, 72% and 89% grain). a Venn diagram of expressed genes. b Functional classification of common expressed genes in three diets. c Expression of diet dependent genes

Table 1

Dietary dependent genes in rumen epithelium

Item

ID

Symbol

Function and description

Mean RPM ± SD

0% grain (mean RPM > 5)

ENSBTAG00000008182

FOSB

FBJ murine osteosarcoma viral oncogene homolog B

11.0 ± 8.6

ENSBTAG00000005182

BOLA-A

Major histocompatibility complex, class 1, A precursor

9.1 ± 3.6

ENSBTAG00000019234

BMP6

Bone morphogenetic protein 6

5.9 ± 4.1

ENSBTAG00000015094

VNN1

Vanin 1

5.2 ± 3.6

72% grain (mean RPM > 5)

ENSBTAG00000009812

Novel

C-X-C motif chemokine 6

5.6 ± 5.0

ENSBTAG00000006214

LOXL2

Lysyl oxidase homolog 2 precursor

5.3 ± 1.6

ENSBTAG00000006367

CTGF

Connective tissue growth factor

5.1 ± 2.1

89% grain (mean RPM > 5)

ENSBTAG00000003668

CXorf57

Chromosome X open reading frame 57

9.3 ± 5.1

ENSBTAG00000009144

Novel

Uncharacterized protein

9.0 ± 6.9

ENSBTAG00000005244

RASL1A

RAS-like, family 11, member A

6.6 ± 4.9

ENSBTAG00000021272

ABCG1

ATP-binding cassette, sub-family G, member1

6.4 ± 3.7

ENSBTAG00000025340

CDHR2

Cadherin-related family member 2

5.8 ± 3.3

ENSBTAG00000010423

LIFR

Leukemia inhibitory factor receptor alpha

5.4 ± 2.3

The differentially expressed genes in the ruminal epithelium between low and high grain diets

The PCA plot revealed that the rumen epithelial transcriptome profiles of heifers fed the low grain (0%) diet were different than those fed high grain diets (72% and 89%) except one outlier from 89% diet (Fig. 2a). When the expression of genes was further compared, 562 genes were differentially expressed (DE). Among them, 432 genes showed higher expression (the fold change (FC) in expression >1.5 with the false discovery rate (FDR) < 0.05) in the ruminal epithelia of heifers fed the high grain diet compared to those fed the low grain diet (Fig. 2b, Additional file 1 worksheet 1) with 130 from 72% grain diet, 191 from 89% grain diet and 111 from both high grain diets (Fig. 2b). These genes were defined as high grain diet enriched (HGD) DE genes. On the other hand, 104 genes had higher expression (the FC in expression >1.5 with the FDR < 0.05) in heifers fed the 0% grain diet comparing to those fed high grain diets with 11 compared with the those fed 72% grain diet, 57 compared with the those fed 89% grain diet, and 36 from both high grain diets (Fig. 2c; Additional file 1). These genes were defined as low grain diet enriched (LGD) DE genes. Although, the functional classification showed no difference in terms of the major functions between the HGD and LGD DE genes (Fig. 2d), the highly expressed HGD DE genes were mainly involved in “Paxillin signaling” and “Integrin signaling” pathways, while highly expressed LGD DE genes were enriched to the canonical pathway of “Complement system” and “Interferon signaling”, based on DAVID and Ingenuity Pathway Analysis (IPA) functional analysis (Additional file1 worksheets 2 and 3).
Figure 2
Fig. 2

Transcriptomic differences between low grain diet (LGD, 0% grain) and high grain diet (HGD, 72% and 89% grain). a Principle component analysis of the total detected genes fed three diets. b Differentially expressed (DE) genes that enriched in HGD. c DE genes that enriched in LGD. d Functional classification of HGD enriched genes and LGD enriched genes

Variation in responses of heifers during dietary transition from 72% to 89% grain diets

When the ruminal pH change for each individual animal was further compared using mean ruminal pH from a previous study [14], mean ruminal pH became lower in all heifers after the first dietary transition, from 0% to 72% grain. However, after second dietary transition, from 72% to 89% grain, the mean ruminal pH of heifers had three patterns: lower, from 6.30±0.09 to 5.87±0.15, P < 0.01 (down group (DG), 5 heifers); similar, from 6.17±0.25 to 6.19±0.25, P = 0.27 (balanced group (BG), 5 heifers); and higher, from 5.84±0.42 to 6.35±0.37, P < 0.05 (up group (UG), 5 heifers) (Fig. 3a). When the acidosis index values (pH·min/kg) were further compared, they were increased in DG animals and decreased in UG animals when the diet transitioned from 72% to 89% grain (Fig. 3b). The DG and UG heifers were then selected for further transcriptome comparison to identify the genes that may be associated with the potential molecular mechanisms behind such varied responses.
Figure 3
Fig. 3

Categorization and Transcriptomic differences between Down group and Up group. a Mean pH values fed three diets. b Acidosis index fed three diets (c) Differentially expressed (DE) genes during diet transition (from 72% to 89% grain) in Down group and UP group. d Expression profile of the 28 conserved DE gene between Down group and UP group

The differentially expressed genes associated with varied responses to ruminal pH change during the second dietary transition

When the transcriptomes were further compared within DG and UG heifers during the second dietary transition period (89% vs. 72%), the gene expression profiles displayed different changes. In total, 67 DE genes (22 up-regulated and 45 down-regulated; FDR < 0.05, FC > 1.5 or < −1.5) were found in DG (Additional file 2 worksheet 1), and 285 DE genes (122 up-regulated and 163 down-regulated; FDR < 0.05, FC > 1.5 or < −1.5) were found in UG (Additional file 2 worksheet 2). Among them, 23 DE genes were commonly detected from both DG and UG, while 44 and 262 DE genes were only found in DG and UG, respectively (Fig. 3c). Most of the common DE genes (n = 18) showed the same change trend in both DG and UG (14 genes were down regulated and 4 gene were up regulated), while expression of 5 DE genes exhibited opposite change trend between the two groups. In particular, the expression of CYP1A1 and CYP1B1 were up regulated in UG (FC = 1.8 and 1.6) but were down regulated in DG (FC = −1.7 and −2.0), whereas the expression of CHAC1, CDSN and KLK10 were up regulated in DG (FC = 2.3, 2.1 and 1.5) but were down regulated in UG (FC = −2.0, −1.4 and −3.3) (Fig. 3d).

Functional analysis of DE genes detected for DG and UG during the second dietary transition

The functional prediction of 67 DE genes in the DG using IPA revealed that their most relevant functions were “Immune response of leukocytes”, “Attraction of phagocytes”, and “Cell movement of leukocytes” (P < 0.05), and were inhibited when fed 89% grain diet (Z-score = −2.2, −2.2 and −2.0, respectively) (Table 2). In addition, the function of “Quantity of Ca2+”, “Fatty acid metabolism”, and “Mobilization of Ca2+” were also inhibited (Z-score < −1.5) for DG after second dietary transition (Table 2). The predicted function of the down-regulated DE genes (using DAVID n = 45) in the DG was related to “Innate immune response” (P < 0.0001) (Table 3).
Table 2

Enriched functions of DE genes in Down and Up group using IPA (89% grain vs. 72% grain)

Enriched functions

p-Value

Z-score

Molecules

Down group

 Immune response of leukocytes

1.14E-02

−2.2

C3, CXCL8, FCGR2B, IL1B, ISG15

 Attraction of phagocytes

7.33E-05

−2.2

C3, CXCL5, CXCL8, IL1B, VCAN

 Attraction of myeloid cells

5.73E-05

−2.2

C3, CXCL5, CXCL8, IL1B, VCAN

 Cell movement of leukocytes

4.17E-02

−2.0

C3, CXCL5, CXCL8, FOS, IL1B, S100A9, VCAN

 Infection of cells

2.35E-03

−2.0

C2, C3, EIF2AK2, FCGR2B

 Stimulation of leukocytes

1.18E-02

−2.0

C3, CXCL5, CXCL8, IL1B

 Quantity of Ca2+

4.61E-03

−2.0

C3, CXCL8, IL1B, S100A9

 Adhesion of neutrophils

2.42E-03

−2.0

C3, CXCL8, IL1B, S100A9

 Apoptosis of macrophages

3.86E-03

−2.0

EIF2AK2, FOS, IL1B, ISG15

 Activation of granulocytes

5.11E-03

−1.9

C3, CXCL8, FCGR2B, IL1B

 Viral Infection

5.73E-05

−1.9

C2, C3, CXCL8, EGR1, EIF2AK2, FCGR2B, IL1B, MX1, S100A9

 Binding of leukocytes

1.39E-02

−1.8

A2M, C3, CXCL8, FCGR2B, IL1B, S100A9

 Fatty acid metabolism

1.78E-03

−1.8

CXCL8, CYP3A4, FABP4, FCGR2B, IL1B

 Apoptosis of leukocytes

1.14E-02

−1.8

CXCL8, EGR1, EIF2AK2, FCGR2B, FOS, IL1B, ISG15

 Cell death of immune cells

1.78E-03

−1.7

C3, CXCL8, EGR1, EIF2AK2, FABP4, FCGR2B, FOS, IL1B, ISG15, MX1, MYH11

 Degranulation of phagocytes

1.57E-04

−1.7

C3, CXCL8, DUSP1, FCGR2B, FOS, IL1B, S100A9

 Adhesion of immune cells

2.77E-02

−1.6

A2M, C3, CXCL8, IL1B, S100A9

 Binding of professional phagocytic cells

5.75E-03

−1.6

A2M, C3, CXCL8, IL1B, S100A9

 Mobilization of Ca2+

1.14E-02

−1.6

A2M, C3, CXCL8, FCGR2B

 Migration of phagocytes

5.28E-03

−1.6

C3, CXCL5, CXCL8, IL1B, S100A9, VCAN

 Migration of granulocytes

2.35E-03

−1.5

C3, CXCL5, CXCL8, IL1B, S100A9

Up group

 Concentration of lipid

8.12E-02

−2.0

ABCA1, CXCL8, CYP1A1, CYP1B1, S100A9

 Proliferation of epithelial cells

9.40E-02

−1.8

CAV1, CDC25B, CXCL8, EDN1, FLT1, FOS, FRS2, GLI1, KRAS, KRT16, MKI67, NR3C1, ODC1, PTHLH, SULF2

 Recruitment of neutrophils

7.13E-02

−1.6

CAV1, CXCL2, CXCL5, CXCL8, CYP1A1, CYP1B1, EDN1, P2RY2

 Quantity of T lymphocytes

1.15E-01

1.5

ACKR4, BIRC5, CCNB2, CD48, FOS, IL6ST, KRAS, LGALS1, MYH11, MYSM1, NR3C1, PTTG1, TACC3, TGFBR1, TSC1, XIAP

 Differentiation of epithelial cells

7.13E-02

1.7

AURKB, CAV1, DSG1, E2F8, EDN1, GNAQ, KRAS, KRT14, KRT16, PLK1, PTHLH, SH3PXD2A

 Transmigration of myeloid cells

7.13E-02

2.0

CXCL5, CXCL8, ITGAV, LGALS1

 Transmigration of phagocytes

9.78E-02

2.0

CXCL5, CXCL8, ITGAV, LGALS1

 Cell cycle progression

7.13E-02

2.2

BIRC5, CDC25B, CKAP2, FOS, KRAS, LGALS4, MAP3K1, PLK1, XIAP

Table 3

Enriched functions of DE genes in the Down and Up group using DAVID (89% grain vs. 72% grain)

GO terms

P value

Molecules

Down Group

Down-regulated genes (n = 45)

 GO:0045087~innate immune response

8.2E-05

EIF2AK2, MX1, S100A9, C2, NLRC5, C1QB

 GO:0042742~defense response to bacterium

1.2E-04

ISG15, CXCL5, TAP

 GO:0006958~complement activation, classical pathway

2.7E-04

C3, C2, C1QB

 GO:0051607~defense response to virus

9.8E-04

OAS1Y, ISG15, MX1, NLRC5

 GO:0045071~negative regulation of viral genome replication

1.2E-03

EIF2AK2, ISG15, MX1

 GO:0009636~response to toxic substance

2.4E-03

EIF2QK2, CYP1A1, CYP1B1

 GO:0009404~toxin metabolic process

1.1E-02

CYP1A1, CYP1B1

 GO:2,000,427~positive regulation of apoptotic cell clearance

1.5E-02

C3, C2

 GO:0006956~complement activation

1.8E-02

C3, C2

 GO:0034340~response to type I interferon

2.2E-02

ISG15, MX1

 GO:0006955~immune response

2.2E-02

OAS1Y, CXCL5, LOC504773

 GO:0070098~chemokine-mediated signaling pathway

4.7E-02

CXCL5, LOC504773

Up-regulated genes (n = 22)

 GO:0035914~skeletal muscle cell differentiation

3.2E-02

EGR1, FOS

Up Group

Down-regulated genes (n = 163)

 GO:0051301~cell division

4.4E-09

SKA3, ASPM, CKS2, PTTG1, CCNB1, BIRC5, CCNB2, AURKA, SPC24, NCAPH, UBE2C, CCNA2, CDCA3, SPDL1, TPX2

 GO:0007080~mitotic metaphase plate congression

6.0E-05

KIF22, KIFC1, CDCA8, KIF2C, CCNB1, SPDL1

 GO:0007018~microtubule-based movement

1.1E-04

KIF22, KIFC1, KIF2C, KIF11, KIF20A, CENPE

 GO:0007059~chromosome segregation

3.1E-04

SKA3, SLC25A5, KIF11, CENPT, BIRC5, HJURP

 GO:0007067~mitotic nuclear division

3.8E-04

PLK1, SKA3, ASPM, CCNA2, PTTG1, CCNB2, NUF2, SPC24

 GO:0031577~spindle checkpoint

6.1E-04

AURKB, BIRC5, SPDL1

 GO:0045143~homologous chromosome segregation

6.1E-04

PLK1, ESPL1, PTTG1

 GO:0000070~mitotic sister chromatid segregation

1.4E-03

PLK1, CDCA8, ESPL1, SPAG5

 GO:0051988~regulation of attachment of spindle microtubules to kinetochore

2.0E-03

RACGAP1, ECT2, SPAG5

 GO:0007052~mitotic spindle organization

2.7E-03

AURKB, CCNB1, NDC80, AURKA

 GO:0000281~mitotic cytokinesis

4.1E-03

PLK1, RACGAP1, KIF20A, CKAP2

 GO:0090307~mitotic spindle assembly

5.3E-03

KIFC1, KIF11, BIRC5, TPX2

 GO:0034501~protein localization to kinetochore

5.4E-03

AURKB, SPDL1, BUB1B

 GO:0010628~positive regulation of gene expression

6.7E-03

MAPK11, CAV1, ACTA2, FN1, ACTG2, VIM

 GO:0000910~cytokinesis

8.0E-03

CIT, KIF20A, BIRC5, ECT2

 GO:0007094~mitotic spindle assembly checkpoint

1.2E-02

PLK1, BUB1B, BUB1

 GO:0071346~cellular response to interferon-gamma

1.9E-02

SLC26A6, GAPDH, LOC504773

 GO:0001578~microtubule bundle formation

2.4E-02

PLK1, KIF20A, TPPP3

 GO:0006096~glycolytic process

2.7E-02

GAPDH, ENO1, ENO2

 GO:0035606~peptidyl-cysteine S-trans-nitrosylation

2.8E-02

GAPDH, S100A9

 GO:0006816~calcium ion transport

3.9E-02

CAV1, ANXA6, CACNA1G

 GO:0045931~positive regulation of mitotic cell cycle

3.9E-02

CCNB1, BIRC5, CDC25B

Up-regulated genes (n = 122)

 GO:0043507~positive regulation of JUN kinase activity

1.4E-02

FZD5, EDN1, EPHA4

 GO:2,001,237~negative regulation of extrinsic apoptotic signaling pathway

2.6E-02

PHIP, ITGAV, TGFBR1

 GO:0009404~toxin metabolic process

2.9E-02

CYP1A1, CYP1B1

 GO:0018406~protein C-linked glycosylation via 2′-alpha-mannosyl-L-tryptophan

2.9E-02

DPY19L3, DPY19L4

 GO:0002904~positive regulation of B cell apoptotic process

2.9E-02

FNIP1, CD24

 GO:0042632~cholesterol homeostasis

3.0E-02

CD24, ABCA1, EPHX2

 GO:0017144~drug metabolic process

3.8E-02

CYP1A1, FMO5

 GO:0045944~positive regulation of transcription from RNA polymerase II promoter

4.1E-02

PHIP, TET2, FZD5, EDN1, CCNT1, FOS, NRIP1, MYSM1, PCGF5

For 285 DE genes in the UG, functions of “Concentration of lipid” and “Proliferation of epithelial cells” were inhibited (Z-score = −2.0, and −1.8, respectively), while the functions of “Cell cycle progression” and “Transmigration of phagocytes” were activated (Z-score = 2.2 and 2.0, respectively) after second dietary transition (Table 2). The function prediction of down-regulated genes (n = 163) in UG using DAVID revealed functions related to activity in cell proliferation, such as “Cell division”, “Mitotic metaphase plate congression”, “Chromosome segregation”, and “mitotic nuclear division” (P < 0.001) (Table 3). Moreover, the enriched function of up-regulated genes (n = 122) in UG were “Positive regulation of JUN kinase activity”, “Toxin metabolic process”, “Cholesterol homeostasis”, and “Drug metabolic process” (P < 0.05) (Table 3).

Further functional pathways analysis using IPA showed eight DE genes in DG were found to be involved in two different pathways including “Role of pattern recognition receptors in recognition of bacteria and viruses” (Z-score = −2.0) and “TREM1 signaling” (Z-score = −2.0) (Table 4). For the UG, the “PTEN signaling” pathway was inhibited (Z-score = −2.0), while the pathways of “NF-κB signaling”, “PDGF signaling”, and “Cell cycle: G2/M DNA damage checkpoint regulation” were activated (Z-score = 2.4, 2.2 and 2.2, respectively). The KEGG pathways of those DE genes enriched by DAVID showed similar results (Table 5). Briefly, eight DE genes in DG were found to be involved in three KEGG pathways including “Staphylococcus aureus infection” and “complement and coagulation cascades” (5 down regulated genes), and “T cell receptor signaling pathway” (3 upregulated genes) (Table 5). In the UG, 25 DE genes were enriched in ten different KEGG pathways. Among them, the downregulated genes were involved in four pathways including “Cell cycle”, “Glycolysis/Gluconeogenesis”, “FoxO signaling pathway”, “p53 signaling pathway”, and “Mineral absorption”; the upregulated genes were mainly related to metabolic pathways including “tryptophan metabolism”, “Retinol metabolism”, “Drug metabolism – cytochrome P450”, and “metabolism of xenobiotics by cytochrome P450” (Table 5).
Table 4

Enriched ingenuity canonical pathways of DE genes in the Down and Up group using IPA (89% grain vs. 72% grain)

Ingenuity canonical pathways

-log (p-Value)

z-score

Molecules

Down group

 Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses

3.42E + 00

−2.0

CXCL8, OAS1, C3, IL1B, C1QB, EIF2AK2

 TREM1 Signaling

2.59E + 00

−2.0

CXCL8, NLRC5, IL1B, FCGR2B

Up group

 PTEN Signaling

6.06E-01

−2.0

TGFBR1, NTRK2, FLT1, KRAS

 Renin-Angiotensin Signaling

9.28E-01

1.6

FOS, MAP3K1, GNAQ, KRAS, MAPK11, FRS2

 IL-8 Signaling

8.71E-01

1.6

RAB11FIP2, CXCL8, FOS, FLT1, ITGAV, KRAS, FRS2

 ERK5 Signaling

9.19E-01

2.0

IL6ST, FOS, GNAQ, KRAS

 JAK/Stat Signaling

8.00E-01

2.0

FOS, GNAQ, KRAS, FRS2

 NF-κB Activation by Viruses

7.90E-01

2.0

MAP3K1, ITGAV, KRAS, FRS2

 VEGF Family Ligand-Receptor Interactions

7.76E-01

2.0

FOS, FLT1, KRAS, FRS2

 Rac Signaling

6.17E-01

2.0

MAP3K1, PIKFYVE, KRAS, FRS2

 PKCθ Signaling in T Lymphocytes

5.46E-01

2.0

FOS, MAP3K1, KRAS, FRS2

 Role of NFAT in Regulation of the Immune Response

4.09E-01

2.0

FOS, GNAQ, KRAS, FRS2

 Cell Cycle: G2/M DNA Damage Checkpoint Regulation

2.15E + 00

2.2

CDC25B, CKS2, CCNB2, PLK1, AURKA, CCNB1

 PDGF Signaling

9.28E-01

2.2

FOS, MAP3K1, CAV1, KRAS, FRS2

 NF-κB Signaling

7.60E-01

2.4

TGFBR1, NTRK2, FLT1, MAP3K1, KRAS, FRS2

Table 5

Enriched KEEG pathways of DE genes in the Down and Up group using DAVID (89% grain vs. 72% grain)

Term

p value

Molecules

Down Group

Down-regulated genes (n = 45)

 bta05150: Staphylococcus aureus infection

1.5E-04

FCGR2B, C3, C2, C1QB

 bta04610: Complement and coagulation cascades

2.9E-04

A2M, C3, C2, C1QB

Up-regulated genes (n = 22)

 bta04660: T cell receptor signaling pathway

9.34E-03

CD3D, CD3G, FOS

Up Group

Down-regulated genes (n = 163)

 bta04110:Cell cycle

1.07E-05

PLK1, CCNA2, ESPL1, CDC20, PTTG1, CCNB1,

 bta00010:Glycolysis/Gluconeogenesis

2.02E-02

GAPDH, LDHB, ENO2, ENO1

 bta04068:FoxO signaling pathway

4.35E-02

PLK1, MAPK11, BNIP3, CCNB1, CCNB2

 bta04115:p53 signaling pathway

4.37E-02

GTSE1, RRM2, CCNB1, CCNB2

 bta04978:Mineral absorption

4.48E-02

SLC26A6, MT1E, MT2A

Up-regulated genes (n = 122)

 bta00380: Tryptophan metabolism

2.77E-02

AOX1, CYP1A1, CYP1B1

 bta00830:Retinol metabolism

3.00E-02

CYP1A1, ADH6, AOX1

 bta00982:Drug metabolism - cytochrome P450

3.73E-02

ADH6, AOX1, FMO5

 bta00980:Metabolism of xenobiotics by cytochrome P450

4.52E-02

CYP1A1, ADH6, CYP1B1

Upstream regulator and network analysis of DE genes detected for DG and UG during the second dietary transition

Four upstream regulators were identified in regulating the DE genes in DG, including Interferon gamma (IFNγ), transmembrane protein 173 (TMEM173), toll like receptor 3 (TLR3), and tumor necrosis factor (TNF) (Fig. 4a). For the DE genes in UG, nine upstream regulators were identified with colony stimulating factor 2 (CSF2) and prostaglandin E receptor 2 (PTGER2) being the hub nodes (Fig. 4b). Two and three networks were enriched (score > 20) in DG and UG, and the related functions were immune response and cell cycle, respectively (Additional file 2 worksheets 3 and 4).
Figure 4
Fig. 4

Upstream regulator analysis of DE genes in Down group (DG) and Up group (UG). a Upstream regulators and target genes in DG. b Upstream regulators and target genes in UG

Differential expression of genes involved in lipid transport, metabolism and intracellular homeostasis regulation in DG and UG

Among above genes identified in DG and UG animals, the expression of 14 lipid transport and 13 fatty acid metabolism related genes (RPM > 1 in 5 heifers of at least one diet) were further investigated in the rumen epithelia (Fig. 5). In DG heifers, the expression of SMCT1 (P < 0.1) was lower while the expression of FABP4 (P < 0.05), HMGCL (P < 0.1), and SREBF2 (P < 0.1) was higher when fed the 89% grain diet than the 72% grain diet. In the UG heifers, the expression of FABP5 (P < 0.05), ABCA2 (P < 0.05), and ABCA7 (P < 0.1), ACAT2 (P < 0.05), HMGCS1 (P < 0.1), HMGCR (P < 0.1), and FDPS (P < 0.1) was lower, whereas the expression of ABCA1 (P < 0.05), and ABCA5 (P < 0.1) was higher when fed the 89% grain diet than the 72% grain diet (Fig. 5a & b). In addition, the expression pattern of genes associated with ion transportation in rumen epithelium was further investigated. The expression of NHE8 (P < 0.1), and MCT4 (P < 0.05) was decreased in DG, while the expression of NHE3 (P < 0.05) was increased in UG after the second dietary transition (Additional file 3: Figure S1). Further validation of the differential expression using reverse transcription quantitative real-time PCR (RT-qPCR) showed the expression of FABP4, ABCA1, NHE3 and ACAT2 was consistent as detected with the RNA-seq data (Additional file 4: Figure S2).
Figure 5
Fig. 5

The patterns of lipid transport and fatty acid metabolism related genes expression during diet transition (from 72% to 89% grain) in DG and UG. a The pattern of lipid transport related genes expression. b The pattern of fatty acid metabolism related genes expression. The legend represents the RPM value scaled by rows, the blue means highly expressed in the 72% grain and yellow means highly expressed in the 89% grain. The data were analyzed by T-test, * indicated P < 0.1, ** indicated P < 0.05, and *** indicated P < 0.01. ↓ means decreased and ↑ means increased when the diet transitioned from 72% grain to 89% grain within each group

In addition, 70 SNPs were detected in 9 genes that relate to lipid transport, metabolism and intracellular homeostasis regulation (Additional file 2 worksheet 5), and association analysis (Fisher’s exact test) of detected single nucleotide polymorphisms (SNPs) indicated that one intronic SNP (g.46834311A > G) in FABP4 was associated with the varied expression in the ruminal epithelial tissue between UG and DG animals (P < 0.1, Additional file 5: Figure S3).

Discussion

The comparative analysis of the whole transcriptome profiles revealed the high grain diet can alter the gene expression in the rumen epithelia. Although gene expression of rumen epithelial tissues in response to high concentrate diets have been studied, most of the previous studies have focused on selected genes using qPCR/microarray [15–17] and/or compared the difference using different groups of animals [11, 18] as well as to study the dairy cows with small numbers of animals [15]. This is the first study that applied RNA-seq based genome-wide transcriptome analysis to study the global gene expression changes in ruminal epithelia under three dietary conditions for the same individual using 15 beef heifers. Furthermore, our effort is the first to identify the varied expression pattern changes in the rumen epithelia of animals having different pH change patterns with the gradual high grain transition from 72% to 89% grain.

It has been known that feeding highly fermentable diets could be associated with the prevalence of ruminal acidosis [6, 19]. When ruminal pH duration time, a widely used diagnostic parameter for the acidosis [20, 21], was analyzed, 8 out of 15 heifers had pH below 5.8 for longer than 5.4 h/d when fed 89% grain, suggesting they may have developed subacute acidosis. In contrast to these heifers, the others had lower ruminal pH duration time when fed the same diet, emphasizing that response variation existed among individuals as previously reported [11, 12], which was reflected by our observation of DG and UG heifers. Together with the mean of continuous ruminal pH, it suggests that each heifer may have developed a different adaptive mechanism, with DG responding less favorably than UG. To understand the mechanism behind such individual variation, the following discussion will be mainly focused on the gene expression pattern variation between in DG and in UG animals.

The most enriched networks of DE genes in DG and UG were related to immune response and cell cycle, respectively, suggesting that the innate immune function might be differently regulated between DG and UG heifers during second dietary transition. Although scarce information is available about the organization of the mucosal immune system in the rumen epithelia, the Langerhans cells were identified in epithelium of bovine forestomach [22] and bovine forestomachs could receive, elaborate and produce signals and mediators of the innate immune response [23]. Here, we found the innate immune response related functions such as “Immune response of leukocytes”, “Attraction of phagocytes”, and “Cell movement of leukocytes” were inhibited in DG but not in UG, suggesting that innate immune response could be one of the mechanisms involved in epithelial adaptation to high grain diet. In addition, the most abundant GO categories enriched by downregulated DE genes in DG were related to immune response. Decreased expression of immune-related genes (including IL-6, IL-10, and interferon) had been reported in ruminal epithelia of lambs when fed concentrate starter [24]. Similar results were found in DG, in particular, the expression of genes encoding complement component 2/3 (C2 / C3), C-X-C chemokine ligand 5 (CXCL5), and NOD-like receptor C5 (NLRC5) were downregulated during the diet transition from 72% to 89% grain. The C2/C3 plays central role in the complement system and contributes to innate immunity [25]. It’s well known that epithelial cells recognized microbial components through means of Toll-like and NOD-like receptors (TLRs and NLRs) [26, 27]. The decreased expression of NLRC5 might limit the recognition of commensal bacteria, which is essential for the development and function of the immune system in the mucosal and peripheral districts [28]. Indeed, the epimural bacterial population differed between acidosis tolerant and susceptible beef cattle and the expression of TLR2 and TLR4 were lower in susceptible cattle [29]. The previous study also showed varied microbial changes associated with the ruminal epithelia for the same animals [14]. The observed potential inhibited innate immune responses further support our previous speculation that host-microbial interactions could play a role in affecting the host adaptation to the high grain diet transition. The lowered innate immune function may negatively affect the rumen function which result in the persistent decrease of ruminal pH during the second transition in DG.

In addition, the enriched pathway of “G2/M DNA damage checkpoint regulation” was activated through downregulation of genes related to cell cycle progression in UG [30–32], including cell division cycle 25 homolog B (CDC25B), cyclin-dependent kinases regulatory subunit 2 (CKS2), cyclin B1/2 (CCNB1/2), polo like kinase 1 (PLK1), and aurora kinase A (AURKA), suggesting the enhanced cell cycle arrest in UG. Furthermore, the enriched KEGG pathway of “p53 signaling” in UG also indicates the higher cell arrest since the activated p53 could activate cell repair and apoptosis procedure to induce cell cycle arrest [33, 34]. Moreover, the pathway of “Metabolism of xenobiotics by cytochrome P450” was enriched in UG heifers. The higher expression of three genes involved in this pathway (alcohol dehydrogenase 6 (ADH6), cytochrome p450 family 1 subfamily A member 1(CYP1A1), and subfamily B member 1 (CYP1B1)) in UG suggest the higher capacity to remove the ruminal toxins in these animals since these genes have been reported to be involved in the metabolism of the some ruminal toxins, such as ethanol and serials of xenobiotics [35, 36]. Taking these together, we speculate that the UG animals may have more activated innate immune responses, cell repair function and toxin removal activities to maintain the ruminal health during the second dietary transition. The upstream regulator analysis further supports our speculation since the identified IFNG, TNF, TLR3, and TMEM173 regulate genes involved in the immune function, while the CSF2 and PTGER2 regulate genes involved in cell proliferation.

Considering the fact that VFAs absorption accounts for up to 50% of the ruminal buffering capacity [7], the increase in intra-epithelial uptake and metabolism of SCFAs could promote the uptake of SCFAs and a stabilization of ruminal pH [8]. We hypothesized that the expression patterns of genes involved in SCFAs absorption and metabolism may also play an important role in individualized animal rumen pH changes when diet was switched from 72% grain to 89% grain. The expression of sodium dependent monocarboxylate transporters 1 (SMCT1) tended to decrease in the DG heifers but no difference was observed in the UG heifers during the second transition period. The SMCT1 and SMCT2 play important roles in SCFA transport in addition to the SCFA−/HCO3− exchanger system [37, 38]. These results suggest that the absorption of SCFA via the SMCT1 in the DG heifers might be decreased, which may partly account for the lower mean pH in DG heifers. In addition to SCFA absorption, SCFA metabolism could also impact on the rumen pH environment. In this study, we focused on the genes involved in butyrate catabolism in the rumen epithelia, which occurs via ketogenesis and beta-oxidation to produce ketone bodies [39, 40]. The expression of ACATs and HMGCS, the genes encoding the rate-limiting enzymes for ketogenesis [41], decreased in UG after the second transition period. However, the expression of HMGCL increased in DG for the same dietary transition period, suggesting that UG heifers have decreased ketogenesis in the epithelia comparing to DG.

In addition, the expression of genes involved in fatty acids and lipids transport and metabolism were also found to be different between DG and UG heifers. For example, the increased expression of FABP4 was observed in DG heifers but not in UG heifers when diet transitioned from 72% to 89% grain. Fatty acid binding proteins (FABPs) usually bind to fatty acids and/or lipids to transport and deliver them to different sites for utilization [42]. It has been reported that increased expression of FABP4 was responsible for the enhanced macrophage lipid accumulation in diabetic patients of human, and this effect was further paralleled with elevated intracellular total cholesterol and triacylglycerol levels [43]. Therefore, the increased expression of FABP4 in the ruminal epithelia of DG heifers suggest the intracellular lipid accumulation, which could be partially supported by the tended increased expression of sterol regulatory element binding protein 2 (SREBF2) in the same group of animals. The SREBF family plays important role in regulating cholesterogenic and lipogenic gene expression in the bovine liver and mammary glands [44, 45], while the isoform SREBF2 preferentially activates cholesterol biosynthesis [46]. Although cholesterol synthesis pathway has not been well characterized in the ruminal epithelia, it’s reported that butyrate can be used as a precursor for cholesterol biosynthesis [46, 47]. Compared with DG, the expression of cholesterol synthesis related genes (HMGCS1, HMGCR and FDPS) was decreased in UG heifers when diet transitioned from 72% to 89% grain. Steele et al. [15] also reported the coordinated downregulation of these genes in the rumen epithelia of dairy cows when fed a 65% grain diet from the first to third week. These suggest that the intracellular cholesterol synthesis could also be one of the mechanisms for the individualized responses when adapted to the high grain diet. It is known that cholesterol is the key component of the plasma membrane and plays crucial role in cellular signal transduction, cell growth, cell polarization, migration, and survival [48]. The cellular content of cholesterol is determined not only by the synthesis but also by the efficiency of influx and efflux [48] transportation of them and their derivatives through ABC transporters [49, 50]. Among them, ABCA1 has been most intensively studied, which could regulate and maintain cellular cholesterol homeostasis through transferring cholesterol to high-density lipoprotein in various types of cells [51, 52]. In this study, the increased expression of ABCA1 was observed in the UG, but not in DG heifers. In addition, the function of “Concentration of lipid” was inhibited in the UG heifers, suggesting the UG might be more efficient in exporting cholesterol out of rumen epithelia and have less stored cholesterol. Taken together, our results suggest that the deposition of cholesterol in addition to cholesterol synthesis in the rumen epithelia is also one of the mechanisms to attribute to the variation observed between DG and UG animals (Additional file 6: Figure S4). Moreover, previous studies showed that increased intracellular cholesterol could inhibit various ion channels in the biological membrane of human [53, 54], and high cholesterol is closely related to membrane permeability and inflammation in human [55]. Therefore, the inhibited function of “Quantity of Ca2+” and innate immunity in the DG animals could be resulted from higher concentration intracellular content of cholesterol during the dietary adaptation from 72% to 89% grain. Further measuring the cholesterol concentration in rumen epithelia is needed to verify the findings obtained from our transcriptome analysis.

Furthermore, this study also identified the genetic variation of genes related to lipid transport and metabolism. Association analysis indicated that 1 SNP (g.46834311A > G) in FABP4 was associated with the varied ruminal pH response between DG and UG heifers, which may be attributed to the changes in the expression of this gene between DG and UG animal after the secondary diet transition. This suggests that the genetic background between “Up” and “Down” heifers could be one of the factors that lead to the observed varied gene expression. To determine the role of genetic variation in regulating the genes involved in different responses to high fermentable dietary transition, future researches on genotyping of large population with varied phenotypic measures (ruminal pH, duration and area of ruminal pH below 5.8) using genome wide association study (GWAS) are needed.

Conclusions

This is the first study that has identified varied molecular mechanisms may contribute to the individual variations in response to the high grain adaptation through analyzing the genome wide transcriptome in the ruminal epithelia of beef heifers. Understanding the whole transcriptome is essential to reveal the molecular events within cells, and to elucidate the mechanisms regulating the adaptive function and physiology of rumen in response to dietary changes. Overall, we suggest that the different genes that controlled immune function, cellular repair function, and intracellular homeostasis (cholesterol) might be the molecular mechanism accounting for individual variation in the response to gradual transition to a high grain diet. As summarized in Fig. 6, feeding high fermentable diet usually results in increased concentration of ruminal VFA and toxins, which may induce the cell damage for both DG and UG. At transcriptomic level, the innate immunity, cell cycle, toxin metabolism, and cholesterol homeostasis was differentially regulated in DG and UG. In the DG, the inhibited innate immune response and accumulated cholesterol in the cells may aggravate the cell damage and finally decrease the cellular stability and homeostasis. However, in the UG, cell cycle arrest and xenobiotic metabolism was activated for repairing and protecting the cell damage, which contributes to increase cellular stability and homeostasis. It is noticeable that, gene expression can be regulated at both transcriptional and post-transcriptional level, and even post-translational modifications can also affect the activity of protein [56, 57]. Future studies on the protein, enzyme and metabolites are needed to comprehensively determine the molecular mechanisms of individualized host response to high grain feeding. Regardless, the transcriptomic information from this study provides the informative clue for further study such as the identified gene networks, upstream regulators (IFNG, TNF, TLR3, CSF2, and PTGER2), and FABP4-A/G SNP could be potential gene and genetic markers for selecting cattle with maintained ruminal pH during high fermentable diet transition.
Figure 6
Fig. 6

Illustrative comparison of the gene networks related to distinct adaptation between Down group and Up group. ↓ means decreased and ↑ means increased when the diet transitioned from 72% grain to 89% grain within each group

Methods

Experimental design and sample collection

Twenty-four Angus-Hereford cross-bred yearling heifers (about 8 months old, weighing 244 kg to 369 kg) from the University of Alberta Kinsella Ranch (Kinsella, AB) were used for this study. Detailed information on the animal study, diets, feeding of diets, and intake has been reported previously [14]. Briefly, heifers were blocked into 1 of 6 blocks based on body weight and, within block, were randomly assigned to 1 of 2 treatments. Each block consisted of 1 heifer assigned to the control (CON, n = 6)) treatment and 3 heifers assigned to the rapid grain adaptation (RA) treatment (n = 18). Unequal weighting of the treatment assignment was used as more variation was expected for the RA treatment than for the CON (Bevans et al., 2005). All the RA heifers we randomly assigned to a gradual grain transition treatment (n = 18) throughout a 29-d experiment period. Briefly, heifers were initially fed a diet containing 0% grain (days 1 to 4) and transitioned to a final diet containing 89% grain using the following intermediate diets: 40% grain (days 5 to 8), 60% grain (days 9 to 12), 72% grain (days 13 to 16), 85% grain (days 17 to 20), and 89% grain (days 21 to 29). Rumen papillae biopsies were collected from the ventral sac of the rumen when the heifers were fed 0% grain on day 3 (the 3rd day of a 4-day feeding period), 72% grain on day 15 (the 3rd day of a 4-day feeding period), and 89% grain on day 26 (the 6th day of a 9-day feeding period). The biopsies were immediately washed with sterile 0.01 M phosphate-buffered saline buffer (pH 6.8) and transferred into RNA-later solution (Invitrogen, Carlsbad, CA). After being kept at 4 °C overnight, all the samples were stored at −20 °C for further analysis. To identify the molecular mechanisms responsible for the varied response to a gradual high grain diet transition, the ruminal papillae samples from RA heifers were used in this study and due to three of animals did not have biopsy samples for all three treatment periods, only fifteen animals (under three dietary conditions) were subjected to transcriptome analysis.

Ruminal pH measurement

Ruminal pH was measured using a stand-alone submersible continuous ruminal pH measurement system following the procedures as reported previously [58]. Briefly, the pH meter was inserted into the rumen through a cannula on day 1 and was removed every four days for standardization. Ruminal pH data were collected every 30 s for each diet treatment (0%, 72% and 89% grain diet), and the daily minimum, mean, and maximum pH values were determined. In addition, the number of episodes that ruminal pH was below the threshold of 5.8, as well as the duration and area of these episodes were calculated. The area was calculated as pH unit × min/day when the ruminal pH was lower than 5.8 and the duration was calculated the time (h/day) that ruminal pH was below 5.8 in the heifers.

The acidosis index, an indicator of the severity of ruminal acidosis, was calculated by the following formula: acidosis index (pH·min/kg) = area that ruminal pH below 5.8 (min)/dry matter intake (kg) [18].

RNA isolation

Ruminal papillae samples were ground while immersed in liquid nitrogen using a frozen mortar and pestle prior to RNA extraction. Total RNA was extracted from 50 mg of the resultant papillae tissue powder using Trizol reagent (Invitrogen, CA, USA) following the manufacturer’s instructions. The quantity and quality of the RNA were determined using Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA) and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA), respectively. The samples with RNA integrity number (RIN) greater than 8.0 were used for library construction.

RNA-seq library construction and sequencing

The RNA-seq libraries were constructed using the TruSeq mRNA Sample Preparation Kit (Illumina, San Diego, CA) according to the manufacturer’s instructions. Total RNA (1.0 μg each) from 45 samples (each with a unique index) were mixed and pooled for transcriptome sequencing (paired ends sequencing, 100 bp) at Génome Québec (Montréal, Canada) using the Illumina HiSeq 2000 system (Illumina) with all 45 samples loaded on a single chip with 9 samples per lane: each lane contains samples from 0%, 72%, and 89% grain diet. All reads were demultiplexed according to their index sequences with CASAVA version 1.8 (Illumina) and reads that did not pass the Illumina chastity filter were discarded.

RNA-seq reads mapping and annotation

RNA-seq reads were aligned to the bovine genome (UMD 3.1) using Tophat 2.0.10 with default parameters [59]. The number of reads mapped to each gene was counted by htseq-count (http://htseq.readthedocs.io/) based on the annotation from the ENSEMBL (http://uswest.ensembl.org/) bovine gene annotation software (v75.30). The expression levels of mRNA in each library was calculated by normalizing reads to reads per million (RPM) by the following formula: RPM = (gene reads number/total mapped reads number per library) × 1,000,000.

Identification of differentially expressed genes

Among the detected genes, those that were detected in at least one sampling point with more than 1 RPM in 15 heifers were considered as expressed genes. Then, differentially expressed (DE) genes were investigated by using bioinformatics tool edgeR [60]. The DE genes affected by diet were identified by comparing any two diets (72% vs. 0%; 89% vs. 0%) (n = 15). After that, the individual DE genes were identified by comparing the two high grain diets (89% vs. 72%) of each group (n = 10). The DE mRNA were identified by false discovery rate (FDR) < 0.05 based on Benjamini and Hochberg multiple testing correction [61] as well as a fold change >1.5 or < −1.5. The expression of four DE genes that related to lipid synthesis (FABP4), lipid transport (ABCA1), Na+/H+ exchange (NHE3), and ketogenesis (ACAT2) was detected using RT-qPCR.

Quantitative real-time PCR (qRT-PCR) analysis

A total of 4 genes were selected to validate the DE genes and the primers used were listed in Additional file 7: Table S1. The RT-qPCR reactions were performed with SYBR Green (Fast SYBR® Green Master Mix; Applied Biosystems) using StepOnePlus™ Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) with the fast cycle and the following program: 20 s pre-denaturalization at 95°C, followed by 40 cycles of 3 s denaturation at 95°C and 30s annealing and extension at 60°C. Gene expression values were normalized to reference gene of β-actin in the same sample. The relative changes in each gene expression were calculated using the 2-ΔΔCT (cycle threshold, CT) method.

Functional analysis

The gene list analysis tool in PANTHER classification system was used for gene ontology (GO) terms analysis of the commonly and highly expressed genes among 0%, 72%, and 89% die, and the ratio was calculated according to the number of ‘hits’ to the terms over the total number of ‘class hits’ [62]. The functional analysis of dietary enriched/DE genes was performed by Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov) [63], and Ingenuity pathway analysis (IPA, Ingenuity Systems, www.ingenuity.com). In details, we used all detected genes in rumen epithelia (11,516 genes) as background of DAVID for functional analysis (Additional file 1), and selected the database of digestive tract tissues (including forestomach, stomach, small intestine, and large intestine) and epithelial/immune cells for IPA analysis in this study. A threshold of P < 0.05 for DAVID and Benjamini and Hochberg multiple testing correction P < 0.05 for IPA, and molecules number > 2 was applied to enrich significant biological functions and pathways. The absolute value of Z-score > 1.5 was used as cutoff for activation or inhibition of biological functions and pathways, and the cutoff for networks was score > 20.

Systematic analysis of lipid transport and fatty acid metabolism related genes in DG and UG during diet transition

The expression of sodium-linked monocarboxylate transporter (SMCT1 and SMCT2), proton-linked monocarboxylate transporter (MCT1, MCT2, MCT3, and MCT4), long chain fatty acid transporter (FABP1, FABP2, FABP3, FABP4, FABP5, FABP6, FABP7, FABP9, and FABP12), cholesterol efflux protein (ABCA1, ABCA2, ABCA3, ABCA4, ABCA5, ABCA6, ABCA7, ABCA9, ABCA10, ABCA12, and ABCA13), bicarbonate transporter (PAT1), anion exchanger (AE2, AE3, AE4, DRA) were analyzed by T-test in rumen epithelia of DG and UG heifers when the diet switched from 72% to 89% grain.

The SCFA metabolism related genes including acyl-CoA synthetases (ACSS1, ACSS2 and ACSS3) that activate the volatile fatty acids, acetyl-CoA acetyl transferases (ACAT1 and ACAT2) and 3-hydroxy, 3-methylglutaryl CoA synthase (HMGCS1) that convert acetyl-CoA to 3-hydroxy, 3-methylglutaryl CoA (HMG-CoA) [47],HMG-CoA lyase (HMGCL) which catalyze the synthesis of ketone bodies acetoacetate and β-hydroxybutyrate [47, 64], cholesterolgenic genes: HMG-CoA reductase (HMGCR), farnesyl-diphosphate farnesyl-transferase 1 (FDFT1), farnesyl diphosphate synthase (FDPS), and lanosterol synthase (LSS), and sterol regulatory element binding transcription factors (SREBF1 and SREBF2) were analyzed by T-test in rumen epithelia of DG and UG heifers during diet transition (from 72% to 89%).

SNP analysis

The 5 UG and 5 DG heifers were chose for SNP analysis. RNA-Seq reads of ruminal papillae fed 0%, 72% and 89% grain diets were first combined for each heifer to increase read coverage. Then the SNP calling was performed using VarScan2 [65]. The minimum based quality of reads was 15, minimum reads depth at a position to call a SNP was 8, and the minimum variant allele frequency threshold was 0.1. The association between alleles and varied ruminal pH response was determined by Fisher exact test, and P < 0.1 was considered as significance difference. The SNP was discarded when less than 2 reads mapped to that location in at least one heifer in this study.

Abbreviations

ABCA : 

Cholesterol efflux protein

ACAT : 

Acetyl-CoA acetyl transferases;

ACSS : 

Acyl-CoA synthetases

ADH6 : 

Alcohol dehydrogenase 6

AE2 : 

Anion exchanger

AURKA : 

Aurora kinase A

C2/C3 : 

Complement component 2/3

CCNB1/2: 

Cyclin B1/2;

CDC25B : 

Cell division cycle 25 homolog B

CDSN : 

Corneodesmosin

CHAC1 : 

ChaC glutathione specific gamma-glutamylcyclotransferase 1

CKS2 : 

Cyclin-dependent kinases regulatory subunit 2

CON: 

Control

CSF2 : 

Colony stimulating factor 2

CXCL5: 

C-X-C chemokine ligand 5

CYP1A1 : 

Cytochrome p450 family 1 subfamily A member 1

CYP1B1 : 

Subfamily B member 1

DAVID: 

Database for Annotation, Visualization and Integrated Discovery

DE: 

Differentially expressed

DG: 

Down group

DRA : 

MHC class II DR alpha chain

FABP : 

Fatty acid binding protein

FC: 

Fold change

FDFT1 : 

Farnesyl-diphosphate farnesyl-transferase 1

FDPS : 

Farnesyl diphosphate synthase

FDR: 

False discovery rate

GO: 

Gene ontology

GWAS: 

Genome wide association study

HGD: 

High grain diet

HMGCL : 

3-methylglutaryl CoA lyase

HMGCR : 

3-methylglutaryl CoA reductase

HMGCS : 

3-hydroxy, 3-methylglutaryl CoA synthase

IFNg : 

Interferon gamma

IL: 

Interleukin

IPA: 

Ingenuity Pathway Analysis

KEGG: 

Kyoto Encyclopedia of Genes and Genomes

KLK10 : 

kallikrein related peptidase 10

LGD: 

Low grain diet

LSS : 

lanosterol synthase

MCT : 

Proton-linked monocarboxylate transporter

NHE : 

Sodium hydrogen exchanger

NLRC5 : 

NLR FAMILY card domain containing 5

NLRC5 : 

NOD-like receptor C5

NLRs: 

NOD-like receptors

PAT1 : 

Bicarbonate transporter

PCA: 

Principal component analysis

PDGF: 

Platelet-derived growth factor

PLK1 : 

Polo like kinase 1

PTGER2 : 

Prostaglandin E receptor 2

RA: 

Rapid grain adaptation

RIN: 

RNA integrity number

RNA: 

Ribonucleic acid

RPM: 

Reads per million

RT-qPCR: 

Reverse transcription quantitative real-time PCR

SCFA: 

Short chain fatty acid

SMCT1 : 

Sodium dependent monocarboxylate transporters

SNP: 

Single nucleotide polymorphisms

SREBF2 : 

Sterol regulatory element binding protein 2

TLR3 : 

Toll like receptor 3

TMEM173 : 

Transmembrane protein

TNF: 

Tumor necrosis factor

TREM1: 

Triggering receptor expressed on myeloid cells 1

UG: 

Up group

Declarations

Acknowledgements

The authors acknowledge the technical assistance by X. Sun for the RNA-seq library construction and the members of Drs. Guan and Oba labs for the sampling.

Funding

This study is supported by Alberta Livestock Meat Agency (ALMA 2011F129R and ALMA2015P008R), NSERC discovery grant for L. L. Guan and Chinese International Postdoctoral Exchange Fellowship Program (20140020) for K. Zhao.

Availability of data and materials

All the sequencing data sets supporting the results in the current study have been deposited in the publicly available NCBI’s Gene Expression Omnibus Database (http://www.ncbi.nlm.nih.gov/geo/) (GSE89874).

Authors’ contributions

KZ did the data analysis and draft the manuscript. LLG, GBP, and MO designed the experiments and revised the manuscript critically for important intellectual content. YHC collected the samples and constructed the libraries. All authors contributed to the writing of the manuscript and approved the submission.

Ethics approval

The animal study was approved by the Livestock Animal Care Committee of the University of Alberta (protocol number OBA077) and all procedures were conducted following the guidelines of Canadian Council on Animal Care.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
Department of Agricultural, Food and Nutritional Sciences, University of Alberta, Edmonton, Canada
(2)
College of Food Engineering and Nutritional Science, Shaanxi Normal University, Xi’an, China
(3)
Department of Animal and Poultry Science, University of Saskatchewan, Saskatoon, Canada

References

  1. Wierenga KT, McAllister TA, Gibb DJ, Chaves AV, Okine EK, Beauchemin KA, Oba M. Evaluation of triticale dried distillers grains with solubles as a substitute for barley grain and barley silage in feedlot finishing diets. J Anim Sci. 2010;88(9):3018–29.View ArticlePubMedGoogle Scholar
  2. Bevans DW, Beauchemin KA, Schwartzkopf-Genswein KS, McKinnon JJ, McAllister TA. Effect of rapid or gradual grain adaptation on subacute acidosis and feed intake by feedlot cattle. J Anim Sci. 2005;83(5):1116–32.View ArticlePubMedGoogle Scholar
  3. Nocek JE. Bovine acidosis: implications on laminitis. J Dairy Sci. 1997;80(5):1005–28.View ArticlePubMedGoogle Scholar
  4. Ametaj BN, Zebli Q, Iqbal S. Nutrition, microbiota, and endotoxin-related diseases in dairy cows. R Bras Zootec. 2010;39:433–44.View ArticleGoogle Scholar
  5. Gressley TF, Hall MB, Armentano LE. Ruminant nutrition symposium: productivity, digestion, and health responses to hindgut acidosis in ruminants. J Anim Sci. 2011;89(4):1120–30.View ArticlePubMedGoogle Scholar
  6. Owens FN, Secrist DS, Hill WJ, Gill DR. Acidosis in cattle: a review. J Anim Sci. 1998;76(1):275–86.View ArticlePubMedGoogle Scholar
  7. Gäbel G, Bestmann M, Martens H. Influences of diet, short-chain fatty-acids, lactate and chloride on bicarbonate movement across the reticulorumen wall of sheep. J Vet Med A. 1991;38(7):523–9.View ArticleGoogle Scholar
  8. Penner GB, Aschenbach JR, Gäbel G, Rackwitz R, Oba M. Epithelial capacity for apical uptake of short chain fatty acids is a key determinant for intraruminal pH and the susceptibility to subacute ruminal acidosis in sheep. J Nutr. 2009;139(9):1714–20.View ArticlePubMedGoogle Scholar
  9. Aschenbach JR, Penner GB, Stumpff F, Gäbel G. Ruminant nutrition symposium: role of fermentation acid absorption in the regulation of ruminal pH. J Anim Sci. 2011;89(4):1092–107.View ArticlePubMedGoogle Scholar
  10. Castillo-Lopez E, Wiese BI, Hendrick S, McKinnon JJ, McAllister TA, Beauchemin KA, Penner GB. Incidence, prevalence, severity, and risk factors for ruminal acidosis in feedlot steers during backgrounding, diet transition, and finishing. J Anim Sci. 2014;92(7):3053–63.View ArticlePubMedGoogle Scholar
  11. Schlau N, Guan LL, Oba M. The relationship between rumen acidosis resistance and expression of genes involved in regulation of intracellular pH and butyrate metabolism of ruminal epithelial cells in steers. J Dairy Sci. 2012;95(10):5866–75.View ArticlePubMedGoogle Scholar
  12. Mohammed R, Stevenson DM, Weimer PJ, Penner GB, Beauchemin KA. Individual animal variability in ruminal bacterial communities and ruminal acidosis in primiparous Holstein cows during the periparturient period. J Dairy Sci. 2012;95(11):6716–30.View ArticlePubMedGoogle Scholar
  13. Humer E, Ghareeb K, Harder H, Mickdam E, Khol-Parisini A, Zebeli Q. Peripartal changes in reticuloruminal pH and temperature in dairy cows differing in the susceptibility to subacute rumen acidosis. J Dairy Sci. 2015;98(12):8788–99.View ArticlePubMedGoogle Scholar
  14. Chen YH, Penner GB, Li MJ, Oba M, Guan LL. Changes in bacterial diversity associated with epithelial tissue in the beef cow rumen during the transition to a high-grain diet. Appl Environ Microbiol. 2011;77(16):5770–81.View ArticlePubMedPubMed CentralGoogle Scholar
  15. Steele MA, Vandervoort G, AlZahal O, Hook SE, Matthews JC, McBride BW. Rumen epithelial adaptation to high-grain diets involves the coordinated regulation of genes involved in cholesterol homeostasis. Physiol Genomics. 2011;43(6):308–16.View ArticlePubMedGoogle Scholar
  16. Steele MA, Dionissopoulos L, AlZahal O, Doelman J, McBride BW. Rumen epithelial adaptation to ruminal acidosis in lactating cattle involves the coordinated expression of insulin-like growth factor-binding proteins and a cholesterolgenic enzyme. J Dairy Sci. 2012;95(1):318–27.View ArticlePubMedGoogle Scholar
  17. Oba M, Mewis JL, Zhining Z. Effects of ruminal doses of sucrose, lactose, and corn starch on ruminal fermentation and expression of genes in ruminal epithelial cells. J Dairy Sci. 2015;98(1):586–94.View ArticlePubMedGoogle Scholar
  18. Penner GB, Taniguchi M, Guan LL, Beauchemin KA, Oba M. Effect of dietary forage to concentrate ratio on volatile fatty acid absorption and the expression of genes related to volatile fatty acid absorption and metabolism in ruminal tissue. J Dairy Sci. 2009;92(6):2767–81.View ArticlePubMedGoogle Scholar
  19. Kleen JL, Hooijer GA, Rehage J, Noordhuizen JPTM. Subacute ruminal acidosis (SARA): a review. J Vet Med A Physiol Pathol Clin Med. 2003;50(8):406–14.View ArticlePubMedGoogle Scholar
  20. Gozho GN, Plaizier JC, Krause DO, Kennedy AD, Wittenberg KM. Subacute ruminal acidosis induces ruminal lipopolysaccharide endotoxin release and triggers an inflammatory response. J Dairy Sci. 2005;88(4):1399–403.View ArticlePubMedGoogle Scholar
  21. Penner GB, Oba M, Gäbel G, Aschenbach JR. A single mild episode of subacute ruminal acidosis does not affect ruminal barrier function in the short term. J Dairy Sci. 2010;93(10):4838–45.View ArticlePubMedGoogle Scholar
  22. Gerneke WH. Langerhans cells in epithelium of bovine forestomach - their role in primary immune response. J S Afr Vet Assoc. 1977;48(3):187–92.PubMedGoogle Scholar
  23. Trevisi E, Amadori M, Riva F, Bertoni G, Bani P. Evaluation of innate immune responses in bovine forestomachs. Res Vet Sci. 2014;96(1):69–78.View ArticlePubMedGoogle Scholar
  24. Liu J, Bian G, Sun D, Zhu W, Mao S. Starter feeding altered ruminal epithelial bacterial communities and some key immune-related genes’ expression before weaning in lambs. J Anim Sci. 2017;95(2):910–21.PubMedGoogle Scholar
  25. Lubbers R, van Essen MF, van Kooten C, Trouw LA. Production of complement components by cells of the immune system. Clin Exp Immunol. 2017;188(2):183–94.View ArticlePubMedGoogle Scholar
  26. Rakoff-Nahoum S, Paglino J, Eslami-Varzaneh F, Edberg S, Medzhitov R. Recognition of commensal microflora by toll-like receptors is required for intestinal homeostasis. Cell. 2004;118(2):229–41.View ArticlePubMedGoogle Scholar
  27. Uehara A, Fujimoto Y, Fukase K, Takada H. Various human epithelial cells express functional toll-like receptors, NOD1 and NOD2 to produce anti-microbial peptides, but not proinflammatory cytokines. Mol Immunol. 2007;44(12):3100–11.View ArticlePubMedGoogle Scholar
  28. Macpherson AJ, Harris NL. Interactions between commensal intestinal bacteria and the immune system. Nat Rev Immunol. 2004;4:478–85.View ArticlePubMedGoogle Scholar
  29. Chen YH, Oba M, Guan LL. Variation of bacterial communities and expression of toll-like receptor genes in the rumen of steers differing in susceptibility to subacute ruminal acidosis. Vet Microbiol. 2012;159(3–4):451–9.View ArticlePubMedGoogle Scholar
  30. Sur S, Agrawal DK. Phosphatases and kinases regulating CDC25 activity in the cell cycle: clinical implications of CDC25 overexpression and potential treatment strategies. Mol Cell Biochem. 2016;416(1):33–46.View ArticlePubMedPubMed CentralGoogle Scholar
  31. Lobjois V, Jullien D, Bouche JP, Ducommun B. The polo-like kinase 1 regulates CDC25B-dependent mitosis entry. Biochim Biophys Acta. 2009;1793(3):462–8.View ArticlePubMedGoogle Scholar
  32. Durlacher CT, Li ZL, Chen XW, He ZX, Zhou SF. An update on the pharmacokinetics and pharmacodynamics of alisertib, a selective aurora kinase a inhibitor. Clin Exp Pharmacol Physiol. 2016;43(6):585–601.View ArticlePubMedGoogle Scholar
  33. Levine AJ. p53, the cellular gatekeeper for growth and division. Cell. 1997;88(3):323–31.View ArticlePubMedGoogle Scholar
  34. Monte M, Benetti R, Buscemi G, Sandy P, Del Sal G, Schneider C. The cell cycle-regulated protein human GTSE-1 controls DNA damage-induced apoptosis by affecting p53 function. J Biol Chem. 2003;278(32):30356–64.View ArticlePubMedGoogle Scholar
  35. Nelson DR, Koymans L, Kamataki T, Stegeman JJ, Feyereisen R, Waxman DJ, Waterman MR, Gotoh O, Coon MJ, Estabrook RW, Gunsalus IC, Nebert DW. P450 superfamily: update on new sequences, gene mapping, accession numbers and nomenclature. Pharmacogenetics. 1996;6(1):1–42.View ArticlePubMedGoogle Scholar
  36. Höög JO, Strömberg P, Hedberg JJ, Griffiths WJ. The mammalian alcohol dehydrogenases interact in several metabolic pathways. Chem Biol Interact. 2003;143:175–81.View ArticlePubMedGoogle Scholar
  37. Miyauchi S, Gopal E, Fei YJ, Ganapathy V. Functional identification of SLC5A8, a tumor suppressor down-regulated in colon cancer, as a Na+-coupled transporter for short-chain fatty acids. J Biol Chem. 2004;279(14):13293–6.View ArticlePubMedGoogle Scholar
  38. Ganapathy V, Thangaraju M, Gopal E, Martin PM, Itagaki S, Miyauchi S, Prasad PD. Sodium-coupled monocarboxylate transporters in normal tissues and in cancer. AAPS J. 2008;10(1):193–9.View ArticlePubMedPubMed CentralGoogle Scholar
  39. Britton R, Krehbiel C. Nutrient metabolism by gut tissues. J Dairy Sci. 1993;76(7):2125–31.View ArticlePubMedGoogle Scholar
  40. Kristensen NB, Harmon DL. Splanchnic metabolism of volatile fatty acids absorbed from the washed reticulorumen of steers. J Anim Sci. 2004;82(7):2033–42.View ArticlePubMedGoogle Scholar
  41. Lane MA, Baldwin RVI, Jesse BW. Developmental changes in ketogenic enzyme gene expression during sheep rumen development. J Anim Sci. 2002;80(6):1538–44.View ArticlePubMedGoogle Scholar
  42. Hertzel AV, Bernlohr DA. The mammalian fatty acid-binding protein multigene family: molecular and genetic insights into function. Trends Endocrinol Metab. 2000;11(5):175–80.View ArticlePubMedGoogle Scholar
  43. Wang XQ, Yang K, He YS, Lu L, Shen WF. Receptor mediated elevation in FABP4 levels by advanced glycation end products induces cholesterol and triacylglycerol accumulation in THP-1 macrophages. Lipids. 2011;46(6):479–86.View ArticlePubMedGoogle Scholar
  44. Harvatine KJ, Bauman DE. SREBP1 and thyroid hormone responsive spot 14 (S14) are involved in the regulation of bovine mammary lipid synthesis during diet-induced milk fat depression and treatment with CLA. J Nutr. 2006;136(10):2468–74.View ArticlePubMedGoogle Scholar
  45. Viturro E, Koenning M, Kroemer A, Schlamberger G, Wiedemann S, Kaske M, Meyer HHD. Cholesterol synthesis in the lactating cow: induced expression of candidate genes. J Steroid Biochem Mol Biol. 2009;115(1–2):62–7.View ArticlePubMedGoogle Scholar
  46. Dempsey ME. Regulation of steroid biosynthesis. Annu Rev Biochem. 1974;43:967–90.View ArticlePubMedGoogle Scholar
  47. Baldwin RL. Use of isolated ruminal epithelial cells in the study of rumen metabolism. J Nutr. 1998;128(2):293S–6S.View ArticlePubMedGoogle Scholar
  48. Simons K, Toomre D. Lipid rafts and signal transduction. Nat Rev Mol Cell Biol. 2000;1(1):31–9.View ArticlePubMedGoogle Scholar
  49. Glavinas H, Krajcsi P, Cserepes J, Sarkadi B. The role of ABC transporters in drug resistance, metabolism and toxicity. Curr Drug Deliv. 2004;1(1):27–42.View ArticlePubMedGoogle Scholar
  50. Sahoo S, Aurich MK, Jonsson JJ, Thiele I. Membrane transporters in a human genonne-scale metabolic knowledgebase and their implications for disease. Front Physiol. 2014;5:91.View ArticlePubMedPubMed CentralGoogle Scholar
  51. Denis M, Bissonnette R, Haidar B, Krimbou L, Bouvier M, Genest J. Expression, regulation, and activity of ABCA1 in human cell lines. Mol Genet Metab. 2003;78(4):265–74.View ArticlePubMedGoogle Scholar
  52. Jessup W, Gelissen IC, Gaus K, Kritharides L. Roles of ATP binding cassette transporters A1 and G1, scavenger receptor BI and membrane lipid domains in cholesterol export from macrophages. Curr Opin Lipidol. 2006;17(3):247–57.View ArticlePubMedGoogle Scholar
  53. Chun YS, Shin S, Kim Y, Cho H, Park MK, Kim TW, Voronov SV, Di Paolo G, Suh BC, Chung S. Cholesterol modulates ion channels via down-regulation of phosphatidylinositol 4,5-bisphosphate. J Neurochem. 2010;112(5):1286–94.View ArticlePubMedGoogle Scholar
  54. Dopico AM, Bukiya AN, Singh AK. Large conductance, calcium- and voltage-gated potassium (BK) channels: regulation by cholesterol. Pharmacol Ther. 2012;135(2):133–50.View ArticlePubMedPubMed CentralGoogle Scholar
  55. Liao JK, Laufs U. Pleiotropic effects of statins. Annu Rev Pharmacol. 2005;45:89–118.View ArticleGoogle Scholar
  56. Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P. Molecular biology of the cell. 5th ed. New York: Garland Science; 2007.Google Scholar
  57. Phillips T. Regulation of transcription and gene expression in eukaryotes. Nat Educ. 2008;1(1):199.Google Scholar
  58. Penner GB, Beauchemin KA, Mutsvangwa T. An evaluation of the accuracy and precision of a stand-alone submersible continuous ruminal pH measurement system. J Dairy Sci. 2006;89(6):2132–40.View ArticlePubMedGoogle Scholar
  59. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.View ArticlePubMedPubMed CentralGoogle Scholar
  60. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.View ArticlePubMedGoogle Scholar
  61. Benjamini Y, Hochberg Y. Controlling the false discovery rate-a practical and powerful approach to multiple testing. J R Statist Soc B. 1995;57(1):289–300.Google Scholar
  62. Mi H, Muruganujan A, Casagrande JT, Thomas PD. Large-scale gene function analysis with the PANTHER classification system. Nat Protoc. 2013;8(8):1551–66.View ArticlePubMedGoogle Scholar
  63. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.View ArticleGoogle Scholar
  64. Rémond D, Ortigues I, Jouany JP. Energy substrates for the rumen epithelium. Proc Nutr Soc. 1995;54(1):95–105.View ArticlePubMedGoogle Scholar
  65. Koboldt DC, Zhang QY, Larson DE, Shen D, McLellan MD, Lin L, Miller CA, Mardis ER, Ding L, Wilson RK. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 2012;22(3):568–76.View ArticlePubMedPubMed CentralGoogle Scholar

Copyright

© The Author(s). 2018

Advertisement

BMC Genomics

ISSN: 1471-2164

Contact us