Introduction
With the advances in 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. Such information helps in the development of marker assisted breeding as well as improving understanding of the molecular mechanisms underlying the target traits.
One of the main sources of human-consumed meat is pork, which represents more than 40% of the meat produced worldwide. The success of pig production is strongly related to improvements in growth and carcass yield. Meat-quality traits are essential for the processing industry and consumer acceptance[1], as a result, these qualitative traits have been the subject of numerous studies in breeding programs. Also, consumer preference for high quality meat has required the identification of new breeding objectives, such as high intramuscular fat (IMF). Skeletal muscle is composed of muscle fiber, fat, connective tissue, blood vessel and nerve tissue. Fat is an essential component. Muscle fat is usually divided into intermuscular fat and IMF. IMF, referred also as marbling, consists of the fat scattered inside a muscle. As an important meat quality trait, IMF content has been extensively studied because of its effect on sensory, flavor, juiciness, tenderness and nutritional quality of meat. In addition, a muscle with an adequate content of IMF is particularly suitable for transformation for dry-cured products[2]. IMF heritability value is estimated to average 0.50, indicating a genetic basis for this trait in pigs[3,4]. However, it is difficult to elucidate the genetic basis of IMF as several biochemical and metabolic processes influence fat deposition in muscles and these processes are determined by a set of interrelated genes and their interaction with environmental factors, including nutrition[5,6]. It has been shown that IMF content varies considerably in diverse breeds. For example, Chinese indigenous breeds have higher IMF contents than the commercial breeds, with the Duroc breed having the highest[7,8].
To date, multiple quantitative trait loci (QTLs) associated with IMF have been detected on pig chromosomes 1, 2, 5, 6, 8, 12, 13 and 17[3] and deposited in the pig QTLdb (http://www.animalgenome.org/cgi-bin/QTLdb/SS/index). A combined linkage scan and GWAS revealed significant pleiotropic regions with effects on both IMF and back-fat tissues[9]. In particular, a QTL on SSC6 harboring the heart fatty-acid binding protein (H-FABP) gene explained 15% to 20% of phenotypic variation in IMF content in different crosses[10–12]. In addition to QTLs, several genes have been characterized for this trait, such as leptin receptor (LEPR)[8], melanocortin 4 receptor (MC4R)[8] and insulin growth factor 2 (IGF2)[13]. Thanks to the high-throughput genotyping PorcineSNP60 BeadChip (Illumina, San Diego, CA, USA), it is feasible to conduct genome-wide association studies (GWAS) for IMF content and fluorescent markers associated with intramuscular fat content.
In the present study, a Duroc pig population, the main breed utilized for terminal sire, was genotyped with PorcineSNP60 v2 BeadChip to identify markers and genes associated with IMF contents.
Materials and methods
Animals
All experimental procedures met the guidelines of the Animal Care and Use Committee, South China Agricultural University, Guangzhou, China (Approval number SCAU#0017) and every effort was taken to minimize animal suffering. Animal and analysis workflow samples from 359 Duroc boars born between 2011 and 2014 were collected from Guangdong Wens Foodstuffs Group Co., Ltd., Yunfu, China.
Phenotyping and estimated breeding values
Pigs were scanned with an Aloka 500V SSD ultrasound machine (Corometrics Medical Systems, USA) to measure their IMF content in longissimus dorsi muscle. The images were collected at 6 to 7 cm off the midline across the tenth to the eleventh ribs, and these images were used to predict IMF content with the BioSoft Toolbox for Swine software (http://www.biotronics-inc.com/Lesson%20Series%20One%20-%20Marbling.pdf).
Genotyping and quality control
Duroc pig genomic DNA was extracted by standard protocols from ear tissue samples, and DNA quality was assessed by ratios of light absorption (A260/280 and A260/230) and electrophoresis. DNA concentration was adjusted to 50 ng·μL-1. All animals were genotyped with the PorcineSNP60 v2 BeadChip which contains 61565 SNP markers across the entire genome[14]. Quality control was performed using PLINK as previously described[15]. Briefly, animals with call rates of>0.99 and SNP with call rates of>0.99, minor allele frequency>0.01, P value>10-6 for the Hardy–Weinberg equilibrium (HWE) test were included. After this filtering procedure, a final set of 41793 informative SNPs from 359 pigs was used in subsequent analyses.
Genome-wide association study
A GWAS was conducted using a linear mixed model implemented in R package GenABEL[16,17]. The model included a random polygenic effect for which the variance-covariance matrix is proportional to genome-wide IBS and includes year and season as fixed effects. The model is Y=m+Xb+ Kw+Sc+Za+ e, where Y is the vector of IMF contents of all genotyped pigs, m is the overall mean, b is the vector of fixed effects including year and season, K is the regression coefficient of individuals IMF contents, and e is the vector of residual errors with e ~ N(0, Ise2) (with I the identity matrix and se2 the residual variance), w is the vector of bodyweight of individuals considered as covariance, c is the vector of SNP effects, a is the vector of random polygenic additive effects calculated as N(0, Gsa2) (with G the genomic kinship matrix calculated from the pedigree and sa2 is the polygenic additive variance), and X, S, and Z are the relative incidence matrix for b, c, and a respectively.
The Bonferroni method was used to determine the genome-wide significance threshold, in which the conventional P-value was divided by the number of tests performed[18]. According to the Bonferroni method, the genome-wide significant (significant) and chromosome-wide significant (suggestive) thresholds were P<0.05/N and P<1/N, respectively, where N is the number of SNPs tested in the analyses. In the present study, the significant and suggestive thresholds were 1.20×10-6 (0.05/41793) and 2.39×10-5 (1/41793), respectively.
Quantile-quantile plot
Given that population stratification greatly impacts GWAS reliability, quantile-quantile (Q-Q) plot analysis is considered an effective way to determine the reliability of the GWAS results. In a Q-Q plot, the horizontal axis represents the expected -log10(P value), and the vertical axis represents the observed -log10(P value). The diagonal line represents y = x. The shaded region shows a 95% confidence interval based on a beta distribution. An overall deviation above the diagonal identity line is generally suggestive of severe population stratification[19]. Deviations from the diagonal 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[20]. The Q-Q plot was constructed using the R software.
Characterization of candidate genes
SNP positions from the most recent Sus scrofa genome (version 10.2) were downloaded from https://www.animalgenome.org/pig. The NCBI annotation version of S. scrofa genome was used to find the genes which were nearest the significant SNPs (https://www.ncbi.nlm.nih.gov/genome/?term=pig). To seek further gene function information, human homologs of these genes were queried in the Ensemble BioMart (http://www.ensembl.org/biomart/martview). To assign significant SNPs positions to previously mapped QTLs in pigs, all pig QTL data was downloaded from http://www.animalgenome.org/cgi-bin/QTLdb/SS/download?file=gbpSS_10.2 (accessed 3 April 2016)[21].
Results
Phenotypes and quality control of genotypes
Prior to GWAS analysis, the distribution of all phenotypes was assessed using the Shapiro test[22]. All phenotypic data conformed to the Gaussian distribution. In total, 18792 markers were excluded as having a low (<1%) minor allele frequency, 489 markers were excluded because of low (<99%) call rate and 840 markers were excluded because they were not in HWE (P<10-6). A final set of 41,973 SNPs was kept for subsequent GWAS analysis of 359 pigs. The number of markers on each chromosome and average distances between two markers after quality control are given in the supplementary materials (Table S1). The average physical distance between two neighboring SNPs on the same chromosome was approximately 73.3 kb, ranging from 54.7 (SSC14) to 191.4 kb (SSC X).
Assessment of population stratification
The Q-Q plot of test statistics in GWAS is shown in Fig. 1a. No overall systematic bias is observed, and the deflation factor l was 1.03, indicating that population stratification had been eliminated.
Fig.1 (a) Quantile–quantile (Q–Q) plot in the left panel shows the observed versus expected -log10 (P-value); (b) Manhattan plot of genome-wide association studies for IMF in male Duroc pigs. -log10 (P-value) of the quantified SNPs were plotted against their genomic positions. Different colors indicate different chromosomes. The solid and dashed lines indicate the 5% genome-wide and chromosome-wide Bonferroni-corrected thresholds, respectively. Chromosome 19 stands for the X chromosome; (c) bar plot showing allelic effect of the most significant SNP on IMF content in pork; (d) in the region plot, -log10 (P-value) of the quantified SNPs were plotted against their genomic positions. The solid and dashed lines indicate the 5% genome-wide and chromosome-wide Bonferroni-corrected thresholds, respectively. The red dot is the most significant SNP. The nearest genes are shown at the bottom of the figure. |
Significant SNPs and candidate genes
The Manhattan plot of GWAS for IMF after genomic control in the tested Duroc pigs is illustrated in Fig. 1b. The allelic effect of the most significant SNP on IMF content in pork is shown in the Fig. 1c. The significantly associated SNPs for this trait are listed in Table 1. Eight significant SNPs were identified. Chromosomes and exact positions based on the S. scrofa genome as well as the genes nearest to significant SNPs are listed in Table 1. Five of eight these SNPs were found in the intergenic regions. The remaining three SNPs were in introns of the protocadherin-7 (PCDH7) gene. Three genes including zinc finger DHHC-type containing 16 (ZDHHC16) pseudogene, LOC102162218 and PCDH7 were proximal to significant SNPs (Fig. 1d).
Tab.1 Significant SNPs and their nearest genes for IMF1content |
| Trait1 | SNP ID | SSC2 | Location3/bp | Adjusted P value | Nearest gene | Gene location4 | Distance5/bp | Homo sapiens homologs |
|---|---|---|---|---|---|---|---|---|
| IMF | ASGA0038200 | 8 | 23147683 | 6.08×10-8 | ZDHHC16 pseudogene | 8:22,792,997..22,794,126 | 353557 | ZDHHC16 |
| H3GA0024419 | 8 | 23477511 | 8.12×10-7 | ZDHHC16 pseudogene | 8:22,792,997..22,794,126 | 683385 | ZDHHC16 | |
| DRGA0008413 | 8 | 23561975 | 8.12×10-7 | LOC102162218 | 8:24,220,615..24,276,003 | - 658640 | NA6 | |
| ASGA0038214 | 8 | 23585035 | 8.12×10-7 | LOC102162218 | 8:24,220,615..24,276,003 | - 635580 | NA | |
| ASGA0089289 | 8 | 24298660 | 8.42×10-7 | LOC102162218 | 8:24,220,615..24,276,003 | 22657 | NA | |
| DRGA0008418 | 8 | 24386346 | 8.42×10-7 | PCDH7 | 8: 24,323,016..24,732,096 | intron | PCDH7 | |
| DRGA0008423 | 8 | 24646964 | 8.42×10-7 | PCDH7 | 8: 24,323,016..24,732,096 | intron | PCDH7 | |
| MARC0038273 | 8 | 24568164 | 9.36×10-7 | PCDH7 | 8: 24,323,016..24,732,096 | intron | PCDH7 |
Note: 1 IMF, intramuscular fat; 2Sus scrofa chromosome; 3 SNP positions in Ensembl; 4 Gene location in Ensembl; 5 The SNP located in the upstream/downstream of the nearest gene; 6 NA, not available. |
Comparison with previously mapped QTLs in pigs
All eight significant SNPs identified in this study span an interval of 1.24 Mb and fall in previously identified QTLs for IMF in Pietrain, Large White, Landrace, and Leicoma pigs (Table S2)[23].
Discussion
The study presented here is a GWAS for IMF content in Duroc pigs. We identified eight novel SNPs that have not been detected in previous studies. Notably, all eight SNPs are located within a 1.24 Mb region containing three candidate genes.
ZDHHC16 is a protein coding gene in the DHHC domain family. The DHHC domain is a protein domain that acts as an enzyme, which adds a palmitoyl chemical group to proteins in order to anchor them to cell membranes[24]. Recent studies indicated that palmitoylation controls lipid metabolism. DHHC17 modulates ClipR-59 plasma membrane binding to regulate Akt signaling and glucose transporter (Glut4) membrane translocation in adipocytes[25]. Recent profiling of adipocytes, major fat storage cells, identified a wide variety of palmitoylated proteins including different transporters and components of lipid and glucose metabolism pathways[26]. One such example is the leptin signaling pathway. Leptin, a peptide hormone produced in adipocytes, regulates metabolism by suppressing feeding and increasing energy consumption. An increase in protein palmitoylation was observed in pancreatic b-cells upon glucose stimulation[27]. Also, a lack of palmitoylation impaired leptin signaling and increased levels of insulin in the serum[28]. In addition, proper expression of DHHC17 can prevent triggering apoptosis of pancreatic b-cell by interleukin 1b during inflammation[29]. Consequently, palmitoylation can impact fat synthesis by regulating the level of leptin and insulin. Another palmitoyltransferase, amyloid precursor protein (App), is one of 20 Drosophila DHHC palmitoyltransferases, transmembrane proteins that add palmitates to cytoplasmic proteins, in order to anchor them to cell membranes. It regulates protocadherin fat signaling in growth by controlling the normal subcellular localization and activity of Dachs[30]. Pseudogenes populate the mammalian genome as remnants of artifactual incorporation of coding mRNAs into transposon pathways[31]. Pseudogenes can regulate paralogous genes expression by means of the RNA interference and piRNA pathway[32–34]. Therefore, ZDHHC16 pseudogene could control lipid metabolism and fat deposition by regulating related genes and is a functionally plausible candidate gene for IMF content in pork.
LOC102162218 is a recently predicted and poorly characterized gene, and no homologous gene has been observed in the human or mouse.
The PCDH7 gene belongs to the protocadherin gene family, a subfamily of the cadherin superfamily. It encodes a protein with an extracellular domain containing 7 cadherin repeats, and is thought to function in cell-cell recognition and adhesion. Alternative splicing yields isoforms with unique cytoplasmic tails[35]. Protocadherin genes directly affect meat-quality traits. protocadherin 19 (PCDH19) is related to the ultimate pH of chicken meat[36]. Moreover, PCDH19 is significantly associated with estimated breeding values for residual feed intake in cattle. Also, the degree of methylation of the protocadherin 15 (PCDH15) precursor promoter was associated with muscle fiber density and drip loss in male three-yellow chickens[37]. Protocadherin genes also affect fat content. According to the International Mouse Phenotype Consortium, protocadherin 18 (PCDH18) knockout mice have a decreased total body fat phenotype, indicating that PCDH18 is associated with fat deposition. Unfortunately, there is no PCDH7 knockout mouse at present. In pigs, a PCDH7 homolog expressed differently in distinct adipocytes, suggesting that the expression level of PCDH7 may affect formation of different lipids[38]. The significant enrichment of rare variants in the protocadherin genes has been observed in a group of extremely obese individuals but not in the general population, indicating an association between rare variants in the protocadherin cluster genes and extreme obesity[39]. Protocadherin genes may also influence IMF by affecting eating. Knocking out the protocadherin gamma subcluster leads to a change in eating behavior in mice[40]. This suggests that rare variants in the PCDH-genes may also have an impact on food intake in pigs and PCDH7 in particular is an interesting candidate gene for IMF content in pork.
Conclusions
In summary, we identified a new genomic region and three genes associated with IMF content in porcine genome. Both ZDHHC16 pseudogene and PCDH7 had never been reported in the previous studies. Identification of the genomic region and putative positional genes associated with lipid metabolism reported here should contribute to the better knowledge of the variation in IMF in Duroc pigs.