Genetic determinism of bone and mineral metabolism in meat-type chickens: A QTL mapping study.

Skeletal integrity in meat-type chickens is affected by many factors including rapid growth rate, nutrition and genetics. To investigate the genetic basis of bone and mineral metabolism, a QTL detection study was conducted in an intercross between two lines of meat-type chickens divergently selected for their high (D +) or low (D -) digestive efficiency. Tibia size (length, diameter, volume) and ash content were determined at 3 weeks of age as well as phosphorus (P) retention and plasma concentration. Heritability of these traits and their genetic correlations with digestive efficiency were estimated. A QTL mapping study was performed using 3379 SNP markers. Tibia size, weight, ash content and breaking strength were highly heritable (0.42 to 0.61). Relative tibia diameter and volume as well as P retention were strongly and positively genetically correlated with digestive efficiency (0.57 to 0.80). A total of 35 QTL were identified (9 for tibia weight, 13 for tibia size, 5 for bone strength, 5 for bone mineralization, 2 for plasma P concentration and 1 for P retention). Six QTL were genome-wide significant, and 3 QTL for tibia relative volume, weight and ash weight on chromosome 6 were fixed, the positive allele coming from the D-line. For two QTL for ash content on chromosome 18 and relative tibia length on chromosome 26, the confidence intervals were small enough to identify potential candidate genes. These findings support the evidence of multiple genetic loci controlling bone and mineral metabolism. The identification of candidate genes may provide new perspectives in the understanding of bone regulation, even beyond avian species.

Background

Growth and feed efficiency are the main criteria used for selection in poultry favouring lean muscle development, especially in fast-growing broilers. According to Havenstein et al. (2003), the body weight of broilers at 42 d of age in a 2001 strain was 4.5 times higher than in 1957 broilers. Efficiency of feed conversion was also largely improved in the 2001 strain compared with the 1957 strain (1.58 vs 2.05 on average from 0 to 42 d of age). However, this large improvement in the genetic potential of chickens may have adverse effects and it is generally recognized that disruption of the dynamics of tissue development, i.e. muscle vs skeleton, contributes to bone weakness. Previous experiments have shown that the cortical bone of fast-growing broilers is highly porous and poorly mineralized (Williams et al., 2000, Shim et al., 2012). Shim et al. (2012) also reported a negative association between tibia breaking strength and growth rate. Skeletal disorders and related leg problems represent one of the major metabolic diseases found in the very efficient modern genotypes, compromising bird welfare and resulting in substantial economic losses. These disorders are generally aggravated by nutritional factors such as inadequate mineral supply and poor feed management. In particular, in recent years the safety margins regarding dietary phosphorus provision, one of the main mineral constituents of bones, have been reduced to tackle with the high cost of feed phosphate and environmental issues. An imbalance between Ca and P supply and mixing errors also contribute to the incidence of these disorders in chicken broilers. Exploring the capacity of the chicken to be more efficient in utilizing P and mineralizing bones, especially through an understanding of the genetic basis therefore represents a valuable approach to counteracting these problems. Information regarding the genetic components of P utilization and bone quality in the chicken is still limited. Quantitative trait loci (QTL) for different bone traits and candidate genes have been identified in adult laying hens (Dunn et al., 2007, Rubin et al., 2007, Johnsson et al., 2015). However, the mechanisms involved in the determinism of bone traits may be different between light adult female hens and heavy male and female growing meat chickens. Divergent selection on digestive efficiency in previous experiments with meat-type chickens using the criterion of apparent metabolizable energy corrected to zero nitrogen (AMEn) induced changes in P utilization and bone characteristics. After 8 generations of divergent selection, birds from the high line had higher levels of P retention (+ 26%), bone breaking strength (+ 32%) and dry tibia weight and ash weight per unit of BW (+ 6.6% and + 8.1%, respectively) than those from the low line (de Verdal et al., 2013).

The aim of this study presented here was therefore to provide new information on the genetic basis of bone-related traits using a F2 intercross between two lines of meat-type chickens divergently selected for their high (D +) or low (D −) digestive efficiency (Mignon-Grasteau et al., 2004). The heritability of tibia size, weight and ash content and of phosphorus retention and blood concentration, and their correlations with digestive efficiency were evaluated. In addition, QTL mapping was performed to identify zones of the genome and candidate genes involved in the determinism of P and bone metabolism.

Materials and methods

Animals

The experiment was conducted according to the guidelines of the French Ministry of Agriculture and European regulations concerning animal experimentation, including authorization N°37–100 from the French Ministry of Agriculture. The experimental unit where birds were kept is registered by the Ministry of Agriculture for animal experimentation under license number C-37-175-1. Measurement of digestive efficiency in individual cages, blood sampling procedures for genotyping and euthanasia procedures by injection of pentobarbital were approved by the Ethics Committee for Animal Experimentation of Val de Loire (00886.02 and 01047.02). This Ethics Committee is registered by the National Committee under number C2EA-19. The personal license number from the French Veterinary Service for this study is 548.

Chicken from the D + and D − lines divergently selected on high or low AMEn, respectively (Mignon-Grasteau et al., 2004), were crossed at the eight generation of selection to produce an F2 design. The F2 generation consisted of 820 animals originating from 6 sires and 60 F1 dams (half from the cross of D + males with D − females, and half from D − males with D + females). Five batches of chicks were produced between January and June 2010.

From hatching to 10 d of age, birds were reared in one pen, and then they were transferred to individual cages. Throughout the experiment, birds were fed a diet similar to the diet containing 55% Rialto wheat used during the selection experiment (Tran et al., 2014). Between 20 and 23 d, they were submitted to a balance trial to determine the coefficient of phosphorus retention (RETP, in %). Total collection of feces was performed (Bourdillon et al., 1990). Excreta were freeze-dried and ground before chemical analyses. Phosphorus concentration in feces was measured by colorimetric analysis. In addition, the efficiency of digestion of energy, lipids, proteins, dry matter and starch was determined by infrared spectrophotometry as proposed by Bastianelli et al. (2010). The results on digestive efficiency can be found in Tran et al. (2014).

At 23 d, birds were weighed (BW) and samples of cardiac blood were collected in EDTA tubes for determination of phosphorus concentration (PP). After centrifugation at 3000 × g for 10 min at 4 °C, plasma was collected and stored at − 20 °C for analysis.

Bone measurements

Birds were slaughtered at 23 d, and right tibias were removed one day later, cleaned of muscle and stored at − 20 °C. Tibia length and diameter were measured with a digital caliper (TL, TD) and the total bone volume (TV) was estimated by considering the tibia to be roughly cylindrical (i.e., (π × TL × TD2)/4).

Tibias were then dried at 100 °C for 12 h and tibia dry matter was weighed (TDM). Tibias were then ashed in a muffle furnace at 550 °C for 12 h, tibia ash weight was recorded (TA) and bone mineral density was assessed according to the ratio of tibia ash to tibia dry matter weight (TA/TDM).

Tibia weight, length, diameter, volume and ash weight were also expressed as relative values by dividing TDM, TL, TD, TV or TA by BW.

Bone breaking strength (BBS) was measured with 3-point loading test using an Instron Universal Testing Instrument (model 5543, Instron SA, Buc, France) according to the method described by Letourneau-Montminy et al. (2008).

Elementary statistics on bone measurements are reported in Table 1.

Table 1

Elementary statistics of bone characteristics at 23 d in F2 birds.

Trait	N	Mean	Standard deviation
TL (mm)	825	64.57	3.30
TD (mm)	826	4.469	0.414
TV (mm3)	825	1028	224
TDM (g)	827	1.357	0.225
TL/BW (mm g− 1)	824	0.145	0.015
TD/BW (mm g− 1)	824	0.010	0.001
TV/BW (mm3 g− 1)	824	2.265	0.321
TDM/BW (%)	824	0.301	0.023
TA/BW (%)	824	0.119	0.010
TA/TDM (%)	827	39.72	1.80
RETP (%)	869	46.89	6.66
BBS (N)	826	88.32	18.91
PP (mg/L)	819	71.09	8.30

TDM, tibia dry matter weight; TDM, tibia dry weight; TL, tibia length; TD, tibia diameter; TV, tibia volume; BW, body weight at 23 d; TA, tibia ash weight; RETP, percentage of retained phosphorus; BBS, bone breaking strength; PP, plasma phosphorus concentration.

Markers and genotyping

All F0, F1 and F2 birds were genotyped with a dedicated Illumina Infinium custom array including 6000 SNP markers chosen for their informativity in our design and distribution across the genome (Tran et al., 2014). The markers presenting deviations from the Hardy-Weinberg balance within families, inconsistent genotyping relative to pedigree or genetic map information or poor quality markers were discarded from the analysis in order to reduce the risk of erroneous results (Tran et al., 2014). Finally, 3379 markers were used. The genetic map was deduced from the physical position of the SNP markers and from the genetic consensus reference map published by Groenen et al. (2009). This set of markers covers 3099.1 cM.

QTL analysis

QTL detection was carried out with the QTLMap software (Gilbert et al., 2008) using a half-sib model (Elsen et al., 1999, Le Roy et al., 1998) with interval mapping based on maximum likelihood estimations (Lander and Botstein, 1989). This model does not make assumptions on the number of QTL alleles segregating in the design. The traits were analyzed separately. Depending on preliminary analysis of variance, data were pre-corrected for fixed effects of batch (4 levels), sex (2 levels), rearing cell (3 levels), cage row (3 levels), slaughter per half-day (2 levels) and experimenter at slaughter (i.e. the person in charge of cutting intestinal segments, 7 levels). QTL analyses were performed by comparing the hypothesis of one QTL (H1) vs no QTL (H0) to test the segregation of a QTL on each linkage group. For chromosome Z, separate analyses were performed for males and females and are presented respectively as Zm and Zf.

The significance threshold at the chromosome-wide level for each trait on each chromosome was calculated from the results of 5000 simulations of performance under the null hypothesis. For the most significant QTLs, 20,000 simulations were made to derive the genome-wide p-value (PG) from the chromosome-wide p-value (PC) using an approximate Bonferroni correction:

PG = 1 − (1 − PC)1/r

where r is the ratio of the length of a specific chromosome to the length of the genome considered for QTL detection as in Tilquin et al. (2005). Confidence intervals for QTLs (95%) were estimated using the LOD drop-off method as proposed by Lander and Botstein (1989).

The significance of the QTL effects within each sire family was tested using a Student test, by assuming an equal distribution of the QTL alleles in the progeny. A QTL effect was retained as significant for Student test p-values < 0.05, and the corresponding sire families were assumed to segregate for this QTL. These familial substitution effects were estimated in families found to segregate significantly for the QTL.

Results

Genetic parameters

Three categories of traits could be distinguished for heritability of the recorded traits (Table 2). The first group included plasma phosphorus concentration, phosphorus retention and relative length of tibia that were not heritable or with a very low heritability estimate (from 0.03 to 0.08). The second group presented moderate but significant heritability estimates (ranging from 0.15 to 0.23) and included relative to body weight measures such as tibia diameter, volume, dry matter weight and ash weight. Finally, the third group showed high values of heritability (between 0.42 and 0.61). This group included raw size and bone weight measures as well as bone mineralization and breaking strength. Bone solidity was positively correlated with bone size (rg = 0.57 to 0.75) but negatively correlated with relative length (rg = − 0.92), indicating that bones that were relatively short compared to the global animal size were more solid. A better solidity was also associated with relatively heavy bones (rg = 0.76) and more mineralized bones (rg = 0.47). Phosphorus retention was positively correlated with relative size and weight of bone, indicating that a better capacity to retain phosphorus was associated with an improved capacity to build bones. In contrast, plasma phosphorus concentration, which results from multiple metabolic pathways, was not correlated with most of the other bone traits.

Table 2

Heritability (± standard errors, on the diagonal) and genetic correlations (± standard errors, above the diagonal) of bone and phosphorus-related traits. Genetic parameter estimates in bold are significantly higher than 0.

Traits	TL		TD		TV		TDM		TL/BW		TD/BW		TV/BW		TDM/BW
TL	0.480	± 0.097	0.496	± 0.105	0.681	± 0.073	0.119	± 0.092	− 0.906	± 0.141	− 0.545	± 0.189	0.085	± 0.213	0.497	± 0.219
TD			0.425	± 0.072	0.978	± 0.010	0.353	± 0.044	− 0.457	± 0.299	0.652	± 0.200	0.981	± 0.029	0.926	± 0.102
TV					0.450	± 0.060	0.295	± 0.105	− 0.646	± 0.202	0.487	± 0.206	0.927	± 0.063	0.929	± 0.181
TDM							0.452	± 0.063	− 0.109	± 0.349	0.705	± 0.075	0.783	± 0.186	0.807	± 0.159
TL/BW									0.030	± 0.023	0.391	± 0.279	− 0.279	± 0.301	− 0.370	± 0.354
TD/BW											0.146	± 0.042	0.757	± 0.101	0.523	± 0.160
TV/BW													0.234	± 0.073	0.907	± 0.065
TDM/BW															0.152	± 0.042

TDM, percentage of tibia dry matter; TL/BW, TDM/BW, TV/BW, TA/BW, TD/BW, tibia length, dry matter weight, volume, ash and diameter relative to body weight at 23 d; TL, TV, TD, tibia length, volume and diameter; TA/TDM, tibia ash relative to tibia dry matter weight; BBS, bone breaking strength; RETP, phosphorus retention; PP, plasma phosphorus concentration.

Several bone traits were correlated with digestive efficiency, which was the criterion for selection of the D + and D − lines. Relative diameter and volume of tibia and the capacity to retain phosphorus were positively correlated with digestive efficiency of energy, which implies that birds with good digestive efficiency also had good digestive efficiency for phosphorus and a greater ability to develop large bones.

Finally, body weight was highly and positively correlated with raw length, diameter and volume of the tibia, which are all indicators of the overall size of birds (Table 3). In contrast, moderate negative correlations were found between the protein digestive coefficient and absolute tibia size.

Table 3

Genetic correlations between bone and phosphorus-related traits (± standard errors) and digestive efficiency and body weight. Genetic parameter estimates in bold are significantly higher than 0.

Traits	AMEn		BW
TL	− 0.286	± 0.199	0.880	± 0.049
TD	− 0.175	± 0.193	0.552	± 0.106
TV	− 0.173	± 0.145	0.701	± 0.066
TDM	− 0.006	± 0.127	0.111	± 0.110
TL/BW	0.322	± 0.489	− 0.970	± 0.061
TD/BW	0.805	± 0.219	− 0.532	± 0.236
TV/BW	0.569	± 0.226	0.166	± 0.254
TDM/BW	0.251	± 0.292	0.468	± 0.287
TA/BW	0.139	± 0.299	0.519	0.258
TA/TDM	− 0.225	± 0.127	0.179	± 0.097
BBS	− 0.254	± 0.189	0.698	± 0.075
PP	− 0.299	± 0.241	0.178	± 0.271
RETP	0.614	± 0.160	− 0.213	± 0.248

TDM, percentage of tibia dry matter; TL/BW, TDM/BW, TV/BW, TA/BW, TD/BW, tibia length, dry matter weight, volume, ash and diameter relative to body weight at 23 d; TL, TV, TD, tibia length, volume and diameter; TA/TDM, tibia ash relative to tibia dry matter weight; BBS, bone breaking strength; RETP, phosphorus retention; PP, plasma phosphorus concentration; AMEn, apparent metabolizable energy corrected to zero nitrogen retention; BW, body weight at 23 d.

Number and position of QTLs

We found a total of 35 QTLs that were significant at the chromosome level: 5 for absolute tibia weight traits, 4 for relative tibia weight traits, 7 for absolute tibia size traits, 6 for relative tibia size traits, 5 for bone solidity, 2 for plasma phosphorus concentration, 1 for phosphorus retention, and 5 for bone mineralization (Table 4, Fig. 1). Among these, 6 were genome-wide significant, 1 on chromosome 1 for relative length of tibia, 3 on chromosome 6 for relative volume, dry matter and ash weight, 1 on chromosome 18 for bone mineralization and 1 on chromosome 26 for relative tibia length.

Table 4

QTLs detected for phosphorus-related traits.

Chromosome	Trait	Position (M)	CIa (M)	Flanking markers of CI		Level of significance		QTL effectb (NFc) Fd
Chromosome-wide	Genome-wide
1	TDM	1.310	1.294–1.328	Gga_rs13865892	Gga_rs15268441	0.002	> 0.200	0.207 (5)
PP	1.940	1.877–1.958	Gga_rs16082461	Gga_rs13887893	0.050	> 0.200	0.214 (4)
TL/BW	4.600	4.523–4.717	GGaluGA059924	Gga_rs15548184	0.0009	0.050	0.215 (5)
TV	4.700	4.531–4.700	Gga_rs13986943	GGaluGA062691	0.046	> 0.200	0.198 (4)
3	TD	0.130	0.079–0.174	Gga_rs14309481	GGaluGA205532	0.025	> 0.200	0.185 (6)
TV	0.130	0.074–0.191	Gga_rs14309317	Gga_rs14316548	0.033	> 0.200	0.185 (6)
TDM	0.130	0.089–0.182	Gga_rs15257971	Gga_rs14316939	0.038	> 0.200	0.210 (5)
BBS	0.130	0.115–0.215	GGaluGA204542	Gga_rs14315147	0.028	> 0.200	0.187 (5)
TDM/BW	0.420	0.359–0.573	Gga_rs14083998	Gga_rs14326932	0.027	> 0.200	0.196 (5)F ±
4	TDM	0.240	0.066–0.359	Gga_rs14420072	Gga_rs10728994	0.050	> 0.200	0.227 (3)
6	TV/BW	0.320	0.278–0.493	Gga_rs10732029	Gga_rs13567432	< 0.0001	0.010	0.261 (5)F ∓
TDM/BW	0.320	0.278–0.387	Gga_rs10732029	Gga_rs16543395	< 0.0001	0.010	0.263 (5)F ∓
TA/BW	0.330	0.282–0.330	Gga_rs10732029	GGaluGA298536	< 0.0001	0.010	0.264 (5)F ∓
7	RETP	0.360	0.222–0.392	Gga_rs13740969	Gga_rs14608624	0.025	> 0.200	0.210 (4)F ∓
11	BBS	0.040	0.016–0.060	Gga_rs13689132	Gga_rs14958789	0.033	> 0.200	0.201 (3)
13	BBS	0.075	0.000–0.127	GGaluGA090214	Gga_rs14991199	0.048	> 0.200	0.272 (2)F0/−
TA/TDM	0.365	0.350–0.415	GGaluGA094679	Gga_rs14061598	0.048	> 0.200	0.214 (3)
15	TA/TDM	0.344	0.321–0.373	GGaluGA109023	GGaluGA109767	0.008	> 0.200	0.191 (5)
18	TA/TDM	0.163	0.143–0.172	Gga_rs15816457	GGaluGA119505	0.001	0.050	0.203 (4)
TD/BW	0.453	0.410–0.513	Gga_rs15042573	Gga_rs16348086	0.004	> 0.200	0.228 (4)F ∓
19	TV	0.060	0.000–0.297	Gga_rs16076471	GGaluGA127011	0.009	> 0.200	0.173 (6)
TD	0.060	0.013–0.128	Gga_rs15043317	Gga_rs13573868	0.010	> 0.200	0.186 (5)
BBS	0.060	0.000–0.313	Gga_rs16076471	Gga_rs14121172	0.003	> 0.200	0.196 (5)
21	TA/TDM	0.380	0.344–0.398	Gga_rs14284716	Gga_rs16074526	0.045	> 0.200	0.374 (1)F0/−
24	PP	0.360	0.343–0.396	GGaluGA192978	Gga_rs14296726	0.044	> 0.200	0.253 (2)F ±
26	TD	0.220	0.000–0.348	Gga_rs14710236	Gga_rs15235289	0.005	> 0.200	0.248 (3)
TV	0.220	0.000–0.356	Gga_rs14710236	Gga_rs13606624	0.020	> 0.200	0.227 (3)
TL/BW	0.250	0.235–0.263	Gga_rs16201749	GGaluGA196721	< 0.001	0.010	0.280 (3)
TD/BW	0.250	0.227–0.375	Gga_rs16201570	Gga_rs14300758	0.049	> 0.200	0.181 (4)
BBS	0.320	0.190–0.375	GGaluGA196172	Gga_rs14300758	0.007	> 0.200	0.203 (4)
TA/TDM	0.350	0.307–0.420	GGaluGA197068	Gga_rs16204548	0.038	> 0.200	0.191 (4)F ±
Zm	TV/BW	0.645	0.459–0.687	Gga_rs14757740	Gga_rs14066443	0.020	> 0.200	0.265 (5)F ∓
Zf	TDM	0.715	0.642–0.715	Gga_rs14754601	Gga_rs16765302	0.035	> 0.200	0.282 (5)
TA/BW	1.405	1.162–1.560	Gga_rs14694865	Gga_rs14765737	0.020	> 0.200	0.312 (5)
TDM/BW	1.425	1.166–1.573	Gga_rs14694865	Gga_rs14765737	0.009	0.120	0.260 (6)

TDM, percentage of tibia dry matter; TL/BW, TDM/BW, TV/BW, TA/BW, TD/BW, tibia length, dry matter weight, volume, ash and diameter relative to body weight at 23 d; TV, TD, tibia volume and diameter; TA/TDM, tibia ash relative to tibia dry matter weight; BBS, bone breaking strength; RETP, phosphorus retention; PP, plasma phosphorus concentration.

QTL effect as a proportion of the phenotypic standard deviation of trait.

Number of F1 sire families heterozygous for the QTL (p < 0.05, Student test).

F, the QTL is fixed in F1 sire families in which it is significant.

1-LOD-drop off confidence interval (lower and upper boundaries, M).

Fig. 1

Map of positions of QTLs for bone and phosphorus-related traits and co-localizations with QTLs for digestive efficiency (AMEn, metabolizable energy corrected to zero nitrogen retention, CDUS, coefficient of digestive use of starch, from Tran et al. (2014) and for the ratio of nitrogen to phosphorus in excreta (N/P, from Mignon-Grasteau et al. (2015).

On chromosomes 6 and Z, 2 and 3 QTLs for relative size and weight of bones were detected at the same position, on chromosomes 1 and 26, 2 and 4 QTLs co-localized for absolute and relative size of bones, and on chromosomes 3 and 19, 4 and 3 QTLs were detected associating absolute size and weight of bones and bone solidity, respectively.

These QTLs were detected in 1–6 of the 6 F1 sire families, and 18 of the 35 QTLs were present in 5 or 6 families. Among these 18 QTLS, 5 QTLs were fixed for relative weight and relative volume, the positive allele coming from the D − line in all families.

The effects of the QTLs were moderate and on average were at 0.227 standard deviation of the traits (ranging from 0.173 to 0.374).

The average confidence interval was 15.6 cM, ranging between 2.8 cM and 40.7 cM. The number of genes detected by annotqtl in the QTL zones were therefore highly variable, from 9 for the QTL or relative length of the tibia on chromosome 26–271 for the QTL of tibia dry matter weight on chromosome 4.

Discussion

Genetic parameters of bone characteristics

Despite the importance of bone-related traits in agriculture and health, very few studies have tried to estimate the genetic parameters of these traits. Despite the high levels of discrepancy between studies depending on the species, age and physiological status of the animals, the three categories of low, medium, and high heritability traits identified in this study were consistent with previous estimates. The heritability estimates for phosphorus retention and plasma concentrations were reported to be low to moderate heritability (Zhang et al., 2003, Ankra-Badu et al., 2010a, Wang et al., 2006) in chickens and quails (Beck et al., 2014). Moderate to high heritability estimates are generally reported for body size and weight, and our estimates of tibia size and weight were on the lower side of the range (0.14 to 0.74, with an average of 0.43) found in the literature for bone size and weight in poultry, cattle, pigs and humans (de Verdal et al., 2013, Albera et al., 2001, Gjerlaug-Enger et al., 2012, Jian et al., 2004, Vallee et al., 2013, Singh and Hjilani, 2008, Gonzalez-Ceron et al., 2015). Finally, our medium to high estimates of heritability of bone strength and ash content were also in the middle of the range reported in the literature, i.e. 0.18 to 0.77 for bone breaking strength in poultry and pigs (de Verdal et al., 2013, Gonzalez-Ceron et al., 2015, Bishop et al., 2000, Peixoto et al., 2010, Storskrubb et al., 2010) and 0.25–0.66 for bone mineral density and ash percentage in bone in man, monkey, and chicken (de Verdal et al., 2013, Peixoto et al., 2010, Liu et al., 2012, Hernandez-de Sosa et al., 2014, Lipkin et al., 2001).

Feed efficiency is one of the main selection criteria in poultry and its correlation with increased muscle mass and reduced fatness has been frequently studied. In contrast, the relationship with bone development has been much less documented. In our study, we found that efficiency was positively correlated with diameter and volume of the tibia, which is consistent with the fact that D + birds have heavier tibias (de Verdal et al., 2013). On the other hand, tibia length was not correlated with efficiency, which is consistent with the poor correlation coefficient between these two traits by de Verdal et al. (2013), but not with the 8% difference observed between the two lines for tibia length.

We also confirmed the absence of correlation or poor correlation between digestive efficiency and bone mineralization and strength that has been reported both in broilers and laying hens (de Verdal et al., 2013, Bordas et al., 1992) and the strong correlation between digestive efficiency and phosphorus bioavailability or retention estimates in (de Verdal et al., 2013, Ankra-Badu et al., 2010a). These results together suggest that selection on digestive efficiency of energy has improved the general capacity to digest, including the absorption of P, leading to greater skeletal development without any deleterious impact on its structure and strength.

Positions of QTLs

The regions identified in this study for tibia size and breaking strength on chromosome 3, breaking strength on chromosome 13 and relative tibia diameter, ash weight and dry matter weight on chromosomes 18, 21 and Z are not mentioned in the literature. Apart from these regions, all the QTLs identified in this study have already been reported in previous studies (Supplementary Table S1). This was for example the case for our largest QTLs on relative volume, dry matter weight and ash weight on chromosome 6, which have already been reported for shank diameter, tibia bone marrow diameter and femur weight in the literature (Sharman et al., 2007, Nadaf et al., 2009). Similarly, the two genome-wide QTLs with a rather low number of candidate genes (i.e., for tibia relative ash weight on chromosome 18 and relative tibia length on chromosome 26) have also previously been identified for various traits concerning dimensions, weight and strength of tibia, humerus and shank in chickens (Schreiweis et al., 2005, Park et al., 2006, Ankra-Badu et al., 2010b). The greater SNP marker density used in our study provided generally smaller confidence intervals (i.e., 2.8 to 2.9 cM on GGA18 and GGA26 vs 11–47 cM in previous studies, and 4.8–21.5 cM vs 13.2–34.5 cM on, GGA6).

As the lines used in the F2 cross had been divergently selected on digestive efficiency and showed marked differences in the anatomy of the digestive tract, we looked for co-localizations between QTL of bone characteristics and feed efficiency, digestive efficiency and anatomy of the digestive tract (Fig. 1). QTL of absolute and relative tibia size and breaking strength co-localize with feed efficiency QTL detected on chromosomes 19 and 26 with the same design (Mignon-Grasteau et al., 2015). Two QTLs of tibia breaking strength and mineralization co-localize with a QTL of coefficient of digestive use of starch on chromosome 26 (Tran et al., 2014). This might indicate a pleiotropic effect of efficiency genes on bone traits such as genes related to the bone morphogenetic protein (see below) involved in both processes.

Co-localization between QTL of bone traits and anatomy of the digestive tract are also present on chromosomes 1, 15, 18, 26 and Z (Tran et al., 2014, Mignon-Grasteau et al., 2015). Among them, we found for example an association between QTLs for tibia dry matter, tibia diameter and duodenum weight on chromosomes 1 and 26. Absorption occurring in the small intestine is an important process involved in the maintenance of phosphorus and calcium homeostasis. The duodenum represents a major site of active calcium and phosphorus transport from the lumen to the blood involving different transporters (Rivoira et al., 2012, Fang et al., 2012). The genes encoding for two of these calcium transporters are present in the QTLs detected in this study, namely TRPV6 which controls the entry of calcium into the enterocyte within the QTL of PP on chromosome 1 and NCX1 which controls the exit of calcium from the enterocyte to the blood on the QTL for relative bone weight on chromosome 3. It can also be hypothesized that selection on digestive efficiency modified the total quantity of phosphorus present in the intestinal lumen. This would in turn modify the passive pathway of phosphorus absorption, which depends on an electrochemical gradient at the level of the tight junctions. This is consistent with the fact that several genes affecting the structure of the tight junctions have previously been identified within QTLs of feed or digestive efficiency with this design (Mignon-Grasteau et al., 2015).

Candidate genes

In most of the QTL regions identified in this study, the confidence interval was too great to search for potential candidate genes. However, the QTLs for ash proportion on chromosome 18 and for relative tibia length on chromosome 26 had small confidence intervals, with 23 and 9 genes respectively, and were genome-wide significant.

Four genes have been identified as functional candidate genes in the QTL for ash percentage on chromosome 18: one gene involved in cell differentiation (meteorin-like METRNL), two genes involved in the cytoskeleton structure and function (Tubulin Folding Cofactor 2 TBCD and Ras analog in brain member 40b Rab40b), and one gene involved in cell proliferation (Forkhead Box K2 FOXK2). METRNL is expressed in osteoblasts both in humans and guinea pigs (Gong, 2007, Jorgensen et al., 2012). Its expression is high during embryogenesis and early post-natal growth and decreases to a low level at adult age (Jorgensen et al., 2012). Moreover, in the presence of TGF-β and BMPII (bone morphogenetic protein II), METRNL expression is depressed in mesenchymal stem cells, precursors of osteoblasts (Sang et al., 2014). Finally, osteoblasts overexpressing METRNL have been found to have fewer mineralization nodes than those with normal expression (Gong, 2007).

The involvement of TBCD in bone mineralization is suggested by detection of a QTL of bone mineral accretion in children. In this study, only two SNP variations were associated with bone mineral accretion, one of them being within an intron of TBCD (Park et al., 2015). Tubulin cofactor D is one of the five cofactors involved in the assembly of microtubules which transport the acidic vesicles to the ruffled border of osteoclasts, where they are used to dissolve the extracellular bone matrix (Itzstein et al., 2011). Rab40b is a member of the Rab family, which is the largest family of small GTPases (Pereira-Leal and Seabra, 2001). This family plays a key role in vesicular trafficking that is necessary for the normal functioning of osteoclasts (Pereira-Leal and Seabra, 2001). Within the Rab family, Rab40 is located in the Golgi apparatus and is probably involved in intra-Golgi trafficking (Stenmark, 2009). Another potential role of Rab40b in bone formation has been suggested by Lee et al. (2007) who showed that it was involved in the activation of the non-canonical Wnt pathway, which in turn affects the development and activation of osteoblasts (Caverzasio, 2009).

FOXK2 was detected in previous study as a positional and functional candidate gene for a QTL of bone mineral density in the Rat genome database (rgd.mcw.edu). Expression of this gene is enhanced in the presence of bone substrate in the mouse and could be involved in bone cell death and the survival process (Purdue et al., 2014).

Four genes were also identified as potential candidate genes for bone length on chromosome 26, all involved in the osteoclast and osteoblast differentiation process. The C4BPA gene (Complement Component 4 binding protein alpha) belongs to the complement family whose role in osteoclast and osteoblast formation was reviewed by Schoengraf (2011). The expression of C4BPA is abnormal in the synovial fluid of persons suffering from osteoarthritis (Wang et al., 2011) and it has been found to interact with BMPR2, one of the main receptors of the bone morphogenetic protein BMP2 (Hassel et al., 2004). However, BMP2 has been shown to be an important gene to ensure normal development of the gastro-intestinal tract, a trait that has been significantly modified by selection on digestive efficiency in our chicken lines (de Verdal et al., 2011).

Two micro-RNA (MIR29C and MIR29B-2) are also functional candidates for this QTL. They are involved in the control of the formation and homeostasis of bone (Liao et al., 2014). The MIR29 family members promote osteoclastogenesis in the mouse through its negative regulation of the RNA of genes such as Cdc42 and Srgap2 which are essential for cytoskeleton organization (Franceschetti et al., 2013). The calcitonin receptor which regulates osteoclast survival and resorption is also a target of MIR29 (Franceschetti et al., 2013).Within this family, MIR29B is a key-regulator of osteoblast differentiation as it downregulates inhibitors of the TGFβ3 and Wnt-β-catenin pathways and upregulates the expression of Runx2, an essential transcription factor for bone formation (Li et al., 2009). Moreover, MIR29B also reduces the synthesis of collagen in the extra-cellular matrix during the osteoblast proliferation phase, which is believed to stimulate collagen fibril maturation and enhance mineral deposition (Li et al., 2009).

The last gene in this region was plexin A2 (PLXN2). Identified as a potential candidate gene in a QTL detection study on human prognathism (Ikuno et al., 2014), this gene both decreases osteoclastic bone resorption and stimulates osteoblastic bone formation (Oh et al., 2012, Hayashi et al., 2012).

Other loci for BBS on chromosomes 13 and 19 have two candidates: osteonectin (SPARC) and heat shock transcription factor family member 5 (HSF5). Osteonectin (SPARC) is an extra-cellular matrix glycoprotein synthesized by osteoblasts that has a role in bone remodeling, cell migration and proliferation, and matrix mineralization. Previous experiments in osteonectin-null mice reported deterioration in bone formation and decreases in osteoblast and osteoclast surface and number leading to severe osteopenia (Delany et al., 2000). Johnsson et al. (2015) recently identify this gene as a candidate for total metaphyseal content using an expression QTL (eQTL) mapping study. There is also evidence that members of the heat shock transcription factor family (HSF) are involved in mechanisms related to bone remodeling. HSF1 is recognized for its role in osteoclast formation suggesting a direct effect of cell stress on bone loss (Delany et al., 2000). Other findings suggest that HSF2 may participate in the modulation of RANK ligand gene expression in osteoblast cells (Chai et al., 2014).

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

All authors contributed to the manuscript preparation and scientific discussion. In addition, AN contributed to the study conception and design, and phenotyping of the digestive tract, CCD and MYB to genotyping of animals, NS to supervision of animal rearing, MC to reproduction of experimental lines and creation of the F2 cross, ELBD to study conception, design and phenotyping and SMG to coordination of the project, study conception and design, selection of experimental lines and genetic analyses. All authors read and approved the final manuscript.