Global liver gene expression differences in Nelore steers with divergent residual feed intake phenotypes

Abstract

Background

Efficiency of feed utilization is important for animal production because it can reduce greenhouse gas emissions and improve industry profitability. However, the genetic basis of feed utilization in livestock remains poorly understood. Recent developments in molecular genetics, such as platforms for genome-wide genotyping and sequencing, provide an opportunity to identify genes and pathways that influence production traits. It is known that transcriptional networks influence feed efficiency-related traits such as growth and energy balance. This study sought to identify differentially expressed genes in animals genetically divergent for Residual Feed Intake (RFI), using RNA sequencing methodology (RNA-seq) to obtain information from genome-wide expression profiles in the liver tissues of Nelore cattle.

Results

Differential gene expression analysis between high Residual Feed Intake (HRFI, inefficient) and low Residual Feed Intake (LRFI, efficient) groups was performed to provide insights into the molecular mechanisms that underlie feed efficiency-related traits in beef cattle. A total of 112 annotated genes were identified as being differentially expressed between animals with divergent RFI phenotypes. These genes are involved in ion transport and metal ion binding; act as membrane or transmembrane proteins; and belong to gene clusters that are likely related to the transport and catalysis of molecules through the cell membrane and essential mechanisms of nutrient absorption. Genes with functions in cellular signaling, growth and proliferation, cell death and survival were also differentially expressed. Among the over-represented pathways were drug or xenobiotic metabolism, complement and coagulation cascades, NRF2-mediated oxidative stress, melatonin degradation and glutathione metabolism.

Conclusions

Our data provide new insights and perspectives on the genetic basis of feed efficiency in cattle. Some previously identified mechanisms were supported and new pathways controlling feed efficiency in Nelore cattle were discovered. We potentially identified genes and pathways that play key roles in hepatic metabolic adaptations to oxidative stress such as those involved in antioxidant mechanisms. These results improve our understanding of the metabolic mechanisms underlying feed efficiency in beef cattle and will help develop strategies for selection towards the desired phenotype.

Background

Feed efficiency-related traits are increasingly being studied because of their importance to the overall profitability of animal production. Moreover, the selection of more efficient animals reduces the land required for feed production, methane emissions and nitrogen excretion resulting from the digestion/metabolic process [1-3]. Heritability estimates for feed efficiency-related traits are moderate in dairy and beef cattle [4-7], including the Nelore breed [8]; however, genetic variation for feed efficiency has not been widely exploited in animal breeding programs because the measurement of this trait is costly [1].

There are several indices that are commonly used to estimate the feed efficiency of growing cattle; one of them being residual feed intake (RFI) which is independent of body weight and weight gain. RFI is used to identify individuals that deviate from their expected level of feed intake given their size and growth rate over at least a 70 day feeding period [3]. Because RFI is not phenotypically dependent on the production traits that are used to estimate expected feed intake, it is possible to compare RFI among individuals that differ in their level of production. This independence has led some researchers to believe that RFI may reflect intrinsic variation in basic metabolic processes [9].

Developments in molecular genetics, specifically high-throughput sequencing methods, offer a unique opportunity to identify genes and pathways that are associated with complex traits and diseases [10]. Current DNA and RNA sequencing methodologies are becoming important tools for unravelling the mechanisms which underlie complex traits, facilitating a new understanding of the genetic regulation of phenotype and allowing for the identification of potential biomarkers for early or more accurate genetic prediction. Gene expression profiling can be applied to identify differentially expressed (DE) genes and isoforms involved in networks that control complex traits, thereby shedding some light on the molecular mechanisms responsible for variation in target traits.

Recent studies have identified putative quantitative trait loci (QTL) for feed efficiency on several chromosomes in Nelore populations [8,11]. However, these studies have largely identified discordant genomic regions, revealing a limitation of genome-wide association studies (GWAS) for identifying loci with significant effects within different subpopulations of the same breed [12]. In this research, two divergent groups of Nelore cattle were selected on their best linear unbiased predictions (BLUP) of additive genetic merit for RFI and classified as either high (HRFI) or low (LRFI). RNA sequencing was used to profile the gene expression of hepatic tissue of 20 sampled animals.

Results

Sequencing throughput, read mapping, and assembly

The RFI phenotypes for this Nelore population were previously used to perform a genome-wide association study (GWAS) and the summary statistics for the population were described [8]. Table 1 presents the BLUP estimates of additive genetic merit, phenotypes, sequencing throughput and mapping statistics for each sample used in this study.

Table 1 Best Linear Unbiased Predictions (BLUP) of additive genetic merit for Residual Feed Intake (RFI), dry matter intake (DMI), average daily gain (ADG), sire, number of reads passing filter and concordant pair alignment rate for each animal within the Low (LRFI, efficient) or High (HRFI, inefficient) groups based on RFI BLUP estimates

After mapping reads with TopHat v2.0.6 [13,14], Cufflinks v2.0.2 [14,15] was used to assemble the transcriptome for each sample. The Cuffmerge utility was then run to create a unique file which contained a parsimonious set of transcripts for these data. The number of detected transcripts that represented potentially new isoforms was very large (~71.44% of the transcripts); nevertheless this was expected considering that almost all genes in mammals undergo alternative splicing [16]. We found a total of 16,962 annotated genes to be expressed in bovine liver; however, 5,707 rare or highly expressed (>1 million reads) genes were not tested in the analysis for differential expression. Lowly expressed genes cannot be statistically tested by the Cuffdiff 2 algorithm while the analysis of highly expressed genes leads to excessive machine memory demands [14,15].

To evaluate sequence quality, we assessed the distribution of transcript abundances for each expressed gene as a box-plot of the log of FPKM values (Additional file 1: Figure S1). Very similar median and quartile values for FPKM estimates were observed for the members of both RFI groups. We also evaluated the expression profiles of selected housekeeping genes Hypoxanthine Phosphoribosyltransferase 1 (HPRT1) and Tyrosine 3-Monooxygenase/Tryptophan 5-Monooxygenase Activation Protein, Zeta (YWHAZ) and found expression patterns for these genes to be similar within each of the treatments. Finally, a principal component analysis (PCA) of FPKM values for all genes indicated that there were sufficient numbers of DE genes to differentiate the RFI groups (Additional file 2: Figure S2).

Genome-wide transcriptome analysis and functional annotation

Differential expression analysis between the HRFI (inefficient) and LRFI (efficient) groups identified 112 DE annotated genes. The sign of the log₂(fold change) was used to partition the DE genes into up- and down-regulated groups with 43 DE genes being down-regulated and 69 up-regulated in the LRFI relative to the HRFI groups (Table 2).

Table 2 Genes found to be differentially expressed in the livers of high and low RFI animals

Six genes that were previously identified in a microarray study that profiled gene expression in the livers of Angus cattle selected for high and low RFI [12] were also identified in this study. The coincident genes included collagen, type I, alpha 1 (COL1A1), glutathione S-transferase M1 (GSTM1), regulator of G-protein signaling 2 (RGS2), ring finger protein 150 (RNF150), solute carrier family 2 (facilitated glucose/fructose transporter), member 5 (SLC2A5) and vimentin (VIM).

Other candidate genes previously described as functioning in the determination of traits related to feed efficiency were also found in this analysis [17-22]. For example, the fatty acid-binding protein 1 (FABP1) also known as liver-type fatty acid-binding protein (L-FABP) was up-regulated in the LRFI group. Uncoupling protein 2 (mitochondrial, proton carrier) (UCP2) and fatty acid desaturase 2 (FADS2) with roles in carbohydrate and/or fatty acid metabolism and mitochondrial function were also found to be DE and up-regulated in the LRFI group.

A joint functional annotation analysis using both the up- and down-regulated genes was performed to avoid the potential loss of pathways in which up-regulated genes down-regulate other DE genes and vice versa. When analyzed using Database for Annotation, Visualization, and Integrated Discovery (DAVID) v6.7 using cattle as the background [23], the identified functional gene clusters were related to signal, glycoprotein, glycosylation, membrane or transmembrane region, integral to membrane, transport, metal ion binding, regulation of transcription, among others (Additional file 3: Table S1).

The top bio functions identified by QIAGEN’s Ingenuity® Pathway Analysis (IPA®, QIAGEN Redwood City, CA www.qiagen.com/ingenuity) were involved in cellular movement, represented by 28 genes, including COL1A1; cytochrome b-245, beta polypeptide (CYBB) and UCP2 and in cell-to-cell signaling and interaction, in which 27 genes were reported as related to this function, including early growth response 1 (EGR1), VIM, and FBJ murine osteosarcoma viral oncogene homolog (FOS). Cellular growth and proliferation (represented by 46 genes including connective tissue growth factor (CTGF); FABP1 and FADS2 - Figure 1) and cellular function and maintenance (represented by 23 genes, including surfactant protein A1 (SFTPA1) and transglutaminase 2 (TGM2)) were also observed. These functions were primarily up-regulated in the LRFI group.

Figure 1

The differentially expressed gene network with functions in cellular growth and proliferation. Genes presented in red are up-regulated in the LRFI phenotype group. Genes presented in green are down-regulated in LRFI animals. The intensity of the colors is related to fold change estimates. Arrows presented in orange, gray and yellow indicate activation, effect not predicted and inconsistency, respectively.

Five KEGG database pathways were found by the DAVID software to be over-represented for genes DE between the divergent RFI groups. These pathways were related to drug or xenobiotic metabolism (BH-adj ≤ 0.44 and BH-adj ≤ 0.27, respectively) complement and coagulation cascades (BH-adj ≤ 0.25) and glutathione (BH-adj ≤ 0.48).

The IPA software reported several other significant canonical pathways involving the 112 DE genes, including complement system (P ≤ 2.16E-05), NRF2-mediated oxidative stress (P ≤ 2.16E-05), melatonin degradation (P ≤ 1.54E-04), glutathione-mediated detoxification (P ≤ 2.08E00), IGF-1 signaling (P ≤ 1.28E-02), TGF-β signaling (P ≤ 7.06E-02), glutathione redox reactions I (P ≤ 8.64E-02) and G-Protein coupled receptor signaling (P ≤ 4.44E-01).

The upstream regulatory analysis performed by IPA predicted regulators based on the consistency of expression direction changes for DE genes within each pathway (Additional file 4: Table S2). The most important regulators identified in this analysis were apolipoprotein E (APOE ) (Figure 2; Additional file 5: Table S3), which was predicted to be inhibited in the LRFI group, endothelin-1 (EDN1) (Figure 3; Additional file 6: Table S4) and arachidonic acid (Figure 4; Additional file 7: Table S5) which were predicted to be activated in the LRFI group. Two additional top upstream regulators were inferred: lipopolysaccharide and lysophosphatidic acid, however, it was not possible to infer their activation or inactivation based upon the DE gene set.

Figure 2

The mechanistic network of the inferred upstream regulator APOE . Genes presented in orange are related to genes up-regulated in the LRFI phenotype group. Genes presented in blue are related to genes down-regulated in LRFI animals. The intensity of the colors is related to fold change estimates. Arrows presented in orange, gray and yellow indicate activation, effect not predicted and inconsistency, respectively.

Figure 3

The differentially expressed gene network of the inferred upstream regulator EDN1 . Genes presented in orange are related to genes up-regulated in the LRFI phenotype group. Genes presented in blue are related to genes down-regulated in LRFI animals. The intensity of the colors is related to fold change estimates. Arrows presented in orange, gray and yellow indicate activation, effect not predicted and inconsistency, respectively.

Figure 4

The differentially expressed gene network of the inferred upstream regulator arachidonic acid. Genes presented in orange are related to genes up-regulated in the LRFI phenotype group. Genes presented in blue are related to genes down-regulated in LRFI animals. The intensity of the colors is related to fold change estimates. Arrows presented in orange, gray and yellow indicate activation, effect not predicted and inconsistency, respectively.

The animals comprising the HRFI and LRFI groups were regrouped based on their phenotypes for the component traits dry matter intake (DMI) and average daily gain (ADG). We performed global DE analyses based on these trait groupings (high vs low DMI and ADG) to provide insights into the molecular mechanisms that underlie RFI in Nelore beef cattle.

In order to generate differentiated groups for these traits we reduced the sample size to 12 (6 high and 6 low) and 8 (4 high and 4 low) animals for DMI and ADG, respectively, however, we consequently lost some statistical power. Of the 58 DE genes for DMI, 35 were also identified for RFI and of the 39 DE genes for ADG 18 were also DE for RFI. While ACE2, AGXT2L1, ARHGEF38, CFD, COL1A1, COL1A2, CYP2B6, EGR1, FABP1, FADS2, FAM115C, FKBP5, GLCE, HDAC10, HOOK1, HOPX, IFI27, LOC100848726, LOC100848941, LOC524810, NUFIP1, PCDH7, PTGER3, PYROXD2, RN28S1, RN5-8S1, RNF150, SFRP2, SFTPA1, SIX1, SLC10A7, SLC5A8, SPTSSB and WDR35 are likely related to RFI by influencing DMI, CACNA2D1, CHPF2, CST3, CYR61, FAM115C, FCGR3A, FKBP5, GPC3, HSPB8, IFI27, ISG15, LOC524810, MSR1, RGS2, RNF150, SFRP2, UGT3A1 and WDR35 influence ADG.

Discussion

The profitability of beef cattle production is based on both input expenses and output prices for the final products, and these can be used to compute a selection index for feed efficiency [1]. Feed has a major impact on the total cost of beef production systems. It is known that feed efficiency traits are heritable and have sufficient genetic variation within populations to facilitate selection [4-8]. The artificial selection of efficient animals would potentially reduce the cost of cattle production; however, selection for this trait is not easy to implement because it is challenging and expensive to measure individual feed intake on large samples of animals. Residual feed intake, a measure of feed efficiency of growing cattle, is a complex trait controlled by different metabolic processes [9].

The integration of multiple sources of genetic information could potentially explain additional genetic variation via the elucidation of the molecular mechanisms controlling important production traits. Gene expression is a key source of variation between individuals and may be used to identify functional candidate genes and pathways that control target traits. Genes that have previously been identified as being DE in a study of liver tissues of Angus cattle selected for RFI [12] were also found in our analysis. These include COL1A1, GSTM1, RGS2, RNF150, SLC2A5 and VIM and suggest that common gene networks underlie RFI regardless of breed genetic background.

Glutathione S-transferase enzymes catalyze the conjugation of glutathione to endogenous compounds such as lipid hydroperoxides and exogenous xenobiotics [24]; the liver is a vital organ for xenobiotic metabolism [25]. The exploration of our genome-wide transcriptome results in DAVID revealed xenobiotic and drug metabolism pathways as being overrepresented and up-regulated in the LRFI group. Chen et al. [12] also found xenobiotic metabolism to be an overrepresented pathway for DE genes, but found this pathway to be down-regulated in the LRFI Angus group in contrast to our findings. Besides GSTM1 and glutathione S-transferase omega 1 (GSTO1) found in our study; other members of the Glutathione S- (GST) family were also reported to be DE by Chen et al. [12]. Genes of the cytochrome P450, family 2, subfamily B, polypeptide 6 (CYP2B6) and UDP glucuronosyltransferase 2 family, polypeptide A3 (UGT2A3) families were also detected in this pathway. The CYP family and UGTs, which are primarily expressed in liver, encode several enzymes with a crucial function on oxidative metabolism of endogenous substrates, including steroids, fatty acids and exogenous molecules [26,27]. These gene families are also likely involved in the NRF2-mediated oxidative stress response pathway which was consistently found to be up-regulated in the LRFI group by the IPA. While glutathione S-transferase functions in the detoxification of products of oxidative stress, cytochrome P450 proteins catalyze reactions involved in drug metabolism and the synthesis of cholesterol, steroids, and other lipids [26,27]. Our findings suggest that inefficient animals have increased oxidative metabolism possibly stimulated by an increased oxidative stress.

The NRF2-regulated signaling pathway plays a role in protecting mitochondria from oxidative stress during fasting and ensures the efficient utilization of fatty acids in mouse liver. A study has shown that Nrf2-knckout mice are predicted to diminish oxidation and increase the accumulation of lipids in liver due to mitochondrial damage [28]. These findings are also pertinent to broilers, which suggest that genes involved in glutathione metabolism may influence feed efficiency due their function in preserving or improving the activity of certain respiratory chain complexes [29].

Besides NRF2-mediated oxidative stress, IPA also pointed to other pathways overrepresented for DE genes, including, melatonin degradation, IGF-1 signaling, G-Protein coupled receptor signaling, and in agreement with the DAVID results, glutathione redox reactions. The IGF-1 signaling pathway was also found by Chen et al. [12], however, while they found the IGF-Binding Protein 3 (IGFBP3) gene to be up-regulated in the LRFI group, we found CTGF and CYR61 genes (cysteine-rich, angiogenic inducer, 61) to be DE in this pathway.

Some of the pathways found in this study, such as IGF-1 signaling have already been reported as functioning in feed efficiency-related traits [30]; however, others are new and may elucidate important unknown mechanisms in Nelore cattle. For example, the involvement of the melatonin degradation pathway in RFI is novel and more studies are necessary to elucidate its role and action in feed efficiency in cattle. Melatonin is responsible for controlling several different biological processes such as a combination of cyclic background and circadian rhythm and also for establishing energy balance and maintaining body weight [31,32]. Its role in energy metabolism and obesity is also recognized [31]; however, the weight-reducing effects of melatonin depend on the actions of several mechanisms, including the circadian clock, energy metabolism and metabolic processes [32]. A functional circadian clock and coordinated metabolic processes are necessary to enhance energy balance and maintenance [32].

The genes of cytochrome P450, families 2 e 4, subfamily B (CYP4B1, CYP2B6) and UDP glucuronosyltransferase 2 and 3 families, polypeptide A (UGT3A1 and UGT2A3), primarily up-regulated in the LRFI group, were also involved in melatonin degradation. Melatonin putatively attenuates oxidative stress by decreasing lipid peroxidation [33]. Peroxidation of lipids produces aldehyde products which induce the activation of hepatic stellate cells [34]; the primary collagen-producing cells within the liver. Collagen genes were consistently observed as being up-regulated in the LRFI group. Furthermore, melatonin interactions with reactive species are effective against oxidative stress by improving the function of the mitochondrial respiratory chain [35]. Melatonin can increase the levels of several antioxidative enzymes, including glutathione peroxidase and glutathione reductase [33]. Our findings consistently predict the activation of functions important to oxidative processes in the inefficient LRFI group.

The RGS2 gene was found to be DE between high and low RFI groups in both Nelore and Angus [12] cattle and may affect feed efficiency via its G protein-coupled signaling activity in different cellular functions including the regulation of body weight and adiposity [36]. RGS2-knockout mice had lower weight than wild-type controls and exhibited reduced fat deposition, decreased serum lipids and leptin levels, resulting in a lean phenotype even when fed the same diet as control animals, however, food intake and energy expenditure were not altered possibly due to altered energy balance and defects in metabolic processes and energy storage [36]. We found RGS2 to be up-regulated in the LRFI and in the low ADG groups in agreement with previous reports [12,36]. Furthermore, also supporting our findings, RGS2 expression has been reported to be up-regulated under conditions of oxidative stress [37].

Many of the enriched functional categories reported by DAVID such as ion transport, metal ion binding, membrane or transmembrane proteins are likely related to the catalysis and transport of substrates through the cell membrane [38]. Transport of substances across cell membranes is required for several vital functions including digestion, absorption of nutrients, cellular signaling, growth, proliferation, cell death and survival which have previously been reported as influencing feed efficiency traits in beef cattle [39]. Some of these biological functions were also found to be enriched for DE genes by the IPA software. Members of the solute carrier group, which are primarily located in the cell membrane (SLC10A7, SLC2A5, SLC41A2, SLC45A3 and SLC5A8), were found to be down-regulated in the HRFI group. The SLC2A5 gene, which facilitates glucose/fructose transport, was found to be the top up-regulated gene in the HRFI group while genes among the most down-regulated in this group were related to lipid catalysis. These results suggest that efficient animals may have an increased ability to absorb glucose, while inefficient individuals overexpress genes related to the catalysis and intracellular transport of fatty acids. This may indicate that the divergent efficiency groups have preferable sources for obtaining the energy required for maintenance.

Feed intake may influence metabolic activity in liver and consequently energy utilization [18]. Kuhla et al. [18] reported a significant down-regulation of FABP1 protein in dairy cows that experienced feed restriction and suggested that this may provide a mechanism for limiting fatty acid oxidation and hepatic triacylglyceride accumulation in the event of negative energy balance. These results are supported by a study in which FABP1 knockout mice demonstrated considerably reduced triacylglyceride levels in liver after fasting [17]. The pattern of FADS2 gene expression is known to regulate the synthesis of polyunsaturated fatty acids. Moreover, FADS2-deficient mice are resistant to obesity and the dysregulation of lipogenesis [20]. This gene may be also important to the peroxidation susceptibility of lipoproteins and their oxidation rate [40] and was up-regulated in the animals with the highest DMI.

The upstream regulatory analysis performed by IPA, which seeks to identify the upstream transcriptional regulatory cascades that are likely to elucidate the observed changes in gene expression [41], may shed some light on the biological activities that occur in the hepatic tissue of animals that are genetically divergent for RFI. This analysis predicted the top upstream regulators to include APOE which was predicted to be inhibited in the LRFI group. The APOE protein functions in lipid transport in liver by assisting in the secretion of very low density lipoprotein (VLDL) [42,43]. Takahashi et al. [43] proposed that serum APOE contents of triglyceride-rich lipoproteins must be controlled by dietary handling in cattle. Wilcox & Heimberg [44] have shown that fasting rats had lower secretions of both VLDL and APOE, therefore having a reduced uptake of VLDL by the liver as compared to fed animals. The inhibition of APOE predicted in the LRFI group may be related to the accumulation of lipoproteins in the liver under conditions of oxidative stress. In a previously performed GWAS study in this population [8], the candidate gene Apolipoprotein A2 (APOA2) which functions to stabilize HDL was detected as being associated with RFI.

EDN1 was also predicted by IPA to be a top upstream regulator of RFI and our results suggest that it is activated in the LRFI group since nine of the eleven DE genes regulated by EDN1 were found to be coactivated. EDN1 was inferred by IPA to be a potential regulator of connective tissue growth factor and early growth response genes such as CTGF and EGR1. Additionally, seven DE genes had expression profiles that were consistent with the activation of arachidonic acid in the LRFI group. These include FABP1 [45], UCP2 [46] and RGS2 [47] which must now be investigated as targets for manipulation through diets containing arachidonic acid. Furthermore, the relative proportion of dietary arachidonic acid to docosahexaenoic acid has been shown to be a determinant of FADS2 expression and consequently influences polyunsaturated fatty acids metabolism in suckling piglets [48].

Despite the fact that genes found to be DE in this study were not detected in the QTL regions found in a previous GWAS study using the same Nelore population [8]; several common biological mechanisms and key drivers were detected. The majority of QTLs identified in the GWAS lies within gene deserts and may affect feed efficiency via regulatory elements that are yet to be identified involved in the modulation of expression of genes. Non-coding functional elements are poorly understood in cattle and can consist of distal enhancers or transcription factor binding sites. The challenge to interpreting the roles of these QTLs lies in the diversity of function of non-coding variants, poor annotation of regulatory elements and potentially unrecognized control mechanisms [49]. However, candidate genes identified in GWAS are known to cause the DE of genes; when an integrated analysis including both GWAS QTLs and RNA-Seq DE genes was performed using IPA, the differentially expressed transcription factors EGR1 and FOS were suggested to be regulating the candidate gene Plasminogen Activator, Urokinase (PLAU) located within a QTL region identified by the GWAS. On the other hand, this gene seems to also coregulate the VIM [50] and CYR61 [51] genes. EGR1 and FOS, found to be up-regulated in the LRFI, are key regulators of genes that are related to cellular growth and differentiation and are also known to be activated in response to oxidative stress [52,53]. Studies targeting the identification of regulatory mutations within the promoters and enhancers/repressors of these genes may be important for understanding the biology of feed efficiency and may have utility for the implementation of genomic selection for feed efficiency in livestock.

Although QTL regions do not have to harbor DE genes, since they can be created by mutations in genes that cause post-translational disruptions affecting the functionality of proteins. The differences in candidate regions/genes found by the GWAS and RNA-Seq may also be explained by the tissue-specific modulation of messenger RNAs (mRNAs). For example, the HRH4 and ADAM12 candidate genes located within a QTL region detected by the GWAS could not be tested for expression differences in this study due to their low expression in liver. This finding does not exclude the implication of the DE for these genes in other tissues on feed efficiency.

Conclusions

We conducted a genome-wide transcriptome profiling study of hepatic tissue from Nelore cattle selected to be genetically divergent for RFI to reveal key metabolic and cell signaling networks. Some previously known mechanisms related to feed efficiency such as xenobiotic metabolism were found; however, new pathways including melatonin degradation were also identified as controlling RFI in Nelore cattle. Overall, our findings demonstrate that changes in gene expression between efficient and inefficient cattle primarily appear to be related to metabolic processes underlying oxidative stress and lipid catabolism. We have potentially identified genes involved in antioxidant mechanisms that play key roles in hepatic metabolic adaptation to oxidative stress. Previous studies have suggested that oxidative stress is increased in inefficient broilers and that this may be related to differences in mitochondrial function [54]. Metabolic response to negative energy balance depends on the availability of fatty acids and ketones as energy sources as well as to the mitochondrial capacity for fatty acid oxidation in tissues with high oxidative energy demands such as liver [55]. The upstream regulators found here guide the future investigation of these molecules to enable the development of intervention strategies such as diet formulation and contribute to the understanding of the physiology and improvement of RFI.

Methods

Animals and sampling

All experimental procedures were approved by the Institutional Animal Care and Use Committee Guidelines of the Brazilian Agricultural Research Corporation – EMBRAPA and were sanctioned by the president, Dr. Rui Machado.

These steers comprised half-sib families produced by the artificial insemination of commercial and purebred Nelore dams, derived from 18 sires representing the main breeding lineages commercialized in Brazil. The 83 calves used in this expression study were allocated to feedlots in Embrapa Southeast Research at about 21 months of age. Within the feedlots, animals were maintained either in individual or collective pens and allowed ad libitum access to feed and water as described by Oliveira et al. [8]. Briefly, animals were fed twice daily, with diets formulated to contain 40% dry matter (DM) in the form of corn silage; crude protein at 13.5% and energy densities of 2.8. The remaining 60% of DM was concentrate, which comprised ground corn, soybean meal, cotton seed, soybean hulls, limestone, mineral mixture, urea and monensin (Rumensin®). Measures of daily feed intake were collected for at least 70 days and body weight was measured every 14 days.

BLUP estimates of genetic merit for RFI were generated for 585 Nelore steers. Liver samples were available for only 83 of the animals which were ranked according to their additive genetic merit for RFI to select 20 animals that were genetically divergent for RFI, as described below. A relationship matrix computed using pedigree information was used in this analysis. Nelore steers that were genetically divergent for RFI (kg/d) were selected based on BLUP estimates of their additive genetic merits produced using the following model:

$$ y = X\beta + Za + \varepsilon $$

Where, y is the vector for average daily feed intake, β is the vector of fixed effects of contemporary group, defined as the combination of season, animal origin and pen type (individual or collective), and partial regressions on age of the animal at entrance to the feedlot, metabolic mid-weight (BW^0.75) and average daily gain, a is the additive genetic merit of the animal for RFI assumed to be normally distributed with E[a] = 0 and Var(a) = $\mathrm{A}{\upsigma}_a^2$ where A is the pedigree numerator relationship matrix, and ε is the vector of residual effects inherent to each observation which was assumed to be normally and independently distributed (0, ${\upsigma}_e^2$), X and Z are design matrices for fixed and random effects, respectively. The model was fit by the MIXED procedure of SAS® software; version 9.3 (SAS Institute Inc.) and selected animals were ranked in the most extreme values for additive genetic merit. Where possible, animals that had common sires were sampled from each end of the BLUP distribution.

Dry matter intake and average daily gain described elsewhere [8] were used to decompose RFI via the regrouping of the animals based on these traits for additional gene expression analyses.

RNA sequencing

Preparation of the mRNA samples for sequencing was performed by ESALQ Genomics Center (Piracicaba, São Paulo, Brazil), using the TruSeq RNA Sample Preparation Kit® (Illumina, San Diego, CA) according to manufacturer’s instructions. Briefly, 100 mg of frozen liver was used to extract RNA using the TRIzol® reagent (Life Technologies, Carlsbad, CA) and 2 μg of total RNA from each liver sample was used for library preparation. The concentration and purity of RNA was measured using NanoDrop™ (Thermos Scientific, Waltham, MA) and then sample integrity was assessed by Bioanalyzer (Agilent, Santa Clara, CA). The mRNA was first enriched from the total RNA by using oligo dT magnetic beads, then the poly(A) RNA was fragmented and cDNA was synthesized. Next, the cDNA underwent end repair, the 3’ ends were adenylated and universal bar-coded adapters were ligated to the cDNA fragments to perform a solid phase PCR to produce the sequencing library. Following library construction, the sequencing library was evaluated and quantified using both an Agilent 2100 Bioanalyzer® and quantitative PCR with the KAPA Library Quantification kit® (KAPA Biosystems, Foster City, CA, USA). Finally, libraries were pooled to perform multiplexing sequencing. Cluster generation and sequencing were performed on the Illumina HiSeq 2000®. Paired-end reads of 2 × 100 bp were produced.

Processing and alignment of sequence reads

Computations were performed on the HPC resources at the University of Missouri Bioinformatics Consortium (UMBC). Low-quality reads were filtered and adapter sequences trimmed using SeqClean software. TopHat v2.0.6 [13,14] was then used to align the reads to the Bos taurus virtual transcriptome internally built by Tophat using the UMD3.1 reference genome. TopHat first extracted the transcript sequences and used Bowtie to align reads to the virtual transcriptome using a provided reference annotation file. The reads that could not be fully mapped to the virtual transcriptome were then mapped to the UMD3.1 reference genome. These reads were converted into genomic mappings and merged with the novel transcriptome mappings and splice junctions in the final output file. A total of 2 mismatches per read were allowed in alignment.

Transcript assembly and quantification

Cufflinks v2.0.2 [15] was initially used to assemble the aligned reads for each sample individually. Cufflinks assembles the aligned reads and provides a parsimonious set of transcripts as a file. Cufflinks also estimates transcript abundances in Fragments Per Kilobase of exon per Million fragments mapped (FPKM), which normalizes transcript expression for transcript length and the total number of sequence reads per sample. The reference annotation supplied to Cufflinks was used to perform a reference annotation-based transcript assembly. The output for each sample included all reference transcripts as well as novel assembled genes and isoforms. Cufflinks assemblies for all samples were then merged using Cuffmerge v2.0.2 which also runs Cuffcompare internally to classify the transcripts. The available annotation file was provided to this analysis to classify the assembled contigs into novel and known transcripts and to maximize the overall quality of the assembly.

Testing for differential expression

Cuffdiff2 software was run to test for DE genes between the RFI groups with geometric normalization used to estimate transcript abundance. Correction for multiple testing (q value) was performed using the Benjamini-Hochberg methodology. Cuffdiff2 calculated the FPKM for each transcript, primary transcript, and gene in each sample. A false discovery rate ≤ 0.05 was adopted to consider a gene as being DE.

Data exploration and visualization was performed using the CummeRbund package [14] implemented in the R programming environment.

Annotation of differentially expressed genes

DAVID v6.7 [23] was used to annotate and interpret the DE gene lists. DAVID software identifies enriched biological themes and gene ontology (GO) terms, clusters functionally related genes and annotation terms for gene lists with EASE scores < 0.1. The Functional Annotation Tool was used to determine the most relevant GO terms within each list of DE genes. The Functional Annotation Clustering algorithm was used to generate a report of related annotation terms and groups of annotation clusters. Finally, DAVID Pathway was used to map the enriched pathways in which DE genes are involved, using the KEGG database.

The IPA (www.qiagen.com/ingenuity) was also used to discover and explore biological processes and the roles of DE genes. The Ingenuity Pathways Knowledge Base comprises relationships such as between genes, mRNAs and proteins to test for significantly overrepresented networks and pathways. We provided the fold changes and q-values of DE among genes from the Cuffdiff analysis to the IPA to perform the statistical analysis for the representation of each network and to visualize the results.

Availability of supporting data

The RNA-seq data sets supporting the results of this study are available in the ENA repository (EMBL-EBI), under accession PRJEB7696.