Advertisement

Mitochondrial protein gene expression and the oxidative phosphorylation pathway associated with feed efficiency and energy balance in dairy cattle

Open ArchivePublished:November 05, 2020DOI:https://doi.org/10.3168/jds.2020-18503

      ABSTRACT

      Feed efficiency and energy balance are important traits underpinning profitability and environmental sustainability in animal production. They are complex traits, and our understanding of their underlying biology is currently limited. One measure of feed efficiency is residual feed intake (RFI), which is the difference between actual and predicted intake. Variation in RFI among individuals is attributable to the metabolic efficiency of energy utilization. High RFI (H_RFI) animals require more energy per unit of weight gain or milk produced compared with low RFI (L_RFI) animals. Energy balance (EB) is a closely related trait calculated very similarly to RFI. Cellular energy metabolism in mitochondria involves mitochondrial protein (MiP) encoded by both nuclear (NuMiP) and mitochondrial (MtMiP) genomes. We hypothesized that MiP genes are differentially expressed (DE) between H_RFI and L_RFI animal groups and similarly between negative and positive EB groups. Our study aimed to characterize MiP gene expression in white blood cells of H_RFI and L_RFI cows using RNA sequencing to identify genes and biological pathways associated with feed efficiency in dairy cattle. We used the top and bottom 14 cows ranked for RFI and EB out of 109 animals as H_RFI and L_RFI, and positive and negative EB groups, respectively. The gene expression counts across all nuclear and mitochondrial genes for animals in each group were used for differential gene expression analyses, weighted gene correlation network analysis, functional enrichment, and identification of hub genes. Out of 244 DE genes between RFI groups, 38 were MiP genes. The DE genes were enriched for the oxidative phosphorylation (OXPHOS) and ribosome pathways. The DE MiP genes were underexpressed in L_RFI (and negative EB) compared with the H_RFI (and positive EB) groups, suggestive of reduced mitochondrial activity in the L_RFI group. None of the MtMiP genes were among the DE MiP genes between the groups, which suggests a non-rate limiting role of MtMiP genes in feed efficiency and warrants further investigation. The role of MiP, particularly the NuMiP and OXPHOS pathways in RFI, was also supported by our gene correlation network analysis and the hub gene identification. We validated the findings in an independent data set. Overall, our study suggested that differences in feed efficiency in dairy cows may be linked to differences in cellular energy demand. This study broadens our knowledge of the biology of feed efficiency in dairy cattle.

      Key words

      INTRODUCTION

      Feed efficiency is well recognized for its high relevance to farm economics, resource sustainability, and climate change (
      • Herrero M.
      • Havlík P.
      • Valin H.
      • Notenbaert A.
      • Rufino M.C.
      • Thornton P.K.
      • Blümmel M.
      • Weiss F.
      • Grace D.
      • Obersteiner M.
      Biomass use, production, feed efficiencies, and greenhouse gas emissions from global livestock systems.
      ;
      • Herrero M.
      • Thornton P.K.
      Livestock and global change: Emerging issues for sustainable food systems.
      ). Feed constitutes a significant recurring expense on-farm (
      • Chamberlain T.
      Understanding the economics of dairy farming Part 1: Income, costs and profit.
      ;
      • Yilmaz H.
      • Gül M.
      • Akkoyun S.
      • Parlakay O.
      • Bilgili M.E.
      • Vurarak Y.
      • Hizli H.
      • Kilicalp N.
      Economic analysis of dairy cattle farms in east Mediterranean region of Turkey.
      ), and selection for feed efficiency is expected to result in animals that have reduced feed or energy requirement for maintenance (
      • Herd R.M.
      • Archer J.A.
      • Arthur P.F.
      Reducing the cost of beef production through genetic improvement in residual feed intake: Opportunity and challenges to application.
      ). Further, selection for feed efficiency is estimated to reduce methane emission by 15% in a decade (
      • de Haas Y.
      • Pryce J.E.
      • Berry D.P.
      • Veerkamp R.F.
      Genetic and genomic solutions to improve feed efficiency and reduce environmental impact of dairy cattle. In Proc. 10th World Congress of Genetics Applied to Livestock Production, Vancouver, Canada.
      ). Feed efficiency is characterized as a complex trait with challenges and costs in measurement of phenotypes (
      • Arthur J.P.F.
      • Herd R.M.
      Residual feed intake in beef cattle.
      ). The definition of feed efficiency differs in growing and lactating animals and is complicated in the latter, considering the rapid mobilization of body reserves following calving and accumulation until next calving (
      • Roche J.R.
      • Friggens N.C.
      • Kay J.K.
      • Fisher M.W.
      • Stafford K.J.
      • Berry D.P.
      Invited review: Body condition score and its association with dairy cow productivity, health, and welfare.
      ). It can be measured as residual feed intake (RFI), which is the difference between actual and predicted DMI (
      • Koch R.M.
      • Swiger L.A.
      • Chambers D.
      • Gregory K.E.
      Efficiency of feed use in beef cattle.
      ;
      • Berry D.P.
      • Crowley J.J.
      Cell Biology Symposium: Genetics of feed efficiency in dairy and beef cattle.
      ). Energy balance (EB) is another indicator trait for feed efficiency closely related to RFI; EB is the difference between energy intake and energy expenditure for lactation, growth, reproduction, and maintenance. Energy balance also shows the status of body reserve mobilization and is considered a good indicator for health, reproduction, and feed management. The EB and health status of dairy cows are also indicated by their blood metabolic profiles after calving.
      Selection for feed efficiency is feasible in dairy cattle and is already included in the breeding objectives of several countries, such as “feed saved” in Australia (
      • Pryce J.E.
      • Wales W.J.
      • de Haas Y.
      • Veerkamp R.F.
      • Hayes B.J.
      Genomic selection for feed efficiency in dairy cattle.
      ,
      • Pryce J.E.
      • Gonzalez-Recio O.
      • Nieuwhof G.
      • Wales W.J.
      • Coffey M.P.
      • Hayes B.J.
      • Goddard M.E.
      Hot topic: Definition and implementation of a breeding value for feed efficiency in dairy cows.
      ). Moreover, selection for feed efficiency demands a better understanding of the biological mechanisms underlying variation to allow for more accurate animal evaluations for the trait in the future.
      Physiologically, variation in RFI is associated with metabolism, feed intake, digestion, activity, and body heat regulation (
      • Herd R.M.
      • Oddy V.H.
      • Richardson E.C.
      Biological basis for variation in residual feed intake in beef cattle. 1. Review of potential mechanisms.
      ;
      • Herd R.M.
      • Arthur P.F.
      Physiological basis for residual feed intake.
      ). The variation in RFI, according to
      • Korver S.
      • van Eekelen E.A.M.
      • Vos H.
      • Nieuwhof G.J.
      • van Arendonk J.A.M.
      Genetic parameters for feed intake and feed efficiency in growing dairy heifers.
      , reflects differences between animals in utilization of metabolizable energy and high RFI animals are likely to have higher ATP production and utilization in their tissues (
      • Del Bianco Benedeti P.
      • Detmann E.
      • Mantovani H.C.
      • Bonilha S.F.M.
      • Serão N.V.L.
      • Lopes D.R.G.
      • Silva W.
      • Newbold C.J.
      • Duarte M.S.
      Nellore bulls (Bos taurus indicus) with high residual feed intake have increased the expression of genes involved in oxidative phosphorylation in rumen epithelium.
      ). Further, efficiency of energy utilization varies between individuals and is closely associated with genetic variation (
      • Nkrumah J.D.
      • Okine E.K.
      • Mathison G.W.
      • Schmid K.
      • Li C.
      • Basarab J.A.
      • Price M.A.
      • Wang Z.
      • Moore S.S.
      Relationships of feedlot feed efficiency, performance, and feeding behavior with metabolic rate, methane production, and energy partitioning in beef cattle.
      ) as well as cattle types (
      • Pfuhl R.
      • Bellmann O.
      • Kühn C.
      • Teuscher F.
      • Ender K.
      • Wegner J.
      Beef versus dairy cattle: A comparison of feed conversion, carcass composition, and meat quality.
      ). Energy balance is a trait positively and closely related to RFI (
      • Hurley A.M.
      • López-Villalobos N.
      • McParland S.
      • Lewis E.
      • Kennedy E.
      • O'Donovan M.
      • Burke J.L.
      • Berry D.P.
      Genetics of alternative definitions of feed efficiency in grazing lactating dairy cows.
      ). This is a complication specific to lactating animals and arises from the interplay of nutrient mobilization and replenishment from body reserves during lactation (
      • Pryce J.E.
      • Wales W.J.
      • de Haas Y.
      • Veerkamp R.F.
      • Hayes B.J.
      Genomic selection for feed efficiency in dairy cattle.
      ). Dairy cattle may differ physiologically from beef cattle in utilization, partitioning, and conversion of nutrients. Therefore, the underlying biology of feed efficiency in lactating dairy cattle requires separate consideration.
      The differences in cellular energy metabolism between high (H_RFI) and low (L_RFI) RFI may provide useful insights into feed efficiency. Cellular energy metabolism occurs in mitochondria, and energy is primarily generated through oxidative phosphorylation (OXPHOS) in complexes of proteins called the electron transport chain (ETC). The proteins of ETC (~130), along with other proteins localized in mitochondria, are referred to as mitochondrial proteins (MiP;
      • Fox T.D.
      Mitochondrial protein synthesis, import, and assembly.
      ) and are vital for mitochondrial function. It has also been proposed that defective proteins in the ETC may lead to suboptimal mitochondrial function and reduce the overall energy efficiency of the animal (
      • Bottje W.G.
      Board Invited Review: Oxidative stress and efficiency: The tightrope act of mitochondria in health and disease.
      ). As such, the variation of metabolic demand across different tissues was reported to be reflected through differential expression of MiP genes (
      • Dorji J.
      • Vander Jagt C.J.
      • Garner J.B.
      • Marett L.C.
      • Mason B.A.
      • Reich C.M.
      • Xiang R.
      • Clark E.L.
      • Cocks B.G.
      • Chamberlain A.J.
      • MacLeod I.M.
      • Daetwyler H.D.
      Expression of mitochondrial protein genes encoded by nuclear and mitochondrial genomes correlate with energy metabolism in dairy cattle.
      ), where a gene is considered as differentially expressed (DE) if the expression in a tissue differed significantly from mean expression across all other tissues. Differential gene expression has been investigated between high and low milk yielders in dairy cattle (
      • Yang J.
      • Jiang J.
      • Liu X.
      • Wang H.
      • Guo G.
      • Zhang Q.
      • Jiang L.
      Differential expression of genes in milk of dairy cattle during lactation.
      ) and between chickens and livestock divergent for feed efficiency (
      • Lassiter K.
      • Ojano-Dirain C.
      • Iqbal M.
      • Pumford N.R.
      • Tinsley N.
      • Lay J.
      • Liyanage R.
      • Wing T.
      • Cooper M.
      • Bottje W.
      Differential expression of mitochondrial and extramitochondrial proteins in lymphocytes of male broilers with low and high feed efficiency.
      ;
      • Kong R.S.G.
      • Liang G.
      • Chen Y.
      • Stothard P.
      • Guan L.L.
      Transcriptome profiling of the rumen epithelium of beef cattle differing in residual feed intake.
      ). Given the importance of energy metabolism in feed efficiency, it is plausible that MiP genes are DE between H_RFI and L_RFI groups, and that a similar pattern would be observed in positive and negative EB groups. Previous studies on the role of MiP gene expression in feed efficiency are restricted to comparatively few selected genes (
      • Lassiter K.
      • Ojano-Dirain C.
      • Iqbal M.
      • Pumford N.R.
      • Tinsley N.
      • Lay J.
      • Liyanage R.
      • Wing T.
      • Cooper M.
      • Bottje W.
      Differential expression of mitochondrial and extramitochondrial proteins in lymphocytes of male broilers with low and high feed efficiency.
      ;
      • Kelly A.K.
      • Waters S.M.
      • McGee M.
      • Fonseca R.G.
      • Carberry C.
      • Kenny D.A.
      mRNA expression of genes regulating oxidative phosphorylation in the muscle of beef cattle divergently ranked on residual feed intake.
      ) and are limited in their scope of inferring meaningful biological pathways.
      More recently, high-throughput RNA sequencing (RNAseq) of entire transcriptomes allows for the identification of biological pathways and key genes behind complex traits and diseases (
      • Salleh S.M.
      • Mazzoni G.
      • Løvendahl P.
      • Kadarmideen H.N.
      Gene co-expression networks from RNA sequencing of dairy cattle identifies genes and pathways affecting feed efficiency.
      ; van Dam et al., 2018;
      • Wang Z.
      • Zhu J.
      • Chen F.
      • Ma L.
      Weighted gene coexpression network analysis identifies key genes and pathways associated with idiopathic pulmonary fibrosis.
      ) using gene correlation network analysis. In recent years, RNAseq has increasingly been used to study feed efficiency in cattle based on transcriptomes of liver (
      • Alexandre P.A.
      • Kogelman L.J.A.
      • Santana M.H.A.
      • Passarelli D.
      • Pulz L.H.
      • Fantinato-Neto P.
      • Silva S.L.
      • Leme P.R.
      • Strefezzi R.F.
      • Coutinho L.L.
      • Ferraz J.B.S.
      • Eler J.P.
      • Kadarmideen H.N.
      • Fukumasu H.
      Liver transcriptomic networks reveal main biological processes associated with feed efficiency in beef cattle.
      ;
      • Salleh M.S.
      • Mazzoni G.
      • Höglund J.K.
      • Olijhoek D.W.
      • Lund P.
      • Løvendahl P.
      • Kadarmideen H.N.
      RNA-Seq transcriptomics and pathway analyses reveal potential regulatory genes and molecular mechanisms in high- and low-residual feed intake in Nordic dairy cattle.
      ), muscle (
      • Zhou N.
      • Lee W.R.
      • Abasht B.
      Messenger RNA sequencing and pathway analysis provide novel insights into the biological basis of chickens' feed efficiency.
      ;
      • Horodyska J.
      • Wimmers K.
      • Reyer H.
      • Trakooljul N.
      • Mullen A.M.
      • Lawlor P.G.
      • Hamill R.M.
      RNA-seq of muscle from pigs divergent in feed efficiency and product quality identifies differences in immune response, growth, and macronutrient and connective tissue metabolism.
      ), rumen epithelium (
      • Kong R.S.G.
      • Liang G.
      • Chen Y.
      • Stothard P.
      • Guan L.L.
      Transcriptome profiling of the rumen epithelium of beef cattle differing in residual feed intake.
      ;
      • Del Bianco Benedeti P.
      • Detmann E.
      • Mantovani H.C.
      • Bonilha S.F.M.
      • Serão N.V.L.
      • Lopes D.R.G.
      • Silva W.
      • Newbold C.J.
      • Duarte M.S.
      Nellore bulls (Bos taurus indicus) with high residual feed intake have increased the expression of genes involved in oxidative phosphorylation in rumen epithelium.
      ;
      • Elolimy A.A.
      • Abdelmegeid M.K.
      • McCann J.C.
      • Shike D.W.
      • Loor J.J.
      Residual feed intake in beef cattle and its association with carcass traits, ruminal solid-fraction bacteria, and epithelium gene expression.
      ), and blood (
      • Khansefid M.
      • Millen C.A.
      • Chen Y.
      • Pryce J.E.
      • Chamberlain A.J.
      • Vander Jagt C.J.
      • Gondro C.
      • Goddard M.E.
      Gene expression analysis of blood, liver, and muscle in cattle divergently selected for high and low residual feed intake 1.
      ). Blood is not a common sample for feed efficiency studies compared with other tissues, presumably due to its lower importance as a component of meat, unlike muscle and liver in beef cattle, as well as the lower mitochondrial activity in blood. On the other hand, lymphocytes, constituting about 30% of white blood cells (WBC), are metabolically oxidative (
      • Kramer P.A.
      • Ravi S.
      • Chacko B.
      • Johnson M.S.
      • Darley-Usmar V.M.
      A review of the mitochondrial and glycolytic metabolism in human platelets and leukocytes: Implications for their use as bioenergetic biomarkers.
      ) and are easily accessible compared with other tissues. Thus, differences in metabolic activity and energy efficiency between efficient and non-efficient groups are likely to be reflected in the WBC transcriptome.
      The main objectives of this study were to profile MiP gene expression in WBC of high and low feed efficiency dairy cattle, to identify key biological pathways and genes underlying feed efficiency in dairy cattle. We hypothesized that MiP genes in the WBC are DE between H_RFI and L_RFI groups as well as positive EB and negative EB groups, and that their biological pathways are related to energy metabolism.

      MATERIALS AND METHODS

      Animals, Residual Feed Intake Calculation, and Ranking

      The Agricultural Research and Extension Animal Ethics Committee (Department of Jobs, Precincts and Regions, Attwood, Victoria, Australia) approved the protocols for this experiment (Application No. 2013-14). From 2013 to 2015, 360 animals, in batches of 40 and 3 batches each year, were run through auto-feeders for 5 wk at the Agriculture Victoria Research farm at Ellinbank, Victoria, Australia. The first week was an adaptation phase, and the remaining 4 wk were for the experiment. Each cow had the DM weight of feed offered and refused (kg), DMI (kg/d), milk yield (kg/d), fat and protein yields (kg/d), BCS, and BW (kg) recorded and used in prediction of RFI using the same model in
      • Pryce J.E.
      • Gonzalez-Recio O.
      • Thornhill J.B.
      • Marett L.C.
      • Wales W.J.
      • Coffey M.P.
      • de Haas Y.
      • Veerkamp R.F.
      • Hayes B.J.
      Short communication: Validation of genomic breeding value predictions for feed intake and feed efficiency traits.
      , as follows:
      DMI = µ + DIM + Batch + Parity + ECM + BW + BCS + RFI,


      where μ = the overall mean DMI/intercept; DMI = average DMI during 28-d experimental period; DIM = DIM at the start of each experiment; Batch was the fixed effect of the experiment (n = 9); Parity is the parity group; ECM = mean ECM yield per day (kg) during the 28-d experimental period; BW was average daily BW, measured using walkover scales (automatic weigh system, model AWS100; DeLaval, Tumba, Sweden). On average, 20 BW measurements were recorded per cow over the 28-d experimental period. Body condition score was assessed twice weekly by 4 assessors, using the 8-point scale described by
      • Earle D.
      A guide to scoring dairy cow condition.
      . A mean BCS of the 4 assessors was recorded per week and averaged over the experimental period. Body condition score (BCS) was included as a covariate in the model to correct for an approximation of body fat content, and RFI is the residual term from the equation.
      Out of 352 animals with RFI phenotypes, we considered 109 that had RNAseq data in this study. Selected animals were ranked based on RFI values, and the top and bottom 14 animals were grouped as H_RFI and L_RFI, respectively, and used as 2 divergent feed efficiency groups for further analysis. The differences in RFI estimates between H_RFI and L_RFI were tested using independent sample t-tests for significance. A principal component analysis was performed to visualize and compare global gene expression of H_RFI and L_RFI group using the ggplot2 (
      • Wickham H.
      ggplot2: Elegant Graphics for Data Analysis..
      ) and ggfortify (
      • Tang Y.
      • Horikoshi M.
      • Li W.
      ggfortify: Unified interface to visualize statistical result of popular R packages.
      ) packages in R version 3.5.2 (
      • R Core Team
      R: A Language and Environment for Statistical Computing..
      ). Similarly, we also predicted EB following
      • de Vries M.J.
      • van der Beek S.
      • Kaal-Lansbergen L.M.
      • Ouweltjes W.
      • Wilmink J.B.
      Modeling of energy balance in early lactation and the effect of energy deficits in early lactation on first detected estrus postpartum in dairy cows.
      , thus:
      DMI = µ + DIM + Batch + Parity + ECM + BW + EB.


      Abbreviations are as described in the equation for RFI estimation, except for exclusion of BCS and inclusion of EB, which is the residual term in this equation.
      The group means of the variables used in the prediction of RFI and EB are provided in Supplemental Table S1 (https://doi.org/10.3168/jds.2020-18503). The gene expression counts of the 14 top- and bottom-ranking animals for RFI and EB groups were analyzed for differential gene expression.
      Further, a subset of animals in our study group had mid-infrared spectral data collected from milk samples, so this was used to predict the profile of 3 metabolites that are good indicators of energy status, BHB, nonesterified fatty acid (NEFA), and BUN, following
      • Luke T.D.W.
      • Rochfort S.
      • Wales W.J.
      • Bonfatti V.
      • Marett L.
      • Pryce J.E.
      Metabolic profiling of early-lactation dairy cows using milk mid-infrared spectra.
      ). There were 9 animals with mid-infrared data in the L_RFI group and 5 in the H_RFI group, and similarly 8 in the negative EB group and 3 in the positive EB group.

      RNAseq Read Alignment, Processing, and Expression Counts

      The blood sampling and processing, RNA extraction, sequencing, and quality controls of the RNAseq reads used in this study have been described by
      • Xiang R.
      • Hayes B.J.
      • Vander Jagt C.J.
      • MacLeod I.M.
      • Khansefid M.
      • Bowman P.J.
      • Yuan Z.
      • Prowse-Wilkins C.P.
      • Reich C.M.
      • Mason B.A.
      • Garner J.B.
      • Marett L.C.
      • Chen Y.
      • Bolormaa S.
      • Daetwyler H.D.
      • Chamberlain A.J.
      • Goddard M.E.
      Genome variants associated with RNA splicing variations in bovine are extensively shared between tissues.
      . In this study, we used the trimmed RNAseq reads that passed quality control (in fastq format) of H_RFI and L_RFI animals. The trimmed high-quality pair-end reads of each library were aligned to Ensembl bovine genome UMD3.1 using STAR version 2.5.3ab (
      • Dobin A.
      • Davis C.A.
      • Schlesinger F.
      • Drenkow J.
      • Zaleski C.
      • Jha S.
      • Batut P.
      • Chaisson M.
      • Gingeras T.R.
      STAR: Ultrafast universal RNA-seq aligner.
      ) and checked for alignment quality using Qualimap2 (
      • Okonechnikov K.
      • Conesa A.
      • García-Alcalde F.
      Qualimap 2: Advanced multi-sample quality control for high-throughput sequencing data.
      ). Only the uniquely mapped reads were used for downstream analyses. A gene count matrix for every sample was generated using the R package featureCounts (
      • Liao Y.
      • Smyth G.K.
      • Shi W.
      featureCounts: An efficient general purpose program for assigning sequence reads to genomic features.
      ).

      Differential Gene Expression

      The R package edgeR (
      • Robinson M.D.
      • McCarthy D.J.
      • Smyth G.K.
      edgeR: A bioconductor package for differential expression analysis of digital gene expression data.
      ) was employed for differential gene expression analysis. Raw gene counts were filtered using edgeR function filterbyExp to remove genes that were not adequately expressed, resulting in 13,469 genes passing this filter. The gene expression counts were normalized and converted to counts per million to correct for variation due to sequencing depth. The edgeR function exactTest was used to identify DE genes between the 2 RFI groups at a threshold of log2 fold change (LFC) > |0.6| and false discovery rate of 5%. The sign of LFC was used to deduce the direction of the gene expression as either overexpressed (positive) or underexpressed (negative). A total of 1,054 MiP genes for cattle that were derived from the human Mitocarta 2.0 (
      • Pagliarini D.J.
      • Calvo S.E.
      • Chang B.
      • Sheth S.A.
      • Vafai S.B.
      • Ong S.E.
      • Walford G.A.
      • Sugiana C.
      • Boneh A.
      • Chen W.K.
      • Hill D.E.
      • Vidal M.
      • Evans J.G.
      • Thorburn D.R.
      • Carr S.A.
      • Mootha V.K.
      A mitochondrial protein compendium elucidates complex I disease biology.
      ;
      • Calvo S.E.
      • Clauser K.R.
      • Mootha V.K.
      MitoCarta2.0: An updated inventory of mammalian mitochondrial proteins.
      ) were considered for the analysis. The DE genes list was further analyzed for functional enrichment and direction of regulation, and gene composition for MiP genes.

      Co-Expression Network Analysis

      To build gene correlation networks and to identify a group of highly correlated gene modules associated with RFI, we used the weighted gene correlation network analysis R package, WGCNA (
      • Langfelder P.
      • Horvath S.
      WGCNA: An R package for weighted correlation network analysis.
      ). The gene expression counts in samples were clustered using hclust (average method) in WGCNA to identify obvious outliers. Outlier samples were removed, and a heatmap of RFI values along with a dendrogram showing a hierarchy of cluster-based gene expression was generated. The construction of an unsigned weighted gene correlation network and module detection were performed using the automatic, one-step function network construction and module detection function blockwiseModule. Module eigengene (ME) represented the average expression level of genes in a module. The module membership is the degree of correlation between genes and module. Gene significance (GS) for each gene was calculated as the correlation between gene expression counts and RFI. The correlation between ME and RFI was calculated to identify modules associated with the trait (r > |0.4|, P-value < 0.05).
      From the modules associated with the RFI, the genes meeting the criteria of module membership > |0.8|, GS > |0.2| and P GS < 0.01 (
      • Wang W.
      • Jiang W.
      • Hou L.
      • Duan H.
      • Wu Y.
      • Xu C.
      • Tan Q.
      • Li S.
      • Zhang D.
      Weighted gene co-expression network analysis of expression data of monozygotic twins identifies specific modules and hub genes related to BMI.
      ;
      • Salleh S.M.
      • Mazzoni G.
      • Løvendahl P.
      • Kadarmideen H.N.
      Gene co-expression networks from RNA sequencing of dairy cattle identifies genes and pathways affecting feed efficiency.
      ;
      • Liu Y.
      • Gu H.-Y.
      • Zhu J.
      • Niu Y.-M.
      • Zhang C.
      • Guo G.-L.
      Identification of hub genes and key pathways associated with bipolar disorder based on weighted gene co-expression network analysis.
      ) were identified as putative hub genes. The putative hub genes were searched for on the interacting genes (STRING) database (
      • Szklarczyk D.
      • Gable A.L.
      • Lyon D.
      • Junge A.
      • Wyder S.
      • Huerta-Cepas J.
      • Simonovic M.
      • Doncheva N.T.
      • Morris J.H.
      • Bork P.
      • Jensen L.J.
      • Mering C.
      STRING v11: Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets.
      ), available online at https://string-db.org/, to construct a protein-protein interaction (PPI) network and screen hub nodes. A putative hub gene with high connectivity in a PPI network plays a critical role in the pathways associated with the module.

      Functional Enrichment Analysis

      The DE genes and genes in correlation network modules significantly correlated with RFI (r > |0.4|; P < 0.05) were subjected to functional enrichment of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (
      • Kanehisa M.
      • Goto S.
      KEGG: Kyoto Encyclopedia of Genes and Genomes.
      ) inbuilt in the online software Database for Annotation, Visualization, and Integrated Discovery (DAVID;
      • Huang D.W.
      • Sherman B.T.
      • Lempicki R.A.
      Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis of large gene lists.
      ,
      • Huang D.W.
      • Sherman B.T.
      • Lempicki R.A.
      Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.
      ), available at https://david.ncifcrf.gov/summary.jsp. The cut-off for considering the pathway as enriched was set at adjusted P < 1.0 × 10−5. The Benjamini-Hochberg correction method in DAVID was used for the calculation of adjusted P-values.

      Validation

      The RFI estimates and blood gene counts from
      • Khansefid M.
      • Millen C.A.
      • Chen Y.
      • Pryce J.E.
      • Chamberlain A.J.
      • Vander Jagt C.J.
      • Gondro C.
      • Goddard M.E.
      Gene expression analysis of blood, liver, and muscle in cattle divergently selected for high and low residual feed intake 1.
      were used to validate the findings. Briefly, the data of
      • Khansefid M.
      • Millen C.A.
      • Chen Y.
      • Pryce J.E.
      • Chamberlain A.J.
      • Vander Jagt C.J.
      • Gondro C.
      • Goddard M.E.
      Gene expression analysis of blood, liver, and muscle in cattle divergently selected for high and low residual feed intake 1.
      were independent of our data set and consisted of RFI phenotypes and gene expression counts from WBC of 19 first-lactation Holstein cows (38 ± 10 DIM). The animals were ranked based on RFI phenotype to identify the top and bottom 8 animals as a H_RFI and L_RFI group, respectively. Their study broadly estimated the correlation of gene counts on the RFI values and identified highly expressed genes that were correlated with RFI using regression models. In the present study, we employed standard differential gene expression and WGCNA approaches for validation of the findings from our primary data set.

      RESULTS

      The mean (SD) RFI values for the H_RFI and L_RFI groups were +1.91 (0.37) and −1.67 (0.28) kg/d, respectively, and differed significantly at P < 0.01. In other words, the L_RFI (i.e., feed-efficient group) consumed on average 3.55 kg/d less feed than the H_RFI group (low feed-efficient group) at the same level of production and maintenance. The principal component analysis plot of overall gene expression showed that principal component 1 explained about 81% of the variation in gene expression between the RFI groups (Figure 1A). Further, animals within an RFI group were more likely to be clustered together (Figure 1B) at the first hierarchical levels, but no clear separation occurred based on overall gene expression. Energy balance was phenotypically highly correlated with RFI (r = +0.99, 109 animals). Except for 2 animals, all the 28 animals in positive and negative EB group were common to the RFI group. The means (SD) of the top 14 positive EB and top 14 negative EB groups were +1.90 (0.38) and −1.64, (0.30) kg/d, respectively, and similar to RFI. The difference between the means of the positive EB and negative EB groups were statistically significant (t = −27.25, P < 0.00001).
      Figure thumbnail gr1
      Figure 1(A) Principal component (PC) analysis plot of global gene expression in white blood cells of 14 high residual feed intake (H_RFI) and 14 low residual feed intake (L_RFI) animals. (B) Hierarchical clustering of gene expression in the 14 H_RFI and 14 L_RFI cows with trait heatmap indicating the intensity of residual feed intake values in a scale of blue (low) to red (high).
      The metabolic profiling of markers for EB (estimated from mid-infrared spectral data) found that the H_RFI group had lower mean (±SD) BHB (0.43 ± 0.06) and NEFA (0.20 ± 0.04) and higher BUN (7.24 ± 0.33) compared with L_RFI (BHB: 0.49 ± 0.08; NEFA: 0.22 ± 0.04; BUN: 6.98 ± 0.52 mmol/L). The positive EB group had lower mean (±SD) BHB (0.45 ± 0.08) and NEFA (0.22 ± 0.04) and higher BUN (7.39 ± 0.20) compared with the negative-EB group (BHB: 0.49 ± 0.08; NEFA: 0.21 ± 0.04; BUN: 7.09 ± 0.41 mmol/L). Although these differences were not significantly different (low number of animals available with data), they were consistent in direction, demonstrating that animals with high feed efficiency (L_RFI and negative EB) also had higher levels of these metabolites (Supplemental Table S2, https://doi.org/10.3168/jds.2020-18503). Altogether, our results support that RFI and EB are highly correlated.

      Differentially Expressed Genes and Functional Enrichment

      Overall, 244 genes were DE between H_RFI and L_RFI groups at LFC ≥ |0.6|) and false discovery rate 5% (Supplemental Table S3, https://doi.org/10.3168/jds.2020-18503). Functionally, DE genes were significantly enriched (number of genes, adjusted P) for OXPHOS (15, 2.1 × 10−7) and ribosome (17, 1.9 × 10−9) pathways. However, 47 DE genes were not identified in DAVID and remained unused for functional enrichment analysis. The number of overexpressed and underexpressed genes among DE genes were 64 and 180, respectively. There were 38 DE MiP genes, and all were underexpressed in the L_RFI group, with LFC ranging from −0.6 to −1.6 (Supplemental Table S4, https://doi.org/10.3168/jds.2020-18503). The DE MiP genes were all composed of NuMiP and, surprisingly, contained no MtMiP genes. Differential gene expression analysis was repeated with only Mt genes, as their mean level of expression was generally much higher than that of NuMiP, which again confirmed that no Mt genes were DE. Further, the expression of Mt genes was less variable (CV 180%) compared with genes from the nuclear genome (CV 300%).
      Similarly, 466 genes were DE between positive EB and negative EB groups. Among these, were 53 NuMiP genes, which were all overexpressed in the positive EB group compared with the negative EB animals, except for 2 genes. The DE genes were enriched for OXPHOS (adjusted P 1.6 × 10−10) and ribosome (adjusted P 7.5 × 10−17) pathways. A similar expression profile between RFI and EB was expected, as these traits were highly correlated and contained almost the same animals in both RFI and EB groups. Thus, the subsequent gene network analysis was conducted only for gene expression in the RFI groups.

      WGCNA

      Overall, 9 correlation network modules out of 28 were significantly related to RFI (r > |0.4|, P < 0.05; Figure 2A). Of these modules, only ME2 and ME3 were significantly enriched for KEGG pathways (Figure 2B, C). Altogether, WGCNA and the functional enrichment analysis supported a strong association between OXPHOS pathway and RFI. Further, WGCNA identified additional pathways associated with RFI traits undetected with differential gene expression analysis.
      Figure thumbnail gr2
      Figure 2Gene expression modules correlated with residual feed intake (RFI) based on weighted gene correlation network analysis in 14 high-RFI and 14 low-RFI cows. (A) Module and number of genes in each module, and the relationship of module with RFI (r), with P-values in parentheses. ME = module eigengene. *Indicates module with a significant relationship (r > |0.4| and P < 0.05) with RFI. The direction of the relationship (correlation) is indicated by colors, where red is positive and blue is negative, and the intensity of color represents the strength of correlation. Functional enrichment of Kyoto Encyclopedia of Genes and Genomes pathways of overall genes module ME2 (B) and ME3 (C) for the main data set.
      Hub genes in a module have high network connectivity and are also highly associated with the corresponding traits. The 173 genes in module ME2 meeting the threshold criteria (module membership > |0.8|, GS > |0.2|, and P GS < 0.01) were identified as putative hub genes (Supplemental Table S5, https://doi.org/10.3168/jds.2020-18503). Networks of hub genes were derived based on PPI network at minimum interaction score of 0.90. The resulting networks of hub genes were associated with KEGG pathways: ribosome, OXPHOS, spliceosome, and proteosome (Figure 3). The hub genes that were MiP genes and associated with OXPHOS pathway were ATP5A1, ATP5J, ATP5J2, NDUFAB1, NDUFA8, and NDUFA10, most of which are noncatalytic proteins in the OXPHOS complexes. For ME3, 74 putative hub genes (Supplemental Table S6, https://doi.org/10.3168/jds.2020-18503), refined based on PPI network, had 6 hub genes for OXPHOS, 5 for ribosomes, and 2 for citrate cycle (Figure 4). The refined hub genes for OXPHOS pathway and citrate cycle were MiP genes (NDUFA5, NDUFA12, NDUFB5, ATP5F1, ATP5C1, SDHD, and PDHB), and the dehydrogenase components were noncatalytic.
      Figure thumbnail gr3
      Figure 3A protein–protein interaction network featuring the network of putative hub genes from genes in module eigengene (ME) 1 (173) with confidence of more than 0.90 from the STRING database (
      • Szklarczyk D.
      • Gable A.L.
      • Lyon D.
      • Junge A.
      • Wyder S.
      • Huerta-Cepas J.
      • Simonovic M.
      • Doncheva N.T.
      • Morris J.H.
      • Bork P.
      • Jensen L.J.
      • Mering C.
      STRING v11: Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets.
      ; https://string-db.org/). Colors indicate different Kyoto Encyclopedia of Genes and Genomes pathways: proteosome (red), spliceosome (dark blue), ribosome (green), and oxidative phosphorylation (yellow).
      Figure thumbnail gr4
      Figure 4A protein–protein interaction network featuring putative hub genes based on genes in module eigengene (ME) 3 (74) with confidence of more than 0.90 in the STRING database (
      • Szklarczyk D.
      • Gable A.L.
      • Lyon D.
      • Junge A.
      • Wyder S.
      • Huerta-Cepas J.
      • Simonovic M.
      • Doncheva N.T.
      • Morris J.H.
      • Bork P.
      • Jensen L.J.
      • Mering C.
      STRING v11: Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets.
      ; https://string-db.org/). Colors indicate different Kyoto Encyclopedia of Genes and Genomes pathways: ribosome (red), oxidative phosphorylation (green), and citrate cycle (blue).

      Validation

      In the independent set of validation cows, the mean (SD) RFI values of H_RFI and L_RFI groups were +1.20 (0.95) and −0.97 (0.71) kg/d, respectively; consequently, as expected, the means of the RFI groups were significantly different [t(14) = 5.20, P < 0.001]. Global gene expression was not distinct but showed greater variability in the H_RFI group compared with the L_RFI group (Supplemental Figure S1, https://doi.org/10.3168/jds.2020-18503). We found 1,695 DE genes between L_RFI and H_RFI groups at LFC ≥ |0.6| and false discovery rate < 0.05, including 151 MiP genes (Supplemental Table S7, https://doi.org/10.3168/jds.2020-18503). The DE genes were significantly enriched for OXPHOS (adjusted P 2.6 × 10−7) and ribosome (adjusted P 5.9 × 10−19) KEGG pathways. All 151 DE MiP genes were NuMiP (18 overexpressed and 133 underexpressed in the L_RFI group compared with the H_RFI group), and none were MtMiP. The 133 underexpressed MiP genes were over-represented (adjusted P) in the following KEGG pathways: OXPHOS (3.4 × 10−28), metabolic pathways (8.6 × 10−13), and ribosome (2.0 × 10−13). We found 31 DE MiP genes in common between the validation and the main data sets (Supplemental Table S8, https://doi.org/10.3168/jds.2020-18503), and they showed concordant underexpression. Altogether, the differential expression of MiP genes between 2 RFI groups, the enrichment of OXPHOS pathway, underexpression of MiP genes in L_RFI group, and lack of representation of MtMiP genes in the validation data set agreed with the findings of the main data set.
      Five co-expression modules were significantly related to the RFI trait in WGCNA analysis: ME1, ME2, ME4, ME5, and ME8. Modules ME1 and ME5 were positively correlated with RFI, whereas ME2, ME4, and ME8 were negatively correlated with RFI (Supplemental Figure S2a, https://doi.org/10.3168/jds.2020-18503). Among the modules, only ME1 was significantly enriched for biological pathways. It was enriched for the ribosome, OXPHOS, spliceosome, metabolic, and proteasome pathways (Supplemental Figure S2b, https://doi.org/10.3168/jds.2020-18503). The top 402 genes were ranked in order of module connectivity and GS to RFI, as well as based on the PPI network-identified hub genes related to ribosome, OXPHOS, and peroxisome (Supplemental Figure S3, https://doi.org/10.3168/jds.2020-18503). The networked hub genes in the OXPHOS pathway that were also MiP genes were ATP5B, COX4I1, UQCRC1, UQCRFS1, SDHA, NDUFS2, NDUFA2, NUDFA8, NDUFA9, and NDUFB10. However, only NDUFA8 overlapped with hub genes in the main data set.

      DISCUSSION

      Residual feed intake and EB were highly correlated and were used to cross-check gene expression patterns in high and low groups. Only 2% of genes were DE from a total of about 13,469 genes that passed the RNAseq quality filters, thus representing only a small fraction of genes. This likely explains the lack of clear clustering when considering expression levels of all genes in the 2 RFI groups in principal component analysis plot and hierarchical cluster (Figure 1A). Within the DE genes, we found that MiP genes were consistently underexpressed in the low-RFI groups, and there was considerable overlap of these DE MiP genes in an independent validation set of cows. Differentially expressed genes were over-represented for OXPHOS and ribosome pathways in both the main and validation data sets. This strongly suggests a role of the OXPHOS pathway in feed efficiency, and this was also evident from the weighted gene co-expression network analysis. However, none of the genes from the mitochondrial genome (MtMiP) were DE or were found in the co-expression module associated with RFI.

      Energy Balance and RFI

      The model for estimation of RFI excluding the changes in BW or BCS gives an estimate of EB (
      • Veerkamp R.F.
      Feed intake and energy balance in lactating animals. In 7th World Congr. Genet. Appl. Livest. Prod. Vol. CD-ROM Communication No. 10-01, Montpellier, France.
      ). Furthermore, it is also very difficult to accurately measure and distinguish between EB and feed efficiency, and no consensus currently exists on the most appropriate mathematical models (e.g.,
      • de Vries M.J.
      • van der Beek S.
      • Kaal-Lansbergen L.M.
      • Ouweltjes W.
      • Wilmink J.B.
      Modeling of energy balance in early lactation and the effect of energy deficits in early lactation on first detected estrus postpartum in dairy cows.
      ;
      • Veerkamp R.F.
      Feed intake and energy balance in lactating animals. In 7th World Congr. Genet. Appl. Livest. Prod. Vol. CD-ROM Communication No. 10-01, Montpellier, France.
      ;
      • Hurley A.M.
      • López-Villalobos N.
      • McParland S.
      • Kennedy E.
      • Lewis E.
      • O'Donovan M.
      • Burke J.L.
      • Berry D.P.
      Inter-relationships among alternative definitions of feed efficiency in grazing lactating dairy cows.
      ). In addition, differences in activity levels between cows have been shown to cause changes in DMI (
      • Olijhoek D.W.
      • Difford G.F.
      • Lund P.
      • Løvendahl P.
      Phenotypic modeling of residual feed intake using physical activity and methane production as energy sinks.
      ). These measures of feed efficiency are highly correlated in mid-lactation cows (0.96), when changes in BW were almost zero (
      • Hurley A.M.
      • López-Villalobos N.
      • McParland S.
      • Kennedy E.
      • Lewis E.
      • O'Donovan M.
      • Burke J.L.
      • Berry D.P.
      Inter-relationships among alternative definitions of feed efficiency in grazing lactating dairy cows.
      ). Changes in BCS primarily occur in early lactation, due to mobilization of body reserves, whereas the 352 animals used in our study estimate RFI were on average approximately 100 DIM and entering mid-lactation, where little change is expected in BCS or BW (
      • Spurlock D.M.
      • Dekkers J.C.M.
      • Fernando R.
      • Koltes D.A.
      • Wolc A.
      Genetic parameters for energy balance, feed efficiency, and related traits in Holstein cattle.
      ;
      • Hurley A.M.
      • López-Villalobos N.
      • McParland S.
      • Kennedy E.
      • Lewis E.
      • O'Donovan M.
      • Burke J.L.
      • Berry D.P.
      Inter-relationships among alternative definitions of feed efficiency in grazing lactating dairy cows.
      ). Therefore, not surprisingly, our study showed a phenotypic correlation close to 1 between RFI and EB. We found similar gene expression patterns and enrichment of pathways in RFI and EB, because almost the same animals were allocated across these 2 study sets. Interestingly, the predicted serum BHB and NEFA levels (metabolite indicators for negative EB) were higher in the more feed-efficient groups (L_RFI and negative EB) compared with the less feed-efficient group (H_RFI and positive EB) but much lower than critical limits (>0.6 mmol/L and 10 mg/dL;
      • Ospina P.A.
      • Nydam D.V.
      • Stokol T.
      • Overton T.R.
      Evaluation of nonesterified fatty acids and β-hydroxybutyrate in transition dairy cattle in the northeastern United States: Critical thresholds for prediction of clinical diseases.
      ) in both groups. It is therefore important to acknowledge that our observations of differential gene expression may be a reflection of genetics underpinning feed efficiency or EB. It is important in the context of genetic selection programs to further explore approaches to better distinguish between RFI and EB traits in dairy cattle (e.g., use of metabolite profiles), because if animals are inadvertently selected for poorer EB, this could affect health and fertility.

      OXPHOS Pathway in Feed Efficiency

      No previous studies have examined MiP genes comprehensively in relation to feed efficiency, although a few genes belonging to MiP or the OXPHOS pathway have been reported in several studies on feed efficiency of meat animal species, mainly beef cattle, pigs, and chicken, involving liver, muscle, and rumen epithelial transcriptomes (
      • Kong R.S.G.
      • Liang G.
      • Chen Y.
      • Stothard P.
      • Guan L.L.
      Transcriptome profiling of the rumen epithelium of beef cattle differing in residual feed intake.
      ;
      • Khansefid M.
      • Millen C.A.
      • Chen Y.
      • Pryce J.E.
      • Chamberlain A.J.
      • Vander Jagt C.J.
      • Gondro C.
      • Goddard M.E.
      Gene expression analysis of blood, liver, and muscle in cattle divergently selected for high and low residual feed intake 1.
      ;
      • Del Bianco Benedeti P.
      • Detmann E.
      • Mantovani H.C.
      • Bonilha S.F.M.
      • Serão N.V.L.
      • Lopes D.R.G.
      • Silva W.
      • Newbold C.J.
      • Duarte M.S.
      Nellore bulls (Bos taurus indicus) with high residual feed intake have increased the expression of genes involved in oxidative phosphorylation in rumen epithelium.
      ), as well as GWA studies (
      • Khansefid M.
      • Millen C.A.
      • Chen Y.
      • Pryce J.E.
      • Chamberlain A.J.
      • Vander Jagt C.J.
      • Gondro C.
      • Goddard M.E.
      Gene expression analysis of blood, liver, and muscle in cattle divergently selected for high and low residual feed intake 1.
      ;
      • Li B.
      • Fang L.
      • Null D.J.
      • Hutchison J.L.
      • Connor E.E.
      • VanRaden P.M.
      • VandeHaar M.J.
      • Tempelman R.J.
      • Weigel K.A.
      • Cole J.B.
      High-density genome-wide association study for residual feed intake in Holstein dairy cattle.
      ). The proteomic analysis of feed efficiency indicated significant differences in the abundance of MiP in the muscle of pigs (
      • Fu L.
      • Xu Y.
      • Hou Y.
      • Qi X.
      • Zhou L.
      • Liu H.
      • Luan Y.
      • Jing L.
      • Miao Y.
      • Zhao S.
      • Liu H.
      • Li X.
      Proteomic analysis indicates that mitochondrial energy metabolism in skeletal muscle tissue is negatively correlated with feed efficiency in pigs.
      ), liver of beef cattle (
      • Baldassini W.A.
      • Bonilha S.F.M.
      • Branco R.H.
      • Vieira J.C.S.
      • Padilha P.M.
      • Lanna D.P.D.
      Proteomic investigation of liver from beef cattle (Bos indicus) divergently ranked on residual feed intake.
      ), and breast muscle in chicken (
      • Kong B.W.
      • Lassiter K.
      • Piekarski-Welsher A.
      • Dridi S.
      • Reverter A.
      • Hudson N.J.
      • Bottje W.G.
      Proteomics of breast muscle tissue associated with the phenotypic expression of feed efficiency within a pedigree male broiler line: I. Highlight on mitochondria.
      ). The gene expression and protein levels of mitochondrial energy metabolism have been linked to feed efficiency in pigs (
      • Vincent A.
      • Louveau I.
      • Gondret F.
      • Trefeu C.
      • Gilbert H.
      • Lefaucheur L.
      Divergent selection for residual feed intake affects the transcriptomic and proteomic profiles of pig skeletal muscle.
      ).
      For dairy cattle, our study is perhaps among the foremost in the investigation of feed efficiency using blood transcriptomes and demonstration of the association of MiP and OXPHOS pathways in feed efficiency. Although blood has low mitochondrial activity compared with muscle and liver tissues (
      • Dorji J.
      • Vander Jagt C.J.
      • Garner J.B.
      • Marett L.C.
      • Mason B.A.
      • Reich C.M.
      • Xiang R.
      • Clark E.L.
      • Cocks B.G.
      • Chamberlain A.J.
      • MacLeod I.M.
      • Daetwyler H.D.
      Expression of mitochondrial protein genes encoded by nuclear and mitochondrial genomes correlate with energy metabolism in dairy cattle.
      ), blood transcriptome in the present study showed distinct MiP gene expression differences between high and low feed efficiency groups and enrichment of the OXPHOS pathways that have not been reported in earlier studies of dairy cattle. These findings altogether suggest that blood is potentially a good tissue for studying gene expression associated with changes in feed efficiency. More importantly, the role of MiP and the OXPHOS pathways in feed efficiency across animal and tissue types is growing more evident.
      Another pathway highly associated with feed efficiency in this study was the ribosome pathway. The enrichment of the ribosome pathway was previously reported in transcriptomes of rumen epithelium (
      • Kong R.S.G.
      • Liang G.
      • Chen Y.
      • Stothard P.
      • Guan L.L.
      Transcriptome profiling of the rumen epithelium of beef cattle differing in residual feed intake.
      ), muscle, and blood (
      • Khansefid M.
      • Millen C.A.
      • Chen Y.
      • Pryce J.E.
      • Chamberlain A.J.
      • Vander Jagt C.J.
      • Gondro C.
      • Goddard M.E.
      Gene expression analysis of blood, liver, and muscle in cattle divergently selected for high and low residual feed intake 1.
      ). Co-regulation of the ribosome with the OXPHOS pathway is expected, considering that protein synthesis is an energy-demanding process, particularly in peptide bonding, where one mole of a polypeptide bond during protein synthesis requires about 4 ATP (
      • MacRae J.C.
      • Lobley G.E.
      Some factors which influence thermal energy losses during the metabolism of ruminants.
      ).

      Direction of Regulation of Mitochondrial Protein Genes

      The MiP genes were underexpressed in the feed-efficient animals. The underexpression of MiP genes in blood WBC and their relation to feed efficiency may in part be explained by the lowered activity of the OXPHOS pathway, which translates into decreased energy production and probably less heat loss (
      • Nkrumah J.D.
      • Okine E.K.
      • Mathison G.W.
      • Schmid K.
      • Li C.
      • Basarab J.A.
      • Price M.A.
      • Wang Z.
      • Moore S.S.
      Relationships of feedlot feed efficiency, performance, and feeding behavior with metabolic rate, methane production, and energy partitioning in beef cattle.
      ). Energy saved from reduced heat loss is believed to be channeled into increased milk production (
      • Goddard M.E.
      • Grainger C.
      A review of the effects of dairy breed on feed conversion efficiency—An opportunity lost?.
      ). Similarly, when using EB, most DE MiP genes were underexpressed in the negative EB group. Given that the calculation of RFI and EB phenotypes is similar, the similarity in strong underexpression of MiP genes is, to some, expected, with most cows overlapping between the 2 phenotypes. The underexpression of MiP genes involved in energy production and use may arise simply from a shortage of energy from feed or from lower consumption of feed, potentially involving a feedback loop to reduce expression.
      In addition, a lowered protein synthesis or turnover, as indicated by underexpression of genes associated with the ribosome pathway in blood of the feed-efficient group, may save energy. Underexpression of both MiP genes and MiP abundance has been reported in feed-efficient pigs (
      • Vincent A.
      • Louveau I.
      • Gondret F.
      • Trefeu C.
      • Gilbert H.
      • Lefaucheur L.
      Divergent selection for residual feed intake affects the transcriptomic and proteomic profiles of pig skeletal muscle.
      ;
      • Fu L.
      • Xu Y.
      • Hou Y.
      • Qi X.
      • Zhou L.
      • Liu H.
      • Luan Y.
      • Jing L.
      • Miao Y.
      • Zhao S.
      • Liu H.
      • Li X.
      Proteomic analysis indicates that mitochondrial energy metabolism in skeletal muscle tissue is negatively correlated with feed efficiency in pigs.
      ). Similarly, less feed-efficient Nellore bulls reportedly had increased expression of genes related to OXPHOS in the rumen epithelium (
      • Del Bianco Benedeti P.
      • Detmann E.
      • Mantovani H.C.
      • Bonilha S.F.M.
      • Serão N.V.L.
      • Lopes D.R.G.
      • Silva W.
      • Newbold C.J.
      • Duarte M.S.
      Nellore bulls (Bos taurus indicus) with high residual feed intake have increased the expression of genes involved in oxidative phosphorylation in rumen epithelium.
      ). Conversely, the transcriptome profile of rumen epithelium had overexpression of OXPHOS genes in feed-efficient beef steers (
      • Kong R.S.G.
      • Liang G.
      • Chen Y.
      • Stothard P.
      • Guan L.L.
      Transcriptome profiling of the rumen epithelium of beef cattle differing in residual feed intake.
      ). The authors justified their findings by suggesting that feed efficiency can be enhanced through better absorption of nutrients from the rumen wall, facilitated by increased morphogenesis and protein turnover in the rumen epithelium. Similarly, the proteome of breast muscle had overexpression of OXPHOS proteins in feed-efficient meat chickens (
      • Kong B.W.
      • Lassiter K.
      • Piekarski-Welsher A.
      • Dridi S.
      • Reverter A.
      • Hudson N.J.
      • Bottje W.G.
      Proteomics of breast muscle tissue associated with the phenotypic expression of feed efficiency within a pedigree male broiler line: I. Highlight on mitochondria.
      ). Considering that growth and weight gain are essential aspects of feed efficiency traits in meat animals, and protein synthesis is one of the mechanisms for muscular growth, this justifies the overexpression of genes associated with OXPHOS and ribosome pathways in the muscle of meat animals.
      Overall, it is currently inconclusive whether the specific direction of the regulation of MiP genes is associated with feed efficiency. However, this has tended to depend on animal species and sampling tissue. Overexpression of MiP genes, along with genes for the ribosome pathway, may increase protein synthesis in muscles and affect weight gain in feed-efficient animals, which could be important for meat species. By contrast, underexpression of MiP genes and lower activity of OXPHOS in dairy cattle suggests lower energy production and utilization for maintenance.

      Key Mitochondrial Protein Genes Related to Feed Efficiency

      Mitochondrial protein genes are encoded by genes from both nuclear and mitochondrial genomes, and OXPHOS and energy metabolism are the critical functions of MiP genes and mitochondria. One of the noteworthy findings of this study was that no MtMiP genes were DE between feed efficiency groups, nor were they among the genes in modules highly correlated with RFI. Further, with exceptions from the studies specifically targeting MtMiP, none of the high-throughput RNAseq studies to date has indicated an association of MtMiPs in feed efficiency (
      • Kong R.S.G.
      • Liang G.
      • Chen Y.
      • Stothard P.
      • Guan L.L.
      Transcriptome profiling of the rumen epithelium of beef cattle differing in residual feed intake.
      ;
      • Khansefid M.
      • Millen C.A.
      • Chen Y.
      • Pryce J.E.
      • Chamberlain A.J.
      • Vander Jagt C.J.
      • Gondro C.
      • Goddard M.E.
      Gene expression analysis of blood, liver, and muscle in cattle divergently selected for high and low residual feed intake 1.
      ). This suggests low variability in MtMiP gene expression between the high and low feed efficient groups, and that MiP variation in feed efficiency is mainly attributable to NuMiP gene expression.
      We found 31 out of 38 DE MiP genes of the primary set in the validation data set. Some of these MiP genes and proteins have been associated with feed efficiency in previous studies (Supplemental Table S9, https://doi.org/10.3168/jds.2020-18503). Some common gene transcripts across studies included COX4I1, ATP5D, and UQCRQ, and differentially abundant MiP relating to feed efficiency were COX4I1, NDUFB11, PRDX2, NME3, and UQCRQ. Overall, the 2 most common MiP and MiP genes across studies were COX4I1 and UQCRQ.
      The COX4I1 gene is located on bovine chromosome 18 between 11,799,175 and 11,807,342 bp (UMD3.1) and encodes a protein that is a component of complex IV of the ETC in mitochondria. Complex IV consists of 13 proteins (
      • Yoshikawa S.
      Beef heart cytochrome c oxidase.
      ), of which 10 encoded by nuclear genes are regulatory, and 3 from the mitochondrial genome constitute the catalytic core in mammals (
      • Kadenbach B.
      • Huttemann M.
      The subunit composition and function of mammalian cytochrome c oxidase.
      ). In humans, a mutation in COX4I1 has been associated with short stature, poor weight gain, and increased chromosomal breaks (
      • Abu-Libdeh B.
      • Douiev L.
      • Amro S.
      • Shahrour M.
      • Ta-Shma A.
      • Miller C.
      • Elpeleg O.
      • Saada A.
      Mutation in the COX4I1 gene is associated with short stature, poor weight gain and increased chromosomal breaks, simulating Fanconi anemia.
      ). We checked 1,646 polymorphic variants in the animals in our study to investigate whether any showed segregation patterns according to our RFI groupings within the COX4I1 region, but no significant pattern was observed. We suggest that the role of COX4I1 and other common genes should be further investigated.
      Highly connected hub genes in a module play essential roles in biological pathways and have been suggested for use as a potential indicator of feed efficiency (
      • Salleh S.M.
      • Mazzoni G.
      • Løvendahl P.
      • Kadarmideen H.N.
      Gene co-expression networks from RNA sequencing of dairy cattle identifies genes and pathways affecting feed efficiency.
      ). In this study, we looked at the top 10% of genes in the modules significantly related to feed efficiency in both the main and validation data sets. We found 62 putative hub genes in common, including 7 OXPHOS and 11 ribosome pathway genes, but only 4 in refined hub genes associated with OXPHOS between the main and the validation data sets within our study. We compared our putative hub genes with hub genes from liver in a published study in Danish Holstein and Jersey cattle (
      • Salleh S.M.
      • Mazzoni G.
      • Løvendahl P.
      • Kadarmideen H.N.
      Gene co-expression networks from RNA sequencing of dairy cattle identifies genes and pathways affecting feed efficiency.
      ). Our primary data set shared only 1 putative hub gene (LPXN) with Holsteins, and no genes were common between the validation set and their study. It was also interesting to note that only 1 hub gene (LCK) was found in common between Holsteins and Jersey in their study. As such, considerable variability exists in hub genes as well as DE genes identified across studies among breeds, species, tissues, and data sets for feed efficiency. Therefore, it appears that there is still much to learn regarding the genes and mutations that underpin differences in feed efficiency, and how genomic selection programs for feed efficiency may affect other traits such as health and fertility. Our results suggest that more efficient animals have lower requirements to generate energy (as measured by MiP gene expression) and thus potentially have better metabolic efficiency of energy utilization compared with less efficient animals.

      CONCLUSIONS

      The findings from our study suggest that mitochondrial protein genes in the blood are differentially expressed between high and low feed efficiency groups of lactating dairy cows. However, all differentially expressed mitochondrial protein genes were from the nuclear genome and none from the mitochondrial genome. The oxidative phosphorylation pathway, which is responsible for energy production, and the ribosome function pathway were associated with feed efficiency. Mitochondrial protein genes were underexpressed in the more feed efficient group, which may suggest a lower metabolic turnover.

      ACKNOWLEDGMENTS

      The authors thank the DairyBio program (a joint venture between Agriculture Victoria, Dairy Australia, and the Gardiner Foundation, Melbourne, Victoria, Australia) for funding. J. D. receives a La Trobe University fee remission scholarship (Bundoora, Australia). The authors are grateful to the staff of Agriculture Victoria at the Ellinbank and Bundoora sites, who were involved in sample collection and animal husbandry. The funding agencies did not have influence on the findings of the study, and the authors have no conflict of interest.

      REFERENCES

        • Abu-Libdeh B.
        • Douiev L.
        • Amro S.
        • Shahrour M.
        • Ta-Shma A.
        • Miller C.
        • Elpeleg O.
        • Saada A.
        Mutation in the COX4I1 gene is associated with short stature, poor weight gain and increased chromosomal breaks, simulating Fanconi anemia.
        Eur. J. Hum. Genet. 2017; 25 (28766551): 1142-1146
        • Alexandre P.A.
        • Kogelman L.J.A.
        • Santana M.H.A.
        • Passarelli D.
        • Pulz L.H.
        • Fantinato-Neto P.
        • Silva S.L.
        • Leme P.R.
        • Strefezzi R.F.
        • Coutinho L.L.
        • Ferraz J.B.S.
        • Eler J.P.
        • Kadarmideen H.N.
        • Fukumasu H.
        Liver transcriptomic networks reveal main biological processes associated with feed efficiency in beef cattle.
        BMC Genomics. 2015; 16 (26678995)1073
        • Arthur J.P.F.
        • Herd R.M.
        Residual feed intake in beef cattle.
        Rev. Bras. Zootec. 2008; 37: 269-279
        • Baldassini W.A.
        • Bonilha S.F.M.
        • Branco R.H.
        • Vieira J.C.S.
        • Padilha P.M.
        • Lanna D.P.D.
        Proteomic investigation of liver from beef cattle (Bos indicus) divergently ranked on residual feed intake.
        Mol. Biol. Rep. 2018; 45 (30178216): 2765-2773
        • Berry D.P.
        • Crowley J.J.
        Cell Biology Symposium: Genetics of feed efficiency in dairy and beef cattle.
        J. Anim. Sci. 2013; 91 (23345557): 1594-1613
        • Bottje W.G.
        Board Invited Review: Oxidative stress and efficiency: The tightrope act of mitochondria in health and disease.
        J. Anim. Sci. 2019; 97 (31247079): 3169-3179
        • Calvo S.E.
        • Clauser K.R.
        • Mootha V.K.
        MitoCarta2.0: An updated inventory of mammalian mitochondrial proteins.
        Nucleic Acids Res. 2016; 44 (26450961): D1251-D1257
        • Chamberlain T.
        Understanding the economics of dairy farming Part 1: Income, costs and profit.
        Livestock (Lond). 2012; 17: 30-33
        • de Haas Y.
        • Pryce J.E.
        • Berry D.P.
        • Veerkamp R.F.
        Genetic and genomic solutions to improve feed efficiency and reduce environmental impact of dairy cattle. In Proc. 10th World Congress of Genetics Applied to Livestock Production, Vancouver, Canada.
        2014
        • de Vries M.J.
        • van der Beek S.
        • Kaal-Lansbergen L.M.
        • Ouweltjes W.
        • Wilmink J.B.
        Modeling of energy balance in early lactation and the effect of energy deficits in early lactation on first detected estrus postpartum in dairy cows.
        J. Dairy Sci. 1999; 82 (10509251): 1927-1934
        • Del Bianco Benedeti P.
        • Detmann E.
        • Mantovani H.C.
        • Bonilha S.F.M.
        • Serão N.V.L.
        • Lopes D.R.G.
        • Silva W.
        • Newbold C.J.
        • Duarte M.S.
        Nellore bulls (Bos taurus indicus) with high residual feed intake have increased the expression of genes involved in oxidative phosphorylation in rumen epithelium.
        Anim. Feed Sci. Technol. 2018; 235: 77-86
        • Dobin A.
        • Davis C.A.
        • Schlesinger F.
        • Drenkow J.
        • Zaleski C.
        • Jha S.
        • Batut P.
        • Chaisson M.
        • Gingeras T.R.
        STAR: Ultrafast universal RNA-seq aligner.
        Bioinformatics. 2013; 29 (23104886): 15-21
        • Dorji J.
        • Vander Jagt C.J.
        • Garner J.B.
        • Marett L.C.
        • Mason B.A.
        • Reich C.M.
        • Xiang R.
        • Clark E.L.
        • Cocks B.G.
        • Chamberlain A.J.
        • MacLeod I.M.
        • Daetwyler H.D.
        Expression of mitochondrial protein genes encoded by nuclear and mitochondrial genomes correlate with energy metabolism in dairy cattle.
        BMC Genomics. 2020; 21: 720
        • Earle D.
        A guide to scoring dairy cow condition.
        J. Agric (Victoria). 1976; 74: 228-231
        • Elolimy A.A.
        • Abdelmegeid M.K.
        • McCann J.C.
        • Shike D.W.
        • Loor J.J.
        Residual feed intake in beef cattle and its association with carcass traits, ruminal solid-fraction bacteria, and epithelium gene expression.
        J. Anim. Sci. Biotechnol. 2018; 9 (30258628): 67
        • Fox T.D.
        Mitochondrial protein synthesis, import, and assembly.
        Genetics. 2012; 192 (23212899): 1203-1234
        • Fu L.
        • Xu Y.
        • Hou Y.
        • Qi X.
        • Zhou L.
        • Liu H.
        • Luan Y.
        • Jing L.
        • Miao Y.
        • Zhao S.
        • Liu H.
        • Li X.
        Proteomic analysis indicates that mitochondrial energy metabolism in skeletal muscle tissue is negatively correlated with feed efficiency in pigs.
        Sci. Rep. 2017; 7 (28345649)45291
        • Goddard M.E.
        • Grainger C.
        A review of the effects of dairy breed on feed conversion efficiency—An opportunity lost?.
        Science Access. 2004; 1: 77-80
        • Herd R.M.
        • Archer J.A.
        • Arthur P.F.
        Reducing the cost of beef production through genetic improvement in residual feed intake: Opportunity and challenges to application.
        J. Anim. Sci. 2003; 81: E9-E17
        • Herd R.M.
        • Arthur P.F.
        Physiological basis for residual feed intake.
        J. Anim. Sci. 2009; 87 (19028857): E64-E71
        • Herd R.M.
        • Oddy V.H.
        • Richardson E.C.
        Biological basis for variation in residual feed intake in beef cattle. 1. Review of potential mechanisms.
        Aust. J. Exp. Agric. 2004; 44: 423-430
        • Herrero M.
        • Havlík P.
        • Valin H.
        • Notenbaert A.
        • Rufino M.C.
        • Thornton P.K.
        • Blümmel M.
        • Weiss F.
        • Grace D.
        • Obersteiner M.
        Biomass use, production, feed efficiencies, and greenhouse gas emissions from global livestock systems.
        Proc. Natl. Acad. Sci. USA. 2013; 110 (24344273): 20888-20893
        • Herrero M.
        • Thornton P.K.
        Livestock and global change: Emerging issues for sustainable food systems.
        Proc. Natl. Acad. Sci. USA. 2013; 110 (24344313): 20878-20881
        • Horodyska J.
        • Wimmers K.
        • Reyer H.
        • Trakooljul N.
        • Mullen A.M.
        • Lawlor P.G.
        • Hamill R.M.
        RNA-seq of muscle from pigs divergent in feed efficiency and product quality identifies differences in immune response, growth, and macronutrient and connective tissue metabolism.
        BMC Genomics. 2018; 19 (30384851): 791
        • Huang D.W.
        • Sherman B.T.
        • Lempicki R.A.
        Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis of large gene lists.
        Nucleic Acids Res. 2009; 37 (19033363): 1-13
        • Huang D.W.
        • Sherman B.T.
        • Lempicki R.A.
        Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.
        Nat. Protoc. 2009; 4 (19131956): 44-57
        • Hurley A.M.
        • López-Villalobos N.
        • McParland S.
        • Kennedy E.
        • Lewis E.
        • O'Donovan M.
        • Burke J.L.
        • Berry D.P.
        Inter-relationships among alternative definitions of feed efficiency in grazing lactating dairy cows.
        J. Dairy Sci. 2016; 99 (26585474): 468-479
        • Hurley A.M.
        • López-Villalobos N.
        • McParland S.
        • Lewis E.
        • Kennedy E.
        • O'Donovan M.
        • Burke J.L.
        • Berry D.P.
        Genetics of alternative definitions of feed efficiency in grazing lactating dairy cows.
        J. Dairy Sci. 2017; 100 (28478005): 5501-5514
        • Kadenbach B.
        • Huttemann M.
        The subunit composition and function of mammalian cytochrome c oxidase.
        Mitochondrion. 2015; 24 (26190566): 64-76
        • Kanehisa M.
        • Goto S.
        KEGG: Kyoto Encyclopedia of Genes and Genomes.
        Nucleic Acids Res. 2000; 28 (10592173): 27-30
        • Kelly A.K.
        • Waters S.M.
        • McGee M.
        • Fonseca R.G.
        • Carberry C.
        • Kenny D.A.
        mRNA expression of genes regulating oxidative phosphorylation in the muscle of beef cattle divergently ranked on residual feed intake.
        Physiol. Genomics. 2011; 43 (20923863): 12-23
        • Khansefid M.
        • Millen C.A.
        • Chen Y.
        • Pryce J.E.
        • Chamberlain A.J.
        • Vander Jagt C.J.
        • Gondro C.
        • Goddard M.E.
        Gene expression analysis of blood, liver, and muscle in cattle divergently selected for high and low residual feed intake 1.
        J. Anim. Sci. 2017; 95 (29293712): 4764-4775
        • Koch R.M.
        • Swiger L.A.
        • Chambers D.
        • Gregory K.E.
        Efficiency of feed use in beef cattle.
        J. Anim. Sci. 1963; 22: 486-494
        • Kong B.W.
        • Lassiter K.
        • Piekarski-Welsher A.
        • Dridi S.
        • Reverter A.
        • Hudson N.J.
        • Bottje W.G.
        Proteomics of breast muscle tissue associated with the phenotypic expression of feed efficiency within a pedigree male broiler line: I. Highlight on mitochondria.
        PLoS One. 2016; 11 (27244447)e0155679
        • Kong R.S.G.
        • Liang G.
        • Chen Y.
        • Stothard P.
        • Guan L.L.
        Transcriptome profiling of the rumen epithelium of beef cattle differing in residual feed intake.
        BMC Genomics. 2016; 17 (27506548): 592
        • Korver S.
        • van Eekelen E.A.M.
        • Vos H.
        • Nieuwhof G.J.
        • van Arendonk J.A.M.
        Genetic parameters for feed intake and feed efficiency in growing dairy heifers.
        Livest. Prod. Sci. 1991; 29: 49-59
        • Kramer P.A.
        • Ravi S.
        • Chacko B.
        • Johnson M.S.
        • Darley-Usmar V.M.
        A review of the mitochondrial and glycolytic metabolism in human platelets and leukocytes: Implications for their use as bioenergetic biomarkers.
        Redox Biol. 2014; 2 (24494194): 206-210
        • Langfelder P.
        • Horvath S.
        WGCNA: An R package for weighted correlation network analysis.
        BMC Bioinformatics. 2008; 9 (19114008): 559
        • Lassiter K.
        • Ojano-Dirain C.
        • Iqbal M.
        • Pumford N.R.
        • Tinsley N.
        • Lay J.
        • Liyanage R.
        • Wing T.
        • Cooper M.
        • Bottje W.
        Differential expression of mitochondrial and extramitochondrial proteins in lymphocytes of male broilers with low and high feed efficiency.
        Poult. Sci. 2006; 85 (17135683): 2251-2259
        • Li B.
        • Fang L.
        • Null D.J.
        • Hutchison J.L.
        • Connor E.E.
        • VanRaden P.M.
        • VandeHaar M.J.
        • Tempelman R.J.
        • Weigel K.A.
        • Cole J.B.
        High-density genome-wide association study for residual feed intake in Holstein dairy cattle.
        J. Dairy Sci. 2019; 102 (31563317): 11067-11080
        • Liao Y.
        • Smyth G.K.
        • Shi W.
        featureCounts: An efficient general purpose program for assigning sequence reads to genomic features.
        Bioinformatics. 2014; 30 (24227677): 923-930
        • Liu Y.
        • Gu H.-Y.
        • Zhu J.
        • Niu Y.-M.
        • Zhang C.
        • Guo G.-L.
        Identification of hub genes and key pathways associated with bipolar disorder based on weighted gene co-expression network analysis.
        Front. Physiol. 2019; 10 (31481902)1081
        • Luke T.D.W.
        • Rochfort S.
        • Wales W.J.
        • Bonfatti V.
        • Marett L.
        • Pryce J.E.
        Metabolic profiling of early-lactation dairy cows using milk mid-infrared spectra.
        J. Dairy Sci. 2019; 102 (30594377): 1747-1760
        • MacRae J.C.
        • Lobley G.E.
        Some factors which influence thermal energy losses during the metabolism of ruminants.
        Livest. Prod. Sci. 1982; 9: 447-456
        • Nkrumah J.D.
        • Okine E.K.
        • Mathison G.W.
        • Schmid K.
        • Li C.
        • Basarab J.A.
        • Price M.A.
        • Wang Z.
        • Moore S.S.
        Relationships of feedlot feed efficiency, performance, and feeding behavior with metabolic rate, methane production, and energy partitioning in beef cattle.
        J. Anim. Sci. 2006; 84 (16361501): 145-153
        • Okonechnikov K.
        • Conesa A.
        • García-Alcalde F.
        Qualimap 2: Advanced multi-sample quality control for high-throughput sequencing data.
        Bioinformatics. 2016; 32 (26428292): 292-294
        • Olijhoek D.W.
        • Difford G.F.
        • Lund P.
        • Løvendahl P.
        Phenotypic modeling of residual feed intake using physical activity and methane production as energy sinks.
        J. Dairy Sci. 2020; 103 (32475658): 6967-6981
        • Ospina P.A.
        • Nydam D.V.
        • Stokol T.
        • Overton T.R.
        Evaluation of nonesterified fatty acids and β-hydroxybutyrate in transition dairy cattle in the northeastern United States: Critical thresholds for prediction of clinical diseases.
        J. Dairy Sci. 2010; 93 (20105526): 546-554
        • Pagliarini D.J.
        • Calvo S.E.
        • Chang B.
        • Sheth S.A.
        • Vafai S.B.
        • Ong S.E.
        • Walford G.A.
        • Sugiana C.
        • Boneh A.
        • Chen W.K.
        • Hill D.E.
        • Vidal M.
        • Evans J.G.
        • Thorburn D.R.
        • Carr S.A.
        • Mootha V.K.
        A mitochondrial protein compendium elucidates complex I disease biology.
        Cell. 2008; 134 (18614015): 112-123
        • Pfuhl R.
        • Bellmann O.
        • Kühn C.
        • Teuscher F.
        • Ender K.
        • Wegner J.
        Beef versus dairy cattle: A comparison of feed conversion, carcass composition, and meat quality.
        Arch. Tierzucht. 2007; 50: 59-70
        • Pryce J.E.
        • Gonzalez-Recio O.
        • Nieuwhof G.
        • Wales W.J.
        • Coffey M.P.
        • Hayes B.J.
        • Goddard M.E.
        Hot topic: Definition and implementation of a breeding value for feed efficiency in dairy cows.
        J. Dairy Sci. 2015; 98 (26254533): 7340-7350
        • Pryce J.E.
        • Gonzalez-Recio O.
        • Thornhill J.B.
        • Marett L.C.
        • Wales W.J.
        • Coffey M.P.
        • de Haas Y.
        • Veerkamp R.F.
        • Hayes B.J.
        Short communication: Validation of genomic breeding value predictions for feed intake and feed efficiency traits.
        J. Dairy Sci. 2014; 97 (24239085): 537-542
        • Pryce J.E.
        • Wales W.J.
        • de Haas Y.
        • Veerkamp R.F.
        • Hayes B.J.
        Genomic selection for feed efficiency in dairy cattle.
        Animal. 2014; 8 (24128704): 1-10
        • R Core Team
        R: A Language and Environment for Statistical Computing..
        R Foundation for Statistical Computing, Vienna, Austria2018
        • Robinson M.D.
        • McCarthy D.J.
        • Smyth G.K.
        edgeR: A bioconductor package for differential expression analysis of digital gene expression data.
        Bioinformatics. 2010; 26 (19910308): 139-140
        • Roche J.R.
        • Friggens N.C.
        • Kay J.K.
        • Fisher M.W.
        • Stafford K.J.
        • Berry D.P.
        Invited review: Body condition score and its association with dairy cow productivity, health, and welfare.
        J. Dairy Sci. 2009; 92 (19923585): 5769-5801
        • Salleh M.S.
        • Mazzoni G.
        • Höglund J.K.
        • Olijhoek D.W.
        • Lund P.
        • Løvendahl P.
        • Kadarmideen H.N.
        RNA-Seq transcriptomics and pathway analyses reveal potential regulatory genes and molecular mechanisms in high- and low-residual feed intake in Nordic dairy cattle.
        BMC Genomics. 2017; 18 (28340555): 258
        • Salleh S.M.
        • Mazzoni G.
        • Løvendahl P.
        • Kadarmideen H.N.
        Gene co-expression networks from RNA sequencing of dairy cattle identifies genes and pathways affecting feed efficiency.
        BMC Bioinformatics. 2018; 19 (30558534): 513
        • Spurlock D.M.
        • Dekkers J.C.M.
        • Fernando R.
        • Koltes D.A.
        • Wolc A.
        Genetic parameters for energy balance, feed efficiency, and related traits in Holstein cattle.
        J. Dairy Sci. 2012; 95 (22916946): 5393-5402
        • Szklarczyk D.
        • Gable A.L.
        • Lyon D.
        • Junge A.
        • Wyder S.
        • Huerta-Cepas J.
        • Simonovic M.
        • Doncheva N.T.
        • Morris J.H.
        • Bork P.
        • Jensen L.J.
        • Mering C.
        STRING v11: Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets.
        Nucleic Acids Res. 2019; 47 (30476243): D607-D613
        • Tang Y.
        • Horikoshi M.
        • Li W.
        ggfortify: Unified interface to visualize statistical result of popular R packages.
        R J. 2016; 8: 474
        • Veerkamp R.F.
        Feed intake and energy balance in lactating animals. In 7th World Congr. Genet. Appl. Livest. Prod. Vol. CD-ROM Communication No. 10-01, Montpellier, France.
        2002
        • Vincent A.
        • Louveau I.
        • Gondret F.
        • Trefeu C.
        • Gilbert H.
        • Lefaucheur L.
        Divergent selection for residual feed intake affects the transcriptomic and proteomic profiles of pig skeletal muscle.
        J. Anim. Sci. 2015; 93 (26115262): 2745-2758
        • Wang W.
        • Jiang W.
        • Hou L.
        • Duan H.
        • Wu Y.
        • Xu C.
        • Tan Q.
        • Li S.
        • Zhang D.
        Weighted gene co-expression network analysis of expression data of monozygotic twins identifies specific modules and hub genes related to BMI.
        BMC Genomics. 2017; 18 (29132311): 872
        • Wang Z.
        • Zhu J.
        • Chen F.
        • Ma L.
        Weighted gene coexpression network analysis identifies key genes and pathways associated with idiopathic pulmonary fibrosis.
        Med. Sci. Monit. 2019; 25 (31177264): 4285-4304
        • Wickham H.
        ggplot2: Elegant Graphics for Data Analysis..
        Springer-Verlag, New York, NY2016
        • Xiang R.
        • Hayes B.J.
        • Vander Jagt C.J.
        • MacLeod I.M.
        • Khansefid M.
        • Bowman P.J.
        • Yuan Z.
        • Prowse-Wilkins C.P.
        • Reich C.M.
        • Mason B.A.
        • Garner J.B.
        • Marett L.C.
        • Chen Y.
        • Bolormaa S.
        • Daetwyler H.D.
        • Chamberlain A.J.
        • Goddard M.E.
        Genome variants associated with RNA splicing variations in bovine are extensively shared between tissues.
        BMC Genomics. 2018; 19 (29973141): 521
        • Yang J.
        • Jiang J.
        • Liu X.
        • Wang H.
        • Guo G.
        • Zhang Q.
        • Jiang L.
        Differential expression of genes in milk of dairy cattle during lactation.
        Anim. Genet. 2016; 47 (26692495): 174-180
        • Yilmaz H.
        • Gül M.
        • Akkoyun S.
        • Parlakay O.
        • Bilgili M.E.
        • Vurarak Y.
        • Hizli H.
        • Kilicalp N.
        Economic analysis of dairy cattle farms in east Mediterranean region of Turkey.
        Rev. Bras. Zootec. 2016; 45: 409-416
        • Yoshikawa S.
        Beef heart cytochrome c oxidase.
        Curr. Opin. Struct. Biol. 1997; 7 (9266181): 574-579
        • Zhou N.
        • Lee W.R.
        • Abasht B.
        Messenger RNA sequencing and pathway analysis provide novel insights into the biological basis of chickens' feed efficiency.
        BMC Genomics. 2015; 16 (25886891): 195