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  <cit> .
heritability estimates for feed efficiency-related traits are moderate in dairy and beef cattle  <cit> , including the nelore breed  <cit> ; however, genetic variation for feed efficiency has not been widely exploited in animal breeding programs because the measurement of this trait is costly  <cit> .
there are several indices that are commonly used to estimate the feed efficiency of growing cattle; one of them being residual feed intake  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  <cit> .
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  <cit> .
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  <cit> .
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  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  for feed efficiency on several chromosomes in nelore populations  <cit> .
however, these studies have largely identified discordant genomic regions, revealing a limitation of genome-wide association studies  for identifying loci with significant effects within different subpopulations of the same breed  <cit> .
in this research, two divergent groups of nelore cattle were selected on their best linear unbiased predictions  of additive genetic merit for rfi and classified as either high  or low .
rna sequencing was used to profile the gene expression of hepatic tissue of  <dig> 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  and the summary statistics for the population were described  <cit> .
table  <dig> presents the blup estimates of additive genetic merit, phenotypes, sequencing throughput and mapping statistics for each sample used in this study.table  <dig>
best linear unbiased predictions  of additive genetic merit for residual feed intake , dry matter intake , average daily gain , sire, number of reads passing filter and concordant pair alignment rate for each animal within the low  or high  groups based on rfi blup estimates
animal_id
phenotype
blup
rfi
dmi
adg
sire
reads passing filter
concordant pair alignment rate
after mapping reads with tophat v <dig> .
<dig>  <cit> , cufflinks v <dig> .
<dig>  <cit>  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 ; nevertheless this was expected considering that almost all genes in mammals undergo alternative splicing  <cit> .
we found a total of  <dig>  annotated genes to be expressed in bovine liver; however,  <dig>  rare or highly expressed  genes were not tested in the analysis for differential expression.
lowly expressed genes cannot be statistically tested by the cuffdiff  <dig> algorithm while the analysis of highly expressed genes leads to excessive machine memory demands  <cit> .
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 .
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  <dig>  and tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, zeta  and found expression patterns for these genes to be similar within each of the treatments.
finally, a principal component analysis  of fpkm values for all genes indicated that there were sufficient numbers of de genes to differentiate the rfi groups .
genome-wide transcriptome analysis and functional annotation
differential expression analysis between the hrfi  and lrfi  groups identified  <dig> de annotated genes.
the sign of the log <dig> was used to partition the de genes into up- and down-regulated groups with  <dig> de genes being down-regulated and  <dig> up-regulated in the lrfi relative to the hrfi groups .table  <dig>
genes found to be differentially expressed in the livers of high and low rfi animals
gene id
locus
mean hrfi, efficient
mean lrfi, inefficient
log
2
*
p-value
q-value
abca3
ace2
acta2
agxt2l1
akr7a3
arhgef38
atp2a2
c1qa
c1qc
c28h10orf57
ca3
cacna2d1
cfd
chpf2
chrne
ckb
col1a1
col1a2
cr2
creld2
cst3
ctgf
cybb
cyp2b6
cyp4b1
cyr61
egr1
eppk1
ero1lb
fabp1
fads2
fam115c
fam174b
fam47e
fbxl14
fcgr3a
fkbp5
folr2
fos
gale
gcsh
glce
gng11
gpc3
gpx3
gstm1
gsto1
hdac10
hebp2
hnf4g
hook1
hopx
hspb8
hyou1
ifi27
ifitm3
irf6
isg15
loc100847320
loc100848726
loc100848941
loc510860
loc524810
loc540627
loc786073
lrrc25
lst1
mir365-2
mknk1
msr1
myom1
npc2
nufip1
pcdh7
pcsk5
pgcp
prune2
ptger3
pyroxd2
rgs2
rn28s1
rn5-8s1
rnase6
rnf150
robo2
s100a11
sall1
sell
sfrp2
sftpa1
siglec12
six1
slc10a7
slc2a5
slc41a2
slc45a3
slc5a8
smad1
sptssb
tcirg1
tgm2
tm4sf5
tmsb10
tmsb4
tnfrsf8
ucp2
uggt1
ugt2a3
ugt3a1
vim
wdr35
wfdc2
*fold change estimates are relative to lrfi, inefficient group.
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  <cit>  were also identified in this study.
the coincident genes included collagen, type i, alpha  <dig> , glutathione s-transferase m <dig> , regulator of g-protein signaling  <dig> , ring finger protein  <dig> , solute carrier family  <dig> , member  <dig>  and vimentin .
other candidate genes previously described as functioning in the determination of traits related to feed efficiency were also found in this analysis  <cit> .
for example, the fatty acid-binding protein  <dig>  also known as liver-type fatty acid-binding protein  was up-regulated in the lrfi group.
uncoupling protein  <dig>   and fatty acid desaturase  <dig>  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  v <dig>  using cattle as the background  <cit> , 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 .
the top bio functions identified by qiagen’s ingenuity® pathway analysis  were involved in cellular movement, represented by  <dig> genes, including col1a1; cytochrome b- <dig>  beta polypeptide  and ucp <dig> and in cell-to-cell signaling and interaction, in which  <dig> genes were reported as related to this function, including early growth response  <dig> , vim, and fbj murine osteosarcoma viral oncogene homolog .
cellular growth and proliferation ; fabp <dig> and fads <dig> - figure 1) and cellular function and maintenance  and transglutaminase  <dig> ) were also observed.
these functions were primarily up-regulated in the lrfi group.figure  <dig>
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  complement and coagulation cascades  and glutathione .
the ipa software reported several other significant canonical pathways involving the  <dig> de genes, including complement system , nrf2-mediated oxidative stress , melatonin degradation , glutathione-mediated detoxification , igf- <dig> signaling , tgf-β signaling , glutathione redox reactions i  and g-protein coupled receptor signaling .
the upstream regulatory analysis performed by ipa predicted regulators based on the consistency of expression direction changes for de genes within each pathway .
the most important regulators identified in this analysis were apolipoprotein e  , which was predicted to be inhibited in the lrfi group, endothelin- <dig>   and arachidonic acid  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  <dig>
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.
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.
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  and average daily gain .
we performed global de analyses based on these trait groupings  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  <dig>  and  <dig>  animals for dmi and adg, respectively, however, we consequently lost some statistical power.
of the  <dig> de genes for dmi,  <dig> were also identified for rfi and of the  <dig> de genes for adg  <dig> were also de for rfi.
while ace <dig>  agxt2l <dig>  arhgef <dig>  cfd, col1a <dig>  col1a <dig>  cyp2b <dig>  egr <dig>  fabp <dig>  fads <dig>  fam115c, fkbp <dig>  glce, hdac <dig>  hook <dig>  hopx, ifi <dig>  loc <dig>  loc <dig>  loc <dig>  nufip <dig>  pcdh <dig>  ptger <dig>  pyroxd <dig>  rn28s <dig>  rn5-8s <dig>  rnf <dig>  sfrp <dig>  sftpa <dig>  six <dig>  slc10a <dig>  slc5a <dig>  sptssb and wdr <dig> are likely related to rfi by influencing dmi, cacna2d <dig>  chpf <dig>  cst <dig>  cyr <dig>  fam115c, fcgr3a, fkbp <dig>  gpc <dig>  hspb <dig>  ifi <dig>  isg <dig>  loc <dig>  msr <dig>  rgs <dig>  rnf <dig>  sfrp <dig>  ugt3a <dig> and wdr <dig> 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  <cit> .
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  <cit> .
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  <cit> .
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  <cit>  were also found in our analysis.
these include col1a <dig>  gstm <dig>  rgs <dig>  rnf <dig>  slc2a <dig> 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  <cit> ; the liver is a vital organ for xenobiotic metabolism  <cit> .
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.
<cit>  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 gstm <dig> and glutathione s-transferase omega  <dig>  found in our study; other members of the glutathione s-  family were also reported to be de by chen et al.
<cit> .
genes of the cytochrome p <dig>  family  <dig>  subfamily b, polypeptide  <dig>  and udp glucuronosyltransferase  <dig> family, polypeptide a <dig>  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  <cit> .
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 p <dig> proteins catalyze reactions involved in drug metabolism and the synthesis of cholesterol, steroids, and other lipids  <cit> .
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  <cit> .
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  <cit> .
besides nrf2-mediated oxidative stress, ipa also pointed to other pathways overrepresented for de genes, including, melatonin degradation, igf- <dig> signaling, g-protein coupled receptor signaling, and in agreement with the david results, glutathione redox reactions.
the igf- <dig> signaling pathway was also found by chen et al.
<cit> , however, while they found the igf-binding protein  <dig>  gene to be up-regulated in the lrfi group, we found ctgf and cyr <dig> genes  to be de in this pathway.
some of the pathways found in this study, such as igf- <dig> signaling have already been reported as functioning in feed efficiency-related traits  <cit> ; 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  <cit> .
its role in energy metabolism and obesity is also recognized  <cit> ; however, the weight-reducing effects of melatonin depend on the actions of several mechanisms, including the circadian clock, energy metabolism and metabolic processes  <cit> .
a functional circadian clock and coordinated metabolic processes are necessary to enhance energy balance and maintenance  <cit> .
the genes of cytochrome p <dig>  families  <dig> e  <dig>  subfamily b  and udp glucuronosyltransferase  <dig> and  <dig> families, polypeptide a , primarily up-regulated in the lrfi group, were also involved in melatonin degradation.
melatonin putatively attenuates oxidative stress by decreasing lipid peroxidation  <cit> .
peroxidation of lipids produces aldehyde products which induce the activation of hepatic stellate cells  <cit> ; 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  <cit> .
melatonin can increase the levels of several antioxidative enzymes, including glutathione peroxidase and glutathione reductase  <cit> .
our findings consistently predict the activation of functions important to oxidative processes in the inefficient lrfi group.
the rgs <dig> gene was found to be de between high and low rfi groups in both nelore and angus  <cit>  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  <cit> .
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  <cit> .
we found rgs <dig> to be up-regulated in the lrfi and in the low adg groups in agreement with previous reports  <cit> .
furthermore, also supporting our findings, rgs <dig> expression has been reported to be up-regulated under conditions of oxidative stress  <cit> .
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  <cit> .
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  <cit> .
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 , were found to be down-regulated in the hrfi group.
the slc2a <dig> 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  <cit> .
kuhla et al.
<cit>  reported a significant down-regulation of fabp <dig> 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 fabp <dig> knockout mice demonstrated considerably reduced triacylglyceride levels in liver after fasting  <cit> .
the pattern of fads <dig> 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  <cit> .
this gene may be also important to the peroxidation susceptibility of lipoproteins and their oxidation rate  <cit>  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  <cit> , 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   <cit> .
takahashi et al.
<cit>  proposed that serum apoe contents of triglyceride-rich lipoproteins must be controlled by dietary handling in cattle.
wilcox & heimberg  <cit>  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  <cit> , the candidate gene apolipoprotein a <dig>  which functions to stabilize hdl was detected as being associated with rfi.
edn <dig> 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 edn <dig> were found to be coactivated.
edn <dig> was inferred by ipa to be a potential regulator of connective tissue growth factor and early growth response genes such as ctgf and egr <dig>  additionally, seven de genes had expression profiles that were consistent with the activation of arachidonic acid in the lrfi group.
these include fabp <dig>  <cit> , ucp <dig>  <cit>  and rgs <dig>  <cit>  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 fads <dig> expression and consequently influences polyunsaturated fatty acids metabolism in suckling piglets  <cit> .
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  <cit> ; 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  <cit> .
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 egr <dig> and fos were suggested to be regulating the candidate gene plasminogen activator, urokinase  located within a qtl region identified by the gwas.
on the other hand, this gene seems to also coregulate the vim  <cit>  and cyr <dig>  <cit>  genes.
egr <dig> 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  <cit> .
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 .
for example, the hrh <dig> and adam <dig> 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  <cit> .
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  <cit> .
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.

 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. heritability estimates for feed efficiency-related traits are moderate in dairy and beef cattle, including the nelore breed; however, genetic variation for feed efficiency has not been widely exploited in animal breeding programs because the measurement of this trait is costly. there are several indices that are commonly used to estimate the feed efficiency of growing cattle; one of them being residual feed intake 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. 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. 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. 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 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 for feed efficiency on several chromosomes in nelore populations. however, these studies have largely identified discordant genomic regions, revealing a limitation of genome-wide association studies for identifying loci with significant effects within different subpopulations of the same breed. in this research, two divergent groups of nelore cattle were selected on their best linear unbiased predictions of additive genetic merit for rfi and classified as either high or low. rna sequencing was used to profile the gene expression of hepatic tissue of sampled animals.  sequencing throughput, read mapping, and assembly the rfi phenotypes for this nelore population were previously used to perform a genome-wide association study and the summary statistics for the population were described. table  presents the blup estimates of additive genetic merit, phenotypes, sequencing throughput and mapping statistics for each sample used in this study.table best linear unbiased predictions of additive genetic merit for residual feed intake, dry matter intake, average daily gain, sire, number of reads passing filter and concordant pair alignment rate for each animal within the low or high groups based on rfi blup estimates animal_id phenotype blup rfi dmi adg sire reads passing filter concordant pair alignment rate after mapping reads with tophat v., cufflinks v. 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; nevertheless this was expected considering that almost all genes in mammals undergo alternative splicing. we found a total of annotated genes to be expressed in bovine liver; however, rare or highly expressed genes were not tested in the analysis for differential expression. lowly expressed genes cannot be statistically tested by the cuffdiff algorithm while the analysis of highly expressed genes leads to excessive machine memory demands. 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. 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 and tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, zeta and found expression patterns for these genes to be similar within each of the treatments. finally, a principal component analysis of fpkm values for all genes indicated that there were sufficient numbers of de genes to differentiate the rfi groups. genome-wide transcriptome analysis and functional annotation differential expression analysis between the hrfi and lrfi groups identified de annotated genes. the sign of the log was used to partition the de genes into up- and down-regulated groups with de genes being down-regulated and up-regulated in the lrfi relative to the hrfi groups.table genes found to be differentially expressed in the livers of high and low rfi animals gene id locus mean hrfi, efficient mean lrfi, inefficient log 2 * p-value q-value abca3 ace2 acta2 agxt2l1 akr7a3 arhgef38 atp2a2 c1qa c1qc c28h10orf57 ca3 cacna2d1 cfd chpf2 chrne ckb col1a1 col1a2 cr2 creld2 cst3 ctgf cybb cyp2b6 cyp4b1 cyr61 egr1 eppk1 ero1lb fabp1 fads2 fam115c fam174b fam47e fbxl14 fcgr3a fkbp5 folr2 fos gale gcsh glce gng11 gpc3 gpx3 gstm1 gsto1 hdac10 hebp2 hnf4g hook1 hopx hspb8 hyou1 ifi27 ifitm3 irf6 isg15 loc100847320 loc100848726 loc100848941 loc510860 loc524810 loc540627 loc786073 lrrc25 lst1 mir365-2 mknk1 msr1 myom1 npc2 nufip1 pcdh7 pcsk5 pgcp prune2 ptger3 pyroxd2 rgs2 rn28s1 rn5-8s1 rnase6 rnf150 robo2 s100a11 sall1 sell sfrp2 sftpa1 siglec12 six1 slc10a7 slc2a5 slc41a2 slc45a3 slc5a8 smad1 sptssb tcirg1 tgm2 tm4sf5 tmsb10 tmsb4 tnfrsf8 ucp2 uggt1 ugt2a3 ugt3a1 vim wdr35 wfdc2 *fold change estimates are relative to lrfi, inefficient group. 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 were also identified in this study. the coincident genes included collagen, type i, alpha, glutathione s-transferase m, regulator of g-protein signaling, ring finger protein, solute carrier family, member and vimentin. other candidate genes previously described as functioning in the determination of traits related to feed efficiency were also found in this analysis. for example, the fatty acid-binding protein also known as liver-type fatty acid-binding protein was up-regulated in the lrfi group. uncoupling protein and fatty acid desaturase 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 v using cattle as the background, 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. the top bio functions identified by qiagen’s ingenuity® pathway analysis were involved in cellular movement, represented by genes, including col1a1; cytochrome b- beta polypeptide and ucp and in cell-to-cell signaling and interaction, in which genes were reported as related to this function, including early growth response, vim, and fbj murine osteosarcoma viral oncogene homolog. cellular growth and proliferation; fabp and fads - figure 1) and cellular function and maintenance and transglutaminase ) were also observed. these functions were primarily up-regulated in the lrfi group.figure 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 complement and coagulation cascades and glutathione. the ipa software reported several other significant canonical pathways involving the de genes, including complement system, nrf2-mediated oxidative stress, melatonin degradation, glutathione-mediated detoxification, igf- signaling, tgf-β signaling, glutathione redox reactions i and g-protein coupled receptor signaling. the upstream regulatory analysis performed by ipa predicted regulators based on the consistency of expression direction changes for de genes within each pathway. the most important regulators identified in this analysis were apolipoprotein e, which was predicted to be inhibited in the lrfi group, endothelin- and arachidonic acid 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 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. 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. 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 and average daily gain. we performed global de analyses based on these trait groupings 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 and animals for dmi and adg, respectively, however, we consequently lost some statistical power. of the de genes for dmi, were also identified for rfi and of the de genes for adg were also de for rfi. while ace agxt2l arhgef cfd, col1a col1a cyp2b egr fabp fads fam115c, fkbp glce, hdac hook hopx, ifi loc loc loc nufip pcdh ptger pyroxd rn28s rn5-8s rnf sfrp sftpa six slc10a slc5a sptssb and wdr are likely related to rfi by influencing dmi, cacna2d chpf cst cyr fam115c, fcgr3a, fkbp gpc hspb ifi isg loc msr rgs rnf sfrp ugt3a and wdr 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. 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. 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. 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 were also found in our analysis. these include col1a gstm rgs rnf slc2a 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; the liver is a vital organ for xenobiotic metabolism. 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. 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 gstm and glutathione s-transferase omega found in our study; other members of the glutathione s- family were also reported to be de by chen et al.. genes of the cytochrome p family subfamily b, polypeptide and udp glucuronosyltransferase family, polypeptide a 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. 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 p proteins catalyze reactions involved in drug metabolism and the synthesis of cholesterol, steroids, and other lipids. 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. 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. besides nrf2-mediated oxidative stress, ipa also pointed to other pathways overrepresented for de genes, including, melatonin degradation, igf- signaling, g-protein coupled receptor signaling, and in agreement with the david results, glutathione redox reactions. the igf- signaling pathway was also found by chen et al., however, while they found the igf-binding protein gene to be up-regulated in the lrfi group, we found ctgf and cyr genes to be de in this pathway. some of the pathways found in this study, such as igf- signaling have already been reported as functioning in feed efficiency-related traits; 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. its role in energy metabolism and obesity is also recognized; however, the weight-reducing effects of melatonin depend on the actions of several mechanisms, including the circadian clock, energy metabolism and metabolic processes. a functional circadian clock and coordinated metabolic processes are necessary to enhance energy balance and maintenance. the genes of cytochrome p families e subfamily b and udp glucuronosyltransferase and families, polypeptide a, primarily up-regulated in the lrfi group, were also involved in melatonin degradation. melatonin putatively attenuates oxidative stress by decreasing lipid peroxidation. peroxidation of lipids produces aldehyde products which induce the activation of hepatic stellate cells; 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. melatonin can increase the levels of several antioxidative enzymes, including glutathione peroxidase and glutathione reductase. our findings consistently predict the activation of functions important to oxidative processes in the inefficient lrfi group. the rgs gene was found to be de between high and low rfi groups in both nelore and angus 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. 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. we found rgs to be up-regulated in the lrfi and in the low adg groups in agreement with previous reports. furthermore, also supporting our findings, rgs expression has been reported to be up-regulated under conditions of oxidative stress. 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. 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. 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, were found to be down-regulated in the hrfi group. the slc2a 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. kuhla et al. reported a significant down-regulation of fabp 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 fabp knockout mice demonstrated considerably reduced triacylglyceride levels in liver after fasting. the pattern of fads 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. this gene may be also important to the peroxidation susceptibility of lipoproteins and their oxidation rate 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, 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. takahashi et al. proposed that serum apoe contents of triglyceride-rich lipoproteins must be controlled by dietary handling in cattle. wilcox & heimberg 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, the candidate gene apolipoprotein a which functions to stabilize hdl was detected as being associated with rfi. edn 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 edn were found to be coactivated. edn was inferred by ipa to be a potential regulator of connective tissue growth factor and early growth response genes such as ctgf and egr additionally, seven de genes had expression profiles that were consistent with the activation of arachidonic acid in the lrfi group. these include fabp, ucp and rgs 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 fads expression and consequently influences polyunsaturated fatty acids metabolism in suckling piglets. 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; 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. 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 egr and fos were suggested to be regulating the candidate gene plasminogen activator, urokinase located within a qtl region identified by the gwas. on the other hand, this gene seems to also coregulate the vim and cyr genes. egr 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. 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. for example, the hrh and adam 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.  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. 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. 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. 
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. 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. we potentially identified genes and pathways that play key roles in hepatic metabolic adaptations to oxidative stress such as those involved in antioxidant mechanisms. among the over-represented pathways were drug or xenobiotic metabolism, complement and coagulation cascades, nrf2-mediated oxidative stress, melatonin degradation and glutathione metabolism. 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. it is known that transcriptional networks influence feed efficiency-related traits such as growth and energy balance. differential gene expression analysis between high residual feed intake and low residual feed intake groups was performed to provide insights into the molecular mechanisms that underlie feed efficiency-related traits in beef cattle. this study sought to identify differentially expressed genes in animals genetically divergent for residual feed intake, using rna sequencing methodology to obtain information from genome-wide expression profiles in the liver tissues of nelore cattle. some previously identified mechanisms were supported and new pathways controlling feed efficiency in nelore cattle were discovered. efficiency of feed utilization is important for animal production because it can reduce greenhouse gas emissions and improve industry profitability. genes with functions in cellular signaling, growth and proliferation, cell death and survival were also differentially expressed.
 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, using rna sequencing methodology to obtain information from genome-wide expression profiles in the liver tissues of nelore cattle.  differential gene expression analysis between high residual feed intake and low residual feed intake groups was performed to provide insights into the molecular mechanisms that underlie feed efficiency-related traits in beef cattle. a total of 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.  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. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords bos indicusrfifeed efficiencytranscriptomicsissue-copyright-statement© the author 2015 