Genome-Wide Association Study of Meat Quality Traits in a White Duroc×Erhualian F2 Intercross and Chinese Sutai Pigs

Thousands of QTLs for meat quality traits have been identified by linkage mapping studies, but most of them lack precise position or replication between populations, which hinder their application in pig breeding programs. To localize QTLs for meat quality traits to precise genomic regions, we performed a genome-wide association (GWA) study using the Illumina PorcineSNP60K Beadchip in two swine populations: 434 Sutai pigs and 933 F2 pigs from a White Duroc×Erhualian intercross. Meat quality traits, including pH, color, drip loss, moisture content, protein content and intramuscular fat content (IMF), marbling and firmness scores in the M. longissimus (LM) and M. semimembranosus (SM) muscles, were recorded on the two populations. In total, 127 chromosome-wide significant SNPs for these traits were identified. Among them, 11 SNPs reached genome-wise significance level, including 1 on SSC3 for pH, 1 on SSC3 and 3 on SSC15 for drip loss, 3 (unmapped) for color a*, and 2 for IMF each on SSC9 and SSCX. Except for 11 unmapped SNPs, 116 significant SNPs fell into 28 genomic regions of approximately 10 Mb or less. Most of these regions corresponded to previously reported QTL regions and spanned smaller intervals than before. The loci on SSC3 and SSC7 appeared to have pleiotropic effects on several related traits. Besides them, a few QTL signals were replicated between the two populations. Further, we identified thirteen new candidate genes for IMF, marbling and firmness, on the basis of their positions, functional annotations and reported expression patterns. The findings will contribute to further identification of the causal mutation underlying these QTLs and future marker-assisted selection in pigs.


Introduction
Meat quality is one of the most important economical traits in farm animals. It is decisive for the suitability of the meat for further processing and storage including retail display. The main attributes of interest are pH, color, firmness, water-holding capacity, fat content and composition, oxidative stability and uniformity [1].
Meat quality homogeneity is a major concern in the pig industry and market, but it is difficult to achieve by traditional selection because most meat quality traits exhibited low to moderate heritabilities [2,3] and measuring them is difficult, expensive, and only possible after slaughter. Fortunately, molecular technologies have played an important role in improving meat quality. Several major genes (such as RYR1, PRKAG3, IGF2) influencing meat quality have been applied in the pig industry, resulting in considerable improvement of meat quality in commercial pig herds [4,5].
In the past decades, quantitative trait loci (QTLs) in livestock have been detected mainly by using linkage mapping method with low-density microsatellite markers across the genome. Thus most QTLs generally span a large chromosomal region (comprising hundreds of genes), from which it is difficult to identify causative genes [6]. To date, 5,024 QTLs for meat quality traits have been deposited in pigQTLdb (http://www.animalgenome.org/cgi-bin/ QTLdb/SS/), but only a handful of causative variants have been identified via QTL fine mapping analysis. During the past few years, the emergence of more cost-effective and high-throughput genotyping platforms, SNP arrays, have rendered association mapping an increasingly popular and powerful approach for QTL mapping in human, animal and plant [7]. In pigs, there is an increasing number of association studies on commercial purebreds or F2 intercross populations to detect SNPs associated with monogenetic [8] and polygenetic traits, such as hematological traits [9,10], T lymphocyte subpopulations [11], body composition and structural soundness traits [12], boar taints [13,14], farrowing traits [15] and meat quality traits [16,17].
White Duroc is a lean-type western pig line and Erhualian is a Chinese fat-type indigenous pig line. They show obvious differences in meat productivity and quality, and are therefore genetically distant from each other. We have previously conducted genetic linkage analyses to detect QTLs for meat quality traits using a White Duroc6Erhualian F2 resource population [18,19].
Here, we carried out GWA analyses in both the F2 population and another population: Sutai pigs. The Sutai pig is a newly developed line which contains 50% Duroc and 50% Chinese Taihu breed (including Erhualian, Meishan and Fengjing strains) and have experienced selective breeding over 18 generations. Because the founder strains of Sutai pigs are close to those of the F2 population, the objectives of this study were not only to identify the precise locations of QTLs for meat quality traits in the two populations, but also to check the consistency of QTL findings across the populations.

Ethics Statement
All procedures involving animals followed the guidelines for the care and use of experimental animals approved by the State Council of the People's Republic of China. The ethics committee of Jiangxi Agricultural University specifically approved this study.

Study Populations and Traits
A three-generation resource population and a Sutai pig population were involved in this study. The former one was created and managed from 2001 to 2006 as described by Ren et al. (2006) [20]. Briefly, two White Duroc sires and 17 Erhualian dams were mated to produce F1 animals, from which nine F1 boars and 59 F1 sows were intercrossed (avoiding full-sib mating) to produce 967 F2 males and 945 F2 females (total n = 1912) in six batches. The Sutai population comprised offspring of four boars and 55 sows. All Sutai piglets were born and raised for 2-3 months at Sutai Pig Breeding Center in Suzhou city, and then they were transferred to a farm in Nanchang city (nearby the farm used for raising the F2) at three different times (July 2, Sep. 3 and Dec. 26,2011). Then they were fed with similar diet (formulated according to age) as that for the F2 animals under a standardized feeding and management regimen, and given free access to water. The F2 and Sutai piglets were weaned at 46 days and 28 days after birth, respectively. The castration was carried out for the F2 boars aged at 90 days and all Sutai piglets aged at 18 days including males and females.
At 24066 days of age, a total of 1030 F2 animals including 549 gilts and 481 barrows and a total of 436 Sutai pigs including 206 gilts and 230 barrows were slaughtered at a commercial abattoir. Meat quality measurements were performed on longissimus muscle (LM) between the 10th-rib and the first lumbar vertebra and semimembranosus muscle (SM) from left-side carcass, as described in detail at elsewhere [18,19,21,22]. The pH values were measured in the LM and SM by a Delta 320 pH Meter (Mettler Toledo, Greifensee, Switzerland) at 45 min and 24 h postmortem. Then, pH drop between the two time points was calculated. Meat color was subjectively assessed according to the color standard (1 = pale; 6 = dark) provided by the US National Pork Producers Council (NPPC) [23], and objectively evaluated using a CM-2600d/2500d Minolta Chroma Meter with parameters L* for lightness, a* for redness and b* for yellowness on the cut surface of the two muscles at 24 h postmortem. Drip loss after 24 h and 48 h storage of the LM and SM were measured using a bag method [24] and an EZ-DripLoss method [25]. Moisture, protein and intramuscular fat (IMF) contents of LM were determined by the routine oven-drying method, a Kjeldahl nitrogen method and an ether extraction method respectively [26]. Subjective marbling score of both muscles and firmness score of the LM were evaluated using NPPC standards [23,27]. For the LM of Sutai, the drip loss was not measured using the bag method and the crude protein content was also not determined. In the study, 933 F2 and 434 Sutai piglets were phenotyped. Descriptive statistics of the phenotype data related to 25 traits are given in Table 1.

Genotyping and Quality Control
Genomic DNA was isolated from ear clip or spleen tissues using a routine phenol/chloroform extraction method, and DNA concentration was diluted to 50 ng/ul. The quality and concentration of genomic DNA fulfilled the requirements for the Illumina Infinium SNP genotyping platform. Genotyping of 62,163 SNPs on the Illumina Porcine 60 K SNP Beadchip was carried out at the Illumina-certified service provider, Beijing Emei Tongde Technology Development Co. Ltd (EMTD). Genotypic data is available on all F2 and Sutai offsprings phenotyped, as well as their parents and/or grandparents. Quality control was carried out using PLINK v1.07 [28] for each population separately. SNP markers were removed if they had genotype-missing rates .0.03 or minor allele frequencies (MAF) ,0.05 or Hardy-Weinberg P, = 10 25 (based on Chi-squared test). Samples were removed on low (,90%) call rate. After quality control, all samples passed the filter and a final set of 39,414 SNPs and 44,532 SNPs was selected for GWA in the F2 and Sutai populations, respectively. The distribution of SNP markers after filtering and marker density on each chromosome are shown in Table S1. Genotype data are deposited in the Dryad repository (http://dx.doi.org/10.5061/ dryad.7 kn7r).

Statistical Analyses
The association analyses were conducted using GenABEL in the R software [29]. SNPs were individually tested for association with all studied traits using a generalized linear mixed model. The model includes a random polygenic effect for which the variancecovariance matrix is proportional to genome-wide identity-by-state (IBS). The model equation is shown below: y~1mzXbzkwzSczZaze where y is the vector of phenotypes of all genotyped and phenotyped F2 or Sutai piglets; m is the overall mean; b is the vector of fixed effects including sex and batch effects; w is the vector of slaughter weight of individuals considered as covariate; c is the vector of SNP effects with Erhualian allele substitute to White Duroc allele; a is the vector of random additive genetic effects with a,N(0, Gs a 2 ), where G is the genomic relationship matrix calculated from the corrected pedigree and s a 2 is the polygenetic additive variance); k is the regression coefficient of slaughter weight and e is the vector of residual errors with e,N(0, Is e 2 ), where I is the identity matrix and s e 2 is the residual variance. X, S and Z are incidence matrices for b, w and c respectively. The herd-year-season effect was contained in the batch effect.
The genome-wide significance threshold was determined by the Bonferroni method, in which the conventional P-value was divided by the number of tests performed [30]. A SNP was considered to have genome-wide significance at P,0.05/N and have chromosome-wide significance at P,1/N, where N is the number of SNPs tested in the analyses. The genome-wide and chromosome-wide significant thresholds were 1.27e-6 (0.05/39414) and 2.54e-5 (1/ 39414) respectively for the F2 population, and were 1.12e-6 (0.05/ 44532) and 2.25e-5 (1/44532) respectively for the Sutai population.
The influence of population stratification was assessed by examining the distribution of test statistics generated from the thousands of association tests and assessing their deviation from the null distribution (that expected under the null hypothesis of no SNP associated with the trait) in a quantile-quantile (Q-Q) plot [6]. In these plots ( Figure 1B and Figure S2), -log 10 P values for each SNP calculated from their observed association statistics (x 2 statistics) were ranked in order from smallest to largest on the yaxis and plotted against the distribution that would be expected under null hypothesis of no association on x-axis. Deviations from the diagonal identity line suggest that either the assumed distribution is incorrect or that the sample contains values arising in some other manner, as by a true association [31]. The Q-Q plot was constructed using R software.
Haplotype or linkage disequilibrium (LD) block analyses were performed for the chromosomal regions with multiple significant SNPs clustered around the peak SNP. The LD blocks were determined using Haploview version 4.2 software with default settings [32].

Population Stratification Assessment
Population stratification for GWAS can lead to false positive results [6]. The Q-Q plots of the test statistics in GWA are shown in Figures 1B and S2. From these plots, it is apparent that there is no clear overall systematic bias in all studied traits. The genomic inflation factors (l) observed in the GWA study were usually less than 1.10, also indicating that no very strong stratification existed.

GWAS Analyses
Both genome-wide significant SNPs and chromosome-wide significant SNPs for the pH, meat color, drip loss, chemical compositions, marbling and firmness are presented in Tables 2, 3, 4, 5. The profiles of the P-values of the tested SNPs for all meat quality traits are shown in Figure 1A and Figure S1. In total, 127 chromosome-wide significant SNPs were identified and among them, 11 showed genome-wise significant association (with underlined P-value in the tables) with different traits: 1 for pH, 4 for meat color, 4 for drip loss and 2 for IMF.
pH values. Five and seven SNPs significantly associated with pH traits were identified in the F2 and Sutai pigs respectively ( Table 2). All the SNPs except for unmapped markers represent five QTL regions on SSC2, 3, 4, 13 and X. The QTL region on SSC3 was common to the two populations. This region from 14. Meat color. We identified 8 and 16 significant SNPs associated with meat color in the F2 population and the Sutai population respectively (Table 3). No common QTL region for the same trait was detected in the two populations. However, SNP ALGA0039930 at 31.27 Mb on SSC7 that was associated with Minolta L* of SM in the F2 population is adjacent to another SNP ALGA0040423 at 37.73 Mb that showed significant association with Minolta a* of SM in the Sutai population. The most significant SNP associated with Minolta a* of both LM and SM in Sutai was the SNP ALGA0060775. This SNP reached genomewide significance level and was located very close to the other two genome-wide significant SNPs ASGA0049740 and MIGA0014909 for the same trait on chromosome 11.
Drip loss. In the F2 population, a total of 12 SNPs were detected to be significantly associated with drip loss of LM and SM after 24 h storage (Table 4). Ten out of these SNPs fall in the region of 1.53 Mb (from 1.31 Mb to 2.84 Mb) on SSC1, and the other two were located at 81.56 Mb and 81.63 Mb on SSC4. As for the Sutai population, there were 30 significant SNPs with effect on drip loss, out of which 3 on SSC15 and 1 on SSC3 reached genome-wide significance level for drip loss of SM after 24 h storage (Fig. 1A). The three most significant SNPs ALGA0086325 (P = 6.74E-08), ALGA0086324 (P = 7.64E-07) and ALGA0110636 (P = 1.03E-06) on SSC15 were in a haplotype block spanning 178 kb (Fig. 1C).
Moisture, protein and IMF contents, marbling and firmness scores. Forty-four SNPs were significantly associated with these traits: 7 for moisture content of LM, 2 for protein content of LM, 32 for IMF of LM, 1 for marbling of SM and 2 for firmness of LM (Table 5). In the F2 population, a 0.46-Mb region from 34.80 Mb and 35.26 Mb on SSC7 contains not only 6 SNPs associated with moisture content of LM, but also 2 SNPs associated with protein content and 2 SNPs associated with firmness of LM. Of the 32 SNPs associated with IMF of LM, 7 were detected in the F2 animals with the most significant SNP MARC0090296 on SSCX (P = 8.92E-07), and 25 in the Sutai pigs with the top SNP ALGA0053636 on SSC9 (P = 1.12E-06). Neither common loci for IMF nor significant SNPs associated with  The White Duroc6Erhualian F2 population and Sutai (ST) population. 2 Description of the traits is available in Table 1. 3 The number of significant SNPs within the QTL regions. 4,5 SNPs position on the Sus Scrofa Build 10.2 assembly. 6 Gene names starting with ENSSSCG represent Ensembl nomenclature while other gene symbols represent HUGO nomenclature. 7,8,9,10 The SNP allele ''A'' frequencies of two F0 Duroc (FA_D), 17 F0 Erhualian (FA_E), the whole F2 population (FA_F2) and Sutai population (FA_ST). 11 Additive effects; positive value indicates that allele ''A'' increased the trait. 12 Genome-wide significant associations are underlined. doi:10.1371/journal.pone.0064047.t002 marbling of LM were found in the two populations. Only one SNP MARC0090739 on SSC13 showed a significant association with marbling of SM in the F2 population.

Discussion
To our knowledge, only one study [16] has applied GWA approach to detect QTL signals for IMF, marbling, meat color and moisture in a Large White6Minzhu F2 population. This article reported that most significant SNPs (except for unmapped SNPs) for these traits were located within a 10.70 Mb region (51.37-61.07 Mb) on SSC12. In this region, we also identified a chromosome-wide significant SNP ALGA0067119 at 58.08 Mb for IMF of LM. The favorable allele (G) that increases IMF derived from Erhualian (Table 5). Whereas our results did not confirm the associations between this region and other phenotypes, and demonstrated that generally more than one genomic region are associated with meat quality traits.

Possible Pleiotropic QTLs
The present results showed that several regions contain multiple significant SNPs associated with different traits. Especially, the SSC7 region from 31.27 Mb to 37.74 Mb harbored SNPs affecting five traits: MARC0069646 for color parameter a*, ALGA0039930 for color parameter L*, MARC0033464 for moisture content, MARC0058766 for protein content and firmness. Our previous QTL mapping study [19] also demonstrated that this region have strong QTL effects on various carcass and meat quality traits measured in the F2 population. So the current GWA result is consistent with the result of linkage analysis. Moreover, the GWA study enhanced the precision of QTL mapping. For example, all 6 significant SNPs associated with moisture content fell into a 0.46 Mb region (34.80-35.26 Mb) on SSC7, much smaller than previously reported QTL interval of 12 cM (approximately 12 Mb).
Additionally, a 3.01 Mb region on SSC3 (from 14.40 Mb to 17.41 Mb) was found to be associated with both pH values (pH 45 min, pH 24 h and pH drop from postmortem 45 min to 24 h) and drip loss. Because the development of drip loss is largely governed by the rate and extent of postmortem pH decline [53], it is likely that there is a common causative variant for these related traits within the region. Similarly, a common SNP DRGA005419 on SSC5 is associated with both L* value of SM and DripEZ_48 h of LM in the Sutai piglets. Combined with the correlation coefficient of 0.45 (significantly greater than zero, P,0.01) between the two traits, it suggests the existence of a pleiotropic QTL simultaneously regulating meat color and drip loss.

GWA QTLs vs. Linkage Mapping QTLs
Previously, a genome-wide significant QTL for IMF was mapped to a region flanked by microsatellite makers SW2456 and S1426 (48-58 cM and 42-103 Mb) on SSCX in the F2 population [19]. This region has a very low recombination rate (average 6 Mb per cM) [54], making it very difficult to fine-map the QTL and to discriminate between multiple QTLs and single QTL by family-based linkage analysis. Fortunately, it is not a big challenge in GWA studies because it can capitalize on all meiotic recombination events in a population, rather than only those occurred currently in the studied families. It is, therefore, not surprising that the association signals for IMF were localized to two distant and small regions on SSCX in this study: one harboring 2 significant SNPs at 45.39 Mb and 46.12 Mb respectively, another harboring 3 significant SNPs from 103.62 Mb to 104.43 Mb. Moreover, the two regions also tended to be associated with marbling score (SNP MARC0090296 at 46.12 Mb with P-value of 2.74E-05 and SNP H3GA0051863 at 104.42 Mb with P-value 1.65E-04, approaching significance level), because IMF and marbling score are highly correlated (r = 0.71, P,0.01).
There are some differences in QTL findings between the present GWA study and the previously published genetic linkage studies using the same F2 population. Several 1% genome-wide significant QTLs that were reported in our previous papers [18,19] failed to replicate in this study, such as two QTLs for IMF of LM on SSC9 and one for color score of SM on SSC11. Such discrepancy maybe due to the following reasons: (1) Linkage analyses were performed under an assumption that the founder lines are fixed for different QTL alleles, whereas GWA analyses systematically investigate SNPs in the entire genome without the constrains of a priori hypotheses; (2) Additive, dominant and even imprinting effects of putative QTL were estimated in the linkage study, whereas only additive effect was tested in this GWA study; (3) We used the QTL linkage mapping procedure suggested by Guo et al. (2008) [55]. This procedure obtained estimates by fitting all identified QTLs as genetic background effects at each step of searching new QTL. In contrast, the linear mixed model was fixed in the GWA study; (4) We used a permutation method and a Bonferroni method to determine the significance thresholds for linkage mapping analysis and GWA analysis respectively. Compared with the permutation method, the Bonferroni correction method operates too conservative, because it assumes the independence of each test even though many of the SNPs are in linkage disequilibrium (LD) and thus correlated with each other. As a result, the Bonferroni power to detect some statistically significant results became relatively weak.

Common and Population-specific QTLs
Repeated detection of a QTL among populations is a way to validate the QTL. Interestingly, the GWA analyses of the Sutai population revealed some genome-wide significant SNPs for IMF and color a* on SSC9 and 11 respectively. They situate within the above-mentioned QTL regions detected in the F2 population. This result thus provided evidence that those genome-wide significant QTLs found in the F2 population are unlikely to be artifacts of linkage analyses.
In addition, several association signals, e.g. those for pH values on SSC3 and for color parameters on SSC7 were repeatly identified in the two populations, validating the existence of these loci. However, many association signals appeared in only one population. This maybe resulted from the differences in environmental background (such as birthplaces, times of weaning and castration, and etc.) and genetic background (because of founder lines, population structure, selection, gene-gene interactions, and etc.) between the two populations.

Candidate Genes
We noticed that the significant SNPs for pH, meat color, drip loss, moisture content and protein content are rarely situated within or near known genes affecting these traits. Only one SNP DRGA0005419 for drip loss on SSC5 was located 436 kb upstream of the ADSL (adenylosuccinate lyase) gene that was found to be possibly associated with drip loss and pH 45 min of LM in Pietrain pigs [56].
In contrast, according to gene biological functions in lipid metabolism, adipocyte and/or muscle development, we identified some candidate genes for IMF, marbling and firmness within 500 kb upstream/downstream of the peak SNPs. In the F2 population, four GWA QTLs for IMF were detected on SSC7, SSC12 and SSCX. The FOS gene is closest to the significant SNP ALGA0043983 on SSC7. This gene encodes a leucine zipper protein that has been implicated as a regulator of signal transduction, cell proliferation and differentiation (e.g. myogenesis) [57]. Furthermore, it was found to be expressed differentially in muscle between the fat type pig breeds (such as Basque and Liangtang) and lean type pig breeds (Large White and Landrace) [58]. Therefore, the FOS gene could be regarded as a prime candidate gene for the QTL. Within the QTL region on SSC12, the MYH1, MYH2 and MYH3 genes that belong to the myosin heavy chain gene family (MYH) have been proposed as candidate genes by , who also identified this QTL in their GWA study. On SSCX, a genome-wide significant SNP MARC0090296 for IMF is located at 46.12 Mb, within the SLC9A7 (Solute carrier family 9 member 7) gene. A promising gene, RGN (also called as SMP30, i.e. regucalcin or senescence marker protein-30), is located 438 kb away from this SNP. Regucalcin plays a multifunctional role as a regulatory protein in intracellular signaling processes in many cell types and is related to lipid metabolism [59]. Regucalcin transgenic rats have been shown to experience hyperlipidemia with increasing age [60]. No apparent candidate genes are located in the vicinity of the SNP ALGA0099852 at 103.62 Mb on SSCX.
In the Sutai population, we found four QTLs for IMF on SSC1, SSC8 and SSC9. A candidate gene for the QTL on SSC1 is the ATG14 (autophagy related 14 homolog) gene that plays an important role in hepatic lipid metabolism [61]. The QTL effect detected on SSC8 could be due to the candidate gene BMPR1B [bone morphogenetic protein (BMP) receptor, type IB], because the ligands of this repceptor is BMPs that can induce commitment of C3H10T1/2 pluripotent stem cells into adipocytes [62,63]. The STEAP4 gene encoding metalloreductase, which is associated with obesity and insulin-resistance in human [64][65][66], is located at 74.98 Mb on SSC9, very close to the strongest association signal (ALGA0053636) detected in the Sutai piglets, and thereby is an excellent positional and biological candidate gene for this QTL. No obvious candidate genes for IMF were found in the distal region (around 152.11 Mb) of SSC9.
The peak SNP MARC0090739 for marbling score is located only 55 kb from the UBASH3A gene, which has a role in immune function and was observed to be differentially methylated in peripheral blood leukocytes between lean and obese adolescents [67]. On SSC7, the peak SNP MARC0058766 (at 34.80 Mb) for the firmness and moisture content, is also significantly associated with the protein content and drip loss of LM in the F2 population. The SNP is located between two candidate genes: the LEM2 gene (at 34.64 Mb) and the HMGA1 gene (at 34.98 Mb). The LEM2 (also called NET25) gene is involved in nuclear structure organization and its mutations cause muscular dystrophies and other disorders [68]. The HMGA1 gene encodes high mobility group AT-hook 1 protein that may play critical role in adipogenesis [69] and serve as a modulator of IGF-I activity [70]. The significant associations between polymorphisms in this gene and backfat thickness as well as drip loss have been reported [71,72].

Conclusions
In summary, this GWA study identified 11 genome-wise significant SNPs and 116 chromosome-wide significant SNPs associated with 25 meat quality traits. Our results narrow down the previously detected QTL intervals, and reveal 7 new QTL positions. At least two QTL regions on SSC3 and SSC7 were found to affect multiple traits and are common to the two populations. However, many QTLs are not conserved across the two populations, reflecting the genetic heterogeneity of these QTLs and the complexity of the genetic basis of pork quality. For some traits including pH values, drip loss and firmness, it is the first time that they are included in a GWA analysis. In the QTL regions, some candidate genes stand out because of their functional annotations, positions and reported expression variation in related tissues. The current findings will contribute to further identification of the causal mutation underlying these QTLs and future improvement of meat quality in pig breeding programmes.