An examination of skeletal muscle and hepatic tissue transcriptomes from beef cattle divergent for residual feed intake 

Abstract

The selection of cattle with enhanced feed efficiency is of paramount importance with regard to reducing feed costs in the beef industry. Global transcriptome profiling was undertaken on liver and skeletal muscle biopsies from Simmental heifers and bulls divergent in residual feed intake (RFI) feed efficiency phenotype, in order to identify genes that may be associated with this trait. We identified 5 genes (adj.p<0.1) to be differentially expressed in skeletal muscle between high and low RFI heifers with all transcripts involved in oxidative phosphorylation and mitochondrial homeostasis. A total of 11 genes (adj.p<0. 1) were differentially expressed in liver tissue between high and low RFI bulls with differentially expressed genes related to amino and nucleotide metabolism as well as endoplasmic reticulum protein processing. No genes were identified as differentially expressed in either heifer liver or bull muscle analyses. Results from this study show a clear effect of gender to the underlying molecular control of RFI in cattle, which may be attributable to differences in the physiological age between heifers and bulls. Despite this we have highlighted a number of genes that may hold potential as molecular biomarkers for RFI cattle. 

Introduction

Global agriculture is currently faced with the ambitious challenge of feeding a rapidly increasing global population, expected to peak at 9.2 billion by 2075¹. This necessary increase in agricultural outputs must also be achieved within the current confines of arable land availability, thus it is essential that animal production systems become more efficient for the continued sustainability of the beef production sector. Animal feed can account for up to 75% of the variable costs in beef production systems; hence reduction of these costs is of paramount importance². Feed efficiency (dietary nutrient utilisation) in beef cattle is a trait of major economic importance³. Indeed, there is significant phenotypic and genetic variation among beef animals in their ability to convert dietary derived nutrients into saleable product^4,5,6. Thus, by improving feed efficiency it is possible to reduce feed intake in cattle while still maintaining growth and skeletal muscle gain and ultimately contributing to beef production profitability and sustainability. 

Residual feed intake (RFI), first described by Koch et al. (1963)⁷ is one such measure of feed efficiency, which can be defined as the difference between an animal’s actual versus its predicted feed intake based on average daily gain (ADG) and metabolic weight⁸. Through measuring an animal’s inherent RFI, feed-efficient animals which have a low-RFI value and consume less than expected as well as feed inefficient animals which have a high-RFI and consume more than expected value may be identified^3,9-15. Indeed, it has been shown that growing beef cattle divergent for RFI can consume up to 20% less feed than their counterparts for the same level of performance^{3,10,11,13-19}. This coupled with a moderate heritability estimate for RFI²⁰ provides a feasible method for effective utilisation of this trait in production systems through genomic selection processes. However, although moderately heritable, a challenge remains to reliably and cost-effectively identify feed-efficient, low-RFI animals and to proliferate their genetics through animal breeding programmes. For example, the primary impediment to genetic progress and adoption of selection strategies based on RFI is both the large-scale logistics and expense of measuring individual animal feed intake and body weight. 

In order to overcome the aforementioned limitations, studies have sought to uncover the underlying biology governing the trait, with the twin-goal for the identification of molecular biomarkers²¹. This not only provides an attractive alternative to direct measurement of dietary intake on large numbers of animals⁶, but it also allows for a better understanding of the biological mechanisms underlying RFI which is paramount for progressing genomic selection. However although results in the literature have reported roles for biological processes including lipogenesis and the immune system toward variation in RFI primarily through transcriptional profiling, there is a distinct lack of commonality of key genes contributing to RFI across these studies. This is undoubtedly due to the multi-faceted nature of RFI as well as the contribution of factors including breed, gender, stage of development as well as production management and dietary intake test period length to the trait²¹. For example a recent study that evaluated the molecular control of cattle divergent in RFI across three contrasting breed types, reported only 5 genes as commonly differentially expressed between high and low-RFI groups²². Furthermore, of these 5 genes only one, SCD, was differentially expressed in only one other molecular based RFI evaluation²³. This lack of commonality, although due to confounding experimental designs represents a major short-fall toward the progress of genomics selection. Equally there is a dearth of information in relation to the potential effect of gender on subsequent variation in RFI, particular under similar rearing conditions. Currently it is not known whether the same molecular mechanisms and genes are contributing to variation in RFI in heifers and bulls destined for beef production. Thus the objective of this study was to evaluate the underlying biology regulating RFI through transcriptome sequencing in Simmental heifers and bulls reared under the same conditions from birth. It is estimated that two thirds of the variation in RFI is due to variation in resting energy expenditure⁵, with both muscle and liver representing important metabolic tissues with skeletal muscle accounting for over 50% body weight¹² and liver accounting for 18-25% of total oxygen body consumption²⁴. Thus our efforts were focused towards examining the transcriptional alterations of these tissues between beef heifers and bulls divergent for high and low-RFI. 

Results

Animal performance

Details of live weight gain, feed intake, and animal performance are presented in Tables 1 for bulls and heifers. Briefly, bulls and heifers had a mean ADG of 1.8 kg and 1.3 kg and dry matter intake (DMI) of 9.3 kg and 9.2 kg respectively during the RFI measurement period. Residual feed intake averaged -0.04 and 0.04 and ranged from -1.92 to 1.04 kg of dry matter(DM)/d and -0.43 to 0.68 of DM/d for bulls and heifers respectively. Bulls and heifers ranked as high RFI consumed 10 and 15 % more than their low RFI counterparts (p<0.05), respectively. Bulls of high and low RFI did not differ (p>0.05) in initial bodyweight (BW), final BW or ADG and likewise for heifers. 

RNAseq read alignment and differential gene expression

RNA sequencing data was successfully created for all samples with approximately 22.4 million sequences per sample generated across all samples (skeletal muscle and hepatic tissue). On average approximately 85% of reads aligned to the bovine reference genome (UMD3.1) across both tissues analysed. A total of 5 genes (adj.p<0.1) were identified to be differentially expressed in skeletal muscle tissue between high and low RFI heifers with all 5 transcripts being up-regulated in the low RFI phenotype (Table 2). Figure 1 provides a visual representation of the correlation matrix pertaining to the differentially expressed genes (DEGs) from heifer skeletal muscle tissue, indicating all DEGs to be highly correlated with each other. A total of 11 genes (adj.p<0.1) were identified to be differentially expressed in hepatic tissue between high and low RFI bulls (Table 3) with 8 transcripts being up-regulated and 3 being down-regulated in the low RFI phenotype. Figure 2 provides a visual representation of the correlation matrix pertaining to the DEGs from bull hepatic tissue, highlighting the most correlated genes. Correlation analyses between DEGs and RFI for each heifer and bull group, identified mostly significant associations. For example DEGs from the heifer skeletal muscle analysis were all significantly negatively associated (p<0.05) with RFI, with the exception of cytb, which tended towards significance (p=0.0622). Similarly within the liver analyses of bulls divergent for RFI, 7 of the 11 DEGs were significantly associated with RFI (p<0.05), with DBP, MANF and GMPPB tending towards significance (p<0.1). ACTA2 expression was not significantly associated with RFI in the liver tissue of bulls (p=0.3078). Correlation results between DEGs and RFI value for each group are presented in Table 4. Among the (approximately) 7,000 genes detected in heifer hepatic tissue, none were found to be differentially expressed. Similarly, no differential expression was detected between RFI phenotypes in the (approximately) 6,000 genes expressed in skeletal muscle tissue from bulls. 

Differentially expressed genes and pathway analysis

Functional analysis of DEGs in hepatic tissue of bulls revealed amino sugar and nucleotide sugar metabolism as well as processes associated with the endoplasmic reticulum to be affected by divergence in RFI (Table 5). Functional analysis of heifer skeletal muscle DEGs indicated a role for oxidative phosphorylation and mitochondrial function towards divergence in RFI (Table 6). All DEG were successfully mapped to a molecular or biological pathway and/or category in the Ingenuity Pathway analysis (IPA) database. DEG were analysed and separated according to their biological function within IPA. The top canonical pathways affected by RFI in heifer skeletal muscle tissue and hepatic tissue of bulls are presented in Table 7. These included enriched pathways related to mitochondrial function and oxidative phosphorylation in heifers and Aldosterone Signalling, GDP-mannose Biosynthesis, Nuclear factor erythroid 2-related factor 2 (NRF2) mediated oxidative stress response, and Eukaryotic Initiation Factor 2 (EIF2) signalling in bulls. 

Discussion

Residual feed intake is a prime feed-efficiency trait to target within beef production systems; this is not only due to the moderate heritability of the trait but also due to its independence from production traits used to calculate it²¹. Furthermore, data from our own group has also shown that RFI is a repeatable trait in beef cattle¹¹. However although well suited toward genetic propagation through genomic selection, results on the underlying genes governing the trait are conflicting across studies^21,25. This is undoubtedly due to the multifaceted nature of the trait as well as the influence of both animal and management factors toward the RFI phenotype. Additional individual experimental parameters may also be contributing to the lack of consistency across results, including for example the measurement of the RFI phenotype and the length of the dietary intake test period²⁶. Furthermore RFI calculated within a specific population, as is the case in RFI transcriptional profiling studies, only reflects natural variation within that specific population, thus the level of divergence may be quite variable across different studies. In addition to these limitations, genomic selection models do not differentiation on the basis of gender, thus if molecular biomarkers are to be successfully employed for a trait it is essential to determine whether biological processes and specific key genes are regulating economically important traits such as RFI across differing gender types. Thus again the objective of this study was to evaluate any potential effects of differing gender to the underlying biological mechanisms regulating variation in RFI in both Simmental heifers and bulls. The animals used in the current study were from a purebred, well-characterised herd, reared as a contemporary group from birth with similar genetics, thus permitting a more equitable comparison of the effect of phenotypic RFI ranking. Our evaluations were focused towards both the liver and skeletal muscle tissue, the metabolic activities of which are both essential for overall body homeostasis and efficiency of an animal with both organs being highly abundant in mitochondria²⁷. Muscle accounts for approximately 50% of body mass and 25% of basal metabolic activity of an animal and plays an important role in resting energy expenditure²⁸. Additionally the liver is a highly oxidative organ accounting for 18-25% of total oxygen body consumption that is responsible for metabolising lipids, proteins and carbohydrates into biologically useful molecules²⁴. We hypothesised that due to the metabolic importance of these organs, that variation in feed efficiency and energy expenditure (measured here using RFI) is likely to be reflected in the transcriptome of tissue from these organs. Understanding the essential biological processes contributing to variation in RFI is critical to elucidating the genetic basis, for this trait. 

Although the RFI values pertaining to the animals used in this study showed clear significant divergence between high and low-RFI groups in both heifers and bulls, we failed to identify any DEGs within the liver tissue of the heifers and the skeletal muscle tissue of the bulls. Similarly, other studies evaluating the molecular control of RFI divergence in crossbred steers reported no DEGs following correction for multiple testing^29-31. Therefore the main observation of the current study is the inconsistent effect of RFI across (i) gender and (ii) tissue for DEG profiles despite a 10% and 15% difference in DMI between high and low RFI heifers and bulls, respectively, with no difference in ADG in the current study. While both genders were of similar age and were reared under the same conditions, it is apparent that both genders were at distinct physiological stages at the time of sample collection. Similarly in a targeted gene expression study, complementary to this current study, we observed a significant effect of gender in the expression of lipogenesis genes within the subcutaneous adipose tissue, this is despite there being no significant difference in subcutaneous fat measurements³². Typically heifers display earlier physiological maturity when compared to bulls of the same age, thus the identification of a gender effect in an energy storage accretion tissue such as adipose tissue was perhaps unsurprising, however the potential effect of physiological age and stage of maturity on metabolic tissues at the same age in the current study was potentially unexpected. Our results indicate that it will be unlikely that the key genes derived from the tissues examined in this study will be accurate predictors of genetic potential for RFI across gender. However this does not mean that other tissues or organs within the body may provide more comparable results and thus reliable biomarkers, based on animal age and not physiological age or stage of development. Despite this, we determined possible molecular mechanisms and biological functions influencing RFI in beef cattle within tissue and gender subgroups. These data highlight a relationship between RFI and the transcriptomic networks involved in mitochondrial function in the skeletal muscle of heifers and evidence for a relationship between RFI and the aldosterone signalling pathway and the NRF2 mediated oxidative stress pathway in hepatic tissue of bulls.

We identified 5 genes that were statistically significantly differentially expressed in skeletal muscle of high and low RFI heifers. Although we did not identify large numbers of DEGs, the level of transcriptional differences are consistent with the findings of other previously published RFI based studies^33-35. Among the 5 DEGs, all 5, COX1 (cytochrome c oxidase subunit 1), ND5 (NADH-dehydrogenase 5), ND6 (NADH-dehydrogenase 6), CYTB (cytochrome b), and COX3 (cytochrome c oxidase subunit 3), were up-regulated in low RFI heifers compared to their high RFI counterparts. Additionally all of these genes were significantly associated with RFI, with the exception of cytb, which only tended towards a significant association, further indicating a role for these genes to RFI variation in the skeletal muscle tissue of heifers. Interestingly all of these genes are components of the electron transport chain (ETC) in the mitochondrial inner membrane and are of major importance to overall energetic efficiency. ND5 and ND6 are subunits of the enzyme complex NADH dehydrogenase (ubiquinone) or complex I of the ETC. CYTB is the main subunit of coenzyme Q: cytochrome c–oxidoreductase or complex III of the ETC. COX1 is the main subunit of cytochrome c oxidase or otherwise known as complex IV of the ETC and COX3 is a transmembrane subunit of this same complex. Previous transcriptomic experiments have reported a relationship between DEG related to the complexes of ETC and RFI^12,36. For example, increased levels of COXII (complex IV) and NADH dehydrogenase subunits have been shown to be associated with efficient animals at both the protein and transcriptome level^12,36-39. However, although genes of the ETC have previously been reported in relation to variation in RFI, this is the first report of an up-regulation of COX1, ND5, ND6, CYTB and COX3 in skeletal muscle tissue of cattle of low vs. high RFI. The lack of commonality of specific key genes underlying processes such as ETC to RFI phenotype may be due to the differences in breed types employed across studies, with Simmental cattle utilised in the current study. For example, in a study examining the effect of various breeds on the molecular control of RFI in liver tissue, Mukiibi et al.²² recorded only 5 genes as commonly differentially expressed across three differing breed types, however an evaluation of the biological processes showed a clear commonality for the underlying biological control of RFI irrespective of breed, but with different key genes dependent on the breed²².

Using gene ontology (GO) enrichment analysis, we identified important processes underlying feed efficiency variation in skeletal muscle of heifers. These included functions related to mitochondrial metabolism, in particular oxidative phosphorylation. Similarly, the top canonical pathways identified by IPA were mitochondrial function and oxidative phosphorylation. Taken together, these analyses highlight a greater capacity for mitochondrial function in the low RFI animals. The relationship between mitochondria and RFI has been addressed previously and mitochondrial dysfunction and oxidative stress have been implicated as contributing to variation in feed efficiency across varying species^12,36,40-43. Mitochondria are highly dynamic organelles that are responsible for 90% of the energy production in the body and are major reactive oxygen species (ROS) regulators⁴⁴. It seems likely that variation in mitochondrial function could contribute to variation in energy utilisation. Furthermore, it has been demonstrated that feed efficient animals exhibit greater capacity to modulate conditions of oxidative stress⁴⁵.

Feed efficient animals have been shown to have a higher activity of all enzymes of the ETC across multiple species including broilers and lambs^40,41,43 and the present study focused on cattle is in agreement with this. The DEG in the current study encode proteins involved in Complexes I, III and IV of the ETC indicating an impaired oxidative phosphorylation system in the skeletal muscle of the less efficient heifers. These results complement previous research suggesting an association between decreased respiration capacity and increased ROS production in less efficient animals^41,46. Moreover, the results of the present study are reinforced by the observation of Kong et al.³⁹ in which the mitoproteome was skewed towards high feed efficiency birds despite no difference in mitochondrial DNA between phenotypes, suggesting an increase in mitochondrial activity in the high feed efficiency phenotype³⁹.

Within the hepatic tissue of bulls we identified 11 DEGs, 8 of which were up-regulated in the low RFI animals. These included; HSPA5 (78 kDa glucose –regulated protein precursor), CRELD2 (cysteine rich with EGF like domains 2), HYOU1 (hypoxia up-regulated protein 1 precursor), MANF (mesencephalic astrocyte derived neurotrophic factor), ACTA2 (actin, alpha2, smooth skeletal muscle, aorta), ENSBTAG00000003151 (DNAJ heatshock protein family (Hsp40) member B11, GMPPB (GDP-mannose pyrophosphorylase B), GNPNAT1 (glucosamine-phosphate N-acetyltransferase 1). Three genes were down-regulated in the low RFI bulls. These included; LIMS2 (LIM zinc finger domain containing 2), GSTT1 (glutathione S-transferase theta 1), DBP (D-box binding PAR bzip transcription factor). HSPA5, CRELD2, HYOU1, ACTA2, DBP, MANF, and GSTT1 have previously been implicated in variation in feed efficiency^{22-24,34,36,45,47-50} and are potential candidate biomarkers for this complex trait. Furthermore, HSPA5, CRELD2, HYOU1 and GSTT1 were all significantly associated and MANF tended towards a significant association with RFI phenotype further implicating the importance of these key genes to RFI, not only to the bulls used in the current study but to other cohorts of cattle divergent for RFI status. Gene ontology analysis of hepatic tissue in bulls identified terms including those related to amino sugar and nucleotide sugar metabolism as well as other related to endoplasmic reticulum protein processing as significantly enriched. Additionally Ingenuity pathway analysis identified Aldosterone Signalling in Epithelial cells, GDP-mannose Biosynthesis, NRF2 mediated Oxidative Stress Response and EIF2 Signalling as the top canonical pathways related to RFI in bull hepatic tissue. The results of the GO enrichment and IPA analysis indicate that oxidative response, protein processing and cell signalling in the liver are likely to be processes that are influencing variation in feed efficiency. 

The aldosterone signalling pathway was identified by IPA as the top canonical pathway due to the up-regulation of the genes HSPA5 and DNAJB11 in the hepatic tissue of low RFI bulls. Aldosterone is secreted by the adrenal glands and has a major role in electrolyte and fluid homeostasis. The aldosterone signalling pathway has been implicated in feed efficiency previously in the spleen of inefficient animals³¹ and, interestingly, a GWAS analysis with cattle identified this pathway to be associated with variation in feed conversion ratio⁵¹. The protein encoded by HSPA5 is a member of the HSP70 family and as this protein interacts with many endoplasmic reticulum (ER) proteins it is likely to be important in monitoring protein transport through the cell⁵². DNAJB11 as a member of the DNA-J family of proteins is involved in the correct folding of proteins⁵³. Specifically DNAJB11 is involved in protein processing and metabolism of proteins and serves as a co-chaperone for HSPA5 in the ER⁵⁰. One of the many functions of this family of proteins is to stabilize new proteins by ensuring correct folding or by helping refold proteins that have already been damaged by cell stress. Both DNAJB11 and HSPA5 have been implicated with feed efficiency previously^24,31,45,54. In agreement with the current study, other work from our group⁵⁵ observed that animals undergoing compensatory growth with concomitant improvements in feed efficiency have a higher hepatic transcript abundance of DNAJB11 and HSPA5. Similarly, an up-regulation of HSPA5 was observed in the hepatic tissue of low RFI animals by Paradis et al.³⁴. Taken together these results are indicative of a greater capacity in controlling cellular function and organisation as well as protein metabolism in more feed efficient animals. 

NRF2 mediated oxidative stress response was also observed to be an overrepresented pathway in the hepatic tissue of high RFI bulls in the current study due to the up-regulation of the genes GSTT1 and DBP. NRF2 is a member of the cap ‘n’ collar basic region leucine zipper (cnc bZip) group of transcription factors⁵⁶. This transcription factor is ubiquitously expressed in tissues but is only activated in response to a range of oxidative and electrophilic stimuli including ROS, antioxidants, glucose induced oxidative damage, heavy metals, and certain disease processes^56-58. This canonical pathway has been associated with feed efficiency in a number of studies^{29,49,50,54,59}. GSTT1 is a member of the glutathione S-transferase family and is involved in metabolism of xenobiotics and in catalysing reactions between the antioxidant glutathione and a host of potentially toxic compounds, highlighting it as an important homeostatic molecule⁶⁰. The glutathione S-transferase family has previously been implicated in feed efficiency in various species^{29,49,50,54,59}. In agreement with the current study Chen et al.⁵⁰ and Lindholm-Perry et al.²⁹ observed an up-regulation of these genes in feed inefficient cattle. DBP is a Protein Coding gene and amongst its cited functions is activation of circadian gene expression. Gene ontology annotations related to this gene include transcription factor activity, sequence-specific DNA binding and transcriptional activator activity and RNA polymerase II core promoter proximal region sequence-specific binding⁶¹. Additionally DBP was also implicated as contributing to variation in RFI phenotype through the network analysis reported by Weber et al.²³. Similarly it has been shown that mice with increased FE have a lower expression of DBP⁶² which is in agreement with the present study. Taken together, our work and that of the aforementioned authors, suggest that less efficient animals are exhibiting an increased oxidative stress, reflected in their increased anti-oxidation activities. 

The differential expression and significant association of HSPA5, CRELD2, HYOU1, GSTT1 and MANF with RFI in the present study are also noteworthy as all six genes have been previously observed as differentially expressed in relation to feed efficiency in cattle^{22,23,45,50,63}. While the biological significance of these genes in relation to feed efficiency remains unclear, due to their consistent presence in the literature, they should not be ruled out as potential biomarkers for this trait. 

Conclusion

The present study contributes to the published knowledge base regarding the transcriptomic regulation of variation in feed efficiency. Our work, in combination with that of others as previously mentioned, highlights common genes underpinning feed efficiency in cattle, as measured by RFI, regardless of breed or genetic background. RNA-seq analysis is an exploratory approach that provides new hypotheses to be further investigated by other complementary approaches including global proteomics and ultimately, potential variation in the genes identified in this study may provide a basis for the selection of candidate biomarkers for the RFI trait and, following appropriate validation, contribute to genomic selection breeding programmes to improve feed efficiency in beef cattle. However, the key message from this work highlights the inconsistency in gene expression profiles across genders of the same genetic merit. While this inconsistency may be explained by the differences in physiological maturity of the two genders, it indicates that extensive further investigation is required before biomarker selection for RFI can be adopted.

Methods

All procedures involving animals in this study were conducted under an experimental licence from the Health Products Regulatory Authority in accordance with the cruelty to Animals Act 1876 and the European Communities (Amendment of Cruelty to Animals Act 1876) Regulation 2002 and 2005.

Animal Model

The animals used in this study were derived from a purebred herd of Simmental cattle originally established to examine various aspects of the biological control of the RFI trait and which has been well characterised to date in the published literature^13-19. In order to generate animals for the current study, the 20 highest (inefficient) and 20 lowest (efficient) ranking cows on phenotypic RFI were subjected to a multiple ovulation and embryo recovery programme and breeding using artificial insemination to bulls with high and low estimated breeding values for RFI (see Crowley et al.²⁰), respectively. There was no crossover of sires used across both cohorts of donor females. Resultant embryos were transferred to crossbred beef heifer recipients. Pregnant heifers were managed under standard protocols and following calving were allowed to suckle their calves for a period of up to 7 days. In order to standardise rearing, calves were then abruptly weaned and were subsequently reared on an electronic calf feeder. Briefly calves were offered milk replacer (MR) (Blossom Easymix; Volac, Co. Cavan, Ireland) and concentrate in pelleted form using an electronic feeding system (Vario; Foster-Tecknik, Engen, Germany), which recorded all feed-related events including intake of both MR and concentrate, drinking speed, as well as number of rewarded (when calves receive milk) and unrewarded (no milk dispensed) visits to the machine. Calves were subsequently weaned at 10 weeks of age and were offered concentrate and hay on a 50:50 dry matter basis until turnout to pasture at approximately six months of age. At approximately 15 months of age all cattle were housed within pens of between 5-7 animals/pen in a slatted floor shed. Cattle were fed once daily (0800 h) and were offered ad libitum concentrate (860 g/kg rolled barley, 60 g/kg soya bean meal, 60 g/kg molasses and 20 g/kg minerals/vitamins) and 3kg grass silage to retain ruminal function. The animals had an acclimatisation period of 14 days to the ad libitum regime and test facilities before the experimental recoding period commenced. Feed intake was recorded daily and BW was recorded twice weekly. The recording period lasted 70 days. Concentrates and silage offered were sampled three times weekly and samples were stored at -20°C pending laboratory analysis. Samples of concentrates and silage were subsequently pooled on a weekly basis for DM determination. Concentrate samples were dried in an oven with forced-air circulation at 98°C for 16 h for DM determination and forage samples dried at 40° for 48 h. 

Statistical Analysis

Average daily live weight gain during the RFI measurement period for each animal was computed as the coefficient of the linear regression of BW (kg) on time (d) by using the GLM procedure of SAS 9.1 (SAS Inst. INC., Cary, NC). Mid-test metabolic BW (MBW) was represented as BW^0.75 35 d before the end of the test which was estimated from the intercept and slope of the regression line. Residual feed intake was calculated for each animal as the difference between actual DMI and expected DMI. Expected DMI was computed for each animal using a multiple regression model, regressing DMI on MBW, ADG and mean lumbar BF change. Animals were ranked within gender and following a power analysis utilising RFI phenotype data previously generated by our group^13-19, a total of 5 animals per group was required for statistical analysis, thus the most extreme animals for high RFI (n=5) and low RFI (n=5) were selected for further analysis. 

Biopsy sample collection

M.longissimus thoracis et lumborum (skeletal muscle) biopsies were harvested as described by Kelly et al.¹² and hepatic tissue was collected by percutaneous punch biopsy as described by McCarthy et al.⁶⁴ from animals deemed high and low RFI under local anaesthetic (5 mL Adrenacaine, Norbrook Laboratories (Ireland) Ltd.) at the end of the RFI measurement period. All surgical instruments used for tissue collection were sterilised and treated with 70% Ethanol and RNaseZap (Ambion, Applera Ireland, Dublin, Ireland). M.longissimus thoracis et lumborumbiopsies were snap frozen in liquid nitrogen directly after collection and hepatic tissue biopsies were washed in sterile phosphate buffered saline and snap frozen in liquid nitrogen. All samples were subsequently stored at -80°C for long-term storage pending further processing. 

RNA isolation and purification  

Total RNA was isolated from 50mg of biopsy samples using QIAzol (Qiagen, UK). Tissue samples were homogenised in 1mL of QIAzol reagent using a rotor-strator tissue lyser (Qiagen, UK) and chloroform (Sigma-Aldrich Ireland, Dublin, Ireland). RNA was subsequently precipitated and purified using the RNeasy plus Universal kit (Qiagen, UK) according to the manufacturer’s guidelines, which included a step to remove any contaminating genomic DNA. The quantity of the RNA isolated was determined by measuring the absorbance at 260nm using a Nanordrop spectrophotometer ND-1000 (Nanodrop Technologies, Wilmington, DE, USA). RNA quality was assessed on the Agilent Bioanalyser 2100 using the RNA 6000 Nano Lab Chip kit (Agilent Technologies Ireland Ltd., Dublin, Ireland). RNA quality was also verified by ensuring all RNA samples had an absorbance (A260/280) of between 1.8 and 2.0 and RIN (RNA integrity number) values of between 8 and 10 were deemed high quality. Any samples that had an (A260/280) absorbance of less than 1.8 were cleaned using Zymo Research RNA clean & concentrator kit (Cambridge Biosciences, UK). High quality RNA samples were selected for cDNA synthesis. 

cDNA library preparation and sequencing 

cDNA libraries were prepared from high quality RNA following the manufacturer’s instructions using the Illumina TruSeq RNA sample prep kit (Illumina, San Diego, CA, USA). For each sample, 1 μg of RNA was used for cDNA library preparation. Resultant cDNA libraries were validated on the Agilent Bioanalyser 2100 using the DNA 1000 Nano Lab Chip kit. cDNA concentration was assessed using Nanordrop spectrophotometer ND-1000 (Nanodrop Technologies, Wilmington, DE, USA) and samples with >25ng/μl were deemed suitable for further analysis. Following quality control procedures, individual RNAseq libraries were pooled based on their respective sample-specific 6bp adaptors and sequenced at 100bp/sequence single-end read using and Illumina HiSeq 2500 sequencer. Approximately 22.4 million sequences per sample were generated.

RNAseq data analysis

FASTQC software (v0.11.5) was used to check the quality of the raw sequencing reads. Input reads were then aligned to the bovine reference genome (UMD3.1) using STAR (v2.5.1). HTSeq (v0.6.1p2) software was used to calculate the number of sequenced fragments overlapping all protein-coding genes from the ENSEMBLv88 annotation of the bovine genome. The number of counts of reads mapping to each annotated gene from HTSeq was than collated into a single file and used for subsequent differential gene expression. The R (v3.4) Bioconductor package, EdgeR (v3.5), which uses a negative binomial distribution model to account for both biological and technical variability, was applied to identify statistically significant DEGs. The analysis was undertaken using moderated tagwise dispersion. An adjusted p-value<0.1 was applied as a threshold to call genes with differential expression levels. The R (v3.4) Corrplot package was used to visualise the correlation matrices pertaining to the DEGs passing multiple correction for each comparison. Gene expression results for DEGs identified through transcriptional profiling were also correlated with RFI value, to determine potential associations between the DEGs and RFI. Correlations were undertaken using the CORR procedure in SAS 9.1 (SAS Inst. INC., Cary, NC).

Gene ontology and pathway analyses

Biological processes, cellular components and molecular functions that were associated with the DEGs were identified using the g:Profiler GO enrichment tool⁶⁵. Fisher’s exact test was used to determine the enrichment of the GO terms. Annotation clusters and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analyses were assessed using the official gene symbols of the DEGs and Bos Taurus as the species where possible. P-values associated with each annotation term inside each cluster are Fisher Exact/EASE Scores. To examine the molecular functions and genetic networks; the RNAseq data was further interrogated using IPA (Ingenuity Systems, Redwood City, CA; http://www.ingenuity.com), a web based software application that enables identification of over-represented biological mechanisms, pathways and functions most relevant to experimental datasets or genes of interest⁶⁶. Data were imported in a flexible format using the gene symbol as the identifier. 

Abbreviations

Declarations

Acknowledgements

The authors wish to acknowledge funding from the Teagasc Walsh Fellowship programme (Project Number: RMIS 6092). Authors also acknowledge assistance from farm and technical staff at Teagasc.

Authors' contributions

DK conceived the animal model and contributed to the statistical analysis and manuscript preparation. CM managed the animal model, conducted laboratory and statistical analysis and prepared the manuscript. KK assisted with molecular and bioinformatics analyses and contributed to manuscript preparation. SMW oversaw the molecular analyses and assisted with manuscript preparation. PC oversaw the bioinformatics analysis. RKP contributed to conception of the study and assisted with manuscript preparation. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests

Data availability

The datasets generated and analysed during the current study are available in the National Centre for Biotechnology Information (NCBI), Gene Expression Omnibus repository, and are accessible through GEO accession number GSE112793 [https://www.ncbi.nlm.nih.gov/geo/] under accession number.

Ethics declarations

All procedures involving animals in this study were conducted under an experimental licence from the Health Products Regulatory Authority in accordance with the cruelty to Animals Act 1876 and the European Communities (Amendment of Cruelty to Animals Act 1876) Regulation 2002 and 2005. All procedures described in this study were carried out in compliance with the ARRIVE guidelines.

Consent for publication

Not applicable

References

1. The state of Food and agriculture (SOFA) - Livestock in the balance. Food and Agriculture Organisation. 2009. http://wwwofao.org/012/i0680e00.htm (2009).
2. Arthur, P F , Archer, J. A. & Herd, R. M. Feed intake and efficiency in beef cattle: overview of recent Australian research and challenges for the future. J. Exp. Agric. 44, 361-369 (2004).
3. Kelly, A. K. et al. Effect of divergence in residual feed intake on feeding behavior, blood metabolic variables, and body composition traits in growing beef heifers. Anim. Sci. 88, 109-123 (2010).
4. Crowley, J. J. et al. Phenotypic and genetic parameters for different measures of feed efficiency in different breeds of Irish performance-tested beef bulls. Anim. Sci. 88, 885-94 (2010).
5. Herd, R. M. & Arthur, P. F. Physiological basis for residual feed intake. Anim. Sci. 87 Suppl 14, E64-71 (2009).
6. Moore, S. S., Mujibi, F. D. & Sherman, E. L. Molecular basis for residual feed intake in beef cattle. J. Anim. Sci. 87 Suppl 14, E41-47 (2009).
7. Koch, R. M., Swiger, L. A., Chambers, D. & Gregory, K. E. Efficiency of feed use in beef cattle. Anim. Sci. 22, 486-494 (1963).
8. Crews, D. H. Jr. Genetics of efficient feed utilization and national cattle evaluation: a review. Mol. Res. 4, 152-165 (2004).
9. Crews, D. H. Jr., Lowerison, M., Caron, N. & Kemp, R. A. Genetic parameters among growth and carcass traits of Canadian Charolais cattle. J. Anim. Sci. 84, 589-597 (2004).
10. Kelly, A. K. et al. Relationship between body measurements, metabolic hormones, metabolites and residual feed intake in performance tested pedigree beef bulls. Sci. 135, 8-16 (2011).
11. Kelly, A. K. et al. Repeatability of feed efficiency, carcass ultrasound, feeding behavior, and blood metabolic variables in finishing heifers divergently selected for residual feed intake. Anim. Sci. 88, 3214-3225 (2010).
12. Kelly, A. K. et al. mRNA expression of genes regulating oxidative phosphorylation in the muscle of beef cattle divergently ranked on residual feed intake. Genomics. 43, 12-23 (2011).
13. Fitzsimons, C., Kenny, D. A., Waters, S. M., Earley, B. & McGee, M. Effects of phenotypic residual feed intake on response to a glucose tolerance test and gene expression in the insulin signaling pathway in longissimus dorsi in beef cattle. Anim. Sci. 92, 4616-4631 (2014).
14. Fitzsimons, C., Kenny, D. A., Fahey, A. G. & McGee, M. Feeding behavior, ruminal fermentation, and performance of pregnant beef cows differing in phenotypic residual feed intake offered grass silage. Anim. Sci. 92, 2170-81 (2014).
15. Fitzsimons, C., Kenny, D. A. & McGee, M. Visceral organ weights, digestion and carcass characteristics of beef bulls differing in residual feed intake offered a high concentrate diet. Animal 8, 949-59 (2014).
16. Lawrence, P., Kenny, D. A., Earley, B., Crews, D. H. Jr. & McGee, M. Grass silage intake, rumen and blood variables, ultrasonic and body measurements, feeding behavior, and activity in pregnant beef heifers differing in phenotypic residual feed intake. Anim. Sci. 89, 3248-3261 (2011).
17. Lawrence, P., Kenny, D. A., Earley, B. & McGee, M. Intake of conserved and grazed grass and performance traits in beef suckler cows differing in phenotypic residual feed intake. Sci. 152, 154-166 (2013).
18. Lawrence, P., Kenny, D. A., Earley, B. & McGee, M. Grazed grass herbage intake and performance of beef heifers with predetermined phenotypic residual feed intake classification. Animal 6, 1648-1661 (2012).
19. Fitzsimons, C., Kenny, D. A., Deighton, M. H., Fahey, A. G. & McGee, M. Methane emissions, body composition, and rumen fermentation traits of beef heifers differing in residual feed intake. Anim .Sci. 91, 5789-5800 (2013).
20. Crowley, J. J. et al. Phenotypic and genetic parameters for different measures of feed efficiency in different breeds of Irish performance-tested beef bulls. Anim. Sci. 88, 885-894 (2010).
21. Fitzsimons, C., McGee, M., Keogh, K., Waters, S. M. & Kenny, D. A. Molecular physiology of feed efficiency in beef cattle in The biology of domestic animals (ed. Scanes, C. G. & Hill, R. A.) 120-163 (CRC Press, 2017).
22. Mukiibi, R. et al. Transcriptome analyses reveal reduced hepatic lipid synthesis and fat accumulation in more efficient beef cattle. Rep. 8, 7303 (2018).
23. Weber, K. L. et al. Identification of Gene Networks for Residual Feed Intake in Angus Cattle Using Genomic Prediction and RNA-seq. PLoS One 11, e0152274 (2016).
24. Connor, E. E. et al. Enhanced mitochondrial complex gene function and reduced liver size may mediate improved feed efficiency of beef cattle during compensatory growth. Integr. Genomics 10, 39-51 (2010).
25. Keogh, K., Waters, S. M. & Kenny, D. A. Invited review: Genomics of feed efficiency and related traits in cattle. British Society of Animal Science. Dublin, Ireland 9-11^th April p. 247 (2018).
26. Culbertson, M. M. et al. Optimum measurement period for evaluating feed intake traits in beef cattle. Anim. Sci. 93, 2482-2487 (2015).
27. Johnson, D. E., Ferrell, C. L. & Jenkins, T. G. The history of energetic efficiency research: where have we been and where are we going? Anim. Sci. 81 Suppl 13, E27-E38 (2003).
28. Henriksson-Larsen, K., Wretling, M. L., Lorentzon, R. & Oberg, L. Do muscle fibre size and fibre angulation correlate in pennated human muscles? J. Appl. Physiol. Occup. Physiol. 64, 68-72 (1992).
29. Lindholm‐Perry, A. K. et al. Differential gene expression in the duodenum, jejunum and ileum among crossbred beef steers with divergent gain and feed intake phenotypes. Genet. 47, 408-427 (2016).
30. Zarek, C. M., Lindholm-Perry, A. K., Kuehn, L. A. & Freetly, H. C. Differential expression of genes related to gain and intake in the liver of beef cattle. BMC Res. Notes 10, 1 (2017).
31. Lindholm-Perry, A. K. et al. Profile of the spleen transcriptome in beef steers with variation in gain and feed intake. Genet. 7,127 (2016).
32. McKenna, C. et al. Residual feed intake phenotype and gender affect the expression of key genes of the lipogenesis pathway in subcutaneous adipose tissue of beef cattle. Anim. Sci. Biotechnol. 9, 84 (2018).
33. Alexandre, P. A. et al. Liver transcriptomic networks reveal main biological processes associated with feed efficiency in beef cattle. BMC Genomics. 16, 1073 (2015).
34. Paradis, F. et al.. Transcriptomic analysis by RNA sequencing reveals that hepatic interferon-induced genes may be associated with feed efficiency in beef heifers. Anim. Sci. 93, 3331-3341 (2015).
35. Salleh, M. S. et al. RNA-Seq transcriptomics and pathway analyses reveal potential regulatory genes and molecular mechanisms in high- and low-residual feed intake in Nordic dairy cattle. BMC Genomics 18, 258 (2017).
36. Tizioto, P. C. et al. Gene expression differences in Longissimus muscle of Nelore steers genetically divergent for residual feed intake. Rep. 6, 39493 (2016).
37. Iqbal, M. et al. Low feed efficient broilers within a single genetic line exhibit higher oxidative stress and protein expression in breast muscle with lower mitochondrial complex activity. Sci. 83, 474-484 (2004).
38. Ojano-Dirain, C., Toyomizu, M., Wing, T., Cooper, M. & Bottje, W.G. Gene expression in breast muscle and duodenum from low and high feed efficient broilers. Sci. 86, 372-381 (2007).
39. Kong, B. W. et al. Proteomics of Breast Muscle Tissue Associated with the Phenotypic Expression of Feed Efficiency within a Pedigree Male Broiler Line: I. Highlight on Mitochondria. PloS One. 11, e0155679 (2016).
40. Iqbal, M. et al. Compromised liver mitochondrial function and complex activity in low feed efficient broilers are associated with higher oxidative stress and differential protein expression. Sci. 84, 933-941 (2005).
41. Bottje, W.G. & Carstens, G. E. Association of mitochondrial function and feed efficiency in poultry and livestock species. Anim. Sci. 87 Suppl 14, E48-E63 (2009).
42. Bottje, W. & Kong, B. W. Cell Biology Symposium: feed efficiency: mitochondrial function to global gene expression. Anim. Sci. 91, 1582-1593 (2013).
43. Sharifabadi, H. R., Zamiri, M. J., Rowghani, E. & Bottje, W. G. Relationship between the activity of mitochondrial respiratory chain complexes and feed efficiency in fat-tailed Ghezel lambs. Anim. Sci. 90, 1807-1815 (2012).
44. Kang, J. & Pervaiz, S. Mitochondria: redox metabolism and dysfunction. Res. Int. 2012, 896751 (2012).
45. Tizioto, P. C. et al. Global liver gene expression differences in Nelore steers with divergent residual feed intake phenotypes. BMC Genomics 16, 242 (2015).
46. Ferguson, M., Mockett, R. J., Shen, Y., Orr, W. C. & Sohal, R. S. Age-associated decline in mitochondrial respiration and electron transport in Drosophila melanogaster. Biochem J. 390, 501-511 (2005).
47. Zhang, R., Zhu, W. & Mao, S. High-concentrate feeding upregulates the expression of inflammation-related genes in the ruminal epithelium of dairy cattle. Anim. Sci. Biotechnol. 7, 42 (2016).
48. Al‐Husseini, W. et al. Expression of candidate genes for residual feed intake in Angus cattle. Anim Genet. 45, 12-19 (2014).
49. Barendse, W. et al. A validated whole-genome association study of efficient food conversion in cattle. Genetics. 176, 1893-1905 (2007).
50. Chen, Y. et al. Global gene expression profiling reveals genes expressed differentially in cattle with high and low residual feed intake. Genet. 42, 475-490 (2011).
51. Santana, M. H. A. et al. Systems genetics investigations for feed intake, feed efficiency and performance in Nellore (Bos indicus) Cattle. Proceedings of the 10^th World Congress of Genetics Applied to Livestock Production. Auckland, 11-18^th February (2018).
52. Calderwood, S. K., Khaleque, M. A., Sawyer, D. B. & Ciocca, D. R. Heat shock proteins in cancer: chaperones of tumorigenesis. Trends Biochem Sci. 31, 164-172 (2006).
53. Dekker, S. L., Kampinga, H. H. & Bergink, S. DNAJs: more than substrate delivery to HSPA. Mol. Biosci. 2, 35 (2015).
54. Kong, B. W. et al. RNA sequencing for global gene expression associated with muscle growth in a single male modern broiler line compared to a foundational Barred Plymouth Rock chicken line. BMC Genomics 18, 82 (2017).
55. Keogh, K., Kenny, D. A., Cormican, P., Kelly, A. K. & Waters, S. M. Effect of dietary restriction and subsequent re-alimentation on the transcriptional profile of hepatic tissue in cattle. BMC Genomics 17, 244 (2016).
56. He, X., Kan, H., Cai, L. & Ma, Q. Nrf2 is critical in defense against high glucose-induced oxidative damage in cardiomyocytes. Mol. Cell Cardiol. 46, 47-58 (2009).
57. Ma, Q. Role of nrf2 in oxidative stress and toxicity. Rev. Pharmacol. Toxicol. 53, 401-426 (2013).
58. Ma, Q. Xenobiotic-activated receptors: from transcription to drug metabolism to disease. Res. Toxicol. 21, 1651-1671 (2008).
59. Zhou, N., Lee, W. R. & Abasht, B. Messenger RNA sequencing and pathway analysis provide novel insights into the biological basis of chickens' feed efficiency. BMC Genomics 16, 195 (2015).
60. Pickett, C. B. & Lu, A. Y. Glutathione S-transferases: gene structure, regulation, and biological function. Rev. Biochem. 58, 743-764 (1989).
61. Khatib, Z. A., Inaba, T., Valentine, M. & Look, A. T. Chromosomal localization and cDNA cloning of the human DBP and TEF genes. Genomics 23, 344-51 (1994).
62. Dankel, S. N. et al. Weight cycling promotes fat gain and altered clock gene expression in adipose tissue in C57BL/6J mice. J. Physiol. Endocrinol. Metab. 306, E210-24 (2014).
63. Kern, R. J. et al. Transcriptome differences in the rumen of beef steers with variation in feed intake and gain. Gene 586, 12-26 (2016).
64. McCarthy, S. D. et al. Negative energy balance and hepatic gene expression patterns in high-yielding dairy cows during the early postpartum period: a global approach. Genomics 42, 188-199 (2010).
65. Reimand, J. et al. g:Profiler-a web server for functional interpretation of gene lists. Nucleic Acids Res. 44(W1), W83-89 (2016).
66. Kramer, A., Green, J., Pollard, J. & Tugendreich, S. Causal analysis approaches in Ingenuity Pathway Analysis. Bioinformatics 30, 523-530 (2014).

Tables

Table 1. Summary of phenotypic data of bulls at end of feed intake trial

¹High RFI is inefficient and low RFI is efficient

²Metabolic BW^0.75, kg is determined as BW^0.75 in the middle of the RFI measurement period which was estimated from the intercept and slope of the regression line after fitting a linear regression line through all metabolic BW (BW^0.75) observations.

³Back fat change is mean of the fat depth at end of intake trial – mean of fat depth at start of intake trial. 

Table 2. Differentially expressed genes in the skeletal muscle of heifers divergent for RFI

Table 3. Differentially expressed genes in the hepatic tissue of bulls divergent for RFI

Table 4. Correlation between DEGs and RFI value in both heifers and bulls

Table 5. Enriched GO functions in hepatic tissue from bulls divergent for RFI. 

Table 6. Enriched GO functions in skeletal muscle from heifers divergent for 

Table 7. Canonical pathways derived from IPA analysis associated with divergence for RFI in bull hepatic tissue and heifer skeletal muscle tissue.