A total of 274 samples from four KNC strains were collected from NIAS. It comprised of 245 KNC H strains (KNCH), 9 KNC S stains (KNCS), 9 KNC R stains (KNCR), and 11 KNC L strains (KNCL). The plumage colors of these strains range from black, black with brown, brown, red-brown, and black. KNC H strain chickens can have black to black and brown plumage and the individuals were classified into seven categories according to the number of body parts it exhibited brown plumage. Plumage color was scored for six specific body parts: head, neck, breast, back, wings, and tail. If the individual only had black feather, it was given a score of zero, however, if an individual showed brown plumage, for every body part it had brown it received 1 point. This classified the individuals into seven categories, ranging from all black to brown in all scored body parts. Blood samples were collected in EDTA tubes and DNA was extracted using Wizard genomic DNA purification kit (Promega, USA) according to the manufacturer’s instruction. The genomic DNA samples were genotyped using the 60K SNP Illumina iSelect chicken array (Illumina Inc., USA).

The 60K SNP Illumina iSelect chicken array contains 57,637 SNPs that are distributed across the chicken genome. SNPs were excluded if it had a missing rate of >5%, a minor allele frequency (MAF) of <0.01, or a Hardy-Weinberg equilibrium (HWE) test p-value of <10^–6 using PLINK 1.07 (Purcell et al., 2007). After the quality control, 53,257 SNPs were retained for further analysis. GWAS analyses on plumage coloration of whole body and the body parts traits were performed using mixed model of GEMMA (v0.93) (Zhou and Stephens, 2012), which accounts for population stratification and sample structure.

Where y is an n-vector of traits (plumage coloration) for n individuals; W = (w₁, w₂, …,w_c) is an n×c matrix of covariates (fixed effects) including a column of 1s; α is a c-vector of the corresponding coefficients including the intercept; x is an n-vector of marker genotypes; β is the effect size of the marker; u is an n-vector of random effects; ε is an n-vector of errors; τ⁻¹ is the variance of the residual errors; λ is the ratio between the two variance components; K is a known n×n relatedness matrix and I_n is an n×n identity matrix. MVN_n denotes the n-dimensional multivariate normal distribution. Relatedness matrix K was calculated as following:

x_i as its ith column representing genotypes of ith SNP, x̄_i as the sample mean, and 1_n as a n×1 vector of 1’s. GEMMA tests the alternative hypothesis H₁:β≠0 against the null hypothesis H₀:β = 0 for each SNP. To correct for multiple hypothesis testing, we obtained adjusted p values by using the Benjamini and Hochberg false discovery rate procedure (Benjamini and Hochberg, 1995), adjusted p-value 0.2 significance level is used. An overview of the test results was shown as a Manhattan plot constructed by the statistical package R. Base pair position of SNP markers were given based on the chicken genome assembly build WASHUC2. Inflation factor was calculated by the R package GenABEL with “median” option (Aulchenko et al., 2007).

We estimated the genetic variance of plumage by using GCTA (Yang et al., 2011). After calculating the genetic relationship matrix (GRM) between all pairs of samples using all the autosomal SNPs, we estimated the genetic component, or heritability, for each trait by REML analysis of an Mixed Linear Model y = Xβ+g_G+ε, where y is a vector of phenotypes, β is a vector of fixed effect such as sex, age with its incidence matrix X, g_G is a vector of aggregate SNP effects as random effect with Var (gG)=AGσG2, and A_G is the GRM estimated from all autosomal SNPs. We defined heritability or the proportion of variance explained by all autosomal SNPs as hG2=σG2/σP2.

Each of the 245 KNC (H strain) was investigated individually for plumage coloration (Figure 1). The predominant plumage color of KNC H strain chickens was black, but 88 out of 245 had brown feathers in addition to the black. This mixing of brown plumage causes the uniformity of Woorimatdag to decrease. Plumage color was investigated in six body-parts: head, neck, breast, back, wings and tail. One point was given for each body part that showed brown plumage (Table 1). Out of the 207 KNC H strain hens, 157 hens only had black plumage color, while 41 hens had brown plumage on the neck and 9 hens had brown plumage on both the head and neck. None of the 38 KNC H strain roosters were pure black.

Roosters and hens, respectively, have ZZ and ZW sex chromosome, which may be the cause of the differential plumage color between sexes. Sex-linked silver locus have been shown to control silver and wild type/gold color and interfere with the coloration of red (Gunnarsson et al., 2007). This result is estimated to be associated with the difference in the color of the hen and rooster. It is possible that the sex-linked plumage coloration is related to the fact that rooster with colorful plumage has an advantage when it comes to mating success (Brawner III et al., 2000).

The genome-wide association study revealed 12 significantly associated SNPs that surpassed the significance level (Figure 2, Table 2). As genomic inflation factor is 0.987, it can be concluded that the GWAS result is not inflated by considering relatedness using GEMMA (Figure 3). Among the significant SNPs, we identified 4 susceptibility SNPs: rs14339964 (Gga3:36327458, p = 4.07 ×10⁻⁹), GGaluGA344987 (Gga3:705798, p = 1.12×10⁻⁶), rs14641648 (Gga8:12987908, p = 2.06×10⁻⁶), and GGaluGA193591 (Gga24:5696828, p = 2.38×10⁻⁶) in the population (Table 2). SNP rs14339964 at Gga3:36327458 is located in an intron region of AKT3 which is known to be regulators of cell signaling in response to insulin and growth factors and involved in a wide variety of biological processes. AKT3 is one of the key genes in the formation of melanoma cells (Tsao et al., 2012). Previous studies reported that through gene-environment interactions pigmentation pathways can contribute to the formation for melanoma and tumours (Gudbjartsson et al., 2008; Ibarrola-Villava et al., 2012). Thus, we indirectly infer that AKT3 mutations may be related to plumage pigmentation. Both SNPs, GGaluGA344987 at Gga3:705798 and rs14641648 at Gga8:12987908, are located in an intergenic region around KRT7 and PAP2 which are associated with pigmentation. PAP2 is another name of LPPR5 which has been found to increase pigmentation (Shan et al., 2009). KRT7 is a member of the keratin gene family and is related with melanocytic tumors (Blum et al., 2010). DDX6 encodes a member of the DEAD box protein family, which has multiple functions including translation suppression and mRNA degradation (Weston and Sommerville, 2006). DDX6 is a previously confirmed gene for vitiligo which is a disease related with pigmentation of skin (Tang et al., 2012). Interestingly, although rs15175679 (Gga20: 8397089, p = 3.91×10⁻⁶) is not significant, the variant exits in gga-mir-668 which is a region that harbors a small RNA. Previous studies of chicken embryogenesis has shown that this small RNA regulates developmental signaling pathways (Shao et al., 2012). The results of GWAS of head plumage, wing plumage, breast plumage, back plumage, neck plumage and tail plumage traits, separately identified the same SNPs: rs14339964 and rs15616451 (near gene: AKT3, ENSGALG00000020136) as the result of GWAS with the whole body trait. Through the concordant result, we infer that quantitative analysis of whole body plumage is not a simple trait (Table 3).

The feather pigmentation related genes including MC1R, TYR, PMEL, MLPH, ASIP, SOX10, and SLC34A2 are well known. However, the related loci of these genes were not found in this study. The chicken 60K SNP chip does not contain SNPs of the MC1R region, and so we could not identify the effects of MC1R in this study. The results of this study are nevertheless meaningful in that novel loci affecting pigmentation at genome-wide level were found. Estimated genetic heritability was 18.2%, but estimated genetic heritability of significant SNPs was 3.1%. The results support a polygenic effect in feather pigmentation. This means previously reported genes MC1R, TYR, PMEL, MLPH, ASIP, SOX10, and SLC34A2 as well as the reported loci in this study are important in plumage coloration. The results may contribute to selecting and breeding of KNC H for plumage color uniformity.