Genome-Wide Association Study Identifies Loci and Candidate Genes for Body Composition and Meat Quality Traits in Beijing-You Chickens

Body composition and meat quality traits are important economic traits of chickens. The development of high-throughput genotyping platforms and relevant statistical methods have enabled genome-wide association studies in chickens. In order to identify molecular markers and candidate genes associated with body composition and meat quality traits, genome-wide association studies were conducted using the Illumina 60 K SNP Beadchip to genotype 724 Beijing-You chickens. For each bird, a total of 16 traits were measured, including carcass weight (CW), eviscerated weight (EW), dressing percentage, breast muscle weight (BrW) and percentage (BrP), thigh muscle weight and percentage, abdominal fat weight and percentage, dry matter and intramuscular fat contents of breast and thigh muscle, ultimate pH, and shear force of the pectoralis major muscle at 100 d of age. The SNPs that were significantly associated with the phenotypic traits were identified using both simple (GLM) and compressed mixed linear (MLM) models. For nine of ten body composition traits studied, SNPs showing genome wide significance (P<2.59E−6) have been identified. A consistent region on chicken (Gallus gallus) chromosome 4 (GGA4), including seven significant SNPs and four candidate genes (LCORL, LAP3, LDB2, TAPT1), were found to be associated with CW and EW. Another 0.65 Mb region on GGA3 for BrW and BrP was identified. After measuring the mRNA content in beast muscle for five genes located in this region, the changes in GJA1 expression were found to be consistent with that of breast muscle weight across development. It is highly possible that GJA1 is a functional gene for breast muscle development in chickens. For meat quality traits, several SNPs reaching suggestive association were identified and possible candidate genes with their functions were discussed.


Introduction
With the advance of high-throughput genotyping platforms, much effort has been spent on identifying molecular markers and genes related to complex traits using genome-wide association studies (GWAS) in several species.
In the field of animal breeding, loci or narrow regions affecting milk production, fertility and growth traits in cattle [1][2][3], body composition, intramuscular fat content, meat color in pigs [4,5] and growth and egg quality in chickens [6,7] have been detected successfully using genome-wide association studies. Such information helps in the development of marker assisted breeding as well as by improving understanding of the molecular mechanisms underlying the target traits.
Body composition and meat quality in broilers are important economic traits. Body composition traits have been analyzed using QTL techniques in F2 crosses between various lines of chickens. A total of 146 QTLs reaching significance have been reported via genome scans based on marker-QTL linkage analyses (http://www.animalgenome.org/cgi-bin/QTLdb/GG/ index,July,2012), which associated with eight body composition traits: carcass weight, dressed percentage, weight and percentage of breast muscle, thigh muscle and abdominal fat. Similar studies on meat quality traits are rare except for 10 QTLs reported for pH value. Application of these QTL results in broiler breeding remains impracticable because of the low precision of mapping.
In the present work, with the aim of identifying potential loci and candidate genes affecting body composition and meat quality traits, 724 chickens from a conservation population of a typical Chinese local breed (Beijing-You chicken) were used in GWAS studies. A total of 16 body composition and meat quality traits were measured or derived.

Ethics Statement
The animal component of this study was conducted in accordance with the Guidelines for the Experimental Animals, established by the Ministry of Science and Technology (Beijing, China). Blood was collected from the brachial vein of the chickens at 80 d by venipuncture, using citrated syringes during a routine health inspection at the experimental station of the Chinese Academy of Agricultural Sciences (CAAS). Animal experiments were approved by the Science Research Department (in charge of animal welfare issue) of the Institute of Animal Sciences, CAAS (Beijing, China).

Birds and Phenotypic Traits
The 728 male Beijing-You chickens were generated from 50 half-sib families. The 50 sires and 250 dams were chosen from conservation stock for this breed and were unrelated. The Beijing-You chickens maintained by the Institute of Animal Sciences, CAAS and have not been subjected to any systematic selection for any trait. All birds were reared in stair-step caging under the same recommended nutritional and environmental conditions. A total of 16 carcass and meat quality traits were measured for the GWAS studies: carcass weight (including feet and head, CW), eviscerated weight (EW), breast weight (BrW), thigh muscle weight (ThW), abdominal fat weight (AbFW), dressed percentage (DP), eviscerated yield, as a percentage of CW (EWP), BrW, ThW and AbFW as percentages of EW (BrP, ThP, AbFP), dry matter and intramuscular fat contents of breast and thigh muscle (DM Br , DM Th , IMF Br and IMF Th ), ultimate pH (pHu), and shear force (SF) of the pectoralis major muscle. The measurements and data were obtained as follows: After a 12-h fast, birds were weighed (LW) and slaughtered at d 100 by standard commercial procedures and CW was recorded. The removable adipose tissues surrounding the proventriculus and gizzard along with those located around the cloaca were weighed as AbFW [8,9]. Then EW was weighed. Carcasses were then dissected into deboned, skinless thighs and breasts for weighing and storage at 220uC until used.
The data for CW and EW were also expressed as percentages on the basis of LW and CW, respectively (DP and EWP). The pHu was the mean of three measurements taken from the left breast muscle after 24-h storage at 4uC (IQ150 pH meter (Hach Company, Loveland, USA). Shear force was determined on breast muscles following Li's method [10] using a universal Warner-Bratzler testing machine MTS Synergie 200 (G-R Manufacturing Company, Manhattan, KS). Dry matter (DM) and IMF in muscle were determined for 20 g samples of the right breast and thigh after removing obvious fat, mincing, and drying in two stages (,12 h each at 65uC then 105uC, as described by Cui et al [11]. The DM content was expressed on the basis of fresh muscle weight. The IMF was measured by Soxhlet extraction with anhydrous diethyl ether and was expressed as a percentage of muscle DM.
One-day-old hatchlings with similar genetic background were reared in the same conditions as above. On each of weeks 4, 8, 10, 12, and 14, six to eight birds of similar weight were randomly selected, stunned, and euthanized using approved procedures. Breast muscle was rapidly dissected, weighed, snap-frozen in liquid nitrogen, and stored at 280uC for mRNA extraction.

Genotyping and Quality Control
Genomic DNA was extracted from the blood samples using phenol-chloroform and was diluted to 50 ng/ml. Each chicken was genotyped using the Illumina 60 K Chicken SNP Beadchip (DNA LandMarks Inc., Saint-Jean-sur-Richelieu, Quebec, Canada). Of the total of 728 chickens, four were excluded because sample call rate was ,95%. Approximately 20% (12,088) of the SNPs were removed for one or more of the following: low call rate (,95%), minor allelic frequency (,0.01) and Hardy-Weinberg equilibrium (HWE) test (p,1E206), or chromosomal location was unknown. After imposing these constraints, 724 individuals and 45,548 SNP markers that distributed on 30 autosomes and the Z chromosome (Table 1) were used for the genome-wide association analyses.

Statistical Analysis
The SNPs that were significantly associated with the phenotypic traits were identified using both simple (GLM, I) and compressed mixed linear (MLM, II) models [12]: where Y is the phenotypic value, X is the genotype (45,548 SNPs), F is the family, and K is the relative kinship matrix; Xa were regarded as fixed effects, while Fb and Km were random effects, and e is the random error. The relative kinship matrix (K) was constructed from 8,006 independent SNP markers, acquired using Plink v1.07 software [13] through all autosomal SNPs, pruned using the indep-pairwise option, with a window size of 25 SNPs, a step of 5 SNPs, and r 2 threshold of 0.2. The analyses were performed using TASSEL 3.0 software [14].
The raw data for some traits (AbFP, AbFW, IMF Br , DM Th , pHu and SF) deviated from normality and Box-Cox or Johnson transformations were applied using Minitab 15 (http://www. minitab.com). Significance thresholds were established from the estimated number of independent SNP markers and LD blocks, defined as a set of contiguous SNPs having pairwise r 2 values .0.40. Using this approach, the estimated total number of independent SNP markers and LD blocks was 19,284. The two threshold P-values were therefore set at 5.19E205 (1/19,284) for suggestive significance, and 2.59E206 (0.05/19,284) for genomewide significance [15].
Quantile-quantile (Q-Q) plots for each trait and Manhattan plots of genome-wide association analyses were produced with R 2.13.2 software (http://www.r-project.org/).

Quantitative Measurement of mRNA
The methods for quantitative measurement of mRNA were referred to Li's methods [10]. Primers for the genes NCOA7, TPD52L1, FABP7, GJA1 and ASF1A were designed (Primer Premier 5.0) from the GenBank sequences (Table 2).

Results
Descriptive statistics of the phenotypic measurements of body composition and meat quality traits in the 724 Beijing-You chickens used for the present GWAS studies are given in Table 3. All non-normal phenotypic data except those for pHu were normalized after the Box-Cox or Johnson transformation ( Table 4).
The results for all SNPs demonstrated to have genome-wide significance (P,2.59E206) or suggestive significance (P,5.19E205) are presented in full in Table S1 and Figure 1. For the traits examined, emphasis is placed on the associations revealed by the compressed MLM analyses because the population structure effect shown under the GLM model could be controlled effectively when the compressed MLM was used ( Figure S1). Of the SNPs revealed by the compressed MLM model, status of SNPs that reached genome-wide significance exposed by the simple model were also indicated.

Loci and Genes for Traits Related to Body Composition
Carcass Weight (CW). As detailed in Table 5, there were 10 SNPs that were significantly associated with CW from the compressed MLM, of which six were of genome-wide significance by analysis with MLM and/or GLM. Seven of the SNPs on chicken (Gallus gallus) chromosome 4 (GGA4) were clustered within a 1.06 Mb region (between 78,475,066 bp and 79,531,679 bp), and are located either within or 19 kb-236 kb away from the nearest known genes: ligand dependent nuclear receptor corepressor-like protein (LCORL), leucineamino peptidase 3 (LAP3), quinoid dihydro pteridine reductase (QDPR); LIM domain-binding protein 2 (LDB2), transmembrane anterior posterior transformation1 (TAPT1). The two SNPs on GGA1 are located within dachshund homolog 1 (Drosophila) (DACH1).
Eviscerated Weight (EW). In the case of EW, nine significant SNPs were identified by compressed MLM, of which five were of genome-wide significance by GLM analysis. The six SNPs on GGA4 are the same as those found for CW but there were slight differences in which of them were of greater significance ( Table 5). The one SNP on GGA1, within DACH1, was also the same as that associated with CW. The remaining two SNPs were on GGA7 and GGA19, neither was associated with any other trait.
Dressed Percentage (DP) and Percentage of Eviscerated Weight (EWP). As seen from Table 5, three SNPs with suggestive significance for DP were identified by the two methods. The two SNPs on GGA4 were in the vicinity (within 92 and  (Table S1).
Thigh muscle Weight (ThW) and Percentage (ThP). Four SNPs, significantly associated with ThW were identified by compressed MLM and occurred on GGA3, GGA11, GGA2 and GGA16; those on GGA3 and GGA11 were of genome-wide significance ( Table 7). The SNP on GGA3 was also associated with LW, though of slightly lesser significance. The significant SNP on GGA11was close (8 kb) to nuclear transport factor 2 (NUTF2) and this locus was also associated with ThP. The other SNP for ThP was of genome-wide significance and was on GGA2, just 2 kb from forkhead box N4 (FOXN4).

Loci and Genes for Meat Quality Traits
The SNPs associated with four traits related to meat quality are provided in Table 8; no SNPs were significantly associated with intramuscular fat content of thigh (IMF Th ) and shear force (SF) of breast muscle.
Dry Matter content in Breast (DM Br ). Five SNPs of significance were identified, located on GGA2 and GGA27. The three SNPs on GGA2 are in the vicinity of or within angiopoietin 1 (ANGPT1). The SNPs on GGA27 have some proximity to FtsJ homolog 3 (E. coli, FTSJ3); SNP Gga_rs16205470 was the only SNP related to meat quality shown to have genome-wide significance.
Intramuscular Fat in Breast (IMF Br ). Two significant SNPs were identified and located on GGA2 and GGA5. The SNP on GGA2 is 120.9 kb away from cholecystokinin (CCK). The SNP on GGA5 is 159.7 kb from Toll interacting protein (TOLLIP).
Dry Matter content in Thigh (DM Th ). Six significant SNPs associated with DM Th were identified and located on GGA1, GGA2 and GGA8. The SNP on GGA2 is 208.4 kb away from suppressor of cytokine signaling 6 (SOCS6). The SNP on GGA8 is 8.5 kb from UMP-CMP kinase (CMPK1).The four SNPs on GGA1 are clustered within a 1.30 Mb region with no annotated genes nearby.
Ultimate pH (pHu). Eleven significant SNPs were identified and located on GGA2, GGA4 and GGA7. The eight SNPs on GGA2 are distributed within a 2.06 Mb region and only one gene (Polypeptide N-acetylgalactosaminyltransferase 1, GALNT1) is in the vicinity. Two SNPs on GGA4 located within protocadherin 19 (PCDH19) or 286.3 kb away from diaphanous homolog 1 (DIAPH1). The SNP on GGA7 located 143.9 kb from the secreted phosphoprotein 2 (SPP2).

Validation of Candidate Genes for BrW and BrP from GWAS by Q-PCR
Expression of candidate genes detected near associated signals for BrW and BrP in the GWAS analysis was further tested by realtime quantitative PCR (Q-PCR) in breast muscle. After measuring the mRNA content of five genes (NCOA7, TPD52L1, FABP7, GJA1, ASF1A) located in a 0.65 Mb region for BrW and BrP, the change in GJA1 expression was found to be consistent with that of the breast muscle weight across development (4, 8, 10, 12 and 14 week) (Figure 2). The correlation between GJA1 mRNA and BrW was moderate (r = 0.554) and significant at the 0.001 level. It is highly possible that GJA1 is a functional gene for breast muscle development in chickens.

Genome-wide Association Analysis
To identify potential loci and candidate genes affecting chicken body composition and meat quality traits, GWAS studies were performed using a conservation population of Beijing-You (BJY) chickens. The BJY chicken is one representative indigenous breed in China [16], a color-feathered, slow-growing chicken and has superior quality of meat products [17,18]. Most of the traits tested showed considerable ranges between maximal and minimal values (Table 3), as would be expected in a population being maintained for the conservation of genetic diversity; this variability would be expected to increase the power of the GWAS.
Two statistical methods, compressed mixed linear model (MLM) and generalized linear model (GLM), were implemented to analyze association between SNPs and phenotypes. Emphasis is placed on the associations revealed by the compressed MLM analyses because population structure effect could be controlled and false positives could be reduced with this approach, as shown in Q-Q plots ( Figure S1). However, because the degree of association might be reduced in MLM [19], so of the SNPs revealed by compressed MLM, status of SNPs that reached genome-wide significance exposed by the simple model were also indicated. In addition, because many SNPs reaching suggestive association in both MLM and GLM located within similar regions to those with genome-wide significance, it is proposed that SNPs with suggestive association also indicate important loci.

Loci and Genes for Traits Related to Body Composition
Breast muscle yield is the most important carcass component in meat-type chickens because of the high premium paid by consumers. Of special interest, one important region (61.83 Mb-68.57 Mb) on GGA3 was identified as being associated with BrW and BrP. Most of the SNPs were located within or near tRNA methyltransferase 11 (TRMT11), nuclear receptor coactivator (NCOA7), tumor protein D53 (TPD52L1), fatty acid binding protein 7 (FABP7) and gap junction protein, alpha 1(GJA1) genes. Of these five genes, only the change of GJA1 (gap junction protein, alpha 1) expression was consistent with that of BrW across development (4, 8, 10, 12 and 14 week) and the correlation between mRNA level and BrW was significant (p,0.001). Thus, it is highly possible that GJA1 is a functional gene for chicken breast muscle development. This result might supply a novel functional gene for breast muscle development and extent the known function of GJA1 which has been found to play a role in skeletal form [20,21]. In practical breeding programs, birds could be selected for breeding stock based on the desired allele in GJA1 for BrW and BrP.
For CW and EW, one consistent region was identified (about 78.47 Mb to 79.53 Mb) on GGA4, which corresponded to QTL regions previously reported [22,23]. There were seven significant SNPs located near four functional genes (LCORL, LAP3, LDB2, and TAPT1). Polymorphisms in LCORL have been detected to associate with human skeletal frame size, linear growth [24][25][26][27] and feed intake and growth of cattle [28]. The LIM domainbinding factor 2 (LDB2) has been associated with chicken body weight (7-12 wk) and average daily gain (6-12 wk) in another GWAS study [15]. The transmembrane protein, TAPT1, involved in transporting molecules across membranes, was speculated [29] to be a downstream effector of HOXC8 that may relate to axial skeletal patterning during development. The amino peptidase LAP3 catalyzes the removal of amino acids and peptides as part of protein maturation and degradation [30]. On GGA1, significant SNPs associated with both CW and EW were found within DACH1 and this gene located near a QTL region known to be related to CW [31]. The gene DACH1 is a target of FGF signaling during limb skeletal development [32]. In addition, two significant SNPs located in GGA3 were found to relate with LW and were also in the QTL region for 40d body weight and 21d body weight, respectively [33,34].
The SNP Gga_rs14402759 on GGA3 was found to associate with both ThW and LW, and this locus was near MYCN. These two associations are consistent with the relatively high correlation (r 2 = 0.75) found between two traits, which possibly reflects their being controlled by common loci.

Loci and Genes for Traits Related to Meat Quality Traits
Three significant SNPs on GGA2 were found to associate with DM Br and were located in the vicinity of or within ANGPT1. Angiopoietin-1 (ANGPT1) is an angiogenesis factor that is also an important modulator of skeletal muscle function [35]. Angiopoietin-1 (ANGPT1) is mainly produced by cardiac, skeletal and smooth muscle cells, and adventitial cells [36]. The additional two SNPs on GGA27 were near FTSJ3 (a putative ortholog of yeast Spb1p) and conditional knockdown revealed that depletion of Figure 1. Manhattan plots showing association of all SNPs with carcass and meat quality traits from compressed MLM. SNPs are plotted on the x-axis according to their position on each chromosome against association with these traits on the y-axis (shown as -log10 p-value). The black dashed line indicates genome-wise significance of suggestive association (p-value = 5.19E205), and the red dashed line shows genomewise 5% significance with a p-value threshold of 2.59E206. Abbreviations: CW, carcass weight; EW, eviscerated weight; DP, dressed percentage; BrW, breast muscle weight; BrP, percentage of breast muscle; ThW, thigh muscle weight; ThP, Percentage of thigh muscle; AbFW, weight of abdominal fat; AbFP, percentage of abdominal fat; DM Br, dry matter content in breast;IMF Br , intramuscular fat in breast; DM Th , dry matter content in thigh; pHu, ultimate pH. doi:10.1371/journal.pone.0061172.g001 Table 5. SNPs with genome-wide and suggestive significance for carcass weight, eviscerated weight and dressed percentage traits. FTSJ3 affects HEK293 cell proliferation and causes pre-rRNA processing defects [37]. For DM Th , four SNPs of suggestive significance on GGA1 were clustered within a 1.3 Mb region, currently lacking any annotated gene. Another SNP on GGA2 was located near SOCS6, which is involved in the regulation of glucose metabolism and plays an important role in regulating insulin action in vivo [38].
A SNP on GGA2 was associated with IMF Br and cholecystokinin (CCK) was found nearby. The gut hormone CCK plays a multiplicity of roles, including those influencing digestion of fat [39].
For ultimate pH, eight SNPs within a 2.06 Mb region on GGA2 were identified although only one known gene (N-acetylgalacto-saminyltransferase 1, GALNT1) has been found near this region. The encoded protein catalyzes the transfer of N-acetylgalactosamine (GalNAc) from UDP-GalNAc to the hydroxyl group of a serine or threonine residues on proteins with O-linked glycosylation, a common post-translational modification. Two significant SNPs were found on GGA4, within PCDH19 and near DIAPH1, respectively, and another on GGA7 near SPP2. Protocadherin-19 (PCDH19) plays a role as an adhesion protein in optic nerve fiber bundling, optic nerve targeting, and/or synapse formation [40]. The gene DIAPH1 encodes a protein that may have a role in the regulation of actin polymerization in hair cells of the inner ear. The SPP2 gene encodes a secreted phosphoprotein that is a member of the cystatin superfamily. Table 6. SNPs with genome-wide and suggestive significance for breast muscle weight and percentage of breast muscle.  SNPs with superscript ''a'' were of genome-wide significance by compressed MLM, those with superscript ''b'' were significant by GLM, and those with superscript ''ab'' were significant by both methods; these SNPs are shown underlined in boldface; 3 The nearest known gene to the genome-wide significant SNPs are shown underlined in boldface. doi:10.1371/journal.pone.0061172.t007 Table 8. SNPs with genome-wide and suggestive significance for meat quality traits. These proteins play important roles in tumorigenesis, stabilization of matrix metalloproteinases, glomerular filtration rate, immunomodulation, and neurodegenerative diseases. A number of candidate loci associated with meat quality and fat traits were identified, but common loci or genes for these traits were rare. These results reinforce the notion of complexity in the genetic basis underlying meat quality and fat deposition; to a certain extent, they might be influenced by epigenetic factors [41]. Figure S1 Quantile-quantile (Q-Q) plots of the GLM (red dots) and compressed models (blue dots) for carcass and meat quality traits. Plotted on the x-axis are the expected p-values under the null hypothesis and on the y-axis are the observed p-values. A: CW, carcass weight; B, EW, eviscerated weight; C, DP, dressed percentage; D, EWP, percentage of eviscerated yield; E, BrW, breast muscle weight; F, BrP, percentage of breast muscle; G, ThW, thigh muscle weight; H, ThP, Percentage of thigh muscle; I, AbFW, weight of abdominal fat; J, AbFP, percentage of abdominal fat; K, DM Br, dry matter content in breast; L, IMF Br , intramuscular fat in breast; M, DM Th , dry matter content in thigh; N, IMF Th , intramuscular fat in thigh; O, pHu, ultimate pH, P, SF, shear force of breast muscle. (TIF)

Supporting Information
Table S1 Information for all significant SNPs related to body composition and meat quality traits. (XLSX) Figure 2. The relative mRNA expression of gap junction protein, alpha 1 gene (a) was consistent with that of the breast muscle weight (b) across development. The mRNA expression is shown as the number of copies (610 5 ) per mg total RNA; Data are means (n = 6). doi:10.1371/journal.pone.0061172.g002