Skip to main content

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

Abstract

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.

Quantitative trait loci (QTL) for traits in animals have been studied for more than 20 years. In chickens, the first QTL associated with FCR and RFI were detected on Gallus gallus (GGA) chromosome 4 by de Koning et al. [10]. Compared to the 594 QTL recorded for FCR across species, only 135 QTL for RFI are recorded in the Animal QTL Database (QTLdb; https://www.animalgenome.org/cgi-bin/QTLdb/GG/index, accessed 07 Jan 2021). Among these, five QTL were identified for RFI in chickens, on GGA3 (51.8–64.9 Mb), GGA10 (19.1–20.3 Mb), GGA12 (16.3–17.3 Mb), GGA19 (0.0–2.0 Mb), and GGAZ (2.1–7.8 Mb) [11, 12].

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].

Methods

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 (MBW0.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:

$$\mathrm{ADFI}=\mu +\mathrm{hatch}+\mathrm{sex}+{\beta }_{1}\mathrm{MWT}+{\beta }_{2}\mathrm{ADG}+{e}_{1},$$
$$\mathrm{ADFI}=\mu +\mathrm{hatch}+\mathrm{sex}+{\beta }_{1}\mathrm{MW}T+{\beta }_{2}\mathrm{ADG}+{\beta }_{3}\mathrm{AbF}+{e}_{2},$$

where \(\mu\) is the intercept, hatch and sex are fixed effects, MWT, ADG, and AbF are as defined above, \({\beta }_{1}\), \({\beta }_{2}\), and \({\beta }_{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:

$$\mathbf{y}=\mathbf{X}\mathbf{b}+{\mathbf{Z}}_{1}\mathbf{a}+{\mathbf{Z}}_{2}\mathbf{m}+{\mathbf{Z}}_{3}\mathbf{c}+\mathbf{e},$$

where \(\mathbf{y}\) is the vector of observations, \(\mathbf{b}\) is the vector of fixed effects, including hatch and sex, \(\mathbf{a}\) is the vector of random direct additive genetic effects, \(\mathbf{m}\) is the vector of random maternal additive genetic effects, \(\mathbf{c}\) is the vector of random common maternal environmental effects, \(\mathbf{e}\) is the vector of random residual effects, and \(\mathbf{X}\), \({\mathbf{Z}}_{1}\), \({\mathbf{Z}}_{2}\), and \({\mathbf{Z}}_{3}\) are incidence matrices for \(\mathbf{b}\), \(\mathbf{a}\), \(\mathbf{m}\), and \(\mathbf{c}\), respectively. The variance–covariance structure assumed for random effects was:

$$\mathrm{Var}\left[\begin{array}{c}\mathbf{a}\\ \mathbf{m}\\ \mathbf{c}\\ \mathbf{e}\end{array}\right]=\left[\begin{array}{cccc}\mathbf{H}{\upsigma }_{\mathrm{a}}^{2}& \mathbf{H}{\upsigma }_{\mathrm{am}}& 0& 0\\ \mathbf{H}{\upsigma }_{\mathrm{am}}& \mathbf{H}{\upsigma }_{\mathrm{m}}^{2}& 0& 0\\ 0& 0& \mathbf{I}{\upsigma }_{\mathrm{c}}^{2}& 0\\ 0& 0& 0& \mathbf{I}{\upsigma }_{\mathrm{e}}^{2}\end{array}\right],$$

where \({\upsigma }_{\mathrm{a}}^{2}\), \({\upsigma }_{\mathrm{m}}^{2}\), \({\upsigma }_{\mathrm{c}}^{2}\), and \({\upsigma }_{\mathrm{e}}^{2}\) are the direct additive genetic, maternal additive genetic, common maternal environmental, and residual error variances, respectively, \({\upsigma }_{\mathrm{am}}\) is the covariance between direct and maternal genetic effects, \(\mathbf{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 \(\mathbf{I}\) is an identity matrix. The covariance between direct and maternal genetic effects was assumed to be zero when convergence problems were encountered. Matrix \(\mathbf{H}\) was constructed as in Legarra et al. [24]:

$$\mathbf{H}=\left[\begin{array}{cc}{\mathbf{A}}_{11}+{\mathbf{A}}_{12}{\mathbf{A}}_{22}^{-1}(\mathbf{G}-{\mathbf{A}}_{22}){\mathbf{A}}_{22}^{-1}{\mathbf{A}}_{21}& {\mathbf{A}}_{12}{\mathbf{A}}_{22}^{-1}\mathbf{G}\\ \mathbf{G}{\mathbf{A}}_{22}^{-1}{\mathbf{A}}_{21}& \mathbf{G}\end{array}\right],$$

where \({\mathbf{A}}_{11}\), \({\mathbf{A}}_{12}\), \({\mathbf{A}}_{21}\), and \({\mathbf{A}}_{22}\) are submatrices extracted from the pedigree-based relationship matrix \(\mathbf{A}\), with indices 1 for non-genotyped dams and 2 for genotyped animals, and \(\mathbf{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 \(\mathbf{H}\) by the pedigree-based relationship matrix \(\mathbf{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://github.com/genetics-statistics/GEMMA/releases) [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:

$$\mathbf{y}=\mathbf{X}\mathbf{b}+\mathbf{S}\beta +\mathbf{u}+\mathbf{e},$$
$$\mathbf{u}\sim \mathrm{MVN}\left(0, \mathbf{G}{\upsigma }_{\mathrm{a}}^{2}\right),$$
$$\mathbf{e}\sim \mathrm{MVN}\left(0,{\mathbf{I}\upsigma }_{\mathrm{e}}^{2}\right),$$

where \(\mathbf{y}\) is a vector of phenotypic values, \(\mathbf{X}\) is the design matrix for fixed effects, including a column of 1 s, \(\mathbf{b}\) is the vector of the corresponding effect estimates, including the intercept, \(\mathbf{S}\) is a vector of genotypes coded major/minor alleles as 0/1 for a given SNP, \(\beta\) is the allele substitution effect for the fitted SNP, \(\mathbf{u}\) is a vector of random animal genetic effects, \(\mathbf{e}\) is a vector of residuals, \({\upsigma }_{\mathrm{a}}^{2}\) is the direct additive genetic variance, \({\upsigma }_{\mathrm{e}}^{2}\) is the residual error variance, \(\mathbf{G}\) is the genomic relationship matrix, \(\mathbf{I}\) is an identity matrix, and \(\mathrm{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-project.org/package=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: \(\mathrm{\%}{V}_{gi}=\frac{2{p}_{i}{q}_{i}{\beta }_{i}^{2}}{{\sigma }_{g}^{2}}\times 100\), where \({p}_{i}\) and \({q}_{i}\) are the allele frequencies for SNP \(i\), \({\beta }_{i}\) is the estimate of the allele substitution effect of SNP \(i\) based on the GWAS model, and \({\sigma }_{g}^{2}\) is the estimate of genetic variance for the trait estimated as described above using the combined genomic-pedigree relationship matrix \(\mathbf{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]: \(\mathbf{y}=\mathbf{X}\mathbf{b}+\mathbf{h}+\mathbf{u}+\mathbf{e}\), where \(\mathbf{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: \(\mathbf{y}=\mathbf{X}\mathbf{b}+\mathbf{Q}\mathbf{q}+\mathbf{u}+\mathbf{e}\), where \(\mathbf{Q}\) is the design matrix for a SNP (additive and dominance) effects, \(\mathbf{q}\) is the additive and dominance effects for the fitted SNP, other components are as defined for the GEMMA model.

Allele frequencies at the most significant SNPs in three generations (generations 5, 6, and 7) were estimated using PLINK (version 1.9) [21].

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://genome-asia.ucsc.edu/cgi-bin/hgGateway?hgsid=472768848_otkBtCHKhHMTV1xrxHuq737iivJ1). 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).

Table 1 Means (± SD) for growth, feed efficiency, and carcass traits for the high and low RFI cohortsa

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 FastQuant 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.

Amplification was performed in a total volume of 10 µL, containing 5 µL of 2 × PCR Master Mix, 0.2 µL of 50 \(\times\) ROX Low Reference Dye, 1.2 µL cDNA (100 ng/µL), 0.3 µLl of each primer, and 3.0 µL ddH2O. The PCR cycle parameters were: 95 °C for 3 min, then 40 amplification cycles at 95 °C for 3 s and 60 °C for 34 s. The following housekeeping genes were used as reference genes: ribosomal protein L32 (RPL32) for breast muscle and thigh muscle [32], ubiquitin B (UBB) for liver tissue [33], and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) for AbF tissue [34]. The comparative CT method [35] was used to determine fold-changes in gene expression calculated as \({2}^{-\mathrm{\Delta \Delta CT}}\). Statistical analyses were performed using analysis of variance (general linear model) and Tukey’s test in SAS 9.4 (SAS Institute, Cary, NC) [36]. Threshold significance values were P < 0.05 (*) or P < 0.01 (**).

Results

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%.

Table 2 Descriptive statistics for the evaluated growth and feed efficiency traits

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.

Table 3 Estimates of variance components and heritabilities for growth and feed efficiency traits using genomic and pedigree relationship matrices

Estimates of genetic (\({\mathrm{r}}_{\mathrm{g}}\)) and phenotypic (\({\mathrm{r}}_{\mathrm{p}}\)) correlations using the combined genomic-pedigree relationship matrix are in Table 4. Estimates based on the pedigree relationship matrix were very similar (\(\mathrm{r}\) = 0.97, [see Additional file 9 Table S9]). Estimates of genetic correlations of ADFI with BW28, BW42, RFI, RFIa, ADG, and AbF were moderate to high, ranging from 0.45 to 0.83, and was lowest (0.26) with FCR (Table 4). Estimates of genetic correlations of BW28 and BW42 with RFI, RFIa, and FCR ranged from -0.09 to 0.15 and did not significantly differ from zero. Estimates of phenotypic correlations of ADG with RFI and RFIa were zero, and the genetic correlations were low (0.17 with RFI and 0.19 with RFIa). Estimates of genetic correlations of AbF with RFI and FCR were 0.48 and 0.37, respectively, whereas that of AbF with RFIa was 0.05. Estimates of genetic and phenotypic correlations of RFI with RFIa and FCR were high and positive (r > 0.71).

Table 4 Estimates of genetic and phenotypic correlations between growth and feed efficiency traits using the relationship matrix that blends genomic and pedigree information

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.

Fig. 1
figure1

Manhattan and quantile–quantile (Q-Q) plots of the GWAS for growth 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. BW28 and BW42 are body weight at 28 d and 42 d of age, and GIF is the genomic inflation factor

Table 5 Information on SNPs associated with growth traits

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 (\(\beta <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.

Fig. 2
figure2

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

Table 6 Information on SNPs associated with feed efficiency traits

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 domain containing (LZIC) genes. The leading variant, AX_172588157, accounted for 4.4% of the genetic variance for RFI.

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–c. For one region on GGA16 (2.34–2.66 Mb), the strongest associated SNPs were AX_101003762 for BW28, within the TRIM39.2 gene (P = 6.81e−10), and AX_172583407 for BW42, within the ZNF692 gene (P = 1.03e−8). LD analysis showed that the levels of LD were moderate for several haplotype blocks. For two regions on GGA1 (91.15–92.61 Mb and 56.78–57.82 Mb), AX_75546765, AX_172588157, and AX_172566874 were identified as leading SNPs in the vicinity of the genes NSUN3, EPHA6, and AGK, for ADFI, RFI, and RFIa, respectively. The two regions also included other genes, including TRIM27.2, receptor for activated C kinase 1 (RACK1), kinesin family member C1 (KIFC1), ADP ribosylation factor like GTPase 13B (ARL13B), syntaxin 19 (STX19), and transmembrane protein 178B (TMEM178B).

Fig. 3
figure3

Association results for the candidate region on GGA16 for BW28 and BW42. a Regional plot of the candidate region on GGA16 (2.34–2.66 Mb) for BW28. b Regional plot of the candidate region on GGA16 (2.34–2.66 Mb) for BW42. In the upper panels, the leading SNPs are highlighted by blue solid circles and those near or within the gene in 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. c Boxplot for BW28 and the genotype at SNP AX_101003762. d Boxplot for BW42 and the genotype at SNP AX_172583407. BW28 and BW42 are body weight at 28 d and 42 days of age

Fig. 4
figure4

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

Distributions of phenotypes for BW28, BW42, ADFI, RFI, and RFIa by genotype of the most significant SNPs are shown in Figs. 3c–d, 4d–f, and Table 7. These data indicate that the extreme phenotypic values corresponded to the homozygous genotypes, whereas the intermediate values corresponded to the heterozygous genotypes. The genotypic averages (and SE) of the five most significant SNPs on growth and feed efficiency traits are in Table 7. SNPs AX_101003762 and AX_172583407 are closely linked SNPs on GGA28. Broilers that were homozygous CC for AX_101003762 and TT for AX_172583407 showed significantly larger BW28, BW42, and ADFI 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.

Table 7 Means (+ SE) of growth and feed efficiency traitsa by genotype of the five most significant SNPs

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.

Table 8 Evolution of frequency of the favorable allele of the most significant SNPs in generations 5 to 7

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, 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.

Fig. 5
figure5

Differentially-expressed genes between the high and low RFI cohorts. Data are expressed as the mean ± SD (n = 15). **P < 0.01

Discussion

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 genomic-pedigree 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 (\({\mathrm{r}}_{\mathrm{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-box-coiled-coil (RBCC)/tripartite motif (TRIM) subfamily of zinc finger proteins that are involved in many biological processes, including cell differentiation [51]. Glycogenin-interacting 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 whole-body 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.9792.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.2792.43 Mb) is an important QTL for feed efficiency. This region was previously identified as a feed efficiency QTL (GGA1: 90.35123.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-tRNAMet to methylate cytosine 34 (C34) at the wobble position [57]. A study on NSUN3-knockout cells showed strongly decreased mt-tRNAMet 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, agouti-related 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 numerous biological processes, including activation of the mammalian target of the rapamycin (mTOR) signaling pathway [72, 73]. On the other hand, AGK is a subunit of the inner membrane 22 (TIM22) complex, which imports and assembles mitochondrial carrier proteins [74, 75]. Inactivation of AGK leads to the reduction of mitochondrial maximal respiration and tricarboxylic acid (TCA) cycle flux in cells [74]. In summary, the NSUN3 and AGK genes play an important role in the regulation of the mitochondrial function.

Mitochondrial function has previously been associated with feed efficiency in poultry and livestock, because mitochondria are responsible for producing approximately 90% of the energy of the cells [76]. In our study, the lower expression of NSUN3 and AGK in breast and thigh muscle of LRFI chickens may result from decreased mitochondrial energy metabolism, which is consistent with previous observations in the skeletal muscle of pigs from LRFI and HRFI lines [77]. It is possible that the significant SNPs identified in this study are in LD with the causative mutations. However, additional strategies are needed in the future for fine mapping the causal mutation.

Conclusions

Heritability estimates for growth and feed efficiency traits in broilers from 28 to 41 days of age were low to moderate, ranging from 0.12 to 0.26. Narrow QTL regions with additive pleiotropic effects were associated with BW28 and BW42 (spanning 0.32 Mb) and RFI (spanning 1.16 Mb). The most likely candidate genes for these QTL, NSUN3, EPHA6, and AGK, are known to be involved in mitochondrial function and behavioral regulation. These results contribute to the identification of candidate genes and SNPs for feed efficiency traits in poultry.

Availability of data and materials

Data sharing is not applicable to this article since no datasets were generated or analyzed during the current study.

References

  1. 1.

    Zampiga M, Flees J, Meluzzi A, Dridi S, Sirri F. Application of omics technologies for a deeper insight into quali-quantitative production traits in broiler chickens: A review. J Anim Sci Biotechnol. 2018;9:61.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  2. 2.

    Crews DH Jr. Genetics of efficient feed utilization and national cattle evaluation: a review. Genet Mol Res. 2005;4:152–65.

    PubMed  Google Scholar 

  3. 3.

    Koch RM, Swiger LA, Chambers D, Gregory KE. Efficiency of feed use in beef cattle. J Anim Sci. 1963;22:486–94.

    Article  Google Scholar 

  4. 4.

    Archer JA, Arthur PF, Herd RM, Parnell PF, Pitchford WS. Optimum postweaning test for measurement of growth rate, feed intake, and feed efficiency in British breed cattle. J Anim Sci. 1997;75:2024–32.

    CAS  PubMed  Article  Google Scholar 

  5. 5.

    Shirali M, Varley PF, Jensen J. Bayesian estimation of direct and correlated responses to selection on linear or ratio expressions of feed efficiency in pigs. Genet Sel Evol. 2018;50:33.

    PubMed  PubMed Central  Article  Google Scholar 

  6. 6.

    Nascimerito ML, Souza A, Chaves AS, Cesar ASM, Tullio RR, Medeiros SR, et al. Feed efficiency indexes and their relationships with carcass, non-carcass and meat quality traits in Nellore steers. Meat Sci. 2016;116:78–85.

    Article  Google Scholar 

  7. 7.

    Faure J, Lefaucheur L, Bonhomme N, Ecolan P, Meteau K, Coustard SM, et al. Consequences of divergent selection for residual feed intake in pigs on muscle energy metabolism and meat quality. Meat Sci. 2013;93:37–45.

    CAS  PubMed  Article  Google Scholar 

  8. 8.

    Mao F, Chen L, Vinsky M, Okine E, Wang Z, Basarab J, et al. Phenotypic and genetic relationships of feed efficiency with growth performance, ultrasound, and carcass merit traits in Angus and Charolais steers. J Anim Sci. 2013;91:2067–76.

    CAS  PubMed  Article  Google Scholar 

  9. 9.

    Hoque MA, Kadowaki H, Shibata T, Oikawa T, Suzuki K. Genetic parameters for measures of residual feed intake and growth traits in seven generations of Duroc pigs. Livest Sci. 2009;121:45–9.

    Article  Google Scholar 

  10. 10.

    de Koning DJ, Windsor D, Hocking PM, Burt DW, Law A, Haley CS, et al. Quantitative trait locus detection in commercial broiler lines using candidate regions. J Anim Sci. 2003;81:1158–65.

    PubMed  Article  Google Scholar 

  11. 11.

    Mignon-Grasteau S, Rideau N, Gabriel I, Chantry-Darmon C, Boscher MY, Sellier N, et al. Detection of QTL controlling feed efficiency and excretion in chickens fed a wheat-based diet. Genet Sel Evol. 2015;47:74.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  12. 12.

    Xu Z, Ji C, Zhang Y, Zhang Z, Nie Q, Xu J, et al. Combination analysis of genome-wide association and transcriptome sequencing of residual feed intake in quality chickens. BMC Genomics. 2016;17:594.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  13. 13.

    Seabury CM, Oldeschulte DL, Saatchi M, Beever JE, Decker JE, Halley YA, et al. Genome-wide association study for feed efficiency and growth traits in U.S. beef cattle. BMC Genomics. 2017;18:386.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  14. 14.

    Lu D, Miller S, Sargolzaei M, Kelly M, Vander Voort G, Caldwell T, et al. Genome-wide association analyses for growth and feed efficiency traits in beef cattle. J Anim Sci. 2013;91:3612–33.

    CAS  PubMed  Article  Google Scholar 

  15. 15.

    Santana MH, Utsunomiya YT, Neves HH, Gomes RC, Garcia JF, Fukumasu H, et al. Genome-wide association analysis of feed intake and residual feed intake in Nellore cattle. BMC Genet. 2014;15:21.

    PubMed  PubMed Central  Article  Google Scholar 

  16. 16.

    Onteru SK, Gorbach DM, Young JM, Garrick DJ, Dekkers JC, Rothschild MF. Whole genome association studies of residual feed intake and related traits in the pig. PLoS One. 2013;8:e61756.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  17. 17.

    Yuan J, Wang K, Yi G, Ma M, Dou T, Sun C, et al. Genome-wide association studies for feed intake and efficiency in two laying periods of chickens. Genet Sel Evol. 2015;47:82.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  18. 18.

    Yuan J, Chen S, Shi F, Wu G, Liu A, Yang N, et al. Genome-wide association study reveals putative role of gga-miR-15a in controlling feed conversion ratio in layer chickens. BMC Genomics. 2017;18:699.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  19. 19.

    Visscher PM, Wray NR, Zhang Q, Sklar P, McCarthy MI, Brown MA, et al. 10 years of GWAS discovery: biology, function, and translation. Am J Hum Genet. 2017;101:5–22.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  20. 20.

    Liu R, Xing S, Wang J, Zheng M, Cui H, Crooijmans R, et al. A new chicken 55K SNP genotyping array. BMC Genomics. 2019;20:410.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  21. 21.

    Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–75.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  22. 22.

    Browning BL, Zhou Y, Browning SR. A one-penny imputed genome from next-generation reference panels. Am J Hum Genet. 2018;103:338–48.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  23. 23.

    Gilmour A, Gogel B, Cullis B, Welham S, Thompson R. ASReml user guide release 4.1 structural specification. Hemel Hempstead: VSN International Ltd; 2015.

    Google Scholar 

  24. 24.

    Legarra A, Aguilar I, Misztal I. A relationship matrix including full pedigree and genomic information. J Dairy Sci. 2009;92:4656–63.

    CAS  PubMed  Article  Google Scholar 

  25. 25.

    VanRaden PM. Efficient methods to compute genomic predictions. J Dairy Sci. 2008;91:4414–23.

    CAS  PubMed  Article  Google Scholar 

  26. 26.

    Zhou X, Stephens M. Genome-wide efficient mixed-model analysis for association studies. Nat Genet. 2012;44:821–4.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  27. 27.

    Gao X, Starmer J, Martin ER. A multiple testing correction method for genetic association studies using correlated single nucleotide polymorphisms. Genet Epidemiol. 2008;32:361–9.

    PubMed  Article  Google Scholar 

  28. 28.

    Aulchenko YS, Ripke S, Isaacs A, van Duijn CM. GenABEL: an R library for genome-wide association analysis. Bioinformatics. 2007;23:1294–6.

    CAS  PubMed  Article  Google Scholar 

  29. 29.

    Al-Mamun HA, Kwan P, Clark SA, Ferdosi MH, Tellam R, Gondro C. Genome-wide association study of body weight in Australian Merino sheep reveals an orthologous region on OAR6 to human and bovine genomic regions affecting height and weight. Genet Sel Evol. 2015;47:66.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  30. 30.

    Belonogova NM, Svishcheva GR, Axenovich TI. FREGAT: an R package for region-based association analysis. Bioinformatics. 2016;32:2392–3.

    CAS  PubMed  Article  Google Scholar 

  31. 31.

    Schifano ED, Epstein MP, Bielak LF, Jhun MA, Kardia SLR, Peyser PA, et al. SNP Set Association Analysis for Familial Data. Genet Epidemiol. 2012;36:797–810.

    PubMed  PubMed Central  Google Scholar 

  32. 32.

    Bages S, Estany J, Tor M, Pena RN. Investigating reference genes for quantitative real-time PCR analysis across four chicken tissues. Gene. 2015;561:82–7.

    CAS  PubMed  Article  Google Scholar 

  33. 33.

    Kim S, Kim T. Selection of optimal internal controls for gene expression profiling of liver disease. Biotechniques. 2003;35:456–60.

    CAS  PubMed  Article  Google Scholar 

  34. 34.

    Tu Y, Su Y, Wang K, Zhang X, Tong H, Gao Y. Gene expression of heart and adipocyte fatty acid-binding protein in chickens by FQ-RT-PCR. Asian-Austral J Anim Sci. 2010;23:987–92.

    CAS  Article  Google Scholar 

  35. 35.

    Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. 2001;25:402–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  36. 36.

    SAS, User’s guide: statistics version 9.4. Cary: SAS Institute; 2016.

  37. 37.

    Bernon DE, Chambers JR. Genetic parameters of unadjusted and adjusted age-constant feed consumption and efficiency of meat type chickens. Poult Sci. 1988;67:1497–504.

    Article  Google Scholar 

  38. 38.

    Abdollahi-Arpanahi R, Morota G, Valente BD, Kranis A, Rosa GJ, Gianola D. Differential contribution of genomic regions to marked genetic variation and prediction of quantitative traits in broiler chickens. Genet Sel Evol. 2016;48:10.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  39. 39.

    Aldridge MN, Vandenplas J, Bergsma R, Calus MPL. Variance estimates are similar using pedigree or genomic relationships with or without the use of metafounders or the algorithm for proven and young animals. J Anim Sci. 2020;98:a019.

    Article  Google Scholar 

  40. 40.

    Harder I, Stamer E, Junge W, Thaller G. Estimation of genetic parameters and breeding values for feed intake and energy balance using pedigree relationships or single-step genomic evaluation in Holstein Friesian cows. J Dairy Sci. 2020;103:2498–513.

    CAS  PubMed  Article  Google Scholar 

  41. 41.

    Bohlouli M, Alijani S, Naderi S, Yin T, Konig S. Prediction accuracies and genetic parameters for test-day traits from genomic and pedigree-based random regression models with or without heat stress interactions. J Dairy Sci. 2019;102:488–502.

    CAS  PubMed  Article  Google Scholar 

  42. 42.

    Powell JE, Visscher PM, Goddard ME. Reconciling the analysis of IBD and IBS in complex trait studies. Nat Rev Genet. 2010;11:800–5.

    CAS  PubMed  Article  Google Scholar 

  43. 43.

    Clément V, Bibé B, Verrier É, Elsen J-M, Manfredi E, Bouix J, et al. Simulation analysis to test the influence of model adequacy and data structure on the estimation of genetic parameters for traits with direct and maternal effects. Genet Sel Evol. 2001;33:369–95.

    PubMed  PubMed Central  Article  Google Scholar 

  44. 44.

    Abdollahi-Arpanahi R, Pakdel A, Nejati-Javaremi A, Moradi Shahrbabak M, Morota G, Valente BD, et al. Dissection of additive genetic variability for quantitative traits in chickens using SNP markers. J Anim Breed Genet. 2014;131:183–93.

    CAS  PubMed  Article  Google Scholar 

  45. 45.

    Schenkel FS, Miller SP, Wilton JW. Genetic parameters and breed differences for feed efficiency, growth, and body composition traits of young beef bulls. Can J Anim Sci. 2004;84:177–85.

    Article  Google Scholar 

  46. 46.

    Ceacero TM, Mercadante ME, Cyrillo JN, Canesin RC, Bonilha SF, de Albuquerque LG. Phenotypic and genetic correlations of feed efficiency traits with growth and carcass traits in Nellore cattle selected for postweaning weight. PLoS One. 2016;11:e0161366.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  47. 47.

    Parsch J, Ellegren H. The evolutionary causes and consequences of sex-biased gene expression. Nat Rev Genet. 2013;14:83–7.

    CAS  Article  Google Scholar 

  48. 48.

    van der Heide EMM, Lourenco DAL, Chen CY, Herring WO, Sapp RL, Moser DW, et al. Sexual dimorphism in livestock species selected for economically important traits. J Anim Sci. 2016;94:3684–92.

    PubMed  Article  Google Scholar 

  49. 49.

    Mebratie W, Madsen P, Hawken R, Jensen J. Multi-trait estimation of genetic parameters for body weight in a commercial broiler chicken population. Livest Sci. 2018;217:15–8.

    Article  Google Scholar 

  50. 50.

    Nikbakht G, Esmailnejad A. Chicken major histocompatibility complex polymorphism and its association with production traits. Immunogenetics. 2015;67:247–52.

    CAS  PubMed  Article  Google Scholar 

  51. 51.

    McNab FW, Rajsbaum R, Stoye JP, O’Garra A. Tripartite-motif proteins and innate immune regulation. Curr Opin Immunol. 2011;23:46–56.

    CAS  PubMed  Article  Google Scholar 

  52. 52.

    Montori-Grau M, Pedreira-Casahuga R, Boyer-Diaz Z, Lassot I, Garcia-Martinez C, Orozco A, et al. GNIP1 E3 ubiquitin ligase is a novel player in regulating glycogen metabolism in skeletal muscle. Metabolism. 2018;83:177–87.

    CAS  PubMed  Article  Google Scholar 

  53. 53.

    Inoue E, Yamauchi J. AMP-activated protein kinase regulates PEPCK gene expression by direct phosphorylation of a novel zinc finger transcription factor. Biochem Biophys Res Commun. 2006;351:793–9.

    CAS  PubMed  Article  Google Scholar 

  54. 54.

    Yang J, Ferreira T, Morris AP, Medland SE. Genetic investigation of ANthropometric Traits (GIANT) consortium, DIAbetes Genetics Replication And Meta-analysis (DIAGRAM) Consortium, et al. Conditional and joint multiple-SNP analysis of GWAS summary statistics identifies additional variants influencing complex traits. Nat Genet. 2012;44:369–375.

  55. 55.

    Shirai T, Inoue E, Ishimi Y, Yamauchi J. AICAR response element binding protein (AREBP), a key modulator of hepatic glucose production regulated by AMPK in vivo. Biochem Biophys Res Commun. 2011;414:287–91.

    CAS  PubMed  Article  Google Scholar 

  56. 56.

    Hansen C, Yi N, Zhang YM, Xu S, Gavora J, Cheng HH. Identification of QTL for production traits in chickens. Anim Biotechnol. 2005;16:67–79.

    CAS  PubMed  Article  Google Scholar 

  57. 57.

    Haag S, Sloan KE, Ranjan N, Warda AS, Kretschmer J, Blessing C, et al. NSUN3 and ABH1 modify the wobble position of mt-tRNAMet to expand codon recognition in mitochondrial translation. EMBO J. 2016;35:2104–19.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  58. 58.

    Nakano S, Suzuki T, Kawarada L, Iwata H, Asano K, Suzuki T. NSUN3 methylase initiates 5-formylcytidine biogenesis in human mitochondrial tRNA(Met). Nat Chem Biol. 2016;12:546–51.

    CAS  PubMed  Article  Google Scholar 

  59. 59.

    Trixl L, Amort T, Wille A, Zinni M, Ebner S, Hechenberger C, et al. RNA cytosine methyltransferase Nsun3 regulates embryonic stem cell differentiation by promoting mitochondrial activity. Cell Mol Life Sci. 2018;75:1483–97.

    CAS  PubMed  Article  Google Scholar 

  60. 60.

    Watson HJ, Yilmaz Z, Thornton LM, Hübel C, Coleman JRI, Gaspar HA, et al. Genome-wide association study identifies eight risk loci and implicates metabo-psychiatric origins for anorexia nervosa. Nat Genet. 2019;51:1207–14.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  61. 61.

    Lindfors C, Nilsson IAK, Garcia-Roves PM, Zuberi AR, Karimi M, Donahue LR, et al. Hypothalamic mitochondrial dysfunction associated with anorexia in the anx/anx mouse. Proc Natl Acad Sci USA. 2011;108:18108–13.

    CAS  PubMed  Article  Google Scholar 

  62. 62.

    Savelieva KV, Rajan I, Baker KB, Vogel P, Jarman W, Allen M, et al. Learning and memory impairment in Eph receptor A6 knockout mice. Neurosci Lett. 2008;438:205–9.

    CAS  PubMed  Article  Google Scholar 

  63. 63.

    Dos Santos FC, Peixoto MG, Fonseca PA, Pires MF, Ventura RV, Rosse ID, et al. Identification of candidate genes for reactivity in Guzerat (Bos indicus) cattle: A genome-wide association study. PLoS One. 2017;12:e0169163.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  64. 64.

    Li C, He J, Chen J, Zhao J, Gu D, Hixson JE, et al. Genome-wide gene-sodium interaction analyses on blood pressure: the genetic epidemiology network of salt-sensitivity study. Hypertension. 2016;68:348–55.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  65. 65.

    Herd RM, Arthur PF. Physiological basis for residual feed intake. J Anim Sci. 2009;87:E64-71.

    CAS  PubMed  Article  Google Scholar 

  66. 66.

    Das G, Yu QL, Hui RY, Reuhl K, Gale NW, Zhou RP. EphA5 and EphA6: regulation of neuronal and spine morphology. Cell Biosci. 2016;6:48.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  67. 67.

    D’Occhio MJ, Baruselli PS, Campanile G. Influence of nutrition, body condition, and metabolic status on reproduction in female beef cattle: a review. Theriogenology. 2019;125:277–84.

    CAS  PubMed  Article  Google Scholar 

  68. 68.

    Reyer H, Hawken R, Murani E, Ponsuksili S, Wimmers K. The genetics of feed conversion efficiency traits in a commercial broiler line. Sci Rep. 2015;5:16387.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  69. 69.

    Waggoner DW, Johnson LB, Mann PC, Morris V, Guastella J, Bajjalieh SM. MuLK, a eukaryotic multi-substrate lipid kinase. J Biol Chem. 2004;279:38228–35.

    CAS  PubMed  Article  Google Scholar 

  70. 70.

    Bektas M, Payne SG, Liu H, Goparaju S, Milstien S, Spiegel S. A novel acylglycerol kinase that produces lysophosphatidic acid modulates cross talk with EGFR in prostate cancer cells. J Cell Biol. 2005;169:801–11.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  71. 71.

    Spiegel S, Milstien S. Functions of the multifaceted family of sphingosine kinases and some close relatives. J Biol Chem. 2007;282:2125–9.

    CAS  PubMed  Article  Google Scholar 

  72. 72.

    Fang Y, Vilella-Bach M, Bachmann R, Flanigan A, Chen J. Phosphatidic acid-mediated mitogenic activation of mTOR signaling. Science. 2001;294:1942–5.

    CAS  PubMed  Article  Google Scholar 

  73. 73.

    Tatsuta T, Scharwey M, Langer T. Mitochondrial lipid trafficking. Trends Cell Biol. 2014;24:44–52.

    CAS  PubMed  Article  Google Scholar 

  74. 74.

    Kang Y, Stroud DA, Baker MJ, De Souza DP, Frazier AE, Liem M, et al. Sengers syndrome-associated mitochondrial acylglycerol kinase is a subunit of the human TIM22 protein import complex. Mol Cell. 2017;67(457–70):e5.

    Google Scholar 

  75. 75.

    Vukotic M, Nolte H, Konig T, Saita S, Ananjew M, Kruger M, et al. Acylglycerol kinase mutated in Sengers syndrome is a subunit of the TIM22 protein translocase in mitochondria. Mol Cell. 2017;67(471–83):e7.

    Google Scholar 

  76. 76.

    Bottje WG, Carstens GE. Association of mitochondrial function and feed efficiency in poultry and livestock species. J Anim Sci. 2009;87:E48-63.

    CAS  PubMed  Article  Google Scholar 

  77. 77.

    Vincent A, Louveau I, Gondret F, Trefeu C, Gilbert H, Lefaucheur L. Divergent selection for residual feed intake affects the transcriptomic and proteomic profiles of pig skeletal muscle. J Anim Sci. 2015;93:2745–58.

    CAS  PubMed  Article  Google Scholar 

Download references

Acknowledgements

The authors thank Dr. Zezhao Wang (Chinese Academy of Agricultural Sciences, CAAS) and Dr. Wengang Zhang (University of Edinburgh) for helpful suggestions when performing the GWAS analyses. The authors thank Dr. Mohamed Shafey Elsharkawy (CAAS) for revising the paper.

Funding

The work was supported by the Agricultural Science and Technology Innovation Program (CAAS-ZDRW202005), the National Nonprofit Institute Research Grant (Y2020PT02), the earmarked fund for the modern agro-industry technology research system (CARS-41), Agricultural Science and Technology Innovation Program (ASTIP-IAS04; ASTIP-IAS-TS-15).

Author information

Affiliations

Authors

Contributions

WL contributed to the study, the analysis of data, and writing of the manuscript. MQZ contributed to analysis and interpretation of data. GPZ, JW, JL, SLW, FRF, DWL, DZ, QHL and LPG contributed to the sample and data collection. YMG contributed to interpretation of data and reviewing of the manuscript. RRL and JW contributed to the design of the study, interpretation of data, and writing of the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Ranran Liu or Jie Wen.

Ethics declarations

Ethics approval and consent to participate

All experimental procedures involving broilers were performed according to the Guidelines for Experimental Animals established by the Ministry of Science and Technology (Beijing, China). Ethical approval on animal survival was given by the animal welfare and ethics committee of the Institute of Animal Sciences (IAS), Chinese Academy of Agricultural Sciences (CAAS, Beijing, China) with the following reference number: IAS2019-44.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Statistics of selection pressure for males and females.

Additional file 2: Table S2.

Statistics of the number of male and female broilers, separately, slaughtered in hatches and generations 5 to 7.

Additional file 3: Table S3.

Distribution of SNPs after quality control.

Additional file 4: Table S4.

Primer sequences for qPCR analysis.

Additional file 5: Table S5.

Phenotypes for individual broilers.

Additional file 6: Table S6.

Descriptive statistics for male and female broilers.

Additional file 7: Table S7.

Estimates of variance components and heritabilities for males and females using the relationship matrix that blends genomic and pedigree information.

Additional file 8: Table S8.

Estimates of variance components and heritabilities for males and females using the pedigree relationship matrix.

Additional file 9: Table S9.

Estimates of genetic and phenotypic correlations among growth and feed efficiency traits based on the pedigree relationship matrix.

Additional file 10: Table S10.

Sequences of SNPs reported in Tables 4 and 5.

Additional file 11: Table S11.

Additive and dominance effects of the five most significant SNPs on growth and feed efficiency traits.

Additional file 12: Table S12.

Means (+ SE) of growth and feed efficiency traits by genotype of the five most significant SNPs in males and females.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Li, W., Zheng, M., Zhao, G. et al. Identification of QTL regions and candidate genes for growth and feed efficiency in broilers. Genet Sel Evol 53, 13 (2021). https://doi.org/10.1186/s12711-021-00608-3

Download citation