A Splice Mutation in the PHKG1 Gene Causes High Glycogen Content and Low Meat Quality in Pig Skeletal Muscle

Glycolytic potential (GP) in skeletal muscle is economically important in the pig industry because of its effect on pork processing yield. We have previously mapped a major quantitative trait loci (QTL) for GP on chromosome 3 in a White Duroc × Erhualian F2 intercross. We herein performed a systems genetic analysis to identify the causal variant underlying the phenotype QTL (pQTL). We first conducted genome-wide association analyses in the F2 intercross and an F19 Sutai pig population. The QTL was then refined to an 180-kb interval based on the 2-LOD drop method. We then performed expression QTL (eQTL) mapping using muscle transcriptome data from 497 F2 animals. Within the QTL interval, only one gene (PHKG1) has a cis-eQTL that was colocolizated with pQTL peaked at the same SNP. The PHKG1 gene encodes a catalytic subunit of the phosphorylase kinase (PhK), which functions in the cascade activation of glycogen breakdown. Deep sequencing of PHKG1 revealed a point mutation (C>A) in a splice acceptor site of intron 9, resulting in a 32-bp deletion in the open reading frame and generating a premature stop codon. The aberrant transcript induces nonsense-mediated decay, leading to lower protein level and weaker enzymatic activity in affected animals. The mutation causes an increase of 43% in GP and a decrease of>20% in water-holding capacity of pork. These effects were consistent across the F2 and Sutai populations, as well as Duroc × (Landrace × Yorkshire) hybrid pigs. The unfavorable allele exists predominantly in Duroc-derived pigs. The findings provide new insights into understanding risk factors affecting glucose metabolism, and would greatly contribute to the genetic improvement of meat quality in Duroc related pigs.


Introduction
In past decades, thousands of quantitative trait loci (QTLs) have been detected for economically important traits in livestock through genetic linkage studies [1]. The recent availability of livestock genome sequences and high-density SNP chips has allowed detection of significant association of nucleotide polymorphisms with complex traits, and effective identification of causal mutations for some monogenic traits [2][3][4]. Despite these progresses, uncovering the quantitative trait genes (QTGs) or nucleotides (QTNs) for complex traits remains a challenging task. Only a handful of QTNs have been convincingly identified in livestock [5][6][7].
More recently, several studies in human, mouse and Drosophila have shown that the integration of phenotypic traits, genetic and gene transcript data enables researchers to identify expression QTL (eQTL), untangle gene-based regulatory networks, infer relationship between gene expression levels and phenotypic traits, and thereby detect novel trait-causing genes [8][9][10]. Therefore, the integrative analysis would play a key role in acquiring the knowledge of mechanisms responsible for livestock complex traits. Indeed, some researchers have identified a number of promising candidate genes for muscle traits and blood lipid traits of pig by integrative analyses of GWAS (or linkage mapping), eQTL and trait-correlated expression [11][12][13][14].
Glycogen storage diseases (GSD) characterized by defects in glycogen metabolism and excess glycogen stored in liver and muscle are multifactorial disorders. Both human and animals suffer from GSD. Agricultural researchers often measure the glycogen content or glycolytic potential [GP = 26 (glucose + glycogen + glucose-6-phosphate) + lactate] in skeletal muscle of farm animals at slaughter. However, instead of the GSD diagnostic indicator, the GP value is mainly used for the prediction of meat quality development during the conversion of muscle to meat, because GP is a determinant of multiple meat quality characteristics, such as pH, color, water holding capacity (drip loss), tenderness and processing yields [15,16].
GSD or extreme GP phenotype can be caused by genetic variation of various enzymes or transporters, which are involved either directly in the synthesis or breakdown of glycogen or in the utilization of its catabolite, glucose-1-phosphate [17]. To date, only one responsible gene, PRKAG3, has been elaboratively evidenced to impact glycogen and GP levels in Hampshire and its related synthetic lines [18,19]. Although more than 10 different QTLs for the trait have been reported (http://cn.animalgenome.org/cgi-bin/QTLdb/ index), QTGs or QTNs underlying these QTLs remain unexplored.
We have previously identified a major QTL for GP and pH values at 42 cM on chromosome 3 in a large scale White Duroc 6 Erhualian F 2 intercross [20]. To decipher the molecular basis of the GP QTL, we herein performed an integrative genetic analysis by utilizing extensive data sets of GP-related traits, high-density genotypes and gene expression profiling from three experimental populations. We show compelling evidence that a splice mutation in the PHKG1 gene is the QTN underlying the major QTL effect on GP and its related traits.

Results
Genome-wide association studies on two experimental populations define the QTL to a ,180-kb region We first performed a genome-wide association study (GWAS) for GP and its components including residual glycogen & glucose (RG), glucose-6-phosphate (G-6-P) and lactate on the White Duroc 6 Erhualian F 2 intercross, in which 877 phenotyped F 2 individuals and its parents/grandparents were genotyped by using Illumina Porcine SNP 60K Beadchips [21]. The quality control filtering and data processing of the GWAS data are described in the Methods. Quantile-quantile plots with genome control l GC values are shown in Figure S1. We found no evidence of systematic inflation of association test results. Consistent with our previous QTL results [20], the GWAS results demonstrated that Sus Scrofa chromosome 3 (SSC3) contained the most genome-wide significant locus (P = 5.85610 222 ) for both GP and RG with the top SNP ss131031160 (pSNP) at 17.09 Mb ( Figure 1). However, this lead SNP was not associated with the glycolysis intermediate product G-6-P and end-product lactate ( Figure S2). These results suggest that the underlying QTG likely influences the conversion between glycogen and glucose (or G-6-P) rather than the conversion between G-6-P to lactate. The pSNP ss131031160 explains 19.6% of the phenotypic variance in RG of the F 2 population.
To validate and fine map the SSC3 QTL for RG, we conducted a second GWAS on 433 Sutai pigs, a Chinese synthetic line derived from a cross between Duroc and Taihu including Erhualian and Meishan after more than twenty-year selection. Besides 62163 SNPs on the Illumina Porcine SNP 60K Beadchip, 53 SNPs derived from DNA sequence comparison between Erhualian and Duroc (see Materials and Methods) within a 1-Mb region surrounding the top pSNP were genotyped for these Sutai pigs. As expected, we confirmed the strong association signal on SSC3 in the Sutai population. The pSNP ss131565361 (P = 9.06610 219 ) at 16.92 Mb ( Figure 1) explains 53.6% of the RG variance. This locus also significantly influenced ultimate pH (measured 24 hours after slaughter) and drip loss ( Figure S2). Based on the LOD drop off 2, the empirical confidence intervals of the QTL in the F 2 and Sutai populations was 920 kb (16.92-17.84 Mb) and 630 kb (16.47-17.10 Mb) respectively ( Figure 2A). Therefore, the most likely QTL interval was their overlapping region of 180 kb . This interval contains 7 annotated genes: GUSB, VKORC1L1, NUPR1L, CHCHD2, PHKG1, SUMF2 and CCT6A (Figure 2A). GWAS results for glycolytic potential and residual glycogen contents in longissimus muscle from the F 2 and Sutai populations. In the Manhattan plots, log 10 (1/P) values of the qualified SNPs were plotted against their genomic positions. ''Chromosome 0'' harbors SNPs that have not yet been mapped to the Sus Scrofa 10.2 assembly. It is most likely that the significant SNPs on chromosome 0 are located on chromosome 3 and have strong linkage disequilibrium with the top SNPs ss131031160 or ss131565361. The red and green dots represent the SNPs that reached 5% genome-wide and suggestive Bonferroni-corrected significances, respectively. The most significant SNPs for residual glycogen are ss131031160 and ss131565361 on SSC3 in the F 2 and Sutai (ST) pigs respectively. GP, glycolytic potential. doi:10.1371/journal.pgen.1004710.g001

Author Summary
Glycogen storage diseases (GSD) are a group of inherited disorders characterized by storage of excess glycogen, which are mainly caused by the abnormality of a particular enzyme essential for releasing glucose from glycogen. GSD-like conditions have been described in a wide variety of species. Pigs are a valuable model for the study of human GSD. Moreover, pigs affected by GSD usually produce inferior pork with a lower ultimate pH (so-called ''acid meat'') and less processing yield due to post-mortem degradation of the excess glycogen. So far, only one causal variant, PRKAG3 R225Q, has been identified for GSD in pigs.
Here we reported a loss-of-function mutation in the PHKG1 gene that causes the deficiency of the glycogen breakdown, consequently leading to GSD and acid meat in Duroc-sired pigs. Eliminating the undesirable mutation from the breeding stock by a diagnostic DNA test will greatly reduce the incidence of GSD and significantly improve pork quality and productivity in the pig. eQTL mapping supports PHKG1 as the most likely responsible gene Using the data of digital gene expression profiles (DGE) tested in longissimus muscle samples from 497 F 2 animals, we identified genome-wide eQTL for the 7 annotated genes in the critical region. We operationally defined cis-eQTL as any association between a gene expression and SNPs within 2 Mb of the gene location. We found 3 cis-eQTL each for PHKG1, SUMF2 and GUSB ( Figure 2B). Intriguingly, the lead SNP for the cis-eQTL affecting PHKG1 expression was identical to the top pSNP (ss131031160) for GP in the GWAS on the F 2 population. The major allele (G) at this SNP was associated with lower residual glycogen and higher expression of PHKG1 (Table 1 and Figure 2C). PHKG1 encodes a catalytic subunit of the phosphorylase kinase (PhK), which can mediate glycogen breakdown. Therefore, the co-localization of the top SNP in both GWAS and eQTL mapping coupled with the biological function highlights PHKG1 as the most likely QTG underlying the major locus on SSC3.
Sequencing the PHKG1 gene identifies a splice mutation as candidate QTN To search for coding variants in the PHKG1 gene, we isolated and sequenced PHKG1 cDNA from muscle samples of 6 F 2 individuals representing three genotypes at the pSNP ss131031160 (2 GG, 2 GA and 2 AA animals). In total, we detected 14 polymorphisms (Table S1). Of note, we found a 32-bp deletion/ insertion (c.del/ins32; Figure 3A&B) polymorphism at exon10 that alters the open-reading frame and causes a premature stop codon, leading to a truncated and nonfunctional protein product.
To test if the 32-bp deletion was also present in the genomic DNA (gDNA), we sequenced the 140 bp genomic region encompassing the deletion using gDNA samples of the above-mentioned 6 F 2 animals. We found only one single nucleotide substitution (C.A) in the region. This SNP (g.8283C.A) lies 5 bp upstream of the c.del/ins32 ( Figure 3C&D), where the pyrimidine (C or T) is highly conserved among vertebrates, including chicken, mouse, dog, cow, sheep and human. The CRA nucleotide transvertion may weaken the strength of the pyrimidine-rich signal near the 39 end of the PHKG1 intron 9, which is known to correlate with splicing efficiency [22]. Using a web-available tool Alternative Splice Site Predictor (ASSP, http://wangcomputing.com/assp/ index.html), we found that the splice site score (1.652) of the mutant sequence was much lower than the wild-type counterpart (3.659), and even lower than the default cutoff value (2.2) for acceptor sites set in the predictor. Therefore, the mutation appears to attenuate or inhibit the role of the constitutive acceptor splice site and this consequently evokes another cryptic acceptor at 32 bp downstream of the mutation, which results in a loss of 32 bp in mature mRNA during PHKG1 transcription ( Figure 3E).
To test if the 32-bp deletion in exon 10 was directly caused by the SNP g.8283C.A, the effect of the variant on splicing was assessed using a minigene splicing assay ( Figure S3). Two minigene expression vectors carrying either the wild-type or the mutant PHKG1 g. 8283C.A segment were constructed and transiently transfected into HeLa cells and 293T cells ( Figure S3A). The minigene transcripts in these transfected cells were analyzed by RT-PCR, using specific primers complementary to exon 9 and exon 10 of the minigene. We found two amplicons of 78 bp and 62 (C) Box plots showing the differences in gene expression level determined by the digital gene expression (DGE) system among three genotypes of the ss131031160 SNP. Different letters above each box plot denote significant differences (P,0.05). doi:10.1371/journal.pgen.1004710.g002 bp from the mutant construct and only one amplicon of 94 bp from the wild-type construct ( Figure S3B). Sanger sequencing analysis revealed that the two truncated amplicons corresponded to two aberrant splicings of the first 16 and 32 nucleotides (nt) at 59 end of exon 10 ( Figure S3C-E). The presence of novel aberrant splicing of 16-nt may be due to that the minigene constructs only contain partial PHKG1 sequence or there is a special splicing factor in the transfected cells. Anyway, the result clearly demonstrates that the variant PHKG1 g. 8283C.A is responsible for the aberrant splicing of 32-nt observed in vivo.

The candidate QTN co-segregates with QTL genotypes of parental boars in the two experimental populations
To obtain further evidence for the causality of the candidate QTN (SNP g.8283C.A), we conducted the concordance test between the SNP genotypes and the QTL genotypes on 9 parental boars (6 F 1 boars and 3 Sutai boars). The QTL genotypes of these boars were determined by the marker-assisted segregation analysis (see Materials and Methods). We found that the SNP genotypes were completely concordant with the QTL genotypes in these boars (Table S2). In contrast, the two top pSNPs ss131031160 and ss131565361 in the GWAS exhibited disconcordance with the QTL genotypes in the 9 parental sires. Concordantly, the association of the g.8283C.A SNP with RG phenotype was strongest among all SNPs genotyped in Sutai pigs, and was at the same significance level as that of the top GWAS SNP ss131031160 in the F 2 population due to their complete linkage disequilibrium (Table S3).
In addition, we re-sequenced a 10-kb segment covering 1 kb upstream of the PHKG1 gene to its 39 end on all 9 sires (6 F 1 and 3 Sutai sires) with deduced QTL status. A total of 142 variants were identified (Table S4). We found that 7 Q chromosomes of these boars shared the same haplotype while 11 q chromosomes corresponded to multiple divergent haplotypes ( Figure S4). Surprisingly, a q haplotype from Sutai boars was nearly identical to the Q haplotypes except for three variants, of which only one (g.8283C.A) showed co-segregation with QTL genotypes across all the 9 sires ( Figure S5). This strongly supports the causality of the g.8283C.A mutation.
Haplotype analysis showed that no other variants within the PHKG1 region were in complete LD with the SNP g.8283C.A across all experimental populations. Based on that, we tried to determine whether other variants can also result in partial QTL effect. When the SNP g.8283C.A was included as a fixed effect in the model, the QTL effects on PHGK1 expression and RG trait vanished in these populations ( Figure S6A). Furthermore, we divided the Sutai pigs into three genotype groups (AA, CA, CC). No eQTL effect was detected within each group ( Figure S6B). It thus suggests that even there is another regulatory SNPs, their effects are likely negligible.
The QTN significantly affects PHKG1 mRNA and protein levels Aberrant mRNA, like the truncated PHKG1 transcript caused by the g.8283C.A SNP, tend to be degraded by a known mechanism of nonsense-mediated mRNA decay (NMD) in the cell. To analyze whether the mutation interferes with PHKG1 mRNA expression levels, qRT-PCR experiments were performed on total RNA isolated from 293T cells transiently transfected by the wildtype and mutant PHKG1 minigenes. The mutant minigene produced nonsense mRNAs at 56% level of normal expression ( Figure  S7). The finding indicates that the PHKG1-32del mRNA bearing the premature termination codon was most likely degraded by NMD.
In fact, electropherogram of RT-PCR products amplified from PHKG1 mRNA in homozygous animals showed that the intensity of the mutant DNA band of 114 bp was much lower than the normal band of 146 bp ( Figure 3A), implying NMD of the mutant PHKG1 mRNA. Accordingly, the PHKG1 cDNA sequence chromatograms illustrated that the fluorescence intensity of wild-type allele (Wt or q) were 2-3 fold stronger than that of the mutant allele (Mt or Q) in animals heterozygous for the g.8283C.A SNP ( Figure 3F). To more accurately estimate the difference in abundance between Wt and Mt transcripts, we designed two pairs of primers: one (Common-59-FP/ RP) for amplification of both transcripts, and the other (Wt-39-FP/ RP) for the specific amplification of Wt transcript ( Figure 4A). After performing qRT-PCR with the two primer sets, we quantified the levels of total mRNA transcripts relative to wild-type transcripts. In heterozygotes, the ratio of (Mt+Wt) to Wt transcripts was 1.4:1 ( Figure 4B), suggesting that the 60% of Mt transcripts were degraded by NMD.
Interestingly, we found that mutant homozygotes (AA or QQ) had PHKG1 Wt, as RT-PCRs using the primers specific for the Wt generated amplicons corresponding to Wt in these mutant homozygotes ( Figure 4C). We further showed that, only about one-eighth of PHKG1 transcripts was Wt ( Figure 4B). The result suggests that the g.8283C.A mutation could strongly decrease but not completely abrogate the original splicing form.
Consistently, the PHKG1 Wt transcript level is highest in muscle samples from CC animals, followed by AC and AA individuals ( Figure 4C). Similar tendency was observed for the PHKG1 Wt protein level in the three genotypes by Western-blot analysis ( Figure 4D). Thus, it is evident that the g.8283C.A mutation affects the expression of PHKG1 at both transcription and translation levels. PHKG1 mRNA level is significantly correlated with the RG phenotype PHKG1 is a critical enzyme in the glycogen metabolism. We expected that the PHKG1 mRNA expression level would be significantly associated with the RG phenotype. However, we did not observe such strong correlation (r = 20.08, P = 0.106; Figure 5B) using the DGE data in the F 2 population. We then made a close examination on the DGE data. We found that the DGE tag for PHKG1 was located in a region upstream of the 32bp deletion, which therefore can not distinguish Wt from Mt transcripts. To correct for the biased effect of the abnormal Mt transcripts, we further used qRT-PCR to measure the Wt level on 117 F 2 individuals and 104 Sutai pigs, of which the number of CC, CA and AA animals at the g.8283C.A locus are nearly identical ( Figure 5A). We then examined the correlation between the qPCR-data and the RG phenotypes. As a result, the PHKG1 Wt (functional) transcript level was significantly correlated with RG phenotype (r#20.4, P,10 25 ) in the two populations ( Figure 5B). Such strong correlation were not detected for any other cisregulated genes within the critical region of the SSC3 QTL ( Figure   Figure 3. Identification of the g.8283C.A splice mutation in the PHKG1 gene. (A) RT-PCR products of the PHKG1 gene on a 2.5% agarose gel, amplified from total RNA from animals with different QTL genotypes (QQ, Qq and qq). Two DNA bands of 146 bp and 114 bp with different intensity were observed. (B) Sequence analysis of the RT-PCR products. A 32-bp deletion in exon 10 was detected in QQ and Qq animals. (C) PCR products of the PHKG1 gene on a 2.5% agarose gel, amplified from genomic DNA from animals with different QTL genotypes (QQ, Qq and qq). (D) Sequence analysis of the PCR products. A point mutation (g. 8283C.A) was identified in QQ and Qq individuals. (E) PHKG1 genomic sequence from exon 9 to exon 10. Intronic sequence is in black and lower case, and exonic sequence in blue and upper case. The normal splice-acceptor and -donor sequences (lower case) and the cryptic splice-donor sequences (upper case) are underlined with two red nucleotide letters in the center. The C.A mutation is denoted as ''M'' in red. (F) Sequence analysis of a heterozygote at the g.8283C.A mutation site. In genomic DNA from a heterozygote, the C and A peaks (the position is indicated by ''M'' and marked by the red bar) were of equal height in sequence of genomic DNA, as expected. Sequence of the cDNA product from the same heterozygote showed a major peak corresponding to the wild-type sequence and a small peak corresponding to the mutant sequence (the +2 position of exon 10 where the nucleotide T belong to the 32-bp deleted portion is marked by the blue bar). doi:10.1371/journal.pgen.1004710.g003 S8), strengthening the causal relationship between PHKG1 and RG.
The PHKG1 g.8283C.A mutation significantly alters the enzymatic activity of the phosphorylase kinase PHKG1 is a catalytic subunit of the phosphorylase kinase (PhK). The g.8283C.A splice mutation in the PHKG1 gene thus likely influence the enzymatic activity of PhK and glycogen phosphorylase. To verify the assumption, we tested the enzymatic activity of PhK in 18 muscle samples_6 per genotype group (CC, AC, and AA). As expected, the AA group showed.6-fold reduction in enzyme activity as compared to the CC and AC groups ( Figure 6). It is reasonable to speculate that reduced PhK enzyme's activity would lead to slow breakdown of glycogen in muscle, which in turn cause excess accumulation of glycogen in the tissue. This is consistent with the SSC3 QTL effect on GP or RG trait.
The PHKG1 g.8283C.A mutation is significantly associated with GP-related meat quality traits in multiple Duroc-derived populations Considering the inter-relationship of GP with other meat quality attributes, we evaluated the effects of the PHKG1 g.8283C.A mutation on all meat quality traits measured in the 864 F 2 individuals from the White Duroc 6 Erhualian intercross, 431 Sutai pigs and 140 three-way hybrid DLY [Duroc 6 (Landrace 6 Yorkshire)] pigs. In the F 2 population, there was a tendency towards lower ultimate pH and higher drip loss in AA pigs with higher GP compared to CA and CC individuals (the mean values of pH 24h in AA, CA and CC groups were 5.70, 5.72 and 5.77, respectively; for drip loss: 1.01%, 0.85% and 0.92%, respectively), although the mean differences did not reach statistical significance, partially due to a low frequency (5%) of the AA genotype in the population.
In Sutai pigs, each genotype group (AA, AC or CC) at the PHKG1 g.8283C.A site has more than 100 individuals. Significant effects of the causal mutation were observed on almost all meat quality traits in this population (Tables 2 and S5). Especially, the RG level in the longissimus muscle (LM) of AA animals was 4 times higher (P = 9.03610 254 ) than those in CC animals. Coincidently, drip loss, rate of pH decline and Minolta a* and b* were higher (P, 0.01) in AA pigs than CC pigs. AA pigs also had lower (P,0.001) intramuscular fat content (IMF) and marbling score compared to CC pigs. In the DLY hybrid commercial population, all 140 individuals were genotyped for the PHKG1 g.8283C.A mutation and 53 surrounding SNPs using the OpenArray platform (see Materials and Methods). Again, we observed the significant effects of the g.8283C.A mutation on RG, pH and drip loss ( Table 2). Moreover, this mutation showed stronger associations with these traits than any other surrounding SNPs ( Figure S9). These results indicate that the A allele has unfavorable effects on these meat quality traits.
To reveal the allele frequency of the PHKG1 g.8283C.A mutation in diverse pig breeds, we genotyped this mutation in a broad  Table 3). The PHKG1 mutant allele was absent in wild boars, suggesting that the mutation likely happened after domestication. However, due to limit number of wild boars examined in this study, we cannot rule out the possibility that it is a standing variation.

Discussion
To date, only a handful of QTG and QTN have been identified for complex traits in farm animals. To the best of our knowledge, this is the first study that has used a system genetics approach including GWAS, eQTL mapping and causality modeling to identify QTG and QTN for complex traits in pigs.
Our study supports the PHKG1 gene as a QTG for GP-related traits based on the following findings: (1) GWAS on the F 2 and Sutai populations enable us to define the major QTL for GP within an interval of 180-kb on SSC3, which contains only 7 genes including PHKG1. (2) Of three positional candidate genes with cis-eQTLs signals, only PHKG1 has the cis-eQTL peak SNP that is identical to the pQTL peak SNP. (3) PHKG1 is a catalytic subunit of the phosphorylase kinase (PhK), which is critical to glycogen degradation. (4) The wild-type transcript level of PHKG1 was significantly correlated to glycogen content in muscle.
Furthermore, our study shows the following evidences for the PHKG1 g.8283C.A mutation as QTN underlying the SSC3 QTL effect: (1) The mutation caused abnormal splicing of PHKG1 mRNA, resulting in a shift in the reading frame with a premature stop codon. (2) The aberrant mRNA transcripts were diminished by NMD, which reduced the PHKG1 protein level and led to PhK enzyme deficiency. (3) The mutation shows the complete concordance between their genotypes and the deduced QTN genotypes across all parental boars. No other variants in the PHKG1 gene exhibited such concordance. (4) This mutation has a consistent and significant effect on GP-related traits across all tested populations.
Here, we illustrate that the QTN (g.8283C.A) is a splice site mutation affecting both transcription and translation levels of the PHKG1 gene, which has at least three implications. One is that cis-eQTL mapping for annotated genes residing the QTL region can not only contribute to identification of regulatory QTNs, but also benefit characterization of protein-altering QTNs including splice and nonsense mutations that affect transcript levels by NMD, just like our identified QTN. Another implication is that if a gene is subjected to alternative splicing, it is necessary to determine  the levels of its different transcripts (instead of only the total amount of these transcripts) by deep RNA sequencing and then examine their individual relationship with phenotypic traits. The third one is that we could identify a cis QTL for PHKG1 protein level that overlaps with the mapped cis-eQTL and pQTL on SSC3 if a big data set on protein abundance was available, and find a strong relationship between protein level and GP-related traits. This information could add a new dimension for the search of QTG. Several studies have recently demonstrated the feasibility of highthroughput proteome quantification and revealed the variation and genetic control of protein abundance in humans and mice [23,24]. Hopefully, whole proteome analysis will be more efficient and routinely added to integrative genomic analysis to accelerate the identification of QTG and QTN for complex traits in livestock. Phosphorylase kinase (PhK) is comprised of four different subunits with a stoichiometry of (abcd) 4 ; a,b, and d are regulatory, while c is catalytic. Each of these subunits has isoforms or splice variants differentially expressed in different tissues. Mutations in 4 PhK subunit genes (PHKA1, PHKA2, PHKB and PHKG2) have been implicated in low PhK activity in liver and/or muscle [25][26][27][28][29]. No variant in the muscle isoform of the PhKcsubunit (PHKG1) has been reported for muscle PhK deficiency [17]. To our knowledge, this study is the first one to confirm the association of the PHKG1 mutation with PhK deficiency, muscle glycogenosis and meat quality traits in pigs.
Pork with high GP (.180 mmol/g) often has low ultimate pH and water holding capacity, and is therefore called ''acid meat''. The R225Q mutation in PRKAG3 gene was the first identified causal mutation for acid meat [19]. This gain-of-function mutation causes an increased glucose uptake and glycogen synthesis in skeletal muscle [30,31]. Its unfavorable allele (225Q) is dominant and specifically present in Hampshire and related breeds [19]. In contrast, the PHKG1 g.8283C.A mutation that we identified is a lossof-function mutation causing the defect in glycogen degradation. This variant occurs predominantly in Duroc or Duroc-crossed pigs. The mutant allele A appears to be partially recessive, since the average RG value of AC heterozygotes approximates to that of CC homozygotes rather than AA homozygotes (Table S3).
Costa et al. [32] have reported that 9.8% of DLY hybrid pigs free of the PRKAG3 225Q allele still show the acid meat phenotype with GP higher than 180 mmol/g. In this study, we found that the PHKG1 QTNs are segregating in DLY hybrid pigs with a frequency of 1.7% (9/540) for the homozygous mutant with GP beyond 180 mmol/g. Therefore, we speculate that a substantial proportion of acid meat is likely caused by the PHKG1 g.8283C.A mutation.
It was observed that the mutation had negative effects on almost all meat quality traits but did not impact any growth traits (e.g. carcass weight; see Table S6) in the F 2 , Sutai and DLY populations. However, Duroc and two Duroc-derived lines carry the mutation at relatively high frequencies ( Table 3). The result may suggests that meat quality has not been the primary breeding goals in these pig populations.
In conclusion, we identified a causal mutation in the PHKG1 gene associated with excess glycogen content in Duroc and its related pigs. The finding highlights the important role of genes encoding PhK subunits in the formation of the GP-related traits. More intriguingly, our finding would be of considerably importance for the pig industry as we can develop a diagnostic DNA test to effectively eliminate the PHKG1 undesirable allele from nucleus herds and consequently improve meat quality.

Ethics statement
All procedures involving animals followed the guideline for the care and use of experimental animals established by the Ministry of Agriculture of China. The ethics committee of Jiangxi Agricultural University specifically approved this study.

Samples and meat quality trait measurement
Three experimental populations were involved in this study: a White Duroc 6 Erhualian F 2 intercross, a Chinese Sutai half-sib population and a commercial DLY hybrid population. The F 2 population was established as described previously [33]. Briefly, two White Duroc boars were mated to 17 Erhualian sows. Nine F 1 boars and 59 F 1 sows were then intercrossed to produce a total of 1912 F 2 animals in 6 batches. Sutai is a Chinese synthetic line that is derived from Chinese Taihu (50%) and Western Duroc (50%) after over 18 generations of artificial selection. The Sutai population comprised offspring of 4 sires and 55 dams. A total of 930 F 2 individuals, 434 Sutai pigs and 540 DLY pigs were used in this study.
The F 2 and Sutai pigs were weaned at day 46 and 28 respectively, and males were castrated at day 90 and 18 respectively. The two populations were raised at an experimental farm of Jiangxi Agricultural University in Nanchang city during their fattening period and were slaughtered for phenotype recording at the age of 24063 days. The DLY pigs grew up at a farm of Xiushui city (about 110 miles away from Nanchang) until the slaughter weight of 90-100kg. All pigs were transported and slaughtered at the same commercial abattoir in Nanchang where the pigs were fasted for 15-20 hours with water available ad libitum. After slaughter, muscle samples were removed within 30 minutes postmortem from the longissimus (LM) and semimembranosus (SM) muscles for RNA isolation and the measurement of meat quality traits. The meat characteristics including the ultimate pH (or pH 24h, measured at 24h postmortem), drip loss, Minolta color parameters (L*, a*, b*), subjective scores of color and marbling, intramuscular fat content (IMF), glycolytic potential (GP), residual glycogen (RG), glucose-6-phosphate (G-6-P) and lactate were measured as described previously [20,34].

Array-based SNP genotyping and GWAS mapping
Genomic DNA was isolated from ear, blood or spleen tissues with a standard phenol/chloroform method. A total of 1020 animals from the F 2 population were genotyped for 62163 SNPs on the Illumina PorcineSNP60 BeadChip [21] according to the standard manufacture's protocol. The 60K Beadchip has 62,163 SNP, of which 54,920 can be mapped to the current pig genome assembly (Sus Scrofa Build 10.2) [35]. The quality control (QC) procedures were carried out by Plink v 1.07 [36]. Briefly, animals with call rate.0.9 and Mendelian error rate ,0.05, and SNP with call rate.0.9, minor allele frequency.0.05, P value.10 25 for the Hardy-Weinberg equilibrium test were included. A final set of 39414 informative SNPs on 930 F 2 pigs were used for subsequent analyses. A total of 434 Sutai pigs were also genotyped using the Illumina PorcineSNP60 BeadChip. To fine map the SSC3 QTL, we increased the marker density in this QTL region with additional 53 SNPs (Table S7), which were identified through comparison of our own whole-genome sequence data from 4 F 0 Erhualian sows and the Duroc reference genome sequence (Sus scrofa Build 10.2). All Sutai pigs were further genotyped for the 53 SNPs by using the TaqMan OpenArray Genotyping System (Life Technologies). After QC filtering as described above, 44560 SNPs were included for further analysis. We also genotyped 140 DLY pigs using our developed OpenArray 53-SNP panel. Out of the 53 SNPs, 28 passing QC were included in association analysis. These genotype data are deposited in the Dryad repository (http://dx. doi.org/10.5061/dryad.7kn7r). The allelic effect of each SNP on phenotypic traits was tested using a general linear mixed model [37]. The model included a random polygenic effect, and the variance-covariance matrix was proportionate to genome-wide identity-by-state [38]. The formula of the model is given in a mathematic expression: Y~mzXbz KwzSczZaze, where Y is the vector of phenotypes; 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. All single-marker GWAS were conducted by GenA-BEL packages [39].
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 [40]. A SNP was considered to have genome-wide significance at P,0.05/N and chromosomewide 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 respectively 1.27e-6 (0.05/39414) and 2.54e-5 (1/39414) for the F 2 population, and 1.12e-6 (0.05/44560) and 2.24e-5 (1/44560) for the Sutai population. The phenotypic variance explained by the top SNPs was estimated by (V reduce -V full )/V reduce , where V full and V reduce are residual variances of models for association analysis with and without SNP term, respectively.
The influence of population stratification was assessed by examining the distribution of test statistics generated from 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 [41]. The Q-Q plot was constructed using R software. Linkage disequilibrium (r 2 ) was estimated for SNPs around the QTL region in the two populations by using PLINK v1.07 [36].

Muscle transcriptome data and eQTL mapping
Total RNA was extracted from longissimus dorsi muscle samples of 497 F 2 animals using Trizol (Invitrogen). RNA quantity and integrity were assessed using a NanoDrop ND-1000 Spectrophotometer (Thermo Fisher Scientific) and a 2100 Bioanalyser (Agilent). Genome-wide transcripts were assayed by digital gene expression (DGE) system and data processing was conducted as described previously [11,42]. In brief, the raw tags were first filtered to produce the clean tag data. For mapping clean tags to reference transcript sets or to the pig reference genome (Sus Scrofa Build 10.2) [35], we created virtual libraries containing all the possible 17-base length sequences of these resources located next to an NlaIII restriction site. The reference transcript sets were downloaded from the database of PEDE (Pig Expression Data Explorer; http://pede. dna.affrc.go.jp/) and pig unigene in NCBI (ftp://ftp.ncbi.nih.gov/ repository/UniGene/Sus_scrofa/). The redundant transcripts overlapped between the two databases were removed from the reference transcript set. For monitoring the mapping events on both strands, virtual sense and antisense tag sequence databases were generated for both full gene and cDNA sequences using in-house Perl scripts. The clean tag sequences were then mapped using SOAP2 [43] allowing up to one mismatches in 21-bp tag sequences. Sense and antisense tag sequences that unsuccessfully mapped to reference transcripts or mapped to multiple genes were filtered. The number of clean tags that uniquely mapped to the reference transcript sequence of each gene was calculated and then normalized to TPM (number of tags mapped to each gene per million clean tags) as expression level of transcript.
The DGE data of 497 F 2 animals were adjusted for gender, batch and kinship using a robust linear regression model [44]. Associations between gene expression levels and the RG phenotype were evaluated with Pearson correlation coefficient by R software. eQTL mapping was performed for the DGE profiles in 497 F 2 animals using mixed linear model implemented by mmscore function of GenABEL in R package. Sex and batch were considered as fixed effects, and the genetic co-variances among samples were also taken into account by fitting kinship matrix derived from whole genome SNP genotypes. Bonferroni correction was applied to adjust the multiple tests. All the above mentioned analyses were carried out with R software.

Resequencing of the PHKG1 gene
The entire cDNA sequence of PHKG1 were determined by RT-PCR (reverse transcriptase-polymerase chain reaction) and RACE (rapid amplication of cDNA ends) analysis. First, RNA was extracted from skeletal muscle of 6 F 2 pigs (2 of each QTL genotype) using Trizol (Invitrogen), then cDNA was synthesized using the PrimerScript RT reagent Kit With gDNA Eraser (Takara) for RT-PCR or the 59-& 39-Full Race Kits (Takara) for RACE. The corresponding primer pairs are listed in Table S8. To get full-length gene sequence, we amplified 14 overlapping segments from genomic DNA, using standard PCR condition and primers reported in Table S9. All PCR products were purified using the QIAquick PCR Purification Kit (Qiagen) and sequenced on both strands using the same primers and BigDye Terminator v3.1 Cycle Sequencing kits (Applied Biosystems) and a 3130 DNA Analyzer (Applied Biosystems). The sequence traces were assembled and analyzed for polymorphisms using the SeqMan program (DNAS-TAR).
Genotyping of the PHKG1 g.8283C.A mutation and association analysis A 140-bp fragment ( Figure 3C) encompassing the g.8283C.A mutation was amplified by PCR from genomic DNA using primers (PHKG1-F1: 59-ATC CCT GTG CTT GCT GGT G-39; PHKG1-R1: 59-CCC GGC GGT ACT GGT AAT-39), digested using enzyme TaqaI and size fractionated by 2% agarose gel elec-trophoresis. The A allele was represented by two fragments of 78 and 62 bp, and the C allele by uncut amplicons.
Association between the PHKG1 g. 8283C.A polymorphism and meat quality traits was calculated using the least square means method of GLM (General Linear Model) procedure in R software. For the F 2 and Sutai populations (n = 930 and 434 respectively), sex and batch was included in the model as fix effects, carcass weight as a covariate. A total of 540 DLY pigs were genotyped for this mutation. The residual glycogen content and other meat quality traits were determined for muscle samples from 140 DLY pigs. The model for DLY included sex and harvest batch as fix effects, and carcass weight as a covariate.
The associations between SNP genotypes (i.e. ss131031160 SNP and PHKG1 g. 8283C.A) and expression levels of three genes (PHKG1, GUSB and SUMF2) were established by One-Way ANOVA analysis in R software. The relationship between the transcript level of PHKG1 and the RG content was calculated using Pearson correlation in R software.
Marker-assisted segregation analysis QTL genotypes of 6 F 1 boars in the F 2 population and 3 Sutai F 0 boars were determined by marker-assisted segregation analysis as described previously [45]. Briefly, a Z-score was calculated for each sire; the score is the log10 of the H 1 /H 0 likelihood ratio where H 1 assumes that the boar is heterozygous at the QTL (Qq), while H 0 postulates that the boar is homozygous QQ or qq. Boars were considered to be Qq when Z.2, QQ or qq when Z ,22, and of undetermined genotype if 22,Z,2.

Quantitative RT-PCR and western blot analysis
Quantitative RT-PCR (qRT-PCR) reactions were performed in a final volume of 10 ml containing 1 ml of 2.5-fold diluted cDNA (corresponding to 20 ng of starting total RNA), 5 ml Power SYBR Green PCR Master Mix (Applied Biosystems), forward and reverse primers (2 pM each) and 3.6 ml free water. PCRs were conducted on an ABI7900HT instrument (Applied Biosystems) under the following cycling conditions: 10 min at 95uC followed by 40 cycles at 95uC for 15 sec and 60uC for 50 sec. Two primer sets Common-5'-FP/RP and Wt-3'-FP/RP ( Figure 4A; Table S10) were used to respectively test total and wild-type transcript (Wt) levels of PHKG1. Beta Actin (ACTB) was included as endogenous controls. Expression of all assays was measured in triplicates and average values of the triplicates were used for the analysis. The quantification of transcripts was performed by the comparative C t (2 -DDCt ) method.
The presence of PHKG1 Wt transcript in animals with three QTL genotypes was assayed by RT-PCR using the primers 59-CAC CCC AAC ATC ATA CAG CT-39 and 59-ACA GAA GCC AGC ACC GTC-39 (RT-3'-RP). PCR products of 698-bp were electrophoresed on 1% agarose gel and visualized by UV illumination ( Figure 4C).
Total protein from pig muscle tissue was extracted using a Protein Extraction Kit (Applygen). Protein samples separated on SDS-PAGE were transferred onto polyvinylidene difluoride membranes and incubated with rabbit anti-PHKG1 (Proteintech) and rabbit anti-mouse b-Actin (as loading control; Beijing Zhong Shan-Golden Bridge Biological Technology) antibodies. Anti-mouse or anti-rabbit secondary antibodies conjugated with horseradish peroxidase were used and visualized using chemiluminescent substrate (Tiangen Biotech).

Glycogen phosphorylase kinase (Phk) activity assay
Phk activity of 18 F 2 animals were determined on their frozen muscle samples using Phk Colorimetric Assay kit (Genmed Scientifics) according to the manufacturer's instructions.

Splicing minigene reporter assay
Generation of minigene constructs. The splicing minigene assay was conducted as previously described [46]. Genomic DNA was isolated from an ear sample of an F 2 pig with AA genotype at the PHKG1 g.8283C.A mutation site. The PHKG1 genomic DNA region of 926 bp from partial exon 8 to partial exon 10 was amplified using forward primer (59-GAC AAG CTT GGA GCA CAG GGG TCA TCA TG-39) and reverse primer (59-AGT GGA TCC ACC CTC GCT CCG TTT CTG TG-39), carrying 59 tails with HindIII and BamHI restriction sites (underlined), respectively. After digestion with HindIII and BamHI enzymes, the PCR products were inserted into the pcDNA3.0 vectors (Invitrogen) with the cytomegalovirus promoter (PCMV). Substitution of adenine (A) to cytosine (C) at the PHKG1 g.8283C.A site was achieved by sitedirected mutagenesis of the mutant expression vector using a QuickChang Site Directed Mutagenesis Kit (Stratagene) with a forward primer (59-GTT CTC AGG CCC GGA CTC GCG CAG GTG ATC G-39) and a reverse primer (59-CAG AGC GAT CAC CTG CGC GAG TCC GGG-39).
Transfection and RT-PCR analysis. The wild-type and mutant minigene constructs were transiently transfected into HeLa cells and 293T cells using Lipofectamine 2000 Transfection Reagent (invitrogen) for 4 h in Opti-MEM (Gibco) medium, according to manufacturer's instructions. Cells were then collected 36 h posttransfection. Total RNA was extracted using TRIzol reagent (Invitrogen). RT-PCR was performed using the PrimerScript RT reagent Kit With gDNA Eraser (Takara), and the forward primer F (59-GGA AGT GCG TCA CTT CAG C-39) and the reverse primer R (59-CGG CGG TAC TGG TAA TAG A-39). RT-PCR products were separated by electrophoresis on 4% agarose gels containing ethidium bromide and visualized by exposure to ultraviolet light. Each DNA band was gel-extracted using the TIANgel Midi Purification Kit (TIANGEN) and sequenced using the Big Dye Terminator cycle sequencing kit (Applied Biosystems) on a 313XL capillary sequencer (Applied Biosystems). The expression levels of the PHKG1 minigenes and b-Actin were measured by quantitative real-time RT-PCR on an 7900HT thermal cycler (Applied Biosystems) as described above. The data were analyzed using the comparative C t (2 -DDCt ) method. Figure S1 Quantile-quantile plot of SNPs after quality control in genome-wide association analyses for meat quality traits. LM, longissimus muscle; SM, semimembranosus muscle; GP, glycogen potential; pH 24, pH measured at 24h postmortem; 24h-drip loss, drip loss of meat after hanging in an EZ-container (KABE Labortechnik) for 24h; ST, Sutai pigs.  Figure S7 The quantification of mRNA levels from wide-type and mutant PHKG1 minigenes using real-time RT-PCR. Each minigene was transfected into 293T cells. Total RNA extracted from each cell group was assayed six times. All data were normalized to b-actin mRNA levels, and then expressed as percentage of the wild-type counterpart. The average and standard deviations of expression values obtained from two cell groups with the same minigene are shown. (TIF) Figure S8 The strength of associations between residual glycogen (RG) content and the mRNA levels of GUSB (r = 0.05, P = 0.280; left panel) and SUMF2 (r = 20.076, P = 0.106; right panel) from the DGE data of 412 F 2 animals. (TIF) Figure S9 Regional association plot for residual glycogen content, 24-h drip loss and pH 36h in longissimus muscle from 140 Duroc 6 (Landrace 6 Yorkshire) pigs. (TIF)