Transcriptome profiling of the rumen epithelium of beef cattle differing in residual feed intake

Background Feed efficient cattle consume less feed and produce less environmental waste than inefficient cattle. Many factors are known to contribute to differences in feed efficiency, however the underlying molecular mechanisms are largely unknown. Our study aimed to understand how host gene expression in the rumen epithelium contributes to differences in residual feed intake (RFI), a measure of feed efficiency, using a transcriptome profiling based approach. Results The rumen epithelial transcriptome from highly efficient (low (L-) RFI, n = 9) and inefficient (high (H-) RFI, n = 9) Hereford x Angus steers was obtained using RNA-sequencing. There were 122 genes differentially expressed between the rumen epithelial tissues of L- and H- RFI steers (p < 0.05) with 85 up-regulated and 37 down-regulated in L-RFI steers. Functional analysis of up-regulated genes revealed their involvement in acetylation, remodeling of adherens junctions, cytoskeletal dynamics, cell migration, and cell turnover. Additionally, a weighted gene co-expression network analysis (WGCNA) identified a significant gene module containing 764 genes that was negatively correlated with RFI (r = −0.5, p = 0.03). Functional analysis revealed significant enrichment of genes involved in modulation of intercellular adhesion through adherens junctions, protein and cell turnover, and cytoskeletal organization that suggest possible increased tissue morphogenesis in the L-RFI steers. Additionally, the L-RFI epithelium had increased expression of genes involved with the mitochondrion, acetylation, and energy generating pathways such as glycolysis, tricarboxylic acid cycle, and oxidative phosphorylation. Further qPCR analysis of steers with different RFI (L-RFI, n = 35; M-RFI, n = 34; H-RFI, n = 35) revealed that the relative mitochondrial genome copy number per cell of the epithelium was positively correlated with RFI (r = 0.21, p = 0.03). Conclusions Our results suggest that the rumen epithelium of L-RFI (efficient) steers may have increased tissue morphogenesis that possibly increases paracellular permeability for the absorption of nutrients and increased energy production to support the energetic demands of increased tissue morphogenesis compared to those of H-RFI (inefficient) animals. Greater expression of mitochondrial genes and lower relative mitochondrial genome copy numbers suggest a greater rate of transcription in the rumen epithelial mitochondria of L-RFI steers. Understanding how host gene expression profiles are associated with RFI could potentially lead to identification of mechanisms behind this trait, which are vital to develop strategies for the improvement of cattle feed efficiency. Electronic supplementary material The online version of this article (doi:10.1186/s12864-016-2935-4) contains supplementary material, which is available to authorized users.


Background
Feed is one of the major costs in beef cattle production, accounting for approximately 70 % of total expenditures [1]. As feed price continues to rise due to factors such as human population growth and decreasing amount of arable land [2], beef cattle feed efficiency has become an increasingly important contributor to reduce the cost of production. Recent efforts in the selection of feed efficient cattle have been focused on residual feed intake (RFI), a measure of feed efficiency that is independent of growth and body weight [3]. The RFI of an individual is calculated as the difference between actual feed intake and expected feed intake based on maintenance and growth requirements [4]. Animals with low (negative) RFI (L-RFI) are considered feed efficient because they consume less feed than expected whereas animals with high (positive) RFI (H-RFI) are feed inefficient due to consuming more feed than expected [4]. Selection for L-RFI animals would not only reduce cost of production, but also decrease the environmental footprint because they produce 28 and 15 % less methane [5] and manure [6], respectively, than H-RFI animals.
Variation in RFI may be due to differences in many biological processes that are influenced by genetic and environmental factors. Physiological processes that have been found to contribute to variation in RFI include those associated with feed intake, digestion, metabolism, activity, and stress [7,8]. Genetic selection and improvement of RFI is possible since RFI is a moderately heritable trait (0.28 to 0.45) [3,9,10], however, the molecular mechanisms underlying RFI are largely unknown. Because variation in RFI is due to differences in many processes, multiple genes from various pathways are likely involved [11][12][13]. Discovering differences at the molecular level may potentially lead to identification of biomarkers for use in the selection of L-RFI (feed efficient) animals.
In ruminant animals, such as cattle, the rumen may play an important role in feed efficiency. It is the main site of microbial feed digestion and fermentation that produces short-chain fatty acids (SCFAs), which are a major energy source for cattle accounting for up to 80 % of total energy requirements [14,15]. SCFAs are mainly absorbed from the rumen [16] and metabolized within the rumen epithelium, liver and peripheral tissues to generate energy for use in maintenance, growth and activity [17,18]. Studies have indicated that SCFAs are absorbed by both simple diffusion [16] and protein-mediated transport through bicarbonate -dependent and -independent mechanisms [19]. Moreover, ketogenic enzymes within the mitochondria of the rumen epithelium play a crucial role in the metabolism of SCFAs, to produce ketone bodies that are an important circulating source of energy throughout the body [20,21].
Differences in the ability of the rumen epithelium to absorb and metabolize SCFAs may be associated with variation in feed efficiency. Therefore, the objective of this study was to investigate differences in gene expression between the rumen epithelial tissues of feed efficient and inefficient beef cattle using RNA-sequencing (RNA-seq) technology and bioinformatics. We hypothesized that the rumen epithelium of L-RFI (efficient) animals would have greater expression of genes involved in absorption and metabolism of energetic substrates such as SCFAs.

Animals and rumen tissue collection
Hereford x Aberdeen Angus hybrid steers (n = 175) were born in the spring of 2006 and raised at the University of Alberta Roy Berg Kinsella Research Station (Alberta, Canada) in accordance with the Canadian Council of Animal Care guidelines [22]. The University of Alberta Animal Care and Use Committee for Livestock approved the experimental protocols (Moore- . Steers were under feedlot conditions while feed intake data was collected using the GrowSafe automated feeding system (GrowSafe Systems Ltd., Airdrie, Alberta, Canada). The finishing diet was fed for 90 days and was composed of 56.7 % barley, 28.3 % oats, 10 % alfalfa pellets, and 5 % feedlot supplement (32 % crude protein beef supplement containing Rumensin (400 mg/kg) and 1.5 % canola oil). After slaughter, a 4-cm 2 piece of rumen tissue was obtained from the central region of the ventral sac and rinsed with sterilized PBS buffer (pH = 6.8) before being placed in a 50 mL tube containing RNAlater solution (Invitrogen, Carlsbad, CA). The samples were then stored at −80°C until further processing.

RNA extraction and sequencing
After thawing, rumen tissue was placed on a petri dish and papillae (~80 mg) were obtained using sterile scissors and scalpels. The total RNA was extracted from the papillae using mirVana kit (Ambion, Austin, TX) according to the manufacturer's instructions. RNA integrity of the total RNA samples was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) and Nanodrop 2000c spectrophotometer (Thermo Scientific, Wilmington, DE) was used to measure the concentration and purity. Only samples with RNA integrity greater than 7 were subjected to RNA-seq library construction.
For each of the 18 samples, a cDNA library was generated from 100 ng of total RNA following the protocol of the TruSeq RNA Sample Preparation v2 kit (Illumina, San Diego, CA). First, total RNA was fragmented and cDNA was synthesized through reverse transcription. The resulting double-stranded cDNA was then subjected to end repair and 3'-end adenylation before ligation of index adapters. Libraries were amplified by 15 cycles of PCR, then validated with the Agilent 2200 TapeStation (Agilent Technologies) and quantified with a Qubit fluorometer using a Qubit dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA). The individual indexed libraries were then pooled and sequenced in five lanes on the Illumina HiSeq 2000 system at the McGill University and Genome Quebec Innovation Centre (Quebec, Canada) to obtain high-quality, 100-bp paired-end reads (Average Phred quality score ≥ 33).

Differential gene expression and principal component analysis
Differential gene expression analysis was performed using the Bioconductor (v3.0) [28] package DESeq2 (v1.6.2) [29] in the R statistical software program (v3.1.2) [30]. Firstly, DESeq2 normalized gene count data using the relative log expression method based on "size factors" to account for RNA-seq library size differences and dispersion estimates were calculated. Then, pairwise comparison of expression was made between the L-RFI and H-RFI group for every gene based on a negative binomial model. Fold changes were obtained along with their associated p-values. The Benjamini Hochberg method (B-H) was used to control the false discovery rate (FDR) by adjusting p-values to correct for multiple testing. A gene was defined as differentially expressed (DE) between the L-RFI and H-RFI group if it had a B-H adjusted p-value (FDR) less than 0.05.
DESeq2 was also used for Principal Component Analysis (PCA) to cluster samples based on gene expression data.
Firstly, count data was normalized based on library size and regularized log transformed. The top 500 genes with the most variance (i.e. Largest standard deviations) were compared between samples to calculate principal components in order to view clustering in two dimensions.
Weighted Gene Co-expression Network Analysis (WGCNA) The WGCNA R software package (v1.41.1) [31] was used to identify modules containing genes that are coexpressed and correlated with the RFI trait. Gene expression data was firstly filtered by keeping genes with at least one read in 50 % of the samples, then normalized by count per million reads (CPM) and Log2 transformed (Log2(CPM + 1)). Removal of the very lowly expressed genes resulted in 14,916 genes for input to WGCNA. Unsigned, weighted correlation network construction and module detection was performed using the automatic one-step function, blockwiseModules. The resulting gene modules were assigned color names by the software and Module-Trait relationships were calculated by Pearson correlation. Modules with statistically significant (p < 0.05) correlations were selected for downstream analysis as potential biologically interesting modules associated with the RFI trait. Additionally, a hierarchical clustering dendrogram of samples and trait heatmap was created using unfiltered, Log2CPM data. Sample clustering was based on Euclidean distance and the traits used in the heatmap were residual feed intake (RFI), dry matter intake (DMI), metabolizable energy intake (MEI), birth weight (BirthWT), weaning weight (WeanWT), metabolic mid-weight (MWT), end weight (EndWT), carcass weight (CWT), and average daily gain (ADG).

Functional gene annotation
Ensembl gene ID lists were converted to gene names and symbols using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) [32,33] and the R/Bioconductor package, biomaRt (v2.22.0) [34]. DAVID and Ingenuity Pathway Analysis (IPA; QIAGEN Redwood City, www.qiagen.com/ingenuity) were used for functional annotation of DE genes and WGCNA gene modules. In DAVID, we examined genes for enrichment in functional categories such as SP_PIR_Keywords and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [35,36]. To correct for multiple testing in DAVID, p-values were adjusted using B-H correction. Biological terms were considered significant if the adjusted p-value was less than 0.05. The genes input to IPA were compared to the manually curated Ingenuity knowledge base to find significantly enriched canonical pathways, and cellular and molecular processes. Statistical significance (p < 0.05) was calculated using the right-tailed Fischer's Exact Test and B-H corrected for multiple testing.
The identified core transcripts (i.e. genes with at least one read in every sample) were input to the Protein Annotation Through Evolutionary Relationship (PANTHER) [37] gene list analysis tool to determine the core functions of the rumen epithelium. Functional classification of genes in categories such as molecular function and biological process was based on Gene Ontology (GO) annotations. or obtained through literatures. The specificity of primers was checked with BLAST (NCBI) and the UCSC In-Silico PCR program [38]. Table 1 lists the housekeeping gene (endogenous control) and the five DE genes selected for qPCR validation along with their GenBank accession number, primer sequences, product length and functional pathway.
All samples were analyzed in triplicate during qPCR using the StepOnePlus Real-Time PCR System (Applied Biosystems). Each reaction contained 10 μL Fast SYBR® Green Master Mix (Applied Biosystems), 1 μL forward primer (20 pmol/μL), 1 μL reverse primer (20 pmol/μL), 7 μL nuclease-free water, and 1 μL cDNA template (2.5 ng/μL). The amplification program consisted of preincubation at 95°C for 20 s followed by 40 cycles of denaturation at 95°C for 3 s and annealing/extension at 64°C for 30 s. Melting curve analysis was performed to confirm single product amplification. The expression of three housekeeping genes (RPS18 (Ribosomal Protein S18), GAPDH (Glyceraldehyde 3-phosphate dehydrogenase), SUZ12 (SUZ12 polycomb repressive complex 2 subunit)) was tested in this study using NormFinder [39] and PROC MIXED procedure of SAS (v9.2; SAS Institute Inc., Cary, NC). We chose SUZ12 because it was the most stably expressed between and within the two groups. For each sample, the expression of each targeted gene was presented as ΔCt, which was obtained by subtracting the average Ct (Cycle threshold) of the SUZ12 housekeeping gene from the average Ct of the target gene. A low ΔCt indicates higher expression and vice versa. An unpaired t-test in SAS (SAS Institute Inc.) was used to analyze whether a gene showed significant difference between RFI groups.

DNA Extraction
Total DNA was extracted and used to quantify the relative mitochondrial genome copy numbers. The DNA was extracted from the rumen papillae (~80 mg) of 104 animals (35 L-RFI, 34 M-RFI, 35 H-RFI) using a bead beating, phenol-chloroform method as described by Li et al., 2012 [40]. Briefly, each sample was placed in a 2-mL microcentrifuge tube containing zirconium beads (0.3 g; 0.1 mm diameter) and washed twice with 1 mL of TN150 buffer (10 mM Tris-HCl [pH 8.0], 150 mM NaCl) by vortexing and centrifugation at 14,600 g for 5 min at 4°C. The tissue was then resuspended in 1 mL of TN150 and the Mini-BeadBeater-8 (BioSpec products, Bartlesville, OK) was used to mechanically disrupt cells at 5000 rpm for 3 min. After bead beating, samples were immediately placed on ice and washed Rev: TGGGCCTGCACACAAGAATG two times by phenol and chloroform-isoamyl alcohol (24:1). The DNA was then precipitated at −20°C overnight using cold 100 % ethanol and 3 M sodium acetate. After precipitation, DNA was washed two times with 70 % ethanol and dissolved in 30 μL nuclease-free water. The DNA concentration and quality was measured using Nanodrop (Thermo Scientific).
qPCR determination of relative mitochondrial DNA copy number per cell The relationship between RFI and the relative copy number of mitochondrial DNA (mtDNA) per cell of the rumen epithelium was performed using correlation analysis. The ratio of mtDNA to nuclear DNA (nDNA) was used to determine the relative copy number of mtDNA per cell. We selected the ND1 (NADH dehydrogenase subunit 1) gene as the mtDNA target (GenBank accession number: NC_006853, Location: 3101-4056) and the single copy, DDX3Y (DEAD box polypeptide 3, Y-linked) gene as the nDNA target (GenBank accession number: AC_000187, Location: 3458510-3468728,). The ND1 primer sequences used to amplify a 200 bp product were 5-GAACCAC TACGACCCGCTACA-3 (forward), 5-GAGTTGGAAGC TCAGCCTGATC-3 (reverse) and the DDX3Y primer sequences used to amplify a 136 bp product were 5-ATCG TGGGCGGAATGAGTGT-3 (forward), 5-CTTGGTGGA AGCGGTTTTGA-3 (reverse) [41]. For qPCR, a standard curve was generated for both the mtDNA and nDNA. First, PCR was performed on pooled DNA to amplify mtDNA and nDNA separately. The PCR mix contained 50 ng/μL DNA along with 10X PCR buffer, 50 mM MgCl 2 , 10 mM dNTP, 20 pmol of each primer, Taq polymerase, and nuclease-free water to a final volume of 50 μL. The PCR program consisted of 1 cycle 94°C for 5 min, 30 cycles of 94°C for 30 s, 62°C for 30 s, 72°C for 30 s and 1 cycle of 72°C for 7 min with hold at 4°C. Each PCR product was then subjected to agarose gel electrophoresis in order to excise DNA bands for purification using the QIAquick Gel Extraction Kit (Qiagen). The concentration of the purified mtDNA and nDNA was measured by Nanodrop (Thermo Scientific) and an eight-point standard curve with 10 fold dilutions was created for each. DNA copy number of the standard curve points was calculated as described by Speicher and Johnson (2012) [42]. Expression of mtDNA and nDNA within 104 samples (35 L-RFI, 34 M-RFI, 35 H-RFI) was then analyzed by qPCR along with their respective standard curve using an annealing temperature of 62°C. DNA copy number within samples was extrapolated from the standard curve and all samples and standard curves were analyzed in triplicate. A one-way ANOVA was performed in SAS (SAS Institute Inc.) to compare RFI groups. Also, a Pearson correlation analysis of relative mitochondrial DNA copy number per cell and RFI was performed in SAS (SAS Institute Inc.) and R [30].

Identification of the rumen epithelial transcriptome
On average there were 42 million high-quality, paired reads generated per sample (mean ± SD = 41,991,418 ± 6,835,817) from RNA-seq and the overall read alignment rate to the bovine reference genome was 85 ± 3.72 %. The total number of genes expressed (genes with at least one read in at least one sample) within the rumen epithelium was 22,338 with an average of 14,677 ± 1062 genes per L-RFI sample and 16,239 ± 1326 per H-RFI sample. There were 11,284 core genes (genes with at least one read in every sample) among samples and of these core genes, 9700 can be annotated by PANTHER. The top molecular functions for the core genes were catalytic activity and binding, which included 30.7 and 30.1 % of the core genes, respectively. In the biological processes category, 32.6 % of the core genes were involved in cellular process and 47.5 % of the genes were associated with metabolic process (Fig. 1).

Comparison of transcriptome profiles between two RFI groups
The PCA plot (Fig. 2) shows that L-RFI and H-RFI gene expression profiles seem to fall into distinct clusters, with those from L-RFI animals grouped together. However, there was a L-RFI outlier sample (271 L) that did not cluster with the rest of the samples (Fig. 2), which may possibly be related to the animal's low ADG and EndWT (Fig. 3). Hierarchical clustering by WGCNA showed that all the L-RFI samples clustered together while the H-RFI samples did not (Fig. 3). When comparing animals with similar ADG and EndWT in the trait heat map, L-RFI individuals had lower RFI, DMI and MEI than H-RFI individuals. The means of the two groups for all traits (Additional file 2) were compared and there was no significant difference in BirthWT, WeanWT, MWT, EndWT, CWT, and ADG. However, the L-RFI group had significantly (p < 0.001) lower RFI, DMI and MEI by 3.82 kg/day, 3.86 kg/day and 4.68 Mcal/day, respectively.

Differential gene expression in rumen epithelium between low-and high-RFI cattle
There were 122 genes DE in the rumen epithelial tissues between L-and H-RFI animals (Additional file 3). Of these DE genes, 85 were up-regulated and 37 downregulated in the L-RFI group. DAVID annotation of the 122 DE genes found no significant KEGG pathways, however the biological term acetylation was significantly enriched (FDR < 0.05). A total of 30 DE genes were involved in the process of acetylation, 26 of which were up-regulated and 4 down-regulated in L-RFI animals ( Table 2). These genes have functions in protein synthesis (i.e. RPL10, RPS15, RPL36), protein degradation (i.e. PSMB6, UBC/UBA52, UBE2V1), cytoskeletal organization (i.e. DSTN, ACTB, TUBB5, TMSB10), energetic substrate generation (i.e. TECR, TPI1), and stress response (i.e. HSPB1, HSF1).
Ingenuity analysis of the 85 differentially up-regulated genes in L-RFI samples showed that the top significantly enriched canonical pathway was remodeling of the epithelial adherens junctions through endocytosis (i.e. DNM2), degradation (i.e. HGS), and cytoskeletal organization (i.e. ACTB, TUBB, TUBA4A) (Table 3 and Additional file 3). The other enriched canonical pathways were involved in functions such as protein synthesis (eIF2 signaling), cytoskeletal dynamics, cell growth, proliferation, apoptosis, and migration (signaling by Rho family GTPases, RhoGDI signaling, ephrin B signaling) ( Table 3). For the 37 DE genes that were down-regulated in L-RFI animals, there were no significant canonical pathways.

qPCR validation of selected DE genes
To show the reliability of the differential gene expression analysis of DESeq2, we further selected five DE genes for validation using qPCR. As shown in Table 4, of these five genes, four (TECR, COX8A, SLC25A39, PKM2) had significantly greater expression (i.e. Lower ΔCt, p < 0.05) in L-RFI while one (TPI1) tended to have significantly greater expression (p < 0.1). These results were in agreement with the differential expression analysis of DESeq2 in that the genes had greater expression in L-RFI rumen epithelium.    Identification of WGCNA gene co-expression modules correlated with RFI WGCNA was conducted to identify gene co-expression modules that are correlated with RFI. We identified a total of eight co-expressed gene modules with three of them significantly correlated with RFI (Fig. 4). The Brown (162 genes) and Green (133 genes) modules had a significant positive correlation with RFI (Brown module: r = 0.54, p = 0.02; Green module: r = 0.51, p = 0.03) while the Turquoise module (764 genes) had a significant negative correlation with RFI (r = −0.50, p = 0.03). Functional analysis of the genes enriched in the Green module did not result in any significantly enriched biological functions or pathways while there was only one significant biological term, which was muscle protein containing four genes in the Brown module. Many significant biological functions were identified for the turquoise module with IPA and DAVID analysis (Tables 5  and 6) and therefore this study focused on functional examination of the 764 genes enriched in this module, which had greater expression in the L-RFI epithelium.

Functional analysis of gene module negatively associated with RFI
When functional analysis was performed using DAVID, 20 biological functions and 6 KEGG pathways were significantly enriched for genes in the turquoise module ( Table 5). The most significant biological functions were acetylation (195 genes) and phosphoprotein (210 genes), which are involved in post-translational modification. Many of the other significant biological functions identified were also involved in protein related processes such as protein synthesis (i.e. ribosomal protein (44 genes), protein biosynthesis (20 genes)), degradation (i.e. proteasome (15 genes), threonine protease (7 genes)), folding (i.e. chaperone (17 genes)), and transport (i.e. transit peptide (35 genes), protein transport (24 genes)). In addition, biological terms associated with energy  . Asterisks indicate modules that are significantly (p < 0.05) correlated with RFI catabolism (i.e. ATP-binding (50 genes)) and anabolism (i.e. mitochondrion (47 genes)) were also identified. KEGG analysis of the 195 genes within the acetylation category indicated significant involvement of the acetylation related genes in ribosome and proteasome pathways. These pathways were also significantly enriched in the turquoise module along with oxidative phosphorylation, adherens junction, regulation of actin cytoskeleton, and spliceosome (Table 5). In the adherens junction pathway, there were 14 genes from the turquoise module that may regulate both the strong and weak adhesion between cells (Fig. 5).
Further functional analysis using IPA showed that the top five most significant molecular and cellular functions for the genes in the turquoise module were protein synthesis, cellular growth and proliferation, posttranslational modification, protein folding and cell death and survival (Table 6). Furthermore, the top canonical pathways were eIF2 signaling, regulation of eIF4 and p70S6K signaling, mTOR signaling and epithelial adherens junctions remodeling and signaling (Table 6).
Additionally, many genes from the turquoise module are involved in energy generating pathways such as glycolysis, tricarboxylic acid (TCA) cycle, and oxidative phosphorylation (Table 7). We also found three DE genes involved in glycolysis and three associated with oxidative phosphorylation, which were up-regulated in the L-RFI epithelium ( Table 7) (Table 7 and Additional file 3).

Relative mtDNA copy number and RFI
The transcriptome analysis results showed up-regulation of genes in the L-RFI epithelium involved in mitochondrial functions such as oxidative phosphorylation. Therefore, we speculated that L-RFI animals would have more relative mtDNA copy numbers (mtDNA/nDNA) in the rumen epithelium than H-RFI animals. To test this hypothesis, we investigated the relationship between RFI and relative mtDNA copy number using a larger population including M-RFI (n = 34) animals along with L-RFI (n = 35) and H-RFI (n = 35) animals. The correlation analysis showed a significant positive relationship (r = 0.21, p = 0.03) between RFI and the relative mtDNA copy number per cell in the rumen epithelium (Fig. 6). The  epithelium of L-RFI animals had 741 ± 35 (mean ± SD) relative mtDNA copy numbers per cell while the epithelium of M-and H-RFI animals had 1104 ± 89 and 1000 ± 80, respectively. There was a significantly greater copy number of mtDNA in the rumen epithelium of M-and H-RFI animals compared to L-RFI animals (p < 0.01). However, no significant difference in mtDNA copy number was observed between M-and H-RFI animals (p > 0.05).

Discussion
Improving cattle feed efficiency could greatly reduce feed expense and increase profit for producers. Our study of steers divergent in feed efficiency as evaluated by RFI showed that L-RFI individuals are more feed efficient compared to H-RFI individuals because they were able to consume 3.82 kg less dry matter per day for similar growth, body weight, and carcass weight. These results coincide with those of others, which found that selection of L-RFI steers with 3.77 kg per day lower feed intake than H-RFI steers may lead to savings of $45.60 per head over a 120-day finishing period and there was no correlation of RFI with growth, body weight and carcass weight [4]. The difference in feed consumption between L-and H-RFI steers for similar production illustrates that differences may exist in the efficiency of feed conversion into energetic substrates, absorption, energy production, and in the distribution of energy for growth, maintenance, and activity requirements. Ruminants utilize the substrates produced by the rumen microbes for their growth. Studies have shown that the rumen microbial community differs between L-and H-RFI steers and may impact feed digestion and fermentation [43]. Although, the rumen microbial community has been associated with divergence in RFI, it is unknown whether there are differences in the rumen tissue that contribute to variation in RFI.
Our study found that 47.5 % of the core genes among the rumen epithelial tissues of beef steers are involved in metabolic processes. Indeed, the rumen epithelium was known to be highly metabolic because it contains a high density of mitochondria in the stratum basale layer [44] and metabolizes a large proportion of absorbed SCFAs [45]. Through hierarchical clustering, we found that the whole transcriptome profiles of the rumen epithelial tissues of H-RFI (inefficient) steers did not cluster together while all L-RFI (efficient) steers did. This indicates that the rumen epithelial transcriptome of H-RFI animals are more variable and there are similarities among the gene expression profiles of L-RFI animals that may be associated with higher feed efficiency. We speculated that there was greater expression of genes involved in absorption of SCFAs from the rumen of L-RFI animals. However, we did not find evidence of differences in the expression of known or candidate SCFA transporters belonging to solute carrier (SLC) gene families such as SLC4A, SLC16A (monocarboxylate transporters), SLC21A, SLC22A, or SLC26A [19] between the L-and H-RFI animals in this study (data not shown). Although we did not observe differences in the expression of genes encoding SCFA transporters, our results do show evidence of differences in the expression of genes that are associated with paracellular permeability, which may impact simple diffusion of SCFAs and other nutrients. For example, there were 14 genes upregulated in the L-RFI epithelium involved in the adherens junction pathway, which is important in modulating cell-cell adhesion [46,47]. These genes encoded the plasma membrane cell adhesion molecules (CAMs), Ecadherin and nectin, as well as intracellular regulators of these molecules such as p120ctn, v-Src, IQGAP1, and ACTB. CAMs require attachment to the actin cytoskeleton in order to mediate cell adhesion [48]. Additionally, the E-cadherin CAM requires interaction with a key regulator, p120ctn, in order to promote strong adhesion between cells [49]. In our study, the L-RFI epithelium had increased expression of the p120ctn regulator, however, there was also increased expression of negative regulators such as v-Src and IQGAP1 that promote weak adhesion. Furthermore, we found several DE genes that were up-regulated in the L-RFI epithelium encoding Dynamin-2 (DNM2), Hepatocyte growth factor-regulated tyrosine kinase substrate (HGS), and cytoskeletal components (Beta-actin (ACTB), Tubulin beta-5 (TUBB5), and Tubulin alpha-4a (TUBA4A)) that are involved in remodeling of epithelial adherens junctions. This suggests there is increased regulation of intercellular adhesion strength in the L-RFI epithelium. We also found that regulation of the actin cytoskeleton was a pathway that was up-regulated in the L-RFI epithelium. The regulation of intercellular adhesion strength and the actin cytoskeleton are vital to many processes such as cell migration, tissue remodeling, and maintenance of tissue integrity [50]. Indeed, we found up-regulated DE genes (RHOG, CFL1, ACTB, MYL12B, GNB1, GNB2, MAPK1) in the L-RFI epithelium involved in several important signaling pathways regulating cell migration including signaling by Rho family GTPases, RhoGDI, and Ephrin B that mediate their effects on cell migration though modulation of intercellular adhesion and regulation of the actin cytoskeleton [51][52][53]. Overall, the increased expression of genes related to modulation of intercellular adhesion strength, cytoskeletal organization, and cell migration signaling pathways suggest there is greater cell mobility and dynamic remodeling within the L-RFI epithelium, which may create gaps between cells resulting in greater paracellular permeability for the absorption of nutrients. Therefore, the L-RFI epithelium may potentially have increased paracellular permeability to maximize nutrient absorption in order to compensate for lower feed intake/substrate availability. Future studies are needed to examine the intercellular spaces, cell migration, and nutrient absorption in the rumen epithelium to determine whether they are truly different between steers divergent in RFI.
In addition to greater cell mobility and dynamic remodeling, further evidence from our transcriptome analyses results support that the rumen epithelium of L-RFI animals may have greater tissue morphogenesis compared to H-RFI animals. Tissue morphogenesis is a process that involves changes to the number, shape, size, position, and gene expression of cells in order to develop and maintain a tissue [54]. For instance, we found upregulation of functions such as cellular growth, proliferation, and death in the L-RFI epithelium, which may indicate increased cell turnover. There was also upregulation of protein synthesis, degradation, folding, transport, and post-translational modification processes that suggest increased protein turnover. Future studies to measure rumen papillae density/surface area or rumen volume are necessary to determine whether the L-RFI epithelium has increased tissue morphogenesis and permeability due to a larger surface area to volume ratio. Greater tissue morphogenesis suggests that the L-RFI epithelium may have greater cell growth and proliferation. This may be due to increased paracellular permeability causing the epithelial cells to be exposed to nutrients at a faster rate.
It is speculated that the difference in the expression of genes involved in tissue morphogenesis could be also associated with the variation in nutrient availability and flux. L-RFI animals had lower feed intake, suggesting that these animals may have lower microbial products available for the animals. However, it has been reported that L-RFI animals tended to have a greater concentration of total SCFAs (p = 0.059) in the rumen fluid with significantly greater concentrations of butyrate (p < 0.001) and valerate (p = 0.006) [43]. High concentrations of nutrients are associated with changes in the rumen epithelium such as extensive sloughing of surface epithelial cells, and increased cell migration and proliferation [55,56]. Although the concentration of ruminal SCFAs was not measured for this study, the lack of difference in the expression of genes associated with SCFA oxidation and ketogenesis between L-and H-RFI suggest that the groups of animals may have similar SCFA levels in the rumen, although low RFI animals eat less. On the other hand, in L-RFI animals, nutrient flux may have been increased to maximize absorption because they have less substrate availability by consuming less feed. Regardless, future studies are needed to measure rate of nutrient absorption as well as the concentration of nutrients in the rumen content and epithelium to define their roles in affecting rumen tissue gene expression and in relation to RFI.
While there was no difference in the expression of SCFA metabolism genes, we found up-regulation of genes in the L-RFI epithelium involved in energy generating pathways such as glycolysis, TCA cycle, and oxidative phosphorylation. The majority of cellular energy, in the form of ATP, is generated through oxidative phosphorylation within the mitochondrial inner membrane [57]. Our results on increased expression of genes encoding mitochondrial proteins involved in oxidative phosphorylation in the L-RFI epithelium suggest that efficient steers may have increased energy production in the rumen epithelium. This may also indicate higher mitochondrial activity and/or greater mitochondrial genome copy numbers in the rumen epithelium of L-RFI steers compared to H-RFI steers. RFI was positively correlated (r = 0.21, p = 0.03), although very weakly, to the mitochondrial genome copy number per cell of the epithelium. Although there were fewer mitochondrial genomes in L-RFI animals, we observed increased expression of mitochondrial genes such as oxidative phosphorylation genes. This suggests that there is greater transcription from less genomes in the mitochondria of L-RFI animals. While we could not directly measure mitochondrial activity in our study, others have reported a more rapid rate of oxidative phosphorylation in the mitochondria of skeletal muscle and the liver in L-RFI beef cattle [58,59]. We speculate that the mitochondria in the L-RFI rumen may also have the same capability, which needs to be further studied. Overall, increased expression of glycolytic and oxidative phosphorylation genes in the L-RFI epithelium suggests increased energy production that provides a resource for the increased tissue morphogenesis occurring to maintain barrier integrity in the L-RFI epithelium, which is an energetically expensive process. Although expensive, we speculate the L-RFI animals have enough energy for other processes throughout the body although they eat less. Future studies on examination of epithelial metabolism and whole body energy balance are needed to support our speculation.
Additionally, our study showed significant enrichment of genes in the L-RFI epithelium associated with acetylation, which is a post-translational modification process. These genes were involved with various processes such as protein synthesis, protein degradation, cytoskeletal organization, metabolism of energetic substrates, and stress response. This suggests that acetylation modulates genes belonging to a wide range of cellular processes. Acetylation is known to be a metabolically sensitive protein modification process [60]. It is known that approximately 35 % of all mitochondrial proteins are acetylated and the majority of these acetylated proteins are involved in energy metabolism [61]. The expression level of genes encoding proteins that have the potential to become acetylated could be an indicator of the level of energy metabolism occurring within animals. Our results on the increased expression of mitochondrial and metabolic genes suggest greater activity of the enzymes and more energy flow in low RFI animals.

Conclusions
The RNA-seq based transcriptome analysis showed increased expression of genes in the L-RFI epithelium involved in modulation of adherens junctions, cytoskeletal organization, and cell migration signaling pathways that may cause large intercellular gaps and result in greater paracellular absorption of nutrients. Also, the L-RFI epithelium showed increased tissue morphogenesis as characterized by increased expression of genes involved in cell and protein turnover, modulation of intercellular adhesion strength, and cell migration, which may be due to increased absorption rate or a larger surface area to volume ratio. There was no difference in SCFA metabolism; however, there was increased expression of genes associated with glycolysis and oxidative phosphorylation in the L-RFI epithelium suggesting greater energy production. The increase in energy production provides a resource for the dynamic processes involved in maintaining tissue barrier integrity. A limitation of this study was that the epithelial tissue used for RNA-seq was obtained from only one location within the rumen. Pooling epithelial samples from multiple locations in the rumen is recommended to provide a better representation of the host