The study protocol and the standard operating procedures (No. 2009-077, C-grade) of Berkshire pigs were reviewed and approved by National Institute of Animal Science’s Institutional Animal Care and Use Committee.

An inbred Berkshire population was used for analysis, and IMF of the Berkshire sample was measured. Seven hundred and four samples were examined. Among them, 367 samples were male, 204 samples were female and the sex of 133 samples was unknown. The IMF of each pig was measured by chemical fat extraction procedures.

The genomic DNAs of pigs were genotyped on the Illumina Porcine 60 K SNP Beadchip. A total of 62,163 SNPs were genotyped. We discarded the markers with low minor allele frequency (<0.05), significant deviation from Hardy-Weinberg equilibrium (p<10⁻³), and low genotype call rate (<95%). Among 62,163 SNPs, 40,191 SNPs passed quality control. The Hardy-Weinberg test failed 3,651 SNPs, 3,304 SNPs failed the genotype missingness test, and 19,829 SNPs failed the minor allele frequency test.

The phenotype (IMF) was standardized to z-scores by subtracting the mean and then dividing by the standard deviation, in each sex group (male, female, unknown) separately. Single trait, univariate linear mixed model was used for the analysis assuming additive effect of SNPs. SNP effects were treated as random effects and sex was added as a covariate. Software GEMMA was used to calculate the genetic relationship matrix of individuals and to test the effects of SNPs by likelihood ratio test [8]. The cutoff for statistical significance of genes was p<0.01.

Gene annotation of significant SNPs was based on the Ensembl Genes 89 database of Sus scrofa genes (Sscrofa 10.2). Significant SNPs were annotated to the nearest genes within a distance of 100 kb. Functional classification of genes was performed on DAVID, an online functional annotation database. Sus scrofa was selected as both species and background option. The cutoff of gene ontology was p<0.05.

The GCTA tool [9] was used to calculate heritability for IMF. We calculated the genetic relationship matrix between all pairs of samples using all the autosomal SNPs. We then estimated the variance of genetic component by restricted maximum likelihood analysis, and heritability by dividing the estimated genetic variance by the total variance measured.

A total of 365 SNPs from all 19 chromosomes were identified as significant SNPs as the result of GWAS. Chromosome 14 contained 53 significant SNPs which was the largest number among all chromosomes. There were 40 and 35 significant SNPs on chromosome 7 and 11 respectively, which contained second and third number of significant SNPs. The statistical significance values of the association between each SNP and IMF calculated as −log10 (p-value) across 18 autosomal chromosomes and chromosome X was plotted in the form of a Manhattan plot (Figure 1).

The 365 significant SNPs found from GWAS were annotated to the nearest genes within 100 kb. Among the 365 significantly identified SNPs, 153 were annotated across 120 genes. There were some SNPs annotated to same genes and none of the significant SNPs on chromosome 8 and 15 had genes within 100 kb distance. Full information of significant SNPs, their chromosome number, position, closest gene, distance from the closest gene, raw p-value is on Supplementary Table S1.

Identified genes were compared with IMF QTL regions base on the Pig QTLdb. Total 25 genes from the 120 significant genes, which is 20.8%, were included in reported pig IMF QTL regions. Seven genes on chromosome 9, 6 genes on both chromosome 2 and 6, and 2 genes on chromosome 4, 7, and 17 each were mapped on QTLs (Figure 2). This suggests a considerable proportion of the genes identified from this study was consistent with the previous QTL studies, and those genes can be considered as genes that are located on the section of the genome having high correlation with IMF of pigs.

Identified genes were classified by their biological function and ontology. Regulation of mitogen-activated protein kinase (MAPK) cascade was the most significant gene ontologies from GOTERM_BP_5. The full result of functionally annotated genes are listed on Supplementary Table S2. Especially, bone morphogenetic protein (BMP)-binding endothelial cell precursor-derived regulator (BMPER), forkhead box protein O1 (FOXO1), ectodysplasin A receptor, ring finger protein 149 (RNF149), cluster of differentiation (CD40), tyrosine-protein phosphatase non-receptor type 1 (PTPN1), SRY (sex determining region Y)-box 9 (SOX9), MYC proto-oncogene (MYC), and macrophage migration inhibitory factor were categorized as genes related to both MAPK cascade and the regulation of MAPK cascade. In addition, growth/differentiation factor 9 and BMP6 were related to the regulation of MAPK pathway. FOXO1, RNF149, PTPN1, MYC were additionally annotated to negative regulation of MAPK cascade and regulation of stress-activated MAPK cascade.

Heritability of IMF was estimated by GCTA. The total variance of the sample was 1.020818 and the genetic variance was 0.526911. The genetic variance was estimated by the variance of genome-wide SNPs. The estimated heritability was 0.516166, approximately 0.52, and the standard deviation of the estimated heritability was 0.061655.

In previous studies, heritability of IMF was estimated at 0.39 [10], 0.44 [11], 0.52 [12], 0.65 [13]. Estimted heritability of IMF in the referred studies had values between 0.39 and 0.65. The heritability estimated from the SNP and phenotype data in this study, 0.52, was in the range of reported estimations and was according with the previous studies.

MAPK cascade was the most significant gene ontology (GO) term from the functional annotation results of significant genes and other GO terms related to MAPK cascade appeared multiple times as well. MAPK pathway regulates various cell functions such as proliferation, differentiation and mitosis [14]. Moreover, MAPK pathway is closely related to the differentiation of preadipocytes to adipocytes [15,16]. Some of the proteins involved in MAPK pathway also regulate adipocyte differentiation. For example, MAPK phosphatase-1 (MKP-1) downregulates the expression of p42/p44 MAPK and plays an important role in adipocyte differentiation [15]. In addition, inhibition of p38MAPK decreases adipocyte differentiation in humans and therefore p38MAPK activation can be seen as a requirement for primary human adipocyte differentiation [16]. Since IMF is determined by the amount of adipocytes, genes related to MAPK pathway could affect IMF by regulating the amount of adipocyte differentiation.

Some of the significant genes related to MAPK cascade or the regulation of MAPK cascade (BMPER, FOXO1, SOX9, PTPN1, CD40) are previously reported to have influence on adipocyte differentiation. BMPER directly interacts with BMPs [17]. Some BMPs activate p38MAPK pathway through the MAPK kinase (MAPKKK) cascade [18] and BMPER could be needed for adipocyte differentiation to activate p38MAPK. Furthermore, BMP4 has an effect on lipid accumulation as well as expression of adipocyte markers [19]. Furthermore, BMP2 and BMP7 induces adipocyte differentiation at low concentrate in C3H10T1/2 cell line [20]. FOXO1 is expressed in the early stages of adipocyte differentiation and acts as a preadipocyte differentiation preventing substituent [21]. Epidermal growth factor (EGF) repeat containing transmembrane protein (pref1) activates MAPK and upregulates SOX9 resulting in inhibition of adipocyte differentiation [22]. CD40 is related to the activation of MAPK [23] and PTPN1 is a negative regulator of CD40 [24]. PTPN1 polymorphisms are reported to be associated with adipocyte related measures such as body fat percentage [25].

The estimated heritability of IMF, 0.52, was relatively high. This means that a substantial part of the total phenotypic variance of IMF is explained by the genetic variance. Here, the variance is that of the population, and thus high heritability suggests high genetic influence in the population on the whole [26]. Heritability is an important parameter for predicting the response to selection [27]. Since the Breeder’s equation is given as R = h² S, where R is the response to selection, S is the selection differential and h is the heritability [28], higher heritability can result in stronger response of selection and effective selection. Therefore, the phenotype information of IMF can be useful information for selecting pigs to breed pigs having high level of IMF.

Pork containing more than 3% IMF tends to have higher sen sory qualities including juiciness, tenderness and taste [29]. As IMF of pork increased from a range of 2.01% to 3.0% to higher than 3%, juiciness, tenderness, and both the intensity and desirability of taste increased. Since the current average of IMF measured from the Berkshire sample was 2.82%, if we increase IMF up to 3% by selection and breeding, we would be able to produce pork with improved juiciness, tenderness and taste.

The tool used for association analysis in this study, GEMMA, adjusts the effect of sex by using sex as a covariate and uses genetic information from the X chromosome in the same way as those from autosomal chromosomes while computing the genetic relationship matrix [8]. However, since females carry two copies of X chromosomes while males carry a single copy, different methods should be used to estimate the genetic relationship for female-female pairs, male-male pairs and female-male pairs respectively in GWAS analyses [9]. Furthermore, among the 704 samples used in this study, the sex of 133, 19% of the sample, was unknown. Also, to balance the allele dosage between sexes, one of the female X chromosome is silenced by random X chromosome inactivation [30]. Therefore, additional information coding which allele was inactivated is needed to adjust GWAS analyses. In this study, information about which allele was inactivated was not provided, and this might together cause inaccuracy in the results from the X chromosome [31]. However, the proportion of significant SNPs on the X chromosome was 2.47% (9 out of 365) which was relatively low. Thus, some part of inaccuracy in the results from the X chromosome may have not affect the overall results of the study.

Owing to the small sample size of animals, the overall signifi cance of the study was low. Small sample size lowers the estimated effect size and consequently lowers the power of the study. The estimated power of the study was only 0.21 [32]. To detect significantly associated SNPs in a study with low power, we had to use liberal statistics and a liberal cutoff (raw p-value and p<0.01). This might cause some significant SNPs to be false positives, but still the SNPs detected in this study can be suggested as candidates for SNPs related to IMF of pig. Besides, we could pick out some SNPs more likely to related to IMF of pig by comparing them with known QTLs or searching their biological pathways. The genes mapped on QTLs or related to MAPK cascade may be stronger candidates for genes that are associated with IMF of pig than others.