Effect of Gender and Residual Feed Intake Phenotype on the Transcriptome of Liver and Skeletal Muscle Tissue in Purebred Simmental Beef Cattle

Abstract

Background: The selection of cattle with enhanced feed efficiency is of paramount importance with regard to reducing feed costs in the beef industry. However the role, if any, of gender to the underlying molecular control of feed efficiency in cattle is not currently known. 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 differentially expressed genes that may be associated with this trait.

Results: We identified 5 genes (p<0.001; false discovery rate (FDR) <0.1) to be differentially expressed in skeletal muscle between high and low RFI heifers with all 5 transcripts being up-regulated in the low RFI phenotype. Among these differentially expressed genes, all transcripts were involved in oxidative phosphorylation and mitochondrial homeostasis. A total of 11 genes (p<0.001;FDR <0.1) were differentially expressed in hepatic tissue between high and low RFI bulls with 8 transcripts being up-regulated and 3 being down-regulated in the low RFI phenotype. These differentially expressed genes were related to oxidative response, protein mediation and cell signalling. No genes were identified as differentially expressed in both heifer liver and bull muscle analyses.

Conclusions: 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. 

Background

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 [1]. 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 [2]. Feed efficiency (dietary nutrient utilisation) in beef cattle is a trait of major economic importance [3]. 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) [7] 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 [8]. 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 [20] 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 animals 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 [21]. This not only provides an attractive alternative to direct measurement of dietary intake on large numbers of animals [6], 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 [21]. 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 [22]. Furthermore, of these 5 genes only one, SCD, was differentially expressed in only one other molecular based RFI evaluation [23]. 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 [5], with both muscle and liver representing important metabolic tissues with skeletal muscle accounting for over 50% body weight [12] and liver accounting for 18–25% of total oxygen body consumption [24]. 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

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 (p < 0.001; false discovery rate (FDR) < 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. 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 (p < 0.001; FDR < 0.1) were identified to be differentially expressed in hepatic tissue between high and low RFI bulls 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. Further details of these DEGs are provided in Tables 3 and 4. 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 5. 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 using both the up- and down-regulated genes was performed using DAVID software and revealed that in hepatic tissue from bulls the DEG were related to metabolism of amino sugars and nucleotides, posttranslational modification, protein turnover and chaperones and protein processing in the endoplasmic reticulum (Table 6). Functional analysis using DAVID in heifer skeletal muscle indicated the DE genes to be involved in oxidative phosphorylation and mitochondrial homeostasis (Table 7). 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 8. 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 [21]. Furthermore, data from our own group has also shown that RFI is a repeatable trait in beef cattle [11]. 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 [26]. 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 [27]. 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 [28]. 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 [24]. 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 [32]. 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, 34, 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 indicated 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 the protein and transcriptome level [12, 36–39] in more efficient animals. 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. [22] 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 [22].

Using gene ontology (GO) enrichment analysis through the DAVID platform, we identified important processes underlying feed efficiency variation in skeletal muscle of heifers. The most enriched cluster contained a plethora of 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 [44]. 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 [45].

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. [39] 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 [39].

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 enrichment using DAVID revealed one enriched cluster in the up-regulated DEG from the low RFI bulls, which included signalling functions, disulphide bonding and peptide signalling. 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 mediation 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 [31] and, interestingly, a GWAS study with cattle identified this pathway to be associated with variation in feed conversion ratio [51]. 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 [52]. DNAJB11 as a member of the DNA-J family of proteins is involved in the correct folding of proteins [53]. Specifically DNAJB11 is involved in protein processing and metabolism of proteins and serves as a co-chaperone for HSPA5 in the ER [50]. 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 [55] 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. [34]. 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 [56]. 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 [60]. 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. [50] and Lindholm-Perry et al. [29] 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 [61]. Additionally DBP was also implicated as contributing to variation in RFI phenotype through the network analysis reported by Weber et al. [23]. Similarly it has been shown that mice with increased FE have a lower expression of DBP [62] 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

Animal Model

The animals used in this study were derived from a purebred herd of Simmental cattle located at Teagasc Grange (Dunsany, Co. Meath, Irealnd) 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 (n=40 in total) 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. [20]), respectively. There was no crossover of sires used across both cohorts of donor females. Following embryo recovery, all cows used were returned to the commercial herd within Teagasc Grange (Dunsany, Co. Meath, Ireland). Resultant embryos were transferred to crossbred beef heifer recipients (n=40), which were purchased from commercial herds specifically for use in this study. In order to overcome potential confounding effects, animals were all purebred from the same genetic pool and were reared under the same conditions as a contemporary group from birth. All animals used were healthy, received regular veterinary interventions as required and were supervised on a daily basis by trained animal care staff. Pregnant heifers were managed under standard protocols and following calving were allowed to suckle their calves for a period of up to 7 days. Following this, all heifers used were returned to the commercial herd within Teagasc Grange (Dunsany, Co. Meath, Ireland). In order to standardise rearing, calves (n=40) 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 and the experimental unit was each individual animal within treatment (20 experimental units in each high and low RFI group). Following completing of the RFI measurement period all heifers and bulls used were slaughtered in an EU licenced abattoir (Eurofarm Foods, Duleek, Co. Meath, Ireland). 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

In order to overcome potential confounding effects, animals were all purebred from the same genetic pool and were reared under the same conditions as a contemporary group from birth. 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 (www.biomath.info/power/ttest.htm) utilising RFI phenotype data previously generated by our group [13-19], a total of 5 animals per group was required for statistical analysis between high RFI and low RFI groups, thus the most extreme animals for high RFI (n=5) and low RFI (n=5) were selected for further analysis, resulting in a total sample size of n=20.

Biopsy sample collection

M.longissimus thoracis et lumborum (skeletal muscle) biopsies were harvested as described by Kelly et al. [12] and hepatic tissue tissue was collected by percutaneous punch biopsy as described by McCarthy et al. [64] 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 to 2.0 and RIN (RNA integrity number) values of between 8 and 10 were deemed high quality. Any samples that had a (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. A FDR of <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).

Pathway analysis

Biological processes, cellular components and molecular functions that were associated with the DEGs were identified using the GO enrichment tool of DAVID (Database for Annotation, Visualization and Integrated Discovery) Bioinformatics Resources v6.8 [65]. 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 [55, 66]. Data were imported in a flexible format using the gene symbol as the identifier.

Abbreviations

Declarations

Ethics approval and consent to participate

Detailed protocols for the study and any licensable procedures involving animals in this study were submitted, reviewed and approved by the Teagasc Animal Ethics Committee and were conducted under an experimental licence issued by the Health Products Regulatory Authority (AE19132/PO11) 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 protocols and procedures employed were sanctioned as adhering to Irish and EU animal welfare and research ethics law. All animals used in this study were bred, reared and owned by Teagasc Grange (Dunany, Co. Meath, Ireland).

Consent for publication

Not applicable

Availability of data and materials

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.

Competing interests

The authors declare that they have no competing interests

Funding

The authors wish to acknowledge funding from the Teagasc Walsh Fellowship programme (Project Number: RMIS 6092).

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.

Acknowledgements

Authors acknowledge assistance from farm and technical staff at Teagasc.

References

1. FAO. The state of Food and agriculture (SOFA) - Livestock in the balance. Food and Agriculture Organisation. 2009. http://wwwofao.org/012/i0680e00.htm. Accessed 30 June 2019.
2. Arthur PF, Archer JA, Herd rm. Feed intake and efficiency in beef cattle: overview of recent Australian research and challenges for the future. Aust J Exp Agric. 2004;44:361–9.
3. Kelly AK, McGee M, Crews DH, Fahey AG, Wylie AR, Kenny DA. Effect of divergence in residual feed intake on feeding behavior, blood metabolic variables, and body composition traits in growing beef heifers. J Anim Sci. 2010;88:109–23.
4. Crowley JJ, McGee M, Kenny DA, Crews DH, Evans RD, Berry DP. Phenotypic and genetic parameters for different measures of feed efficiency in different breeds of Irish performance-tested beef bulls. J Anim Sci. 2010;88:885–94.
5. Herd RM, Arthur PF. Physiological basis for residual feed intake. J Anim Sci. 2009;87(Suppl 14):E64–71.
6. Moore SS, Mujibi FD, Sherman EL. Molecular basis for residual feed intake in beef cattle. J Anim Sci. 2009;87(Suppl 14):E41–7.
7. Koch RM, Swiger LA, Chambers D, Gregory KE. Efficiency of feed use in beef cattle. J Anim Sci. 1963;22:486–94.
8. Crews DH Jr. Genetics of efficient feed utilization and national cattle evaluation: a review. Genet Mol Res. 2004;4:152–65.
9. Crews DH Jr, Lowerison M, Caron N, Kemp RA. Genetic parameters among growth and carcass traits of Canadian Charolais cattle. Can J Anim Sci. 2004;84:589–97.
10. Kelly AK, McGee M, Crews DH Jr, Lynch CO, Wylie AR, Evans RD, et al. Relationship between body measurements, metabolic hormones, metabolites and residual feed intake in performance tested pedigree beef bulls. Livest Sci. 2011;135:8–16.
11. Kelly AK, McGee M, Crews DH Jr, Sweeney T, Boland TM, Kenny DA. Repeatability of feed efficiency, carcass ultrasound, feeding behavior, and blood metabolic variables in finishing heifers divergently selected for residual feed intake. J Anim Sci. 2010;88:3214–25.
12. Kelly AK, Waters SM, McGee M, Fonseca RG, Carberry C, Kenny DA. mRNA expression of genes regulating oxidative phosphorylation in the muscle of beef cattle divergently ranked on residual feed intake. Physio Genomics. 2011;43:12–23.
13. Fitzsimons C, Kenny DA, Waters SM, 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. J Anim Sci. 2014;92:4616–31.
14. Fitzsimons C, Kenny DA, Fahey AG, McGee M. Feeding behavior, ruminal fermentation, and performance of pregnant beef cows differing in phenotypic residual feed intake offered grass silage. J Anim Sci. 2014;92:2170–81.
15. Fitzsimons C, Kenny DA, McGee M. Visceral organ weights, digestion and carcass characteristics of beef bulls differing in residual feed intake offered a high concentrate diet. Animal. 2014;8:949–59.
16. Lawrence P, Kenny DA, Earley B, Crews DH 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. J Anim Sci. 2011;89:3248–61.
17. Lawrence P, Kenny DA, Earley B, McGee M. Intake of conserved and grazed grass and performance traits in beef suckler cows differing in phenotypic residual feed intake. Livest Sci. 2013;152:154–66.
18. Lawrence P, Kenny DA, Earley B. McGee M.Grazed grass herbage intake and performance of beef heifers with predetermined phenotypic residual feed intake classification. Animal. 2012;6:1648–61.
19. Fitzsimons C, Kenny DA, Deighton MH, Fahey AG, McGee M. Methane emissions, body composition, and rumen fermentation traits of beef heifers differing in residual feed intake. J Anim Sci. 2013;91:5789–800.
20. Crowley JJ, McGee M, Kenny DA, Crews DH Jr, Evans RD, Berry DP. Phenotypic and genetic parameters for different measures of feed efficiency in different breeds of Irish performance-tested beef bulls. J Anim Sci. 2010;88:885–94.
21. Fitzsimons C, McGee M, Keogh K, Waters SM, Kenny DA. Molecular physiology of feed efficiency in beef cattle. In: Scanes CG, Hill RA, editors. The Biology of Domestic Animals. Boca Raton: CRC Press; 2017.
22. Mukiibi R, Vinsky M, Keogh K, Fitzsimmons C, Stothard P, Waters S, et al. Transcriptome analyses reveal reduced hepatic lipid synthesis and fat accumulation in more efficient beef cattle. Sci Rep. 2018;8:7303.
23. Weber KL, Welly BT, Van Eenennaam AL, Young AE, Porto-Neto LR, Reverter A. Identification of Gene Networks for Residual Feed Intake in Angus Cattle Using Genomic Prediction and RNA-sEq. PLoS One. 2016;11:e0152274.
24. Connor EE, Kahl S, Elsasser TH, Parker JS, Li RW, Van Tassell CP, et al. Enhanced mitochondrial complex gene function and reduced liver size may mediate improved feed efficiency of beef cattle during compensatory growth. Funct Integr Genomics. 2010;10:39–51.
25. Keogh K, Waters SM, Kenny DA. Invited review: Genomics of feed efficiency and related traits in cattle. British Society of Animal Science annual meeting, Dublin, 9-11th April 2018. p. 247.
26. Culbertson MM, Speidel SE, Peel RK, Cockrum RR, Thomas MG, Enns RM. Optimum measurement period for evaluating feed intake traits in beef cattle. J Anim Sci. 2015;93:2482–7.
27. Johnson DE, Ferrell CL, Jenkins TG. The history of energetic efficiency research: where have we been and where are we going? J Anim Sci. 2003;81(Suppl 13):E27–38.
28. Henriksson-Larsen K, Wretling ML, Lorentzon R, Oberg L. Do muscle fibre size and fibre angulation correlate in pennated human muscles? Eur J Appl Physiol Occup Physiol. 1992;64:68–72.
29. Lindholm-Perry AK, Butler AR, Kern RJ, Hill R, Kuehn LA, Wells JE, et al. Differential gene expression in the duodenum, jejunum and ileum among crossbred beef steers with divergent gain and feed intake phenotypes. Anim Genet. 2016;47:408–27.
30. Zarek CM, Lindholm-Perry AK, Kuehn LA, Freetly HC. Differential expression of genes related to gain and intake in the liver of beef cattle. BMC Res Notes. 2017;10:1.
31. Lindholm-Perry AK, Kern RJ, Keel BN, Snelling WM, Kuehn LA, Freetly HC. Profile of the spleen transcriptome in beef steers with variation in gain and feed intake. Front Genet. 2016;7:127.
32. McKenna C, Porter RK, Keogh K, Waters SM, McGee M, Kenny DA. Residual feed intake phenotype and gender affect the expression of key genes of the lipogenesis pathway in subcutaneous adipose tissue of beef cattle. J Anim Sci Biotechnol. 2018;9:84.
33. Alexandre PA, Kogelman LJA, Santana MG, Passarelli D, Pulz LH, Fantinato-Neto P, et al. Liver transcriptomic networks reveal main biological processes associated with feed efficiency in beef cattle. BMC Genom. 2015;16:1073.
34. Paradis F, Yue A, Grant JR, Stothard P, Basarab JA, Fitzsimmons C. Transcriptomic analysis by RNA sequencing reveals that hepatic interferon-induced genes may be associated with feed efficiency in beef heifers. J Anim Sci. 2015;93:3331–41.
35. Salleh MS, Mazzoni G, Hoglund JK, Olijhoek DW, Lund P, Lovendahl P, 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 Genom. 2017;18:258.
36. Tizioto PC, Coutinho LL, Oliveira PSN, Cesar ASM, Diniz WJS, Lima AO, et al. Gene expression differences in Longissimus muscle of Nelore steers genetically divergent for residual feed intake. Sci Rep. 2016;6:39493.
37. Iqbal M, Pumford NR, Tang ZX, Lassiter K, Wing T, Cooper 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. Poult Sci. 2004;83:474–84.
38. Ojano-Dirain C, Toyomizu M, Wing T, Cooper M, Bottje WG. Gene expression in breast muscle and duodenum from low and high feed efficient broilers. Poult Sci. 2007;86:372–81.
39. Kong BW, Lassiter K, Piekarski-Welsher A, Dridi S, Reverter-Gomez A, Hudson NJ. 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. 2016;11:e0155679.
40. Iqbal M, Pumford NR, Tang ZX, Lassiter K, Ojano-Dirain C, Wing T, et al. Compromised liver mitochondrial function and complex activity in low feed efficient broilers are associated with higher oxidative stress and differential protein expression. Poult Sci. 2005;84:933–41.
41. Bottje WG, Carstens GE. Association of mitochondrial function and feed efficiency in poultry and livestock species. J Anim Sci. 2009;87(Suppl 14):E48–63.
42. Bottje W, Kong BW. Cell Biology Symposium: feed efficiency: mitochondrial function to global gene expression. J Anim Sci. 2013;91:1582–1593.
43. Sharifabadi HR, Zamiri MJ, Rowghani E, Bottje WG. Relationship between the activity of mitochondrial respiratory chain complexes and feed efficiency in fat-tailed Ghezel lambs. J Anim Sci. 2012;90:1807–15.
44. Kang J, Pervaiz S. Mitochondria: redox metabolism and dysfunction. Biochem Research International. 2012;896751.
45. Tizioto PC, Coutinho LL, Decker JA, Schanbel RD, Rosa KO, Oliveira PSN, et al. Global liver gene expression differences in Nelore steers with divergent residual feed intake phenotypes. BMC Genom. 2015;16:242.
46. Ferguson M, Mockett RJ, Shen Y, Orr WC, Sohal RS. Age-associated decline in mitochondrial respiration and electron transport in Drosophila melanogaster. Biochem J. 2005;390:501–11.
47. Zhang R, Zhu W, Mao S. High-concentrate feeding upregulates the expression of inflammation-related genes in the ruminal epithelium of dairy cattle. J Anim Sci Biotechnol. 2016;7:42.
48. Al-Husseini W, Gondro C, Quinn K, Herd RM, Gibson JP, Chen Y. Expression of candidate genes for residual feed intake in Angus cattle. Anim Genet. 2014;45:12–9.
49. Barendse W, Reverter A, Bunch RJ, Harrison BE, Barris W, Thomas MB. A validated whole-genome association study of efficient food conversion in cattle. Genetics. 2007;176:1893–905.
50. Chen Y, Gondro C, Quinn K, Herd, Parnell PF, Vaneslow B. Global gene expression profiling reveals genes expressed differentially in cattle with high and low residual feed intake. Anim Genet. 2011;42:475–90.
51. Santana MHA, Kadarmideen HN, Pant SD, Alexandre PA, Oliveira GA, Gomes RC, et al. Systems genetics investigations for feed intake, feed efficiency and performance in Nellore (Bos indicus) Cattle. Proceedings of the 10th World Congress of Genetics Applied to Livestock Production. Auckland, 11-18th February 2018.
52. Calderwood SK, Khaleque MA, Sawyer DB, Ciocca DR. Heat shock proteins in cancer: chaperones of tumorigenesis. Trends Biochem Sci. 2006;31:164–72.
53. Dekker SL, Kampinga HH, Bergink S. DNAJs: more than substrate delivery to HSPA. Front Mol Biosci. 2015;2:35.
54. Kong BW, Hudson N, Seo D, Lee S, Khatri B, Lassiter K, 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. 2017;18:82.
55. Keogh K, Kenny DA, Cormican P, Kelly AK, Waters SM. Effect of dietary restriction and subsequent re-alimentation on the transcriptional profile of hepatic tissue in cattle. BMC Genom. 2016;17:244.
56. He X, Kan H, Cai L, Ma Q. Nrf2 is critical in defense against high glucose-induced oxidative damage in cardiomyocytes. J Mol Cell Cardiol. 2009;46:47–58.
57. Ma Q. Role of nrf2 in oxidative stress and toxicity. Annu Rev Pharmacol Toxicol. 2013;53:401–26.
58. Ma Q. Xenobiotic-activated receptors: from transcription to drug metabolism to disease. Chem Res Toxicol. 2008;21:1651–71.
59. Zhou N, Lee WR, Abasht B. Messenger RNA sequencing and pathway analysis provide novel insights into the biological basis of chickens' feed efficiency. BMC Genom. 2015;16:195.
60. Pickett CB, Lu AY. Glutathione S. -transferases: gene structure, regulation, and biological function. Annu Rev Biochem. 1989;58:743–64.
61. Khatib ZA, Inaba T, Valentine M, Look AT. Chromosomal localization and cDNA cloning of the human DBP and TEF genes. Genomics. 1994;23:344–51.
62. Dankel SN, Degerud EM, Borkowski K, Fjaere E, Midtbo LK, Haugen C, et al. Weight cycling promotes fat gain and altered clock gene expression in adipose tissue in C57BL/6J mice. Am J Physiol Endocrinol Metab. 2014;306:E210-24.
63. Kern RJ, Lindholm-Perry AK, Freetly HC, Snelling WM, Kern JW, Keele JW, et al. Transcriptome differences in the rumen of beef steers with variation in feed intake and gain. Gene. 2016;586:12–26.
64. McCarthy SD, Waters SM, Kenny DA, Diskin MG, Fitzpatrick R, Patton J, et al. Negative energy balance and hepatic gene expression patterns in high-yielding dairy cows during the early postpartum period: a global approach. Physiol Genomics. 2010;42:188–99.
65. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:R14.
66. Keogh K, Kenny DA, Cormican P, McCabe MS, Kelly AK, Waters SM. Effect of dietary restriction and subsequent re-alimentation on the transcriptional profile of bovine skeletal muscle. Plos one. 2016;11:e0149373.