Identification of QTL regions and candidate genes for growth and feed efficiency in broilers

Background Feed accounts for about 70% of the total cost of poultry meat production. Residual feed intake (RFI) has become the preferred measure of feed efficiency because it is phenotypically independent of growth rate and body weight. In this study, our aim was to estimate genetic parameters and identify quantitative trait loci (QTL) for feed efficiency in 3314 purebred broilers using a genome-wide association study. Broilers were genotyped using a custom 55 K single nucleotide polymorphism (SNP) array. Results Estimates of genomic heritability for seven growth and feed efficiency traits, including body weight at 28 days of age (BW28), BW42, average daily feed intake (ADFI), RFI, and RFI adjusted for weight of abdominal fat (RFIa), ranged from 0.12 to 0.26. Eleven genome-wide significant SNPs and 15 suggestively significant SNPs were detected, of which 19 clustered around two genomic regions. A region on chromosome 16 (2.34–2.66 Mb) was associated with both BW28 and BW42, and the most significant SNP in this region, AX_101003762, accounted for 7.6% of the genetic variance of BW28. The other region, on chromosome 1 (91.27–92.43 Mb) was associated with RFI and ADFI, and contains the NSUN3 and EPHA6 as candidate genes. The most significant SNP in this region, AX_172588157, accounted for 4.4% of the genetic variance of RFI. In addition, a genomic region containing the gene AGK on chromosome 1 was found to be associated with RFIa. The NSUN3 and AGK genes were found to be differentially expressed in breast muscle, thigh muscle, and abdominal fat between male broilers with high and low RFI. Conclusions We identified QTL regions for BW28 and BW42 (spanning 0.32 Mb) and RFI (spanning 1.16 Mb). The NSUN3, EPHA6, and AGK were identified as the most likely candidate genes for these QTL. These genes are involved in mitochondrial function and behavioral regulation. These results contribute to the identification of candidate genes and variants for growth and feed efficiency in poultry.


Background
Feed efficiency is the most important economic trait in poultry and livestock industries because up to 70% of the total production cost is due to feed [1]. Improvement of feed efficiency increases profitability for producers and reduces the environmental footprint of production. Feed efficiency is generally defined as the ability of an animal to convert feed nutrients into body mass or other useful products. The traditional measure of feed efficiency, feed conversion ratio (FCR), has been widely used in animal breeding programs. However, selection for FCR has not necessarily resulted in an improvement in the efficiency of feed utilization due to its strong correlation with body weight and growth rate [2]. An alternative measure of feed efficiency that has gained popularity in farm animal production is residual feed intake (RFI), which was first used by Koch et al. [3] for beef cattle. RFI is generally defined as the difference between the actual and the predicted feed intake during the measurement period. The predicted feed intake is based on the expected or average requirement of an animal for the maintenance of body weight and production [4] and is typically derived as the regression of ADFI on metabolic body weight and average daily gain (ADG). Compared with FCR, one of the main advantages of selection on RFI is reduction of feed intake without jeopardizing production traits, such as body weight and growth rate [5].
RFI-efficient animals tend to have leaner carcasses with lower subcutaneous fat thickness [6,7]. Inclusion of backfat thickness in the model to account for difference in fat versus lean deposition when calculating RFI has been studied in cattle [8] and pigs [9]. In the current study, RFI adjusted for abdominal fat weight (RFIa) was calculated as the difference between the observed and the predicted ADFI using multiple regression of ADFI on metabolic body weight at mid-test (MWT), average daily gain (ADG), and abdominal fat (AbF) weight.
With the advent of high-density SNP genotyping arrays, genome-wide association studies (GWAS) have become a powerful tool for the detection of QTL in farm animals. GWAS for RFI have been performed in cattle [13][14][15], pigs [16], layers [17,18], and meat-type chickens [12]. However, to date only a limited number of QTL for RFI have been identified in livestock. One reason is that RFI appears to be regulated by many genes, each with a small effect, and thus a relatively large sample size is required for the detection of these QTL [19].
Compared with cattle and pigs, generation of data on complex traits in a large population of chickens is relatively easy because of operational convenience, low unit cost of animals, short duration of the study, etc. The objectives of this study were to estimate genetic parameters, and identify significant loci and genes that affect feed efficiency traits in 3314 broilers, based on genotyping with a customized 55 K chicken single nucleotide polymorphism (SNP) array [20].

Experimental birds
The chickens used in this study were obtained from the pure line B of fast-growing chickens that has been selected for seven generations for increased body weight and growth rate traits by Foshan Gaoming Xinguang Agricultural and Animal Industrials Co., Ltd. (Foshan, China). The statistics of selection pressure in male and female selection candidates for generations 2 to 7 are in Additional file 1: Table S1. Body weight and feed intake were measured during the growth phase (28 to 41 days of age) on 2000 male and 1365 female chickens from generations 5 to 7 across 11 hatches. Among these, 1137 males and all females were selected according to hatch number and slaughtered at 42 days of age to obtain phenotypes on carcass traits, including abdominal fat (AbF) weight (see Additional file 2: Table S2). During the test period from 28 to 41 days of age, broilers were housed in individual cages with free access to feed and water. The corn-soybean meal diet contained 12.35 MJ/kg metabolic energy and 178 g/kg crude protein. Blood samples for DNA were obtained at 40 days of age via wing vein punctures using citrated syringes during a routine health inspection.

Phenotypes
The eight traits that were either measured or calculated consisted of two growth traits, five feed intake and efficiency traits, and AbF weight. Growth traits were body weight at 28 days (d) of age (BW28) and BW42. Feed intake and efficiency traits were average daily feed intake (ADFI), residual feed intake (RFI), residual feed intake adjusted for weight of abdominal fat (RFIa), average daily gain (ADG), and feed conversion ratio (FCR). Total feed intake for each broiler for the total test period was used to calculate ADFI. Average daily gain (ADG) was calculated from BW28 and BW42. FCR was obtained from total feed intake divided by total weight gain from 28 to 42 days of age. Metabolic weight at mid-test (MWT) was calculated for each bird as the average of BW28 and BW42 (MBW) raised to the power of 0.75 (MBW 0.75 ). After slaughter at 42 days of age, abdominal fat was removed and weighed for each individual to obtain AbF. RFI was estimated as the difference between the observed and the predicted ADFI with (RFIa) or without accounting for AbF. RFI and RFIa were derived as the residuals from the following two models: where µ is the intercept, hatch and sex are fixed effects, MWT, ADG, and AbF are as defined above, β 1 , β 2 , and β 3 are partial regression coefficients, and e is the residual.
Quality control was applied to the phenotypes of all traits and data on 35 birds with phenotypes that deviated by more than three standard deviations from the mean were removed.

Genotyping, imputation, and quality control
Genomic DNA was extracted from blood samples using the phenol-chloroform method. Genotyping was conducted with a custom chicken 55 K SNP array (Beijing Compass Biotechnology Co., Ltd., Beijing, China), which is designed based on the Gallus_gallus-5.0 assembly and includes 52,060 SNPs [20]. The physical positions of these SNPs were liftovered to the GRCg6a assembly using the UCSC liftOver tool (20 May 2019). For quality control of the genotypic data, we used the PLINK (version 1.9) software [21]. Sixteen birds with a sample call rate lower than 90% were excluded. In total, 10,809 SNPs were removed with a call rate lower than 90% or a minor allelic frequency (MAF) lower than 0.05. Eighty-seven SNPs that were not assigned to chromosome or linkage group were also excluded. Missing alleles were imputed using the Beagle 5.0 software [22]. Finally, 41,164 SNPs and 1972 males and 1342 females passed the quality control criteria (see Additional file 3 Table S3).

Estimation of heritability and genetic correlations
Univariate and bivariate animal models were fitted by restricted maximum likelihood (REML), using the ASReml v4.1 software package [23]. In preliminary analyses, the fixed effects of hatch and sex were significant (P < 0.01) based on Wald F statistics and were, therefore, included in the subsequent analyses. The following model with maternal genetic and maternal environmental effects was used: where y is the vector of observations, b is the vector of fixed effects, including hatch and sex, a is the vector of random direct additive genetic effects, m is the vector of random maternal additive genetic effects, c is the vector of random common maternal environmental effects, e is the vector of random residual effects, and X , Z 1 , Z 2 , and Z 3 are incidence matrices for b , a , m , and c , respectively. The variance-covariance structure assumed for random effects was: where σ 2 a , σ 2 m , σ 2 c , and σ 2 e are the direct additive genetic, maternal additive genetic, common maternal environmental, and residual error variances, respectively, σ am is the covariance between direct and maternal genetic effects, H is a relationship matrix that combines genomic and pedigree relationships because 1304 dams (including 675 dams from generation 4) had only pedigree information, and I is an identity matrix. The covariance between direct and maternal genetic effects was assumed to be zero when convergence problems were encountered. Matrix H was constructed as in Legarra et al. [24]: where A 11 , A 12 , A 21 , and A 22 are submatrices extracted from the pedigree-based relationship matrix A, with indices 1 for non-genotyped dams and 2 for genotyped animals, and G is the genomic relationship matrix based on the SNPs following the first method of VanRaden [25]. The same model was also used for pedigree-based analyses by replacing H by the pedigree-based relationship matrix A . The pedigree used in the analysis consisted of 4865 birds from generations 3 to 7.

Genome-wide association study
A GWAS was performed for each trait using a univariate linear mixed model implemented in GEMMA version 0.98.1 software (https ://githu b.com/genet ics-stati stics /GEMMA /relea ses) [26], one SNP at a time. The fixed effects of hatch and sex were included in the model for all traits, except for RFI and RFIa. The hatch and sex were included as fixed effects when the phenotypic values of RFI and RFIa were estimated. The fixed effects were used to construct a design matrix using the model.matrix() function in R (version 3.6.0) and included in the GEMMA software similar to a covariate. The statistical model used was: where y is a vector of phenotypic values, X is the design matrix for fixed effects, including a column of 1 s, b is the vector of the corresponding effect estimates, including the intercept, S is a vector of genotypes coded major/ minor alleles as 0/1 for a given SNP, β is the allele substitution effect for the fitted SNP, u is a vector of random animal genetic effects, e is a vector of residuals, σ 2 a is the direct additive genetic variance, σ 2 e is the residual error variance, G is the genomic relationship matrix, I is an identity matrix, and MVN is the multivariate normal distribution.
The Wald test was used as a criterion to identify SNPs that were significantly associated with the investigated traits. Genome-wide significance was assessed using the simpleM method [27]. This method first derives the composite linkage disequilibrium (LD) correlation matrix from the SNP genotypes, then calculates the eigenvalues of the component LD matrix by principal component analysis, and finally infers the effective number of independent tests as the number of principal components that jointly contribute 99.5% of the variation in SNP genotypes. Manhattan and quantile-quantile (Q-Q) plots were constructed for each trait by the qqman package (http://cran.r-proje ct.org/packa ge=qqman ) in R (version 3.6.0). The genomic inflation factor (GIF) was calculated by the GenABEL R package [28]. The percentage of genetic variance that was explained by each significant SNP was calculated following Al-Mamun et al. [29] as: where p i and q i are the allele frequencies for SNP i , β i is the estimate of the allele substitution effect of SNP i based on the GWAS model, and σ 2 g is the estimate of genetic variance for the trait estimated as described above using the combined genomic-pedigree relationship matrix H.

Further analysis of QTL regions
Significant and adjacent SNPs for the analyzed traits were further investigated by region-based association analyses using the fast family-based sequence kernel association test (FFBSKAT) program (SNP effects fitted as random effects) implemented in the R package FREGAT [30]. Region-based association analysis is an efficient approach to identify causal SNPs. The fixed effects of hatch and sex were added to the FFBSKAT, except for RFI and RFIa. The linear mixed model was given following Schifano et al. [31]: y = Xb + h + u + e , where h is the vector of random SNP effects and the other parameters are the same as for the GEMMA model. The LD among SNPs in these regions was estimated using PLINK (version 1.9) [21]. Boxplots of the phenotype distribution by genotype of SNPs were produced by the ggplot2 package in R (version 3.6.0).
Additive and dominance effects of SNPs on the analyzed traits were estimated using the ASReml v4.1 software package [23]. To estimate the additive effect, a covariate was fitted with values 1, 0, and -1 for the homozygous genotype for the major allele, the heterozygous genotype, and the homozygous genotype for the minor allele, respectively. To estimate the dominance effect, a covariate was fitted with value 0 for the homozygous genotypes and 1 for the heterozygous genotype. The model was: y = Xb + Qq + u + e , where Q is the design matrix for a SNP (additive and dominance) effects, q is the additive and dominance effects for the fitted SNP, other components are as defined for the GEMMA model.

Analysis of the expression of candidate genes by qPCR
The genes that were located closest to genome-wide significant and suggestive SNPs were identified based on the UCSC annotation of the GRCg6a genome version (http:// genom e-asia.ucsc.edu/cgi-bin/hgGat eway?hgsid =47276 8848_otkBt CHKhH MTV1x rxHuq 737ii vJ1). Expression of these candidate genes (NSUN3, EPHA6, and AGK) was assessed in the relevant tissues using qPCR analyses. One hundred and seventy-five males from generation 6 with growth and feed efficiency between 28 and 41 d of age were slaughtered at 42 d of age to measure carcass traits. Tissue samples from birds with the highest (n = 15) and lowest (n = 15) RFI phenotypes were collected (Table 1).
Total RNA was isolated from liver, breast muscle, thigh muscle, and AbF tissues using the TRIzol reagent (Life Technologies, Inc., Carlsbad, CA). First-strand cDNA was synthesized from 2 µg total RNA using the Fast-Quant RT Kit (with gDNAase) (TIANGEN BIOTECH, Beijing, China). Power SYBR ® Fast qPCR Master Mix (KAPA, Wilmington, MA) was used to analyze mRNA expression of the selected genes. qPCR analysis was performed with an ABI Q7 Flex Real-time detection system (Applied Biosystems, Foster City, CA). Specific primers [see Additional file 4 Table S4] were designed using the Primer Premier 6.0 software based on chicken coding sequences and were subsequently synthesized by the Beijing Genomics Institute (Beijing, China). The desired PCR product size was set to 100 to 250 bp and the best primer pair among the output was selected.

Descriptive statistics of growth and feed efficiency traits
The descriptive statistics in Table 2 show the distribution of growth and feed efficiency traits of the 3314 broilers (1972 males and 1342 females) used in the GWAS. Phenotypes for individual broilers [see Additional file 5: Table S5] and descriptive statistics for male and female broilers are shown separately [see Additional file 6: Table S6]. The average body weight at 28 days of age (BW28) was 1.08 kg and reached 2.21 kg at 42 days of age (BW42). Chickens consumed an average of 151 g/ days of feed and gained 80.6 g/days in weight during the growth phase (28 to 41 days of age). Values of RFI and RFIa ranged from -19.3 to 19.0 g/days and from -18.0 to 18.1 g/days, respectively. The coefficients of variation for the five growth and feed efficiency traits ranged from 7.8 to 19.5%.

Genetic parameters of growth and feed efficiency traits
Estimates of variance components and heritabilities using genomic and pedigree relationship matrices are in Table 3. Estimates of genomic heritabilities for two growth traits (BW28 and BW42) and the five feed efficiency traits (ADFI, RFI, RFIa, ADG, and FCR) were low to moderate, ranging from 0.12 to 0.26. Genetic correlations between direct and maternal effects for BW28, BW42, ADFI, RFI, RFIa, and AbF were moderate to high, ranging from -0.77 to -0.45. Estimates of heritabilities and genetic correlations between direct and maternal effects based on pedigree information were higher than genomic estimates. Estimates of maternal genetic and maternal environmental effects on growth and feed efficiency traits were low (0.00 to 0.09). Estimates of variance components and heritabilities based on genomic and pedigree relationship matrices are shown separately for male and female broilers in Additional file 7: Table S7 and Additional file 8: Table S8. Estimates of additive genetic variances and heritabilities of growth and feed efficiency traits were lower for males than for females.
Estimates of genetic ( r g ) and phenotypic ( r p ) correlations using the combined genomic-pedigree relationship matrix are in Table 4. Estimates based on the pedigree relationship matrix were very similar ( r = 0.97, [see Additional file 9 Table S9]). Estimates of genetic correlations

Genome-wide association study of growth and feed efficiency traits
The number of independent tests was estimated to be 26,217. Based on this, genome-wide and suggestive significance thresholds were set equal to 1.91*10 -6 (0.05/26,217) and 3.81*10 -5 (1/26,217), respectively.

Growth traits
Manhattan and quantile-quantile (Q-Q) plots for BW28 and BW42 are shown in Fig. 1. Detailed information on SNPs significantly associated with BW28 and BW42 is in Table 5 and sequence information regarding these SNPs is in Additional file 10: Table S10. For BW28, seven genome-wide significant SNPs were identified on GGA16, and five suggestively significant SNPs were detected on GGA4 and 16. The genomic inflation factor (GIF) was 1.02 for BW28, which suggests that the population structure was well controlled. Eleven of the genome-wide significant SNPs were within a 316 kb region on GGA16 (2.34-2.66 Mb). The minor allele of the most significant SNP in this region, AX_101003762, had a negative effect estimate ( β < 0 ) and accounted for 7.6% of the genetic variance for BW28. These SNPs on GGA16 were located either within genes or 0.03 to 47.43 kb away from the nearest genes, which include zinc finger protein 692 (ZNF692), tripartite motif containing 39.2 (TRIM39.2), and tripartite motif containing 27.2 (TRIM27.2).
For BW42, three genome-wide significant SNPs and five suggestively significant SNPs were located on GGA1, 4, 10, and 16. The GIF was 1.00, which indicates that this association analysis was hardly affected, if at all, by population stratification. Three genome-wide significant SNPs and two suggestively significant SNPs were located within a 206.0 kb region (GGA16: 2.34-2.55 Mb) and four SNPs were the same as those found for BW28. The leading SNP for BW42, AX_172583407, accounted for 5.0% of the genetic variance.

Feed efficiency traits
Manhattan and quantile-quantile (Q-Q) plots for the three feed efficiency traits are shown in Fig. 2. Detailed information on SNPs significantly associated with the feed efficiency traits is in Table 6 and sequence information regarding these SNPs is in Additional file 10: Table S10.
For ADFI, one SNP exceeded the suggestive significance threshold (GGA1: SNP AX_75546765 at 91.27 Mb). This SNP was located near the NOP2/Sun RNA methyltransferase family member 3 (NSUN3) gene and accounted for 3.4% of the genetic variance for ADFI.
For RFI, one genome-wide significant SNP and three suggestively significant SNPs were detected, two of which were located within a 457.9 kb region on GGA1 (91.97-92.43 Mb). This region was located about 700.6 kb downstream of the SNP AX_75546765 that was associated with ADFI. These SNPs were located within or near the EPH receptor A6 (EPHA6), NSF attachment protein gamma (NAPG), and leucine zipper and CTNNBIP1

Table 4 Estimates of genetic and phenotypic correlations between growth and feed efficiency traits using the relationship matrix that blends genomic and pedigree information
a Upper diagonal is genetic correlation, and lower diagonal is phenotypic correlation BW28, body weight at 28 d of age; BW42, body weight at 42 d of age; ADFI, average daily feed intake; RFI, residual feed intake; RFIa, residual feed intake adjusted for weight of abdominal fat; ADG, average daily gain; FCR, feed conversion ratio; AbF, weight of abdominal fat   2 Manhattan and quantile-quantile (Q-Q) plots of the GWAS for feed efficiency traits. Each dot is a SNP in the dataset. The horizontal red and blue lines indicate the genome-wide significant (P value = 1.90e−6) and suggestive thresholds (P value = 3.80e−5), respectively. ADFI, RFI, and RFIa are average daily feed intake, residual feed intake, and residual feed intake adjusted for weight of abdominal fat, respectively, and GIF is the genomic inflation factor For RFIa, one genome-wide suggestively significant SNP (GGA1: 57.28 Mb) was located in the first intron of the acylglycerol kinase (AGK) gene. This SNP accounted for 6.3% of the genetic variance.

Region-based association test, linkage disequilibrium, and allele frequency analysis
The two region-based association plots for GGA16 and GGA1, which included multiple genes in the associated LD region, are shown in Figs. 3a-b and 4a- than those that were homozygous TT and CC (P < 0.05), respectively. ADFI and RFI were significantly lower for individuals that were homozygous CC for AX_75546765 and CC for AX_172588157 than those that were homozygous AA and TT (P < 0.01), respectively. The homozygous AA individuals for AX_172566874 had lower RFIa, and higher BW42 than those with homozygous GG (P < 0.05). In addition, these five SNPs had significant additive effects for multiple traits (P < 0.05) and did not exhibit significant dominance effects (P > 0.05, [see Additional file 11: Table S11]). The effects of the five leading SNPs also differed between male and female broilers [see Additional file 12: Table S12]. Estimates of the additive effects of AX_101003762 on BW28, of AX_75546765 on ADFI, and of AX_172566874 on RFIa were larger for females than for males, while those of AX_172583407 on BW42 and of AX_172588157 on RFI were smaller for females than for males. Except for SNP AX_172583407, the dominance effects of these SNPs were larger for females than for males. These results suggest that the five leading SNPs have additive and pleiotropic effects on growth and feed efficiency traits that are sex-specific.
The evolution in the frequency of the favorable allele of the most significant SNPs over three generations (generation 5, 6, and 7) is shown in Table 8. The frequency of the favorable C allele of SNP AX_101003762, which was associated with BW28, increased from 0.71 to 0.78 under selection pressure for increased body weight and growth rate. The frequency of the favorable T allele of SNP AX_172583407, which was associated with BW42, increased from 0.72 to 0.78. The frequency of the favorable A allele of SNP AX_172566874, which was associated with RFIa, increased from 0.43 to 0.49. The allele frequencies of the other two SNPs were only slightly different between generations.

Expression of candidate genes in high and low RFI males
Expression of the three nearest genes (NSUN3, EPHA6, and AGK) within the QTL regions for ADFI, RFI, and RFIa were evaluated by qPCR analysis in males with high and low RFI phenotypes (HRFI and LRFI, Fig. 5). In breast and thigh muscle, the relative expressions of NSUN3 and AGK were significantly higher in HRFI than in LRFI males (P < 0.01). In abdominal fat (AbF) tissue, a b c d e f Fig. 4 Association results for the candidate region on GGA1 for ADFI, RFI, and RFIa. a Regional plot of the candidate region on GGA1 (91.15-92.61 Mb) for ADFI. b Regional plot of the candidate region on GGA1 (91.15-92.61 Mb) for RFI. c Regional plot of the candidate region on GGA1 (56.78-57.82 Mb) for RFIa. In the upper panels, the leading SNPs are highlighted by blue solid circles and those near or within the gene by red color. Different levels of linkage disequilibrium between the leading SNP and surrounding SNPs are indicated in different colors. P-values are based on analyses in FFBSKAT. d Boxplot for ADFI and genotype at SNP AX_75546765. e Boxplot for RFI and genotype at SNP AX_172588157. f Boxplot for RFIa and genotype at SNP AX_172566874. ADFI, RFI, and RFIa are average daily feed intake, residual feed intake, and residual feed intake adjusted for weight of abdominal fat, respectively the expressions of NSUN3 and AGK were significantly lower in HRFI than in LRFI individuals (P < 0.01). No significant differences between the HRFI and LRFI cohorts were found for the expressions of NSUN3 and AGK in liver. There was no detectable expression of EPHA6 in the liver, breast muscle, thigh muscle, or AbF tissues.

Estimates of genetic parameters
The genomic heritability estimates for growth and feed efficiency traits were low to moderate (0.12 to 0.26) and lower than those reported by Bernon and Chambers [37] and Xu et al. [12], who reported pedigree-based estimates that ranged from 0.22 to 0.44 in fast-growing broilers from 28 to 42 days of age, and from 0.22 to 0.56 in medium-growing purebred broilers from 44 to 83 days of age. However, our estimates are consistent with those of Abdollahi-Arpanahi et al. [38], who reported a genomic heritability estimate of 0.23 for BW35 in broilers. Yuan et al. [18] reported low genomic heritability estimates for ADFI (0.15), RFI (0.17), and FCR (0.21) in layer chickens. Heritability estimates for BW28, BW42, ADFI, RFI, RFIa, ADG, and AbF based on the combined genomicpedigree relationship matrix were lower than those based on the pedigree relationship matrix. Similar findings were previously reported for pigs [39] and dairy cattle [40,41]. Several factors can contribute to these differences. First, due to ascertainment bias in the SNP array, particular types of variants, such as rare variants, variants with low MAF, copy number variants, and structural variants may be absent, which can reduce genomic heritabilities. Second, genomic heritabilities and pedigree-based heritabilities apply to different base populations. The pedigree relationship matrix is based on the identity-by-descent (IBD) with respect to a base population consisting of the founders of the pedigree. The genomic relationship matrix is constructed based on allele frequencies in the  current population rather than those in the base population because genotypes for previous generations are unknown [25]. The use of the current population results in smaller estimates of additive genetic variance since the current population is expected to be more inbred than the base population [42]. Third, pedigree-based analyses resulted in lower estimates for maternal genetic and maternal environmental variances than analyses using the combined genomic-pedigree relationship matrix for several traits, which may affect estimates of the additive genetic variance [43] and lead to lower genomic heritability estimates. Estimates of genetic correlations between the eight traits were not significantly different based on these two relationship matrices, which is consistent with the result of Abdollahi-Arpanahi et al. [44] in broilers. The estimate of the genetic correlation between RFI and RFIa was strong (0.83). Similar results for RFI and RFI adjusted for backfat thickness (RFIb) in young beef bulls have been reported by Schenkel et al. [45]. The estimate of heritability was slightly higher for RFI than for RFIa, which is consistent with the results of Ceacero et al. [46], who found slightly higher estimates of heritability for RFI (0.24) than for RFIb (0.20). Our results also indicated that RFI has a moderate genetic correlation with AbF ( r g = 0.51), whereas that of RFIa with AbF is close to zero. This suggests that selection for RFI will also result in changes in AbF, while selection for RFIa would change feed intake without affecting AbF. RFIa is usually not a desired trait for a fast-growing broiler breeding program because less abdominal fat weight is preferred. Inclusion of other traits in the model for calculating RFI could be studied in the future.
Sexual dimorphism, largely caused by differences in gene expression [47], has been investigated for economically important traits in livestock species [48]. In the current study, this dimorphism was observed for growth and feed efficiency traits, with males having significantly higher BW28, BW42, ADFI, and ADG than females. In addition, estimates of heritability for growth and feed efficiency traits were lower for males than for females, which is consistent with the results of Mebratie et al. [49] and of van der Heide et al. [48]. Interestingly, the additive effects of the five leading SNPs on growth and feed efficiency traits were also sex-biased, which suggests that the genes located near these SNPs might be differentially expressed between male and female broilers.

Genome-wide association study of growth and feed efficiency traits Loci and genes for growth traits
A 316.0-kb genomic region on GGA16 was associated with both BW28 and BW42. Previous studies have reported QTL on GGA16 for BW9, BW56, and BW84 in medium-growth broilers [11] and in Iranian indigenous chickens [50]. In the absence of other factors, positive selection will result in a subtle evolution of the favorable allele frequencies over generations. In agreement with this, the frequency of the favorable allele of the most genome-wide significant SNPs (AX_101003762 and AX_172583407) associated with BW28 and BW42 continued to increase from generation 5 to generation 7, because selection for these traits was maintained as the main emphasis. These results substantiate evidence that this region on GGA16 (2.34-2.66 Mb) is associated with BW28 and BW42.
The SNPs, AX_101003762 and AX_172583407, were found to be within the TRIM39.2 and ZNF692 genes, respectively. TRIM39.2 is a member of the RING-B-boxcoiled-coil (RBCC)/tripartite motif (TRIM) subfamily of zinc finger proteins that are involved in many biological processes, including cell differentiation [51]. Glycogenininteracting protein 1 (GNIP1) is an isoform of GNIP/ TRIM7 and a tripartite motif (TRIM) protein. In vivo overexpression of GNIP1 in mouse muscle activates the protein kinase B-glycogen synthase kinase-3-glycogen synthase cascade and subsequent glycogen synthesis in muscle, leading to decreased blood glucose levels, lactate levels and mouse body weight, without affecting wholebody insulin or glucose tolerance [52]. ZNF692 was first identified as a transcription factor associated with AMPK signaling [53]. In humans, ZNF652, a member of the zinc finger protein gene family is significantly associated with height [54]. Shirai et al. [55] suggested that ZNF692 is a key modulator of hepatic glucose production regulated by AMPK in vivo.

Loci and genes for ADFI and RFI
An important region on GGA1 (GGA1: 91.97-92.43 Mb) was associated with RFI and a suggestive significant SNP (GGA1: 91.27 Mb) was associated with ADFI. This finding indicates that the 1.16-Mb region (GGA1: 91.27-92.43 Mb) is an important QTL for feed efficiency. This region was previously identified as a feed efficiency QTL (GGA1: 90.35-123.03 Mb) in a meat-type X egg-type resource chicken population by Hansen et al. [56]. This region contains two candidate genes, namely NSUN3 and EPHA6.
The NSUN3 gene is a member of the Nol1/Nop2/Sun domain (NSUN) family, which is closely related to mitochondrial function and is localized to the mitochondrial (mt) matrix, where it interacts with mt-tRNA Met to methylate cytosine 34 (C34) at the wobble position [57]. A study on NSUN3-knockout cells showed strongly decreased mt-tRNA Met methylation and formylation, as well as a reduction in mitochondrial translation, protein synthesis, and reduced oxygen consumption, leading to deficient mitochondrial activity [58,59]. Recently, Watson et al. [60] showed that NSUN3 is significantly associated with anorexia nervosa in human. Common symptoms of anorexia nervosa are related to functional impairments of the mitochondrial oxidative phosphorylation system complex I, including poor feeding, neurodegeneration, and muscle weakness [61].
EPHA6 affects physical activity, and thus energy expenditure. EPHA6 knockout mice display behavioral deficits associated with learning and memory disabilities [62]. Dos Santos et al. [63] found that EPHA6 was significantly associated with behavioral reactivity in cattle, which is a temperament trait. In humans, EPHA6 is strongly associated with blood pressure [64]. In addition, Herd and Arthur [65] reported that physical activity accounted for 9% of the phenotypic variation for RFI in beef cattle. As previously reported, the brain plays critical roles in the regulation of feeding behavior and body energy homeostasis [1]. EPHA6 plays key roles in a variety of biological functions, such as brain development and behavioral regulation. It is highly expressed in most regions of the adult mouse brain, including the hypothalamus [66], which plays a key role in modulating feed intake through the regulation of metabolic hormones and their receptors, such as leptin, neuropeptide Y, agoutirelated protein, IGF1, and ghrelin [1,67].

Loci and genes for RFIa
Only one suggestively significant SNP (GGA1: 57.28 Mb) was identified for RFIa. This SNP was located in the first intron of the AGK gene. Previous studies found a region on GGA1 (57.0-58.0 Mb) to be significantly associated with FCR in a commercial broiler line and found another SNP, GGaluGA019865, located in the first intron of AGK, to be a leading SNP [68]. AGK, also known as multi-substrate lipid kinase (MULK), is abundantly expressed in muscle, heart, kidney, and brain [69,70]. AGK has lipid kinase-dependent and kinase-independent functions within the mitochondria. On the one hand, by acting as a lipid kinase, it phosphorylates both monoacylglycerol and diacylglycerol to form lysophosphatidic acid (LPA) and phosphatidic acid (PA), respectively [70]. LPA plays an important role in the synthesis of phospholipids and cell proliferation by transactivating the epidermal growth factor (EGF) tyrosine kinase receptor [71]. PA acts both as an essential molecule for mitochondrial ultrastructure and function, and as a second messenger that regulates