Bivariate Genome-Wide Association Analysis of the Growth and Intake Components of Feed Efficiency

Single nucleotide polymorphisms (SNPs) associated with average daily gain (ADG) and dry matter intake (DMI), two major components of feed efficiency in cattle, were identified in a genome-wide association study (GWAS). Uni- and multi-SNP models were used to describe feed efficiency in a training data set and the results were confirmed in a validation data set. Results from the univariate and bivariate analyses of ADG and DMI, adjusted by the feedlot beef steer maintenance requirements, were compared. The bivariate uni-SNP analysis identified (P-value <0.0001) 11 SNPs, meanwhile the univariate analyses of ADG and DMI identified 8 and 9 SNPs, respectively. Among the six SNPs confirmed in the validation data set, five SNPs were mapped to KDELC2, PHOX2A, and TMEM40. Findings from the uni-SNP models were used to develop highly accurate predictive multi-SNP models in the training data set. Despite the substantially smaller size of the validation data set, the training multi-SNP models had slightly lower predictive ability when applied to the validation data set. Six Gene Ontology molecular functions related to ion transport activity were enriched (P-value <0.001) among the genes associated with the detected SNPs. The findings from this study demonstrate the complementary value of the uni- and multi-SNP models, and univariate and bivariate GWAS analyses. The identified SNPs can be used for genome-enabled improvement of feed efficiency in feedlot beef cattle, and can aid in the design of empirical studies to further confirm the associations.


Introduction
Optimization of feed efficiency in livestock production demands the consideration of the system inputs and outputs. In feedlot beef cattle enterprises, feed consumption dominates the input (and costs) and weight gain dominates the output (and return). Average daily gain (ADG) per animal, computed as the difference between final and initial trial weight divided by the number of days in the trial, is a frequently used indicator of weight gain. The cost of feed represents 62% to 84% of the total costs in a beef cattle production unit [1]. Dry matter intake (DMI) per day and animal is a frequently used indicator of feed consumption. In addition, 70 to 75% of the total energy feed intake is spent on maintenance functions (e.g. body temperature, digestion) in beef cattle [2]. Metabolic body weight (MBW) per animal, computed as BW 0.73 , is an accepted indicator of maintenance requirements.
Genomic improvement of feed efficiency in beef cattle relies on the identification of genomic variants (single nucleotide polymorphisms or SNPs) associated with feed efficiency components. A genome-wide association study (GWAS) can be used to identify SNPs to be included in genome-enabled selection decisions. The study of feed efficiency requires the consideration of output (ADG) and input (DMI) indicators, adjusted for maintenance requirements (MBW).
The majority of the SNPs reported to be associated with feed efficiency were identified from the analysis of each component (ADG or DMI) separately, or functions thereof such as residual feed intake and residual average daily gain [3][4][5]. On one hand, the analysis of feed efficiency components separately may fail to exploit the covariation between the components and consequently loose statistical precision to detect SNPs. On the other hand, the feed efficiency functions adjust either component by the other, thus imposing the selection of one component as the response and the assumption that the other component is an explanatory variable measured without error. Furthermore, the analysis of these functions fails to consider the uncertainty of the adjusted values. Bivariate analysis can augment the statistical precision to detect SNPs associated with both feed efficiency components. This gain stems from the consideration of covariation between the components that can augment the SNP signal relative to the noise or error [6][7][8][9][10][11][12]. No bivariate GWAS of feed efficiency components in beef cattle has been reported. The objectives of this study were: 1) to identify and characterize SNPs associated with feed efficiency components in a feedlot beef cattle population using bivariate analysis; 2) to compare the results from bivariate and univariate analyses; 3) to evaluate the results from uni-and multi-SNP models identified in a training data set on a validation data set; and 4) to enhance the interpretation of the results using functional genomic analyses and network visualization. Single nucleotide polymorphisms that exhibited favorable associations with both feed efficiency components, or that exhibited a favorable association with either component while minimizing a disfavorable trend on the other component were highlighted. These SNPs are well-suited for genome-enabled selection programs to improve feed efficiency and for follow-up empirical confirmation.

Ethics statement
All procedures were conducted following the guidelines recommended in the Guide for the Care and Use of Agricultural Animals in Agricultural Research and Teaching [13] with the approval of the University of Illinois Institutional Animal Care and Use Committee. The respective owners of the animals granted permission for their use in this study.

Beef cattle steers studied
Measurements were collected from 1,321 feedlot steers from five ranches in Montana between 2005 and 2008. The combination of ranch, harvest group, and harvest year resulted on 27 contemporary groups (CGs). Pedigree and breed information from 3,331 animals [14] were used to define the breed composition of each steer and to infer the relationship matrix. Steers pertain to one of five breed compositions: purebred Angus (AN), 3/4 Angus (3/ 4AN), crossbred Angus and Simmental (ANSM), 3/4 Simmental (3/4SM), or purebred Simmental (SM). The trial lasted an average (6 standard deviation) of 165 (616) days. Each steer received one of the twelve diets [15]. The diets were further grouped into five dietary treatments according to the main ingredient, total net energy, and non-degradable fiber (Table 1, [15]).

Measurements
Two feed efficiency components, ADG and DMI, were analyzed. Individual steer ADG (kg) was the difference between adjusted final weight (FW) and initial weight (IW) in the trial divided by days in the trial. The FW was estimated by dividing the individual hot carcass weight by the average dressing percentage of the harvest group. Individual daily DMI (kg) was measured using the GrowSafe automated feeding system (GrowSafe Systems Ltd., Airdrie, Alberta, Canada). Individual MBW was calculated using the estimated BW mid-trial. The age of the steer at mid-trial (mA; days) was also recorded. The average (6 standard deviation)  ADG, daily DMI, IW, FW, MBW, and mA

Genotyping and quality control
The DNA was extracted from blood samples using the salting out method [16]. The SNP genotypes were obtained from IlluminaH BovineSNP50 BeadChips v1 and v2 platforms (Illumina Inc., San Diego, CA) that include 54,001 and 54,609 SNPs, respectively. Quality control was performed in two steps on the 52,340 SNPs present in both platforms. In the first quality control step, SNPs not assigned to chromosomes, according to the Bos_taurus_UMD_3.1 assembly [17], and having a GenCall score ,0.2 (suggesting unreliable genotype [18]) were filtered. From this step, 519 and 16 SNPs were excluded. The software PLINK v.1.07 [19] was used to perform the second quality control step. In this step, SNPs and steers were removed when not meeting either one of the following criteria: missing steer per SNP ,20%) [20]; Hardy-Weinberg equilibrium test P-value .0.00001 [21]; missing SNP per steer ,10% [22]; and minor allele frequency .5% [22]. After the second quality control step, 264 SNPs, 1,202 SNPs, 9 steers, and 9,811 SNPs were not considered for further analysis applying the previous criteria, respectively. The final data set included 1,312 steers, 40,528 SNPs, and a total genotyping rate of 99.55%.

Statistical analyses
Uni-SNP model, univariate and multivariate analyses. Univariate (ADG or DMI) and bivariate (ADG and DMI) analyses of uni-SNP mixed-effect models were used to detect SNPs associated with feed efficiency. The uni-SNP model for the univariate analysis (Equation 1; Eq.1) was: where Y ijklmn denoted the observed ADG or DMI, m denoted the overall mean, SNP i denoted the fixed effect of an individual SNP genotype, B j denoted the fixed effect of breed (5 levels), D k denoted the fixed effect of diet (5 levels, Table 1), CG l denoted the random contemporary group effect (27 levels) that has a Normal distribution (0, s 2 CG ), b 1 denoted the fixed effect regression coefficient for the covariate mA, b 2 denoted the fixed effect regression coefficient for the covariate MBW, a ijklmn denoted the random animal polygenic effect that has a Normal distribution (0, As 2 a ) where A denoted the additive relationship matrix, and e ijklmn denoted the random normal distributed error (0, s 2 e ). The corresponding bivariate analysis (Eq. 2) was: where Y ADG and Y DMI denoted the vectors of observed ADG and DMI, respectively; X ADG and X DMI denoted the incidence matrices for the fixed effects for ADG and DMI, respectively; b ADG and b DMI denoted the vectors of solutions associated with X ADG and X DMI , respectively; Zu ADG and Zu DMI denoted the incidence matrices for the random contemporary groups for ADG and DMI, respectively; u ADG and u DMI denoted the vectors of solutions associated with Zu ADG and Zu DMI ; respectively; Za ADG and Za DMI denoted the incidence matrices for the random animal polygenic effects for ADG and DMI, respectively; a ADG and a DMI denoted the vectors of solutions associated with Za ADG and Za DMI , respectively; and e ADG and e DMI denoted the vectors of random errors associated with Y ADG and Y DMI , respectively; assuming random effects distributed as multivariate Normal that had mean equal to zero and covariance matrix:  where s 2 aADG and s 2 aDMI denoted the random animal polygenic variance for ADG and DMI, respectively; s aADG,DMI denoted the random animal polygenic covariance between ADG and DMI; s 2 CG ADG and s 2 CG DMI denoted the random contemporary group variance for ADG and DMI, respectively; I denoted the identity matrix; s 2 eADG and s 2 eDMI denoted the random error variance for ADG and DMI, respectively; and s eADG,DMI denoted the random error covariance between ADG and DMI. The models used in the univariate and bivariate analyses included the same explanatory variables.
The GWAS was implemented using Qxpak v.5.05 [23] and SNPs exhibiting associations with the feed efficiency components at an unadjusted P-value ,0.0001 were deemed significant. The additive and dominance effects were estimated for SNPs on autosomal chromosomes, and the additive effect was estimated for SNPs located on chromosome X. The additive effect estimate was computed relative to the less frequent (minor) allele among the steers studied. The additive effect estimate was defined as the change on the feed efficiency component per additional minor allele in the SNP genotype. The dominance effect estimate was defined as the difference on the feed efficiency component between the heterozygous and homozygous steers. Model assumptions including independence of residuals, homogeneity of variance, and normality were evaluated.
Multi-SNP model selection. A multi-SNP model was developed for the univariate and bivariate analyses. Starting with the SNPs detected at P-value ,0.001 in the uni-SNP models and all other explanatory variables equal, a stepwise feature selection approach was used. The final multi-SNP model included the SNPs that entered (were added to the model) and stayed (were kept in the model after consideration of all other SNPs in the model) in the multi-SNP model at P-value ,0.0001.
Training and validation of the uni-and multi-SNP results. The SNPs were first identified using complementary models and analyses on a training data sets. Subsequently, the findings were evaluated on a validation data set that included a separate group of steers. Training and validation data sets were generated from the records that passed the quality control based on sire family [4,24,25]. Steers were randomly assigned to either the training set (976 steers; 75%) or the validation set (336 steers; 25%; Table 2). Data partitioning ensured that each sire was represented in only one of the data sets to minimize potential confounding between SNP and individual associations [4,24].
The SNPs detected (P-value ,0.0001) in the training data set using the uni-SNP model and univariate and bivariate analyses were validated at P-value ,0.05 [24]. The trend (sign) of the genetic estimates was also compared between the training and validation data sets. For the multi-SNP models, the SNPs were validated based on the change in the model prediction accuracy, termed model adequacy (MA), between training and validation data sets. For the univariate analyses of ADG and DMI, the square root of the mean square error (RMSE) was used to indicate the difference between the observed and predicted values and thereof, model inadequacy. The change in MA for the univariate analyses of ADG and DMI (Eq. 3) was: where RMSE V is the RMSE from the validation data set and RSME T is the RMSE from the training data set. For the bivariate analysis, model inadequacy was the average of three root mean (co)variance terms: RMSE for ADG, RMSE for DMI, and the root means square covariance (RMSC) between ADG and DMI. The change in MA was computed as for the univariate analysis. Linkage disequilibrium. Some of the SNPs detected could be an artifact of the dependency between SNPs that exhibit high linkage disequilibrium (LD). This situation is the result of the average probe spacing (49.4 kb) of the platform [26] and the large number of SNPs tested. Statistical dependencies between significant SNP pairs located less than 500 kb apart that could suggest LD was assessed using the standard r 2 statistic [27] in PLINK. The LD extent in cattle is estimated to be 500 kb [28].
Genetic parameters. The genetic parameters of ADG and DMI were estimated to assess the potential amount of genetic variability that could be associated to SNPs. Heritability and genetic and phenotypic correlations between ADG and DMI were estimated using an animal model and univariate and bivariate analyses implemented in WOMBAT [29,30]. The explanatory variables included in the animal models encompassed those described in Eq. 4.
Functional and gene network analyses. The detected SNPs were mapped to harboring or proximal (within 2 kb of the 59 untranslated region or 0.5 kb of the 39 untranslated region of a gene) genes in the Bos_taurus_UMD_3.1 assembly. The SNP mapping and gene information was obtained from the National Center for Biotechnology Information, SNP and Gene databases [31].
Functional analysis of the genes corresponding to the SNPs detected (P-value ,0.01) in the bivariate analysis offered insights into the categories enriched among the genes. The consideration of genes from the bivariate analysis was motivated by the goal of identifying functional categories among genes that could have pleiotropic effects on both feed efficiency components. Genes farther upstream and downstream from the detected SNPs were not included in the functional analysis because the number of spurious (false positive) genes added to the functional analysis could have overwhelmed the fewer real (true positive) loci, potentially biasing the results. The enrichment of Gene Ontology (GO) FAT categories and KEGG pathways among the genes was investigated using Fisher's exact test in DAVID [32,33]. The GO FAT categories are a subset of the broadest GO terms, filtered to minimize overshadowing of more specific terms due to repetition of more general categories. Functional annotation charts were considered significant at P-value ,0.001 using the Bos Taurus genome as background.
Gene networks associated with feed efficiency were visualized using the genes affiliated to the enriched functional categories. The network was visualized using the BisoGenet plug-in Cytoscape [34], with default settings. Identified (or target) genes and intermediate connecting genes from the NCBI database genes were represented by nodes. The final pathway included target genes separated by at most two intermediate genes. Edges denoted known relationships between genes.

General results
The heritability estimates of ADG and DMI were 0.14 and 0.25, respectively. These heritability estimates confirmed the opportunity for genome-based improvement of these feed efficiency components. The phenotypic and genetic correlations between ADG and DMI were 0.52 and 0.18, respectively. These estimates were consistent or slightly lower than in previous reports [35,36]. The positive genetic correlation supports the hypothesis that SNP alleles that have positive association with ADG and negative association with DMI could be identified and could assist with genome-based improvement feed efficiency in beef cattle.
A summary of the number of significant SNPs (and corresponding genes) detected by the univariate and bivariate analyses is presented in Table 3. Among the 28 SNPs detected (P-value ,0.0001), 19 SNPs overlapped between the univariate and bivariate analyses. The bivariate analysis detected the highest number of SNPs (11 SNPs) followed by the univariate analyses of DMI (9 SNPs) and ADG (8 SNPs). Similar associations between SNPs and ADG or DMI have been previously reported [37,24]. The partial overlap of SNPs confirmed the complementary information offered by the univariate and bivariate analyses. For certain SNPs, the bivariate analyses could gain precision through the consideration of covariation between ADG and DMI relative to univariate analyses. For other SNPs, univariate analyses benefited from lower noise of each trait studied separately, relative to the bivariate analysis. Associations with SNPs were identified on 10 chromosomes (Tables 4 and 5 The potential relationship between feed efficiency and the genes harboring or in the proximity of the detected SNPs was investigated. The connection between the feed efficiency components and genes was based on gene annotation information available at the National Center for Biotechnology Information, Gene database [31]. This information was complemented with literature review where relevant.

Univariate uni-SNP analysis of ADG
Previous studies reported genomic regions associated with ADG on BTA 2,[4][5][6][7]9,11,[14][15][16][17][18][19][20]22,23,26, and 28 [38][39][40][41][42][43][44][45][46][47][48][49][50][51][52]. The SNPs associated with ADG in the present study (Table 4) were mapped to chromosomes previously linked to ADG, with the exception of rs41629972 located on BTA13. This SNP is located approximately 33 kb upstream Kruppel-like factor 6 (KLF6) and is within the 500 kb LD span reported in cattle [28,29]. The zinc finger protein encoded by this gene has been associated with cell proliferation, differentiation, signal transduction, and cell death [31, 53,54]. KLF6 regulates genes in the transforming growth factors b signaling pathway [55]. Transforming growth factors b superfamily members affect both muscle development and postnatal skeletal muscle mass [31, 56,57]. Every additional T allele was associated with 0.04 kg higher ADG relative to the C allele. Also associated with ADG but not mapped to a gene, rs41565199 was mapped on BTA14 within LD reach (464 kb downstream) of zinc finger and homeobox 2 (ZHX2). This QTL   region has been linked to ADG in Japanese Black (Wagyu) cattle [47]. Homozygous CC steers exhibited the highest ADG. The remaining six SNPs associated with ADG were mapped to known genes. Steers heterozygous for rs109934193 on BTA2and rs110787048 on BTA4, had higher ADG than the average homozygous steer. These SNPs were mapped to NCK-associated protein 5 (NCKAP5) and dipeptidyl-peptidase 6 (DPP6), respectively. The three SNPs on BTA15 associated with ADG were: rs41620774 mapped to engulfment and cell motility/CED-12 domain containing 1 (ELMOD1), rs108964818 mapped to Lys-Asp-Glu-Leu containing 2 (KDELC2), and rs41768978 mapped to paired-like homeobox 2a (PHOX2A). Both rs41620774 and rs108964818 were mapped less than 1 Mb apart, and for both SNPs, steers homozygous for the minor allele (CC and TT, respectively) had the lowest ADG. The protein encoded by ELMOD1 has a GTPase-activator function on small G proteins of the arf family [58]. These proteins have a central role in the organization of the secretory and endocytic pathways [59]. Mapped 4 Mb from a QTL previously associated with ADG [45], steers heterozygous for rs41768978 exhibited 0.060 kg higher ADG than the average homozygous steer. This SNP was mapped to the intronic region of PHOX2A, a gene associated with respiratory rhythm (and thus biochemical energy) and autonomic nervous system development [31,60]. Lastly rs42342964, mapped to the PAK1 interacting protein 1 (PAK1IP1) on BTA23, was associated with ADG. The p21-activated protein kinase-interacting protein 1 encoded by this gene has been associated with cell proliferation and signal transduction [31,60].
Four additional SNPs associated with DMI were not mapped to genes. Every additional C allele on rs41663978 was associated with lower daily DMI (20.2260.05 kg) relative to the A allele. This SNP was mapped on BTA6, approximately 5 Mb from a QTL previously associated with DMI [45]. In addition, rs41663978 is located within 500 kb from ADAM metallopeptidase with thrombospondin type 1 motif, 3 (ADAMTS3), groupspecific component vitamin D binding protein (DBP), and neuropeptide FF receptor 2 (NPFFR2). These genes have been associated with protein processing [60], leanness [31], and obesity [31], respectively. Similarly, rs41632270 was mapped to a QTL region on BTA13 associated with DMI [63] and within 500 kb of several genes including kinesin family member 16B (KIF16B), Nacetylneuraminic acid phosphatase (NANP), otoraplin (OTOR), phosphoribosylaminoimidazole carboxylase pseudogene (PAICSP), and small nuclear ribonucleoprotein polypeptide B (SNRPB2). In particular, NANP muschlparticipates on the amino sugar and nucleotide sugar metabolism [31]. Heterozygous steers for rs41632270 exhibited the lowest daily DMI (20.4060.10 kg) relative to the average homozygous steer.
Although rs42128656 was associated with DMI, this SNP was not mapped to a known gene. This SNP was 600 kb from other two SNPs (rs43291568 and rs43291603) on BTA15 also associated with DMI. This pair of SNPs were mapped within 27 kb of each other, however the pairwise LD among these SNPs was low (r 2 = 0.052). These two SNPs are in the intronic region of the coxsackie virus and adenovirus receptor-like membrane protein (CLMP). This gene encodes for a type I transmembrane protein of the CTX family and transmemrane proteins have activating or suppressing roles on cell growth [65].
Among the SNPs associated with DMI and mapped to known genes, rs108942504 was found on a gene on BTA22 that encodes a structural protein, the transmembrane protein 40 (TMEM40). Similarly, rs41588990 was mapped to the intronic region of CCR4-NOT transcription complex, subunit 6-like (CNOT6) on BTA6. This gene plays a role in the deadenylation of mRNAs in the cytoplasm and deadenylation has been associated with cell growth [66]. The results from the univariate GWAS offer a first glimpse of SNPs that could be used in genomic improvement of feed efficiency. However, univariate analyses may miss detecting SNPs because these analyses do not exploit the correlation between ADG and DMI. This situation could be especially detrimental in scenarios with limited data size, limited effect of the SNP or limited disequilibrium between the SNP and the QTL influencing the feed efficiency components. Bivariate uni-and multi-SNP analyses shed additional light on these SNPs.

Bivariate uni-SNP analysis of ADG and DMI
The SNPs simultaneously associated with ADG and DMI (Pvalue ,0.0001) were presented in Table 5. Many of these SNPs that have a pleiotrophic association with both feed efficiency components were also found in either univariate analysis. The SNPs that overlapped between the univariate and bivariate analyses were summarized in Table 6. Four of the eleven SNPs detected in the bivariate analysis were also detected in the univariate ADG analysis and other five SNPs were detected in univariate DMI analysis ( Table 6). The two additional SNPs detected by the bivariate analysis but not detected by the univariate analyses were rs41722387 and rs110522962. These results emphasize the need to consider the results of multivariate and univariate GWAS to precisely detect and characterize SNPs associated with feed efficiency. Among the SNPs uncovered by the uni-SNP bivariate analysis, rs41722387 was mapped approximately 450 kb from rs41565199, a SNP associated with ADG on BTA14 (Table 4). Despite being physically close, the LD between these SNPs was low (r 2 = 0.046). Furthermore, rs41565199 was mapped to pseudo metalloendopeptidase (OMA1) that inhibits growth and approximately 400 kb upstream from hyaluronan synthase 2 (HAS2) that mediates cellular growth [31]. Homozygous TT steers for rs41565199 exhibited the highest feed efficiency due to higher ADG and lower DMI relative to steers that had other genotypes. The other SNP detected solely by the bivariate analysis was rs110522962. Mapped on BTA17, this SNP was approximately 4 Mb from a QTL region previously associated with ADG [48]. Homozygous TT steers exhibited the highest ADG and lowest daily DMI relative to steers that had other genotypes.
Among the remainder nine SNPs detected by the bivariate uni-SNP analyses, rs42128656, rs111010038, and rs108942504 were also associated with DMI in the univariate analyses. Six additional SNPs were detected by the bivariate analyses and had the same signs for both feed efficiency components. Positive associations of the same allele with both feed efficiency components are not always undesirable because a significant increase in ADG could compensate for a less significant increase in DMI. For example, rs108964818 (that maps to KDELC2) had positive associations with ADG and DMI albeit a much higher additive estimate for ADG than DMI. The SNPs detected by the one univariate and the bivariate analyses (Table 6) were also detected at a less stringent threshold (P-value ,0.01; data not shown), by the other univariate analysis.

Univariate and multivariate multi-SNP analyses
The polygenic nature of ADG and DMI can be adequately described with a multi-SNP model. In turn, the multi-SNP function can be used to predict feed efficiency or in genomeenabled selection programs to improve feed efficiency. Findings from the un-SNP analyses were used to develop a multi-SNP predictive equation. A multi-SNP model was developed using stepwise selection and 53, 58, and 84 SNPs (P-value ,0.001) identified by the uni-SNP univariate ADG, DMI and bivariate analyses. The final multi-SNP ADG, DMI, and the bivariate models included nine, eight, and seven SNPs, respectively (P-value ,0.0001; Table 7). These SNPs encompassed 21 unique SNPs on 10 genes. Among these, 11 SNPs were detected by the uni-SNP analyses.
The multi-SNP model selection for the bivariate analysis of ADG and DMI resulted in three SNPs not detected by the uni-SNP bivariate model. Among these, rs109945988 was associated in the multi-SNP ADG model, and rs41600811 and rs42459305 represent new associations. In particular, rs41600811 was mapped 100 kb downstream the cell adhesion molecule with homology to L1CAM (close homolog of L1; CHL1) on BTA22. The cell-cell adhesion function of this gene enables cells to assemble into organized tissues. The remaining SNP, rs42459305, was mapped to a gene dense region on BTA29 with several predicted loci and located less than 3 kb downstream from the olfactory receptor, family 8, subfamily G, member 5 (OR8G5) and could impact feed consumption.

Validation
The SNPs detected on the training data set were evaluated on the validation data set. Findings from the uni-SNP models were confirmed at P-value ,0.05 in the validation data set. This threshold was used for two reasons. First, this validation constitutes the second of the two-phase approach. The SNPs have already been detected in the training data set at a P-value ,0.0001. Second, a limited number of SNPs required validation. The SNPs (and analyses) validated were: rs108942504 (univariate DMI); rs41629972 and rs108964818 (univariate ADG); and rs108964818, rs41768978, and rs110522962 (bivariate).
For the multi-SNP models, validation was assessed by the change in the MA between the validation and training data sets, relative to the validation set (Table 7). Overall, the MA in the small validation data set was comparable to that in the larger training data set used to detect the SNPs. This result confirmed that the SNPs detected have a high likelihood to be replicable in additional populations. For the multi-SNP univariate DMI (ADG) analysis, the RMSE only increased 7.21% (11.67%) despite the fact that the validation data set was 300% smaller than the training data set. The higher loss in MA for the bivariate multi-SNP analysis on the validation data set (19.4%) may be due the higher parameterization of the model and lower precision of each estimate relative to univariate analyses.

Functional analyses and gene networks visualization
The 236 genes corresponding to the SNPs detected at P-value ,0.01 by the uni-SNP bivariate analysis were considered for functional analysis. The P-value threshold was selected because of the high number of SNPs detected by the univariate analyses at Pvalue ,0.01 that were also detected in the bivariate analysis.
Seven functional categories were enriched (P-value ,0.001) among the genes corresponding to the detected SNPs (Table 8). The most significant categories included the GO molecular functions of cation channel activity and metal ion transmembrane transporter activity. Both categories encompassed 10 genes. Affiliated to these two GO categories was transient receptor  potential channel 2 (TRPC2). This gene was reported to be associated with several behavioral responses [60] and could be related to consumption and energy maintenance requirements. This gene corresponded to rs41603221 that was detected (P-value = 0.0006) in the bivariate analysis. Furthermore, this SNP was mapped 4 Mb from a QTL reported to be associate with ADG on BTA 15 [45]. Ion channel activity was associated with maintenance of normal gradient on plasma membranes, participation in cellular de-and re-polarization, neurotransmitter release, immune function, insulin secretion, and active transport mechanisms required for the digestion and absorption of nutrients [67,68,69].
A comprehensive network of the genes affiliated to the enriched molecular functions is reconstructed (Figure 1). In this network, the highest numbers of connections were displayed by the target gene FGF2 and the intermediate gene ubiquitin C (UBC). These genes could have driver or hub role on feed efficiency components. In addition, FGF2 was implicated on smooth muscle cell differentiation and signaling [31,70].

Conclusions
Single nucleotide polymorphisms associated with the feed efficiency components ADG and DMI in feedlot beef steers were identified using uni-SNP and multi-SNP models and univariate and bivariate analyses. The complementary set of SNPs detected by the univariate and bivariate analyses confirmed the value of considering both GWAS approaches. For certain SNPs, the bivariate analyses could gain precision through the consideration of covariation between ADG and DMI relative to univariate analyses. For other SNPs, univariate analyses could benefit from lower noise of each trait studied separately, relative to the bivariate analysis. Genomic loci that had favorable associations with ADG and DMI simultaneously, or favorable associations with either trait with minimum detrimental association with the other trait, while accounting for the body maintenance requirements, were identified. The validation of models and SNPs suggest that the findings could be replicable. Functional analysis and gene network visualization facilitated the interpretation of the association between SNPs mapping to genes that have ion channel-related molecular function and feed efficiency components. Results from this study can be used for genome-enabled improvement of feed efficiency in feedlot beef cattle, to support further empirical confirmation of the associations, and as proof of concept of the value of complementary association analyses.