Genomic Regions Associated with Feed Efficiency Indicator Traits in an Experimental Nellore Cattle Population

The objective of this study was to identify genomic regions and metabolic pathways associated with dry matter intake, average daily gain, feed efficiency and residual feed intake in an experimental Nellore cattle population. The high-density SNP chip (Illumina High-Density Bovine BeadChip, 777k) was used to genotype the animals. The SNP markers effects and their variances were estimated using the single-step genome wide association method. The (co)variance components were estimated by Bayesian inference. The chromosome segments that are responsible for more than 1.0% of additive genetic variance were selected to explore and determine possible quantitative trait loci. The bovine genome Map Viewer was used to identify genes. In total, 51 genomic regions were identified for all analyzed traits. The heritability estimated for feed efficiency was low magnitude (0.13±0.06). For average daily gain, dry matter intake and residual feed intake, heritability was moderate to high (0.43±0.05; 0.47±0.05, 0.18±0.05, respectively). A total of 8, 17, 14 and 12 windows that are responsible for more than 1% of the additive genetic variance for dry matter intake, average daily gain, feed efficiency and residual feed intake, respectively, were identified. Candidate genes GOLIM4, RFX6, CACNG7, CACNG6, CAPN8, CAPN2, AKT2, GPRC6A, and GPR45 were associated with feed efficiency traits. It was expected that the response to selection would be higher for residual feed intake than for feed efficiency. Genomic regions harboring possible QTL for feed efficiency indicator traits were identified. Candidate genes identified are involved in energy use, metabolism protein, ion transport, transmembrane transport, the olfactory system, the immune system, secretion and cellular activity. The identification of these regions and their respective candidate genes should contribute to the formation of a genetic basis in Nellore cattle for feed efficiency indicator traits, and these results would support the selection for these traits.


Introduction
Costs associated with animal feeding can be up to 50% of the total cost in beef cattle production systems, and there is growing interest concern in adopting strategies to reduce these costs [1]. So improving feed efficiency would increase both sustainability and profitability in the beef cattle industry. Selection for feed efficiency indicator traits in beef cattle can reduce production cost, decrease use of natural resources, and reduce impacts on the environment, optimizing production efficiency [2].
Feed efficiency (Kg gain/ Kg feed intake) showed a moderate positive correlation with weight gain and mature weight [3]. Thus, selection to improve feed efficiency could increase mature weight and increase energy maintenance requirements [4]. Also, residual feed intake (RFI), as proposed by Koch et al. [5], defined as the difference between actual feed intake and the feed intake required to meet maintenance requirements and growth, offers advantages over G:F since there is no correlation with performance traits [2,6,7,8,9].
Heritability estimates for RFI have been reported as moderate by many authors, ranging from 0.30 to 0.45 [3,10,11,12,13], suggesting that this trait can be improved by selection. However, RFI is difficult and expensive to measure, which often limits its implementation a selection criteria in beef cattle breeding programs. Recently, genome wide association studies (GWAS), using a high-density genotyping array, have been applied aiming to discover genomic regions associated with feed efficiency traits [11,14,15,16]. It is important to highlight that most of these studies have been implemented with taurine breeds (Bos taurus), and there are few studies for indicine breeds (Bos indicus) under tropical conditions. Therefore, in order to identify genomic regions associated with feed efficiency indicator traits, as well as to elucidate the genetic basis of them, it is important to encourage genomic studies with zebu animals, since zebu breeds are prevalent in herds under tropical and subtropical conditions. The objective of this study was to identify genomic regions and metabolic pathways associated with dry matter intake (DMI), average daily gain (ADG), feed efficiency (G:F) and residual feed intake (RFI) in an experimental Nellore cattle population.

Material and Methods Data
This study was approved by ethics committee of the Faculty of Agrarian and Veterinary Sciences, Sao Paulo State University (UNESP).
The data set used in this study was provided by the APTA Beef Cattle Center-Institute of Animal Science (IZ), Sertãozinho, São Paulo, Brazil. Phenotypic information is animals born from 2004 to 2012. These animals belong to three experimental lines of Nellore cattle, which have been selected since 1978 for yearling weight: selection line (NeS) is a closed herd selected for higher yearling weight; traditional line (NeT) is submitted to the same selection criterion as NeS but, eventually, receives animals from NeS; and a control line (NeC) selected for average of yearling weight [17]. The analyzed data was obtained by feed efficiency tests performed from 2005 to 2013, consisting of 541 males and 355 females. Some of these animals (n = 683) were restricted to individual troughs, which offered daily feed and refusal was controlled, while the others (n = 213) were held in two collective pens equipped with the GrowSafe 1 feeding system.
After weaning, animals were kept in the test for during 83±15 days, preceded by 28 days of adaptation, for evaluation the feed intake and average daily gain (ADG). Animals were weighed every 14 days after 12 hours of fasting. The diet was formulated with 67% of total digestible nutrients (TDN) and 13% of crude protein (CP), allowing ADG of 1.0 kg/day.

Traits
In order to ensure ad libitum feed intake the food supply was adjusted daily, allowing refusals varying from 5-10% of offer. The following feed intake data was not considered in the analyses: for days when animals were handled outside of the facilities for a number hours, during failure (GrowSafe) and when no refusals were found. Diet dry matter percentage was determined from weekly samples of offer and refusals. The ADG in each test was considered as the linear regression coefficient of body weight (BW) on days in test (DIT): where, y i = BW in i th observation; α = intercept of regression equation corresponding to the initial BW; β = linear regression coefficient corresponding to ADG; DIT i = days in test for i th observation; and ε i = random error associated with each observation. Metabolic weight (BW 0.75 ) was calculated as: Ã 0:75 , with α and β assuming the values obtained by the equation described above. Feed efficiency (G:F) was calculated as the ratio of ADG to DMI. The residual feed intake (RFI) was considered as error of linear regression equation of dry matter intake on average daily gain and metabolic weight within each contemporary group (CG: sex, year of birth, and pen), as shown below: where, βT, βTA, and βTB are regression coefficients of classificatory variable test group (TG) and of interactions between TG and covariates ADG and BW 0.75 , respectively; and ε is RFI. The descriptive statistics for DMI, ADG, G:F and RFI are presented in Table 1.

DNA extraction
The extraction of DNA from blood samples was performed using the DNeasy Blood & Tissue Kit (Qiagen). The DNA purification was performed using a column containing silica fragments (column purification). Firstly, DNA binds to the membrane of the extraction column and then it was washed until it has high purity. At the end of the protocol, the DNA was eluted with buffer AE (blood). The amount of scanning and quality of the material obtained was taken with the use of a spectrophotometer apparatus (NanoDrop 1000, Thermo Scientific, USA, 2008). Quality was measured by absorption ratio A260 / A280. A ratio of less than 1.8 suggests contamination from protein, and more than 2.0 suggests RNA contamination.

Genotyped Animals
Animals were genotyped using the high-density SNP chip (Illumina High-Density Bovine-BeadChip 777,000 SNP). Markers with minor allele frequency (MAF) and call rate higher than

Quantitative genetic analyses
The contemporary groups (CGs) included animals born on the same farm in the same year, and from the same management group as yearlings. The CGs with fewer than 3 records were eliminated from the analyses to maintain variability in the CGs. Records exceeding 3 standard deviations above or below the mean of each CG were excluded, to avoid the inclusion of possible measurement error or outliers. The model included the random additive animal effect, the fixed effects of CG, calving month, age of animal at beginning of the test (linear effect), and the dam age as co-variable (linear and quadratic effect). The (co)variance components were estimated using the single step genomic BLUP (ssGBLUP), under Bayesian inference [18]. The ssGBLUP is a modified version of the animal model (BLUP) with additive relationship matrix A -1 replaced by H -1 [19]: where A 22 is a numerator relationship matrix for genotyped animals and G is the genomic relationship matrix created as described by VanRaden et al. (2009) [20]: where Z is the gene matrix containing allele frequency adjustment; D is the matrix that have the SNP weight (initially D = I); and, q is a weighting / standardization factor. According to Vitezica et al. (2011) [21], such factors can be obtained by ensuring that the G average diagonal is next to A 22 The model can be represented by the following matrix equation: where y is the observations vector; β is the vector of fixed effects; a is the additive direct vector; X is known as incidence matrix; Z is the incidence genetic random effects additive direct matrix (the β vector associated with the y vector); e is the residual effect vector. The priori distributions of vectors y, a and e were given by: where H is the relationship coefficients matrix among animals obtained from the single-step analyzes (single-step); R is the residual variance matrix; I is the identity matrix; G is the genetic additive variance matrix and is the Kronecker product. An inverted qui-square distribution was used for the prior values of the direct and residual genetic variances. A uniform distribution was used for the priori for the fixed effects a uniform distribution. A total of 3,342 animals were considered in the pedigree file. Analyses were performed using GIBBS2F90 [19,22]. The a posteriori estimates were obtained using the application POSTGIBBSF90 [22]. The analyses were originated from chain lengths of 1,000,000 interactions, where the first 200,000 interactions were discarded. For parameter estimation, the samples were stored at each 100 cycles, building samples with 80,000 information. The data convergence was verified with the interactions versus evaluation graphic of sampled values and using the criteria proposed by Geweke; Heidelberger, Welch; and Raftery, Lewis [23][24][25] through analysis package Bayesian Output Analysis (BOA) in the software R 2.9.0 (The R Development Core Team, 2009).

Genome-wide association analysis
Genome-wide association analysis for each trait was performed using the single-step GWAS (ssGWAS) methodology [26]. The same linear animal model used to estimate the (co)variance components was applied. The effects were decomposed in genotyped (a g ) and ungenotyped (a n ) animals, as describe by Wang et al. [26], considering the effect of genotyped animals as: where Z is a matrix that relates genotypes of each locus and u is a vector of marker effects, and the variance of animal effects was assumed as: where D is a diagonal matrix of weights for variances of markers (D = I for GBLUP), σ 2 u is the genetic additive variance captured by each SNP marker when no weights are present, and G Ã is the weighted genomic relationship matrix.
The ratio of covariance of genetic effects (a g ) and SNPs (u) is: sequentially: where λ is a variance ratio or a normalizing constant. According to VanRaden et al. [20], where M is the number of SNP and p i is the allele frequency of the second allele in the i th SNP. According to Stranden e Garrick [27], the markers effects can be described by: The estimated SNP effects can be used to estimate the variance of each individual SNP effect [28] and apply a different weighting for each marker, such as: The following iterative process described by Wang et al. [26] was used considering D to estimate the SNP effects: The effects of markers were obtained by 2 iterations from step 2 to 7. The percentage of genetic variance that is responsible for i th region was calculated as described by Wang et al. [26]: where a i is the genetic value of the i th region that consists of contiguous 10 consecutive SNPs, σ 2 a is the total genetic variance, Z j is the vector of gene content of the j th SNP for all individual, and û j is the marker effect of the j th within the i th region. The results were presented by the proportion of variance expressed by each window of 10 SNPs. In addition, the genes located at ±500 Kb of each window were considered.

Search for genes
The chromosome segments that are responsible for more than 1.0% of additive genetic variance were selected to explore and determine possible quantitative trait loci. The bovine genome Map Viewer was used for identification of genes, available at "National Center for Biotechnology Information" (NCBI -http://www.ncbi.nlm.nih.gov) [29] in UMD3.1 version bovine genome and Ensembl Genome Browser (http://www.ensemble.org/index.html) [30]. The classification of genes for biological function, identification of metabolic pathways and enrichment of genes was performed on the website "The Database for Annotation, Visualization and Integrated Discovery (DAVID) v. 6.7" (http://david.abcc.ncifcrf.gov/) [31] and GeneCards (http:// www.genecards.org/) [32].

Genetic parameter estimates
The (co)variance components and heritability estimates for DMI, ADG, G:F and RFI are showed in Table 2. The criteria used to diagnose the chain convergence indicated convergence of all estimated parameters. Thus, the burn-in period considered was sufficient to reach the convergence in all parameter estimates. The posterior marginal distributions of heritability estimates for feed efficiency indicator traits showed accurate values, tending to normal distribution. The symmetric distributions of central tendency statistics of analyzed traits indicated that the analyses are reliable.
The heritability estimated for G:F showed low magnitude. Grion et al. [12] and Ceacero et al. [33] found low heritability for the same trait (0.17 ± 0.07, 0.14 ± 0.06, respectively). For DMI and ADG, the estimated heritability was moderate to high. Baldi et al. [34] in a study with Nellore animals and modeling the weights with random regression models, found lower heritability estimates for average daily gain during the performance test (0.21). Similar heritability estimates (0.47) to those obtained in the present study for DMI was reported by Bolormaa et al. [11] with taurine and zebu animals. The heritability obtained for RFI showed moderate magnitude. Recently, Grion et al. [12], working with Nellore animals obtained a higher RFI heritability estimate (0.33± 0.10). Bolormaa et al. [11] also reported higher heritability estimates for RFI (0.36), in a study with nine herds of Bos Taurus, Bos indicus, and crossbreeds. Also, Silva et al. [35] reported similar heritability estimates for ADG (0.39±0.08) and DMI (0.43±0.08) in a study using the same experimental population. The results of this study pointed out that there is genetic variability in selecting for feed efficiency indicator traits in Nellore cattle. Thus, it is important to known whether there are more genes involved to better understand the genetic architecture of these traits.
The known genes found in the regions that accounted for more than 1.0% of additive genetic variance are presented in tables according to the studied trait. The results indicated a total of 8, 17, 14 and 12 different windows with known genes that are responsible for more than 1.0% of the genetic variance for DMI, ADG, G:F, and RFI, respectively.

Genomic regions
For DMI, eight genomic regions that are responsible for more than 1.0% of the additive genetic variance were found (Table 3; Fig 1). The window that is responsible for the most part of additive genetic variance for DMI was located in chromosome BTA4 where one candidate gene, called NXPH1, was found associated with the DMI. The RFX6, GPRC6A, FAM162B, KPNA5 and ZUFSP genes were identified in the window located in chromosome BTA9 at the 34 Mb position that is responsible for 3.40% of the additive genetic variance. The RFX6 gene is related to the regulation of transcription. Results of a study in rats have suggested that this gene acts on the differentiation of cells in insulin production [36].
In the associated window located in BTA15 at the 46 Mb position the OR2D2, OR2D3 and OR10A4 genes associated with DMI were found, which have been reported to play roles of olfactory receptors, coupled to G proteins, to recognize and mediate the olfactory signals in humans and rats [37,38]. The ZNF214, ZNF215 and NLRP14 genes were found in the same chromosome. The NLRP14 gene encodes a protein with activity related to the immune system of the animal [39]. Three genomic regions located on chromosome BTA18 were found associated with DMI. Among the genes found in those regions, the CACNG7 and CACNG6 genes are associated with calcium channels [40].
The results of gene enrichment analysis and functional clusters showed that genes associated to DMI are involved in functions related to ion transport (Table A in S1 Tables). The metabolic pathways (Table 4) showed that those significant genes (p-value < 0.05) encode for a protein domain known as a zinc finger that acts as a structural element in proteins. Zinc is essential for several biochemical and cellular signaling pathways, participating in DNA and RNA synthesis and division and cell activation. Also, Zinc is essential for immune response [41].
For ADG, seventeen genomic regions located in 14 different chromosomes are responsible for more than 1.0% of the additive genetic variance (Table 5; Fig 2). In the window located in BTA1 at position 75 Mb, two candidate genes associated with ADG (MB21D2 and FGF12) were found. The OR52J3 and OR51A7 genes found in a window in BTA15 encode olfactory  receptors. The CAPN8, CAPN2 and TP53BP2 genes are located in BTA16, where the first two encode subunits of the calpain enzyme. Calpain is an enzyme related to the tenderness of beef after slaughter [42]. Ten genes were located in the window, in BTA18 at position 49 Mb, that are responsible for the greatest proportion of additive variance (6.07%) of ADG. Among these genes, three of them (MAP3K10, CNTD2 and AKT2) encode kinases protein, which belongs to the largest family of proteins in eukaryotes. These genes also play an important part in intracellular communication, regulation and signal transduction, and catalyzed the phosphorylation of proteins by ATP transfer [43]. The enrichment analyses for ADG revealed functional clusters related to the catabolism of macromolecules, transcription, protein kinase and binding nucleotides ( Table B in S1 Tables). The gene enrichment and metabolic pathways showed that those significant genes (p <0.05) are related to endocytosis (extracellular transport) and merger of myoblasts, which are precursor cells of muscle fibers, and formation of syncytia (multinucleated cells) ( Table 6).
Several genes were found in the 14 genomic regions that are responsible for more than 1% of the additive genetic variance of G:F (Table 7, Fig 3). The GOLIM4 and SERPINI1 genes were found in the window located in BTA1 at position 100 Mb. The GOLIM4 gene encodes an integral membrane protein localized in the Golgi apparatus, which is the main organelle in the secretory pathway in eukaryotic cells [44]. In chromosome BTA4 at position 26 Mb, the candidate genes PRPS1L1 and HDAC9 were found. The HDAC9 gene is a histone deacetylase enzyme, which is related the transcription factor. The position 78 Mb (BTA7) made an important contribution to the additive genetic variance for G:F (7.74%), but no candidate gene was found in this genomic region. The RASEF and FRMD3 genes were found in the window located in BTA8 at position 77 Mb. The RASEF gene encodes the GTP binding protein and binding calcium ions involved in the regulation of membrane traffic.
The associated candidate genes NREP, YTHDC2, and KCNN2 were found at the window at 2 Mb (BTA10). The KCNN2 gene is related to the activity in calcium/potassium channels. The Ca/K channels are present in the cytosol modulate tissue concentrations regulating the potential membrane and when present in the liver it is believed that they play a role in the response to metabolic stress [45]. Nine genes were found in the window located at the position 1.7 Mb in BTA11 (ACOXL, BUB1, TPC3, NPHP1, MALL, MRPS5, ZNF514, ZNF2 and PROM2). At the position 9 Mb in the same chromosome (BTA11), seven other candidate genes were listed: MRPS9, GPR45, TGFBRAP1, FHL2, TACR1, POLE4 and HK2. The GPR45 gene belongs to a family of receptors present in G proteins, and some studies in humans and mice have demonstrated that the GPR gene family is responsible for and is functional in the brain [46]. The SF3B3, COG4, FUK, ST3GAL2, and AARS genes were identified in BTA18, close to the 1700 kb region. The COG4 gene is responsible for in the Golgi apparatus, performing functions related to secretion [44]. In the same chromosome at position 60 Mb two genes were identified: ZNF677 and ZNF729. Both of them encode type Zinc Fingers binding proteins, which is characterized by the coordination and stabilization of zinc ions in several processes of ion exchange [16].
Seven candidate genes for G:F were found in BTA20 at the position 71 Mb: CEP72, SLC9A3, EXOC3, PDCD6, SDHA, CCDC127, and LRRC14B. In humans the SLC9A3 gene has the function of pH regulation, eliminating the acids produced by the metabolism and has proton antiporter activity, and solute carrier family 9 is involved in the exchange of sodium ions and protons, because comprises Na + /H + exchanger proteins [47]. The window located at 5.6 Mb (BTA21) presented five candidate genes GABRG3, VIMP, CHSY1, LRRK1, and ALDH1A3. BTA22 presented six associated candidate genes FRMD4B, LMOD3, ARL6IP5, UBA3, TMF1, and EOGT. In BTA27 the associated genomic region did not show any gene because it is an intergenic region. The results of gene enrichment and functional analysis reported clusters related to nucleotides, Golgi apparatus, protein transport and acetylation (Table C in S1 Tables). The metabolic pathways linked to carbohydrate metabolism, such as six-carbon polysaccharides fructose and mannose. These polysaccharides are intermediates in the glucose degradation process and glycolytic pathway, which are the main energy source for all cell types from mammals, being responsible for the ATP supply in aerobic and anaerobic conditions [48] (Table 8).  A total of 12 SNP windows distributed in BTA1, BTA4, BTA7, BTA8, BTA10, BTA18, BTA21, and BTA24 are responsible for more than 1% of the additive variance for RFI (Fig 4). Mujibi et al. [49], in a study of crossbreed beef cattle, reported 11 SNP windows associated with RFI. Santana et al. [16] found two SNP windows (located in BTA8 and BTA21) associated with RFI in a Nellore cattle population using two different densities of SNP markers.
Many candidate genes identified in regions associated with RFI are coding proteins involved in the ion transportation system. This system can consume more than 10% of the total energy used by ruminants. Thus, animals with a reduced energy expenditure on this of this system can redirect energy to spend on other processes, which may influence feed efficiency [16]. In the present study, the window located at 121 Mb (BTA1) is responsible for 5.67% of the additive genetic variance for RFI. This window harbored the ZIC1 e ZIC4 genes, which are related to ion transport, since codified protein subunits which depend on zinc ions for stability. Other genes are also related to ion transportation system, such as GPR98, KCNV2 (acting on potassium channels) and ZDHHC7.
A total of eight candidate genes was found in BTA4, of which three of them are linked to the perception of food: TAS2R3 and TAS2R4, which are related to taste through taste receptors, and OR9A, which encodes odor receptors. In another region of BTA4 (118 Mb) six candidates genes were found. The LMBR1 gene, located at 118 Mb (BTA4), has the function to encode membrane proteins. In addition, other candidate genes encodes proteins with biological functions associated with many cellular processes such as transcription (POLR3G and TBPL2), cellular secretion (KTN1 and GOLIM4), and transport across membranes (TMEM178B and  Table 9). The results of enrichment pathway analysis indicated that genes associated with RFI (p <0.05) are related to sensory receptors that operate in food taste perception and receptors coupled to G proteins (Table 10). The enrichment by functional cluster for RFI showed that those genes are linked to cell membranes that are related to the thermodynamic equilibrium of cells. (Table D in S1 Tables). The most common mechanism to maintain the thermodynamic equilibrium of cells is ion exchange activity, maintaining cell differentiation potential through the sodium-potassium pump. For adequate pump operation, in order to maintain thermodynamic equilibrium, energy is needed (ATP to keep the sodium-potassium pump equilibrium). This process requires approximately 25% of the basal energy expenditure of an individual [50].
The GOLIM4 and SERPINI1 genes, located in BTA1, have been linked to more than one trait, such as G:F and RFI. It was observed that few SNP windows located nearby in the same chromosome were associated with more than one trait like in BTA8, where two SNP windows were associated with RFI and G:F. Seven genes related the zinc finger protein domain located in BTA18 were found in several nearby SNP windows associated with G:F and DMI. These results could be due to pleitropic effects, which means that the expression of different traits could be influenced by the same set of genes which acts in a coordinated manner to contribute to feed efficiency.
It is important to highlight that the results obtained in this study are also supported by previous studies. Rolf et al. [14] working with an Angus population, reported genomic regions associated with DMI, RFI and ADG close to those obtained in this study for the same traits. Bolormaa et al. [11] reported seven genomic regions near (at a maximum distance of 3 Mb) those identified in this study for RFI. Recently, Karisa et al. [51] reported a candidate gene (CYP2B) associated with RFI in BTA18 (49 Mb), at the same position that the present study identified candidate genes for ADG. Recently, in a study with a Nellore cattle population, De Oliveira et al. [52] reported differently located candidate genes for feed efficiency traits than found in this study (HRH4, ALDH7A1, APOA2, LIN7C, CXADR, ADAM12 and MAP7). However, the genes described in that study have similar gene ontology (immune system, energy and ion metabolism) to the genes reported in the present study.
The large number of genomic regions associated with feed efficiency traits obtained in this study should support a better understanding the genetic and physiological mechanisms that determine growth, feed intake and feed efficiency in zebu animals. The results demonstrate the probability that these traits have their expression controlled by many QTL with small individual effects, confirming their polygenic nature. The identification of relevant genes might be a difficult task, since the additive genetic variance contribution from each region or SNP window for many traits was lower than expected. Thus, strategies such as genomic selection that take into consideration the variability among all markers might be a more adequate alternative to improve these traits.
In recent years there has been growing concern about the contribution of the beef industry to climate change. Livestock in particular has been identified as a major contributor to global warming according to the FAO report [53]. Climate change and variability impacts on livestock productivity, and especially on the economic and political behavior of international markets, constituting the main current threat to beef exports in many countries. The ability to become a supplier country of reliable and safe food to the world should be combined with environmental sustainability. The results obtained in this study show that it is possible to improve beef cattle feed efficiency, trough selection using genomic information, and reach additional  benefits to the environment by reducing greenhouse gas emissions. Finally, the data set used in this study belongs to a beef cattle research farm, which provides selected sires to commercial herds in many regions of Brazil. Therefore, the information found in this study will contribute to the selection for animals with better feed efficiency and increased environmental and social sustainability.

Conclusion
The results of this study pointed out that selection for feed efficiency indicator traits is feasible in Nellore cattle under tropical conditions. It is expected that the response to selection would be higher for RFI than for G:F. Several genomic regions harboring possible QTL for feed efficiency indicator traits were identified. The candidate genes identified are involved in energy and protein metabolism, ion transport, transmembrane transport, the olfactory system, the immune system, secretion (Golgi apparatus) and cellular activity (cell multiplication). The identification of these regions and their respective candidate genes should contribute to the formation of a genetic basis for Nellore feed efficiency indicator traits, and these results would support the selection for these traits.
Supporting Information S1 Tables. Table A-Functional cluster enrichment analysis via DAVID database for DMI.