In total, 83 F₀ (15 sires and 68 dams) and 595 F₁ KNCs were used for this study. They were composed of 68 half-sib families and 15 full-sib families. In order to consist of families, three sires were mated with 14 or 15 dams in each line, and 595 progenies were produced by artificial insemination. All the bird specimens were provided from National Institute of Animal Science (NIAS), Republic of Korea, and raised under the same environmental conditions until 20 weeks of age. The blood samples from each bird were collected for DNA extraction [9]. This animal experiment strictly followed “The Guide for Care and Use of Laboratory Animals” published by the Institutional Animal Care and Use Committee of the NIAS (2012-C-037) in Korea.

In total, 171 DNA markers (132 microsatellites and 39 single nucleotide polymorphisms [SNPs]) across the 26 autosomes were used in this study. The microsatellite markers were genotypes using Genetic Analyzer 3130xl (Applied Biosystems, Foster City, CA, USA) and SNP markers were genotyped by polymerase chain reaction–restriction fragment length polymorphism and Fluidigm Genotyping Technology (Fluidigm, South San Francisco, CA, USA). The order of genetic markers and genetic map distances were determined using CRI-MAP software version 2.4 [10]. The total map length of autosomal chromosome was 3,019.3 cM and the average genetic map distance between markers was 17.5 cM [11].

The thigh and breast meats of each bird were used in this study. Lipid extraction for fatty acid analysis was performed by the Folch method [12]. The fatty acid composition (%) was analyzed using a gas chromatography (HP 7890A, Agilent Technologies, Wilmington, DE, USA).

Prior to conduct QTL mapping, normal distribution of the phenotype data were validated. If putative outliers were detected, they were excluded based on normality ascertainment using Ryan-Joiner (RJ) method implemented in Minitab software (Minitab, State College, PA, USA). The RJ score ≥0.99 was employed for the ascertainment of normality. Several fatty acid composition related traits showed significant deviation from normality. Therefore, they were transformed with natural logarithm or square root to eliminate their skewness. A general linear model was carried out using the Minitab software (version 14, Minitab Inc., USA) to identify factors affecting phenotypic variation. Estimation of heritabilities and genetic correlations were conducted by ASReml-R package [13].

QTL analysis was performed using two different approaches: half-sib QTL analysis and full-sib QTL analysis. Phenotypic records used for half-sib QTL mapping were pre-adjusted for fixed effects of sex, batch, line, and carcass weight. A half-sib QTL analysis was conducted with GridQTL program using paternal half-sib families [14]. The conditional probability of alternative QTL genotypes of given marker genotypes was computed at 1-cM intervals within and across the half-sib pedigrees. The pre-adjusted phenotypic data of offspring were then regressed onto the probability of QTL genotype to access the significant effect of allele substitution using the following equation:

Where, y_ij was the phenotypic record of the jth bird, originating from the ith sire; m_i was the average effect for the ith half-sib family; a_i was the allelic substitution effect for a putative QTL within the ith sire; x_ij was the conditional probability for jth individual that could transmit a given allele from the ith sire; e_ij was the residual. Chromosome-wide significance thresholds (i.e., 1% for significant threshold, 5% for suggestive threshold) were established by 1,000 permutations [15].

A variance component QTL mapping was performed with SOLAR program using full-sib families [16]. First, the identical-by-decent (IBD) matrix was computed based on the markers, the genetic linkage map, and the family information at each putative position along the chromosomes. The Z chromosome was not included for IBD matrix construction. The computed IBD matrix was then used to model a QTL using the following full linear mixed model:

Where, y was a vector of the phenotypic data for fatty acid composition related traits; b was a vector of fixed effects including sex, line, batch, and carcass weight; u was a vector of random polygenic effects; q was a vector of random QTL effect; e was a vector of residual; The mean and variance for random polygenic effects was defined as: u~N(0, Aσ_a²), where A was the numerator relationships matrix computed from the resource F1 pedigree in this study and σ_a² was the additive variance component. The mean and variance for additive QTL effects were defined as: q~N(0, Gσ_q²), where G was the IBD matrix and σ_q² was the QTL effect variance component, e was a vector of residual with a distribution of N(0, Iσ_e²), where I was the identity matrix and σ_e² was the residual variance component. X, Z, W were incidence matrices for b, u, and q, respectively. The reduced linear mixed model used for the null hypothesis was:

Model 3 used the same variable definitions as in model 2. Likelihood ratio (LR) for a given genomic position was computed by comparing log likelihood of the two models. Under null hypothesis, the LR value asymptotically followed x²_df=1-distribution. Chromosome-wide significance thresholds (i.e., 1% for significant threshold, 5% for suggestive threshold) were established by a numerical method as described previously [17]. Unless otherwise stated, the p-values are nominal. The 1-LOD (logarithm of odds) drop method was employed to estimate the confidence intervals for the identified QTL by half-sib and full-sib QTL analyses [18].

Positional candidate genes in the significant QTL regions were screened using Ensembl-BioMart (http://asia.ensembl.org). The candidate genes were selected based on the gene name and previous scholarly literature using National center for Biotechnology Information (NCBI; https://www.ncbi.nlm.nih.gov) and Google Scholar (scholar.google.co.kr) database.

The results of heritability estimations for the fatty acid composition of the KNC meat are shown in Table 1. The heritabilities of each trait were estimated to be in the range of 0.025 to 0.560. All estimated values were moderate to high except thighs C18:2 and 20:1 and breast C20:1. The estimation results of genetic and phenotypic correlations between each fatty acid trait are described in Supplementary Table S1. A wide range of correlations were observed between fatty acids; phenotypic correlation showed a range of −0.903 to 0.916, and genotypic correlation showed a range of −0.968 to 0.981. There were high negative correlations between monounsaturated fatty acid and polyunsaturated fatty acid contents. These negative correlations were also observed in the estimation of genetic and phenotypic correlation results for thigh and breast meats (Table 2).

It is important to know the heritability of each group for each fatty acid trait, and their correlation structures, to enable the selection of desirable animals; meat fatty acid composition has relevance to meat quality. Therefore, the heritability estimation results of this study can be used for further study on the meat quality of KNC, especially its unique flavor and taste. However, further elucidation of the genetic and phenotypic relationships between each fatty acid trait, between saturated and unsaturated fatty acids, and between fatty acids and other carcass traits, is needed. The fatty acid traits are not a single trait; therefore, it is difficult to apply those in real-life breeding programs [8].

Table 3 shows the results of half-sib and full-sib QTL analyses. Using half-sib analysis, 30 QTL regions were identified in the KNC thigh and breast meat, but the 12 QTL regions of thigh meat and 18 regions of breast meat did not match each other. A total of nine QTL regions were detected at a 1% chromosome-wide significant level (p_significant<0.01); C16:1, C18:3, and C20:1 in the thigh meat and C16:0, C16:1, C18:0, C18:2, C18:3, and C22:6 in the breast meat. In total, 21 other QTL regions were detected using a 5% chromosome-wide suggestive levels (p_suggestive<0.05).

Full-sib QTL analysis identified a total of seven QTL regions in fatty acid composition-related positions. These QTL regions did also not match between the thigh and breast meat. In the GGA2, two QTLs were investigated for breast C16:0 and C22:6, with 5% chromosome-side suggestive levels (p_suggestive<0.05). Chromosome 7 was related to C20:4 in thigh meat (p_suggestive <0.05), and C16:1 (p_suggestive<0.05) and C18:1 (p_significant<0.01) in breast meat, and chromosome 10 had two suggestive QTL regions for C16:0 and C18:0 in the thigh meat (p_suggestive<0.05).

The overall result shows that in total, 16 autosomes were related with the fatty acid composition of KNC meat. Chromosome 7 had the largest number of QTLs (5 QTLs), while four QTL regions were investigated in chromosome 2, 3, and 12, respectively. In chromosome 2, there were coincident regions which affected breast C16:0 in the half and full-sib analyses. In total, 10 fatty acid traits had significant QTL regions (p_significant<0.01) across nine chromosomes.

There were differences in the QTL regions between thigh and breast meat, one of the reason is the level of fat in the meat [8]. Generally, there are differences in the general composition and cholesterol content of each part of the animal, which also affected the results of this study. The reason that the half-sib and full-sib analysis results showed different QTL regions is dependent on which genetic information was used for this study; the half-sib analysis used only maternal genetic information, while the full-sib analysis used information from both parents when the identical-by-descent proportion was calculated [19].

Table 3 also presented positional candidate genes which related the fatty acid composition of this study. These candidate genes were shown to be directly/indirectly involved in lipid and fat metabolism in previous studies. Adiponectin is encoded by the adiponectin, C1Q and collagen domain containing (ADIPOQ) gene, ADIPOR2 (Adiponectin receptor 2) is one of the ADIPOQ binding receptors. Adiponectin is adipocyte-specific protein, which has an important role in insulin resistance; it is currently being studied in obesity, diabetes, and arteriosclerosis research [20]. Lipoprotein aids in the transportation of hydrophobic lipid molecules, and they are classified by density (high-density lipoprotein [HDL], low-density lipoprotein [LDL], intermediate-density lipoprotein, very-low-density lipoprotein, and chylomicrons). LDL receptors mediate the endocytosis of cholesterol-rich LDLs. LDL receptor-related proteins conduct a wide range of physiological processes, including the regulation of lipid metabolism [21]. HDL binding proteins re regulated by the level of cellular cholesterol, and may play a role in protecting cells from the over-accumulation of cholesterol [22]. Fatty acid binding proteins (FABP) are a family of lipid-binding proteins (FABP1 to FABP9 genes), and are distributed by tissue specificity. Their main role is the regulation of fatty acid uptake and intracellular transport [23]. Elongation of very long-chain fatty acids proteins (ELOVL) encode the enzyme elongase, which catalyzes the fatty acid elongation process. The ELOVL family (ELOVL1 to ELOVL7 genes) have a characteristic fatty acid substrate specificity and tissue-specific gene expression [24,25]. Fatty acyl-CoA reductase (FAR) catalyzes the reduction of fatty acyl-CoA to fatty alcohol, and FAR1 and FAR2 encode the reductase enzyme, which is needed for the first step of the wax biosynthesis. Mannosyl-glycoprotein N-acetylglucosaminyltransferases are a subset of glycosyltransferases, which catalase glucosidic bond formations using sugars that are offered from nucleoside phosphates or lipid phosphates [26]. Glucosaminyl (N-Acetyl) transferase (GCNT) also includes glycosyltransferase, and GCNT2 was reported as a glycosphingolipid biosynthesis-related gene in chickens [27]. The P-type ATPase catalyzes the translocation of phospholipids of cell membranes. The phospholipid transporting ATPase 11C is active at the plasma membrane as a major type IV P-ATPase flippase, and transports phosphatidylserine from the outer to inner leaflet of the membrane lipid bilayer of most cells [28].

Based on the role of these candidate genes associated with lipid/fat metabolism, some of the genes have already been reported in beef and pork, in which they were found to affect meat quality traits. As a representative gene, ADIPOQ was found to significantly affect traits such as marbling, back fat thickness, carcass weight, and fatty acid oxidation [29–32]. FABP4 was also found to have an effect on the intramuscular fat composition, marbling, and fatty acid composition [33–36].

Consumers in modern society are increasingly interested in healthy livestock products and high-quality meat. Therefore, there is a growing need for genetic studies on the meat fatty acid composition to maintain and improve these qualities. Fatty acid composition is closely related with lipid and fat metabolism, so it is necessary to understand the synthesis process and the genetic relationships between fatty acids and lipids, and between different fatty acids in meat. Research is also needed to thoroughly elucidate the function of the genes within the investigated QTL region that were unknown in this study.

This study provides insights into the fatty acid composition of chicken meat because there is lack of genetic research in this area. The estimation of heritability and the corresponding correlations observed in this study can be used for the genetic improvement of meat fatty acid composition in the future.

Taken together, the results of this study can be provide a baseline for the causative genes that affect the fatty acid composition of chickens, especially regarding meat flavor. Flavor is one of the most important traits that reflects the consumer’s preferences in Korea. Moreover, techniques similar to this study are expected to be used for breeding strategies of KNC as well as the production of new commercial chicken breeds using indigenous chickens, in the future.