Identification of genomic regions and candidate genes for chicken meat ultimate pH by combined detection of selection signatures and QTL

In chicken as in pigs, muscle ultimate pH (pHu) is a major factor of variation of both meat quality and processing ability. Normal values of pHu of broiler breast meat are approximately 5.8 to 5.9. The further the pHu deviates from this value, the more quality defects occur. Meat with pHu values greater than 6.1 is classified as dark, firm, and dry (DFD)-like meat while meat with pHu values less than 5.7 is classified as acid meat which is often referred to as pale, soft, and exudative (PSE)-like meat in broiler [1]. These quality defects lead to economic losses induced by reduced water holding capacity, increased cooking loss, decreased tenderness, and reduced emulsification capacity in the case of acid meat. In case of DFD-like meat, atypical color, after-flavor, dry and sticky texture, and decreased product shelf life are observed [2]. The near perfect genetic correlation (− 0.97) between pHu and glycogen content of breast muscle highlights a common genetic background between the two traits [3]. Glycogen is the major storage form for carbohydrates particularly in the liver and skeletal muscle [4]. In chicken, selection for increased growth rate and breast muscle mass has been associated with reduced glycogen storage [3, 5, 6]. Moreover, recent studies have reported reduced muscle glycogen content or elevated ultimate pH in breast muscles affected by degenerative disorders such as white striping and wooden breast [7, 8]. The understanding of the biological determinism of glycogen content and ultimate pH is thus a key point for the control of muscle integrity and meat quality in chicken.

It is now well established that chicken breast meat ultimate pH has a moderate to high heritability, as shown by the heritability values of 0.3 to 0.5 found in different genetic lines [3, 9–12]. However, in contrast to pigs where numerous QTLs (Animal QTLdb) and a major gene [13] have been identified for pHu, the genetic architecture of this trait remains poorly understood in chicken.

In the present study, we took advantage of a unique resource population constituted from two lines divergently selected for pHu, to identify the genomic regions underlying the genetic variability of this trait and to search for candidate genes. As described in Alnahhas et al. [14], the two lines have been divergently selected since 2009 with the breeding value of the pHu of the pectoralis major muscle (PM-pHu) as selection criterion. This population originated from a grandparental fast-growing line of broiler chicken, selected for a balance between growth and reproduction traits. As the measure of PM-pHu requires sacrificing the birds, a sib-selection was applied. As shown in Fig. 1, the selection process has been quite efficient. After six generations of selection, mean PM-pHu was estimated at 5.67 in the pHu- line while it was equal to 6.16 in the pHu + line (p-value < 0.0001). This means that 61% of the breast meat in the pHu- line could be classified as acid or PSE-like (PM-pHu < 5.7) and 63% of breast meat in the pHu + line as DFD (PM-pHu > 6.1). Significant changes (p-value < 0.0001) were also observed in the thigh since the pHu of the Sartorius muscle was on average 6.20 and 6.53 in the pHu- and pHu + lines, respectively. As expected, by comparison to the PM muscle of the pHu + line, the PM muscle of the pHu– line was characterized by a higher Glycolytic potential [14]. In order to identify the genomic regions potentially affected by this divergent selection, a genome-wide scan for loci with outstanding genetic differentiation between the two lines was performed on the first five generations of selection. Although the divergent selection focused on a specific trait, several other traits diverged between the two lines. This was the case, as expected, for several parameters of meat quality (such as color, drip loss, texture) but also for the percentage of meat (thigh plus breast) which was higher in the pHu + line than in the pHu- line (14). Consequently, it was not possible to establish a direct relationship between particular selection signatures and the phenotypic expression of PM-pHu. Therefore, genome-wide QTL mapping was conducted on more than 550 birds from the sixth generation of selection to identify genomic regions with direct impact on the ultimate pH of the breast and the thigh meat, two traits which are known to be strongly genetically correlated [14].

Chicken is a good model for studying genomic regions under selection. This is mainly due to the availability of valuable animal populations, such as divergent lines, and increasing genomic resources. The effects of single trait divergent selection have already been studied in this species with the 60 K SNP array. Johansson et al. [18] analyzed the genome of two divergent lines after 50 generations of selection for body weight and concluded that from 50 up to over 100 regions were under selection. Zhang et al. [19] considered a shorter term experiment, i.e. they looked for selection signatures after 11 generations of divergent selection for the percentage of abdominal fat and identified more than 50 regions within each of the two lines. In our study, even fewer generations (5) of divergent selection were considered, which should limit the detection of false positive selection signals due to genetic drift. We also controlled the risk of false positives by using statistical tests that explicitly estimate the amount of drift in sampled populations, and detect regions where genetic diversity is significantly inconsistent with a neutral model including this amount of drift. This resulted in the detection of 63 genomic regions under selection, at a false discovery rate (FDR) level of 5%. The inclusion of reference genotypes from G0 increased the detection power of this analysis. Indeed, the key point in hapFLK (resp. FLK) is to quantify, for each locus, the haplotype (resp. allele) frequency variation from an ancestral founder population to each sampled population (see the methods for more details). In most studies this ancestral population is unobserved, so ancestral frequencies are estimated from the frequencies in observed populations, but having direct access to these ancestral frequencies necessarily improves the approach.

Another original aspect of our study was to consider genotypes at successive generations, which allowed characterizing the dynamics of selection in regions underlying the genetic progress. Out of the 10 most significant hapFLK and FLK signals, with a p-value below 1e-6, four (named hapFLK-1b, hapFLK-1c, hapFLK-2c, hapFLK-5a in Table 1) were detected as soon as the first generation (G1) and were recovered at each of the four following generations (G2 to G5), and another one (hapFLK-9e) was detected at three generations (G2, G4 and G5). Four others (hapFLK-1f, hapFLK-2a, hapFLK-26a in Table 1 and FLK-13b in Table 2) corresponded to a later onset of selection, as they were detected only at G5. As shown in Table 1, various profiles were observed for the less significant selection signals, most of them being detected at only one generation and some others at 2 or 3 generations of the selection process. The dynamics of selection in candidate regions could also be visualized by plotting haplotype frequencies along generations in the two lines (Fig. 3, Additional file 2: Figure S2 and Additional file 3: Figure S3).

From a methodological perspective, our study also highlighted the interest of using haplotype instead of single SNP information for detecting selection signatures. Indeed, 53 candidate regions under selection were detected with hapFLK, with p-values as small as 1e-16 (Table 1), while only 14 were detected with FLK, with p-values down to 1e-8. The increased power provided by haplotypes has been demonstrated by Fariello et al. [16] using computer simulations, and has been confirmed since then by the outcomes of several applied studies in sheep [20], chicken [21] or cattle [22]. In their computer simulations, Fariello et al. [16] considered selection at a single SNP; in this case the higher power of hapFLK was related to the fact that the causal SNP itself was generally not genotyped, so that allele frequencies at this SNP (and their variation between populations) was better quantified using local haplotypes than other individual SNPs. However, as noted by the authors, another strong advantage of hapFLK is that selection can also act directly on a combination of linked alleles, due to epistatic effects; this situation is clearly easier to capture using haplotypes, because allele frequencies at each of the linked variants may only show minor variations.

Nevertheless, among the 14 regions detected by FLK, 10 were not detected by hapFLK (Table 2). This is a bit unexpected, because at the short time scale of our experiment recombination is limited, so the increase in frequency of one allele should be associated to that of a larger haplotype. In order to check that the FLK signal in these regions was not spurious, we considered the most significant SNP in each of the ten regions, and evaluated how likely it is that the allele frequency trajectories observed at this locus in the two lines, are reached due to drift. For this purpose, we used forward simulations without selection, but with population sizes mimicking our experiment (see the Methods for more details). When simulations were performed using a standard Wright-Fisher model, p-values remained small but generally increased, and some of them were no more significant (Additional file 6: Table S1). While the FLK and simulation approaches are both based on the allele frequency variation from generation 0 to 5, FLK p-values are computed by assuming that this variation follows a normal distribution. The accuracy of this approximation is known to decrease when the allele frequency in the ancestral population is close to 0 or 1, which is generally the case for the 10 SNPs considered here. This likely explains the difference between the two approaches. However, when considering a more realistic simulation model accounting for the exact number of male and female reproducers at each generation of the experiment, most p-values were found again very small. Besides, it should be noted that 7 of the 10 regions of Table 2 include several significant FLK SNPs or / and suggestive GWAS hits (BF > =5). Overall, we therefore believe that most of these 10 regions are true selective signals. The fact that they were not detected with hapFLK might be related to the difficulty of reconstructing haplotypes in regions where linkage disequilibrium is weak. Indeed, the average r2 within G0, for markers distant by less than 500 kb, was 0.442 (sd 0.005) for the 10 regions of Table 2, versus 0.467 (sd 0.0008) genome-wide and 0.512 (sd 0.004) for the 53 regions of Table 1. This suggests that combining hapFLK with FLK might still be useful in some cases, such as low SNP density or high recombination rate.

As a result of the divergent selection for PM-pHu, genomic regions showing selection signatures were enriched in QTL controlling this trait. On average, the Bayes Factor resulting from our association study in the whole population constituted from the two divergent lines was significantly larger for SNPs located within the candidate regions under selection of Tables 1 and 2, than for other SNPs in the genome (p-value of 0.0018 for the regions of Table 1 and 0.0078 for the regions of Table 2, with a two-sided Student test). Consistent figures were also found for the BF values issued from within-line analysis (p-value of 2.453e-06 for the regions of Table 1 and 0.04187 for the regions of Table 2). However, among the 24 QTL regions for PM-pHu trait of Table 3 (with BF > 20), only 7 corresponded to a selection signature (hapFLK-1c, hapFLK-1e, hapFLK-2a, hapFLK-5a, hapFLK-12b, hapFLK-14a, and FLK-8a). This is likely related to the polygenic architecture of PM-pHu, implying that QTL regions for this trait cannot all be efficiently selected in only 5 generations (at least not enough to be detected by our approach). Conversely, among the 63 regions under selection listed in Tables 1 and 2, only 7 corresponded to one of the 24 QTL regions for PM-pHu, and 8 others (including regions hapFLK-1b and hapFLK-2c discussed above) showed substantial evidence of association with PM-pHu (BF ≥ 5, data not shown). Thus, 48 regions with little evidence for association with PM-pHu showed significant selection signatures, and only 5% of those are expected to be false positives (see the methods). These regions likely correspond to QTL with relatively small effects, leading to moderate selection signatures. Indeed, as discussed previously, the use of individuals from G0 and, in the case of hapFLK, of haplotype information, allows detecting selection events with high power. In contrast, the power of the GWAS study is largely driven by single SNP allele frequency differences between the two lines, and moderate differences should not allow detecting association with PM-pHu. Finally, we cannot exclude the fact that traits or functions have been modified by the selection process although they are not directly related to muscle energy metabolism and meat ultimate pH.

The study of the two divergent lines allowed a better understanding of the genetic architecture of chicken meat ultimate pH and revealed several outstanding regions with highly significant signatures of selection and/or direct association with PM-pHu and SART-pHu phenotypes. In order to identify putative candidate genes, these regions were examined in light of the results of the transcriptomic analysis performed on samples from the pectoralis major muscle of the pHu + and pHu- lines [23] and of the knowledge of genes involved in human glycogen storage diseases (GSD) [24]. The most significant PM-pHu QTL identified in the current study (PM-1a) was located on GGA1 (BF = 300). It co-localized with a SART-QTL (named SART-1a) and with the most significant signature of selection detected by hapFLK at all generations of selection (named hapFLK-1c). The SNP exhibiting the strongest relation with PM-pHu (i.e. Gga_rs13842050) was segregating between the two lines, 86% of the pHu- birds being homozygous AA and 77% of the pHu + birds homozygous GG at the 6th generation of selection. It is likely that this SNP is in linkage disequilibrium with one or several mutations that affect PM-pHu. As shown in Additional file 7: Table S2, several differentially expressed (DE) genes were evidenced in this region. PPP1R3A, located ~ 506 kb away from the most significant SNP, is of special interest as it codes for a muscle-specific regulatory subunit of protein phosphatase 1 (PP1). By steering the catalytic subunit (PP1c) to glycogen, it promotes dephosphorylation of glycogen synthase (GS) and glycogen phosphorylase (GP) and thereby glycogen synthesis [24]. Knockout mice lacking PPP1R3A exhibit a 90% reduction in muscle glycogen [25]. A few years ago, a mutation was identified in human and was shown to be the first prevalent mutation known to impair glycogen synthesis and to decrease glycogen levels in human skeletal muscle [26]. PPP1R3A was significantly over-expressed in the pHu- line [23] which is consistent with the higher muscle glycogen content observed in this line [14]. Another gene, SLC37A4, involved in human GSD was present in the top list of the regions of interest (Additional file 7: Table S2). It is located ~ 507 kb away from one of the most significant SNP detected for PM-pHu on GGA24 (BF = 110). It corresponds to glucose-6-phosphate translocase (also called G6PT) which transports the glucose 6-phosphate from the cytoplasm to the endoplasmic reticulum, where it works together with the glucose 6-phosphatase to break down the sugar molecule and release free glucose that can leave the cell. Mutations in SLC37A4 are estimated to account for about 20% of the GSD type I (or Von Gierke disease) characterized by the accumulation of glycogen in some organs and tissues [27]. As indicated in Additional 7 Table S2, several other DE genes were identified in the regions of interest. Two of them (RHOC and LOC107052650) were part of a set of 20 biomarkers evidenced as pertinent predictors of the variation of PM-pHu between the two lines [23]. This set also comprised VTI1B, SLC2A1 and CAV3 which were within signatures of selection or QTL regions detected at a lower level of significance in the current study (hapFLK-5e for VTI1B, QTL PM-21a for SLC2A1, and QTL PM-12c for CAV3). SLC2A1, also called GLUT1, is coding for Glucose transporter-1, which is considered to play a key role in maintaining basal glucose transport in most chicken cell types as in mammals [28]. Further investigation by eQTL detection should be envisaged for the most promising of these positional and expressional candidate genes. It is also worthwhile to note that two other genes responsible for human GSD [24] were present in the regions of interest even if they were not differentially expressed. This is the case for GAA (located within the QTL PM-14a) coding for lysosomal alpha-glucosidase which is essential for the degradation of glycogen to glucose in lysosomes and responsible for GSD type II or Pompe disease. It is also the case for PHKA1 located ~ 20 kb away from the most significant SNP (BF = 257) we identified for SART-pHu on GGA4 (QTL SART-4a). The PHKA1 gene encodes the muscle alpha regulatory subunit of phosphorylase kinase which catalyzes phosphorylation of glycogen phosphorylase and thus promotes glycogen degradation. Mutations in PHKA1 cause GSD type VIII which is usually a mild myopathy with slight elevation of plasma creatine kinase concentration and muscle glycogen content.

The divergent selection conducted on PM-pHu allowed the creation of a unique resource population. This population, as shown in the present study, is highly useful to understand the genetic control of not only the selection criterion but more generally that of muscle glycogen storage and metabolism in chicken. Sixty-three genomic regions showed significant signatures of positive selection and were enriched in QTL for PM-pHu. Still, several QTL with strong effect have not yet responded to selection, suggesting that the divergence between lines might be further increased. A few regions of major interest, with significant selection signatures all along the selection process and/or highly associated with PM-pHu or SART-pHu phenotypes, were evidenced for the first time in chicken. They suggested several candidate genes, directly involved in glycogen synthesis and degradation or in the balance between these two systems, which necessitate further investigation in chicken.