Genome-wide Association Study of Chicken Plumage Pigmentation.

To increase plumage color uniformity and understand the genetic background of Korean chickens, we performed a genome-wide association study of different plumage color in Korean native chickens. We analyzed 60K SNP chips on 279 chickens with GEMMA methods for GWAS and estimated the genetic heritability for plumage color. The estimated heritability suggests that plumage coloration is a polygenic trait. We found new loci associated with feather pigmentation at the genome-wide level and from the results infer that there are additional genetic effect for plumage color. The results will be used for selecting and breeding chicken for plumage color uniformity.

INTRODUCTION

Poultry industry, in the process of rapid industrialization, developed commercial chicken strains from a small number of breeds. To increase the productivity of native chickens, they were bred for economic traits. Although this process resulted in higher productivity, at the same time it decreased genetic diversity (Tadano et al., 2007). In recent years, it has become increasingly important to protect national endemic genetic resources and use local breeds to create commercial strains that can adapt to the changing environment.

In Korea, the National Institute of Animal Science (NIAS) has been studying the process of indigenization of foreign breeds in to Korea and methods to restore Korean native chicken breeds. Korean native chickens (KNCs) as defined by NIAS in 2008 are chickens that have been bred true for at least seven generations. The commercial KNC called Woorimatdag (WR CC) was developed by crossing three native chicken breeds (Heo et al., 2011). Woorimatdag has contributed to the industrialization of KNCs because of its rapid growth and the texture of the meat in comparison to the native chickens (Park, 2010). However, the use KNC H strain in the paternal line to create Woorimatdag has led to the decrease in plumage uniformity. Unlike typical white broilers, KNCs usually have colored feathers and various pigmentation patterns. Plumage color is an important factor that is used by consumers to distinguish between KNC strains. Although plumage color is easily observed, the genetics behind the feather coloration is governed by both qualitative and quantitative features (Klungland and Vage, 2000). In chickens, mutations in MC1R and TYR genes have been shown to be associated with feather pigmentation (Kerje et al., 2003; Liu et al., 2010). However, there is a lack of research on the genetics of plumage coloration in Korean chicken at the genome-wide scale. The purpose of this study is to characterize the genetic polymorphism underlying different plumage color using the chicken 60K SNP chip through GWAS (genome-wide association study) and to increase plumage color uniformity of Woorimatdag. The results will also be used for selecting and breeding KNC H strain.

MATERIALS AND METHODS

Sample and phenotype collection, and genotyping

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).

Genome-wide association test

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 = (w1, w2, …,w) 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; τ−1 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 is an n×n identity matrix. MVN denotes the n-dimensional multivariate normal distribution. Relatedness matrix K was calculated as following:

x as its ith column representing genotypes of ith SNP, x̄ as the sample mean, and 1 as a n×1 vector of 1’s. GEMMA tests the alternative hypothesis H1:β≠0 against the null hypothesis H0:β = 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).

Estimating genetic variance

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β+gG+ε, where y is a vector of phenotypes, β is a vector of fixed effect such as sex, age with its incidence matrix X, gG is a vector of aggregate SNP effects as random effect with , and AG is the GRM estimated from all autosomal SNPs. We defined heritability or the proportion of variance explained by all autosomal SNPs as .

RESULTS AND DISCUSSION

Plumage color of KNC H strain

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.

Figure 1.

Example of plumage coloring in chicken. (a) (b) female, (c) (d) male.

Table 1.

Plumage color pattern of KNC H strain

Feather color						Point	Number of chickens
Head	Neck	Back	Breast	Wings	Tail
Male
0	1	0	0	0	0	1	3
0	1	1	0	0	0	2	2
1	1	0	0	1	0	3	9
1	1	1	0	0	0	3	1
1	1	1	0	1	0	4	14
1	1	1	0	1	0	4	4
1	1	1	1	1	0	5	5
Female						Shank color
0	0	0	0	0	0	0	157
0	1	0	0	0	0	0	41
1	1	0	0	0	0	0	9

Feather color: 0 = black, 1 = black+brown.

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 SNPs associated with feather pigmentation

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−9), GGaluGA344987 (Gga3:705798, p = 1.12×10−6), rs14641648 (Gga8:12987908, p = 2.06×10−6), and GGaluGA193591 (Gga24:5696828, p = 2.38×10−6) 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−6) 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).

Figure 2.

Manhattan plot of GWAS result. GWAS for the integrated phenotypes using Illumina chicken 60K SNP BeadChip of 274 samples. The x-axis of the Manhattan plot shows the genomic position, the y-axis represents the log10 base transformed p-values, Benjamini and Hochberg false discovery rate procedure (Benjamini and Hochberg, 1995), the dashed line show significance level of adjusted p<0.2.

Table 2.

Top SNPs associated with plumage coloration

rs	CHR	Position	Min /Maj1	Freq	Beta	SE	p value	Q value*	gG SNP**	Gene	Location
rs15304667	1	70248953	G/A	0.0623	1.221	0.230	2.66E-07	0.014	0.143	STK38L	Intron
rs15408789	1	125252900	G/A	0.115	0.548	0.116	3.58E-06	0.186	0.112	AP1S2	Intergenic
GGaluGA172731	2	149522175	G/A	0.421	1.180	0.244	2.43E-06	0.126	0.575	ENSGALG000 00018081	Intergenic
rs14339964	3	36327458	A/C	0.210	1.186	0.194	4.70E-09	0.000	0.394	AKT3	Intron
GGaluGA239670	3	110847381	G/A	0.206	1.366	0.241	4.50E-08	0.002	0.447	TFAP2B	Intergenic
rs15616451	4	75015502	A/G	0.053	0.597	0.125	3.30E-06	0.171	0.060	ENSGALG000 00020136	Intergenic
rs16445392	4	85858579	A/C	0.025	0.908	0.186	1.85E-06	0.096	0.044	MXD4	Intron
rs15790835	6	19031740	A/C	0.132	1.087	0.227	2.89E-06	0.150	0.249	ENSGALG000 00005969	Intergenic
rs14641648	8	12987908	A/G	0.115	1.165	0.239	2.06E-06	0.107	0.237	PAP2	Intron
rs15047928	19	5092926	A/G	0.329	1.325	0.251	2.87E-07	0.015	0.585	FAM211A	Intron
GGaluGA193591	24	5696828	G/A	0.121	1.162	0.240	2.38E-06	0.123	0.247	DDX6	Intron
GGaluGA344987	E22C19 W28	705798	G/A	0.287	1.110	0.222	1.12E-06	0.058	0.454	KRT7	Intergenic

Minor allele/Major allele.

Adjusted p value.

Estimated variance.

Figure 3.

Quantile-Quantile plot of GWAS result Inflation factor (lambda) = 0.9877269.

Table 3.

Top SNPs associated with plumage coloration by each part

Plumage part	rs ID	CHR	Position	Min /Maj1	Freq	Beta	SE	p value	Q value*	Gene	Location
Head	rs14398623	3	97933258	A/G	0.115	−0.149	0.030	1.56E-06	0.081	RNF144A	Intron
	rs318020030	7	19055437	A/G	0.350	0.220	0.045	1.79E-06	0.093	PPP1R9B	Intron
	rs16669242	9	14516013	A/C	0.296	0.235	0.041	2.74E-08	0.001	FGF12	Intergenic
	GGaluGA105119	14	14228307	G/A	0.321	0.154	0.032	2.76E-06	0.144	ENSGALG000 00023628	Intergenic
	rs14112979	18	7210201	G/A	0.341	0.219	0.042	3.26E-07	0.017	HELZ	Intron
Wing	rs13982792	1	1.83E+08	A/G	0.463	0.153	0.033	3.48E-06	0.181	XPO4	Intron
	rs14131527	2	4895894	A/G	0.346	0.114	0.022	4.14E-07	0.022	XIRP1	Intron
	rs14188826	2	59639362	A/G	0.058	0.163	0.034	2.70E-06	0.140	PRL	Intergenic
	rs14339964	3	36327458	A/C	0.210	0.146	0.029	1.19E-06	0.062	AKT3	Intron
	rs14404313	3	103537929	C/A	0.146	0.203	0.036	3.75E-08	0.002	NT5C1B	Intergenic
	GGaluGA110134	15	9605683	G/A	0.333	0.175	0.037	3.52E-06	0.184	MSI1	Intron
	rs15027075	17	10503821	A/G	0.357	0.111	0.023	2.70E-06	0.141	PBX3	Intron
Breast	rs14316836	3	7925459	G/A	0.293	0.062	0.013	2.49E-06	0.129	LCLAT1	Intron
Back	rs14401050	3	100155987	A/G	0.204	0.174	0.036	2.23E-06	0.116	E2F6	Intergenic
	rs15616451	4	75015502	A/G	0.053	0.118	0.024	1.71E-06	0.089	ENSGALG000 00020136	Intergenic
	GGaluGA287070	5	49185847	G/A	0.393	0.198	0.038	4.53E-07	0.024	VRK1	Intergenic
	GGaluGA095084	13	11601056	G/A	0.082	0.175	0.033	1.81E-07	0.009	SGCD	Intron
	rs15022353	15	7826821	A/G	0.216	0.178	0.029	1.84E-09	9.56E-05	TTC28	Intron

Minor allele/Major allele.

Adjusted p value.

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.