A Major Locus for Quantitatively Measured Shank Skin Color Traits in Korean Native Chicken.

Shank skin color of Korean native chicken (KNC) shows large color variations. It varies from white, yellow, green, bluish or grey to black, whilst in the majority of European breeds the shanks are typically yellow-colored. Three shank skin color-related traits (i.e., lightness [L*], redness [a*], and yellowness [b*]) were measured by a spectrophotometer in 585 progeny from 68 nuclear families in the KNC resource population. We performed genome scan linkage analysis to identify loci that affect quantitatively measured shank skin color traits in KNC. All these birds were genotyped with 167 DNA markers located throughout the 26 autosomes. The SOLAR program was used to conduct multipoint variance-component quantitative trait locus (QTL) analyses. We detected a major QTL that affects b* value (logarithm of odds [LOD] = 47.5, p = 1.60×10-49) on GGA24 (GGA for Gallus gallus). At the same location, we also detected a QTL that influences a* value (LOD = 14.2, p = 6.14×10-16). Additionally, beta-carotene dioxygenase 2 (BCDO2), the obvious positional candidate gene under the linkage peaks on GGA24, was investigated by the two association tests: i.e., measured genotype association (MGA) and quantitative transmission disequilibrium test (QTDT). Significant associations were detected between BCDO2 g.9367 A>C and a* (PMGA = 1.69×10-28; PQTDT = 2.40×10-25). The strongest associations were between BCDO2 g.9367 A>C and b* (PMGA = 3.56×10-66; PQTDT = 1.68×10-65). However, linkage analyses conditional on the single nucleotide polymorphism indicated that other functional variants should exist. Taken together, we demonstrate for the first time the linkage and association between the BCDO2 locus on GGA24 and quantitatively measured shank skin color traits in KNC.

INTRODUCTION

Genetic studies of shank skin colors have long history. In 1935, Knox investigated the qualitative inheritance of shank skin colors in chickens (Knox, 1935). He summarized loci and their alleles that affect shank skin colors, including the W locus (W allele for white shank: w allele for yellow shank) and the E locus (E allele for the extension of melanic pigmentation for plumage and shank skin colors: e allele for non-extension of black plumage and shank skin colors). Smyth (1990) also reviewed the Mendelian inheritance of pigmentation in shank skin colors in chickens: The various shank skin colors are influenced by combination of specific genes that influence carotenoid and melanin pigmentations, unidentified polygenic modifiers and environmental factors such as diets. The carotenoid pigmentation is affected by carotenoid in feed, while the melanin pigmentation is affected by the melanocytes in the dermal and epidermal tissue. In common with Knox, the autosomal dominant W locus is known to be responsible for the dominant white (W) and recessive yellow (w) shank skin colors, which are related to the carotenoid pigmentation. For the melanin pigmentation, the sex-linked Id locus has causal relationship with dermal melanization. The inhibitor of dermal melanin (Id) is dominant to dermal melanin (id). The extended black (E) locus affects epidermal melanization, which is associated with dark shank skin colors. In addition, the fibromelanosis (Fm) locus has also been described as the causative locus for hyperpigmentation in dermal pigment.

Recently, causal genes responsible for Mendelian shank skin color variations were identified in chickens. The beta-carotene dioxygenase 2 (BCDO2) gene on GGA24 (GGA for Gallus gallus) was identified as a causal gene of white/yellow shank that corresponds to the W locus (Eriksson et al., 2008). The endothelin 3 (EDN3) gene on GGA20 was reported as the causal gene for fibromelanosis (Fm) in black Silkie chicken (Dorshorst et al., 2010; 2011). In addition, the E locus is known as the melanocortin 1-receptor (MC1R) gene on GGA11 which is associated with dominant extended black to the recessive yellow in coat/plumage colors in mammals and avians (Andersson et al., 2003; Kerje et al., 2003; Klungland and Våge, 2003; Ling et al., 2003; Rees, 2003). However, the phenotype collection of these studies have used a non-objective scoring of shank skin colors, which makes the results between experimentations and between different research groups difficult to compare.

Korean native chicken (KNC) is an important indigenous species in Korea (Choi et al., 2015). In particular, Korean consumers prefer these KNC meats because of their unique texture and flavor (Choe et al., 2010; Jung et al., 2011). The conservation project for KNC was started in 1994 by the National Institute of Animal Science (NIAS). As a result, five lines of KNC have been established and they are classified mainly based on plumage colors [Gray-Brown (G), Black (L), Red-Brown (R), White (W) and Yellow-Brown (Y); DAD-IS; http://dad.fao.org/]. When these lines were established, the selection criteria for shank skin color was charcoal gray or dark green colors so as to differentiate from commercial broilers which have yellow shanks. However, the shank skin colors in all lines have shown large color variations, even after more than 20 generations of selection (Figure 1). The shank skin color is one of the important economic traits because it influences the choice of consumer (Sirri et al., 2010). Although the price of KNC is approximately twice that of broilers in Korea, consumer acceptance of KNC has steadily increased due to KNC’s taste and unique texture of meat together with superior nutritional contents (Jeon et al., 2010). Since general consumers in Korea consider dark shank skin as a major characteristic of native chicken, they prefer chickens with dark shank skin. From a perspective of basic science, variation in skin color can provide a unique opportunity to describe the relationship between phenotypic and genotypic diversity (Barsh, 2003).

Figure 1

The phenotypic variations of shank skin color in Korean native chicken. The range and median values of L*, a*, and b* are presented.

To date, quantitative trait loci (QTL) analysis have not been focused on the identification of the genomic regions affecting shank skin color traits in chickens. Therefore, the aim of this study is to identify the loci and associated positional candidate genes for quantitatively measured shank skin color traits in KNC.

MATERIALS AND METHODS

Animals

There are five lines classified by plumage color in KNC (FAO DAD-IS, http://dad.fao.org/). Within each line, three KNC sires were mated to 4~5 dams to generate a two-generation resource pedigree that consisted of 83 founders (i.e., generation 0, G0; 15 sires and 68 dams) and 585 progeny (i.e., generation 1, G1; 282 males and 303 females). The 585 G1 progeny were classified into the 5 lines [i.e., 110 Gray-Brown (G), 88 Black (L), 135 Red-Brown (R), 122 White (W) and 130 Yellow-Brown (Y)] based on their plumage colors. This two-generation resource population consisted of 68 nuclear families that ranged in size from 3 to 20 birds. All animals were raised under standard indoor conditions at the experimental facilities of National Institute of Animal Science (NIAS), Cheonan, Republic of Korea. The Animal Care Committee at the NIAS approved all experimental procedures.

Phenotypic measurements

After the sacrifice, shank skin color data were collected from the 585 G1 birds using a spectrophotometer (Minolta CM-3500d, Osaka, Japan) and expressed as lightness (L*), redness (a*), yellowness (b*). The phenotypic measurement was performed by the same person to minimize possible variations during the measurement.

Genetic marker data

DNA isolation was conducted by a manual method (Miller et al., 1988). All experimental animals were genotyped with 171 informative DNA markers (132 microsatellites and 39 single nucleotide polymorphism [SNPs]) across the 26 autosomes and Z chromosome, selected from the Ark (http://www.thearkdb.org/arkdb/) and ENSEMBL (http://www.ensembl.org/Gallus_gallus/) databases. The polymerase chain reaction (PCR) amplifications of 132 microsatellite markers were performed with fluorescently labeled primers (with forward primer labeled with one of the four fluorescent dyes – FAM, NED, PET, VIC). In addition, 39 SNP markers were selected from our candidate genes SNPs data that was generated by direct sequencing and next generation sequencing technology (Seo et al., 2015). They were genotyped using PCR-restriction fragment length polymorphism and Fluidigm Genotyping Technology (Fluidigm, San Francisco, CA, USA). A list of markers used in this linkage analysis can be found in Supplementary Table 1. The map order and genetic map distance were determined using the CRIMAP program version 2.4 (Green et al., 1990). The total genetic map length was 3,019.3 cM. The average distance between makers was 17.5 cM. This genetic linkage map showed similar genetic distances and same genetic marker orders compared with the chicken consensus linkage map (Groenen et al., 2000). The sex-average linkage map was used for further QTL analysis.

Statistical and quantitative genetic analyses

Before the estimation of heritability and QTL analysis, we obtained descriptive statistics and identified sources (i.e., sex, batch, line and carcass weight) that affect shank skin color data using Minitab program (Minitab Inc., State College, PA, USA). The Ryan-Joiner (RJ) score, a correlation coefficient indicating degree of normality with a perfect score being 1, was used to assess the residual normality. A residual distribution with the RJ score over 0.99 was regarded as a normal distribution in this study.

The following quantitative genetic analyses were conducted in each of the three shank skin color traits using the SOLAR program version 6.6.2 (Almasy and Blangero, 2010). To compute heritabilities (h) of each trait, we used a maximum likelihood (ML)-based methods implemented in the SOLAR program to estimate variance components for random additive polygenic effect and random residual environment effect after considering the fixed effects of covariates (i.e., sex, batch, line and carcass weight).

Quantitative trait linkage analyses

We carried out genome scans to identify QTL for shank skin color traits by using a multipoint variance component linkage analysis approach. The maximum likelihood based SOLAR program was used to estimate the genetic variance due to the genomic region around a specific genetic marker (Almasy and Blangero, 2010). This method is based on calculating the expected genetic co-variances between relatives in the experimental pedigree as a function of the identical-by-decent relationships at a given DNA marker locus assumed to be closely linked to a QTL. In addition, additive polygenic effect due to correlation among family members was incorporated in the multipoint variance component linkage analysis model. The hypothesis of linkage was evaluated by comparing the likelihood of a reduced model in which the variance attributable to the QTL for shank skin color traits equals zero (i.e., no linkage) to a full model in which the variance due to the QTL affecting the trait of interests was estimated (i.e., linkage). The difference between the 2 log10 likelihoods generates a logarithm of odds (LOD) score that is correspondent to the classical LOD score of linkage analysis. Multipoint LOD scores were then computed at 1-cM intervals along each genomic position in this study. Whenever, the variance component linkage analysis detected evidence of potential linkage signals in more than one genomic location, we applied the oligogenic QTL approach: i) a locus with the highest LOD score was incorporated in the variance component model, ii) the genome scan for a second locus was repeated, conditional on the QTL detected in the previous scan. Linkage analysis of Z chromosome was conducted by using half-sib QTL analysis option in GridQTL (http://www.gridqtl.org.uk). Only males were used for the linkage analysis of Z chromosome. Multiple testing issues in QTL scans were addressed by calculation of highly significant, significant and suggestive empirical thresholds using a numerical method (Piepho, 2001). Genome-wide thresholds for significant linkage (α = 0.05), and highly significant linkage (α = 0.001) were employed. Chromosome-wide 5% significance levels were used as the threshold for suggestive linkage. The p-values are nominal unless otherwise stated. For each QTL with significant and suggestive linkage thresholds, we applied the 1-LOD drop method to estimate support intervals for QTL.

Measured genotype association, QTDT, and conditional linkage analysis

We conducted a quantitative-trait association analysis under a highly significant linkage peak on GGA24 using a measured genotype approach. This approach was used on the basis of mixed effects model approach to test genotype-specific differences in the means of quantitative phenotypes while accounting for the additive polygenic effect. The association test model includes the fixed effects (i.e., sex, batch, line, SNP marker genotypes, and carcass weight) together with the random additive polygenic effect. In addition, the quantitative transmission disequilibrium test (QTDT), an extended version of the transmission disequilibrium test (TDT) to quantitative traits, was conducted to further evaluate the association between the SNP marker and the phenotype data using the nuclear families since the QTDT was developed to access association in the presence of linkage signal without increasing type I error rate (Abecasis et al., 2000). The qtld command in the SOLAR program was used for measured genotype association analysis and QTDT (Havill et al., 2005). In order to estimate genotypic values and their standard error of each SNP genotype, the predict command in ASReml-R was used (Gilmour et al., 1995). The level of significance for an association was set at p<0.01.

To evaluate linkage conditional on the identified association, we applied a combined linkage and association analysis by simultaneously incorporating the variance component of linkage and the genotypic-specific means of the measured genotype association. The conditional linkage analysis were conducted using the SOLAR program (Almasy and Blangero, 2004).

RESULTS AND DISCUSSION

Shank skin color is composite and complex phenotype which can be represented by quantitatively measured lightness (L*), redness (a*), and yellowness (b*). Different types of chickens, domestic or indigenous chickens, have a diverse shank skin color phenotypes. Recently, the shank skin color traits have been used to reveal their origin of domestication or genetic relationship in chickens. However, the previous genetic studies on shank skin colors were performed based on discrete skin color phenotypes (Ericksson et al., 2008; Siwek et al., 2013; Li et al., 2014). Moreover, QTL studies of shank phenotypes have focused largely on shank diameter, length, and weight (http://www.animalgenome.org/cgi-bin/QTLdb/GG/ontrait?trait_ID=2295) Therefore, this study is the first to identify QTL affecting variation in quantitatively measured shank skin color traits in chicken. The genetic basis of shank skin color can have significant impacts on basic skin color biology and domestic processes. In addition, the shank skin color can often have economic values when used as an indicator for indigenous chicken breed recognition since most of native chicken have pigmented shank color. All five lines of KNC are currently regarded as having gray shank and foot color (DAD-IS: http://dad.fao.org/). But, in fact, they have diverse shank skin color variation (Figure 1). This large variation in shank skin colors was used as the fundamental source of this QTL study.

Mean (standard deviation, number of birds) lightness (L*) over all animals was 54.04 (4.78, 583) and ranged between 36.58 and 68.56. Mean (standard deviation, number of birds) redness (a*) over all animals was 0.131 (1.06, 583) and ranged between −2.60 and 4.48, whereas yellowness (b*) ranged between −2.20 and 24.11 (mean±standard deviation; 8.34±6.53, 583 birds) (Figure 1). The lightness (L*) exhibited a high level of residual kurtosis (1.78), which may affect the test statistic (i.e., LOD score). Thus, we applied inverse normalization for this trait. The quantitative genetic analysis detected significant (p< 0.001) heritabilities for all three shank skin color traits measured in this study. The heritability estimates for L* (0.87), a* (0.71) and b* (0.70) were high.

To identify genomic regions that influence shank skin color traits in KNC, we conducted genome scans by using a multipoint variance component linkage analysis approach (Table 1). For L*, only one suggestive QTL on GGA2 was detected at 325 cM (LOD = 1.62, p = 0.004). For a* and b*, two highly significant QTL and two suggestive on GGA 4, 7, and 24 were revealed, respectively. The linkage analysis identified a highly significant QTL that affected b* (LOD = 47.53, p = 1.60×10−49) on GGA24 (Figure 2A). At the same location, a QTL that influence a* was mapped (LOD = 14.20, p = 6.14×10−16) (Figure 2B). The 1-LOD drop support interval of the QTL for a* and b* was located between LEI0069 and AGL9 (Figure 2). On GGA4 and GGA7, we detected suggestive level of QTL for b*. The LOD score of QTL on GGA4 was 2.12 (p = 0.0018), whereas the LOD score of QTL on GGA7 was 2.54 (p = 0.000626).

Table 1

Quantitative trait loci for shank skin color traits in Korean native chicken

Trait	Location (GGA@cM)	Logarithm of odds	p-value	Flanking markers
L*	2@325	1.62*	0.00441	ADL0146–MCW0157
a*	24@130	14.20***	6.14×10−16	LEI0069–ALG9
b*	4@0	2.12*	0.00178	ADL0255–ADL0203
b*	7@94	2.54*	0.000626	ROS0019–ADL0169
b*	24@130	47.53***	1.60×10−49	LEI0069–ALG9

Suggestive chromosome-wise significance;

significant experiment-wise threshold;

highly significant experiment-wise threshold.

Figure 2

The marginal and conditional logarithm of odds (LOD) scores for b* (A) and a* (B) in Korean native chicken (KNC). The conditional LOD score means that beta-carotene dioxygenase 2 (BCDO2) g.9367 A>C variant was added when the LOD score was computed. The marker map with distances between genetic markers in Kosambi centimorgans (cM) is shown on the x-axis. The y-axis represents LOD scores. The horizontal lines represents 5% genome-wide significance.

To search for additional genomic regions affecting b*, we conducted a multi-locus linkage analysis, which took into account the effect of the QTL on GGA 24, by fixing this position in the model and then rescreening the genome. In the 2nd linkage pass, which was conditional on the QTL on GGA 24 that was identified in the initial genome screen, evidence of linkage on GGA7 remained (LOD = 2.43, Table 2). It is important to note that we could detect new evidence of QTL on GGA8 (37 cM; LOD = 3.10) after considering the major QTL for b* on GGA24 (Table 2). In the 1st pass, the LOD score of this QTL was 1.06. This new QTL on GGA8 was located between markers MCW0275 and MCW0160. For L* and a*, the multi-locus linkage analysis could not identify additional linkage signals.

Table 2

Results from initial and conditional genome scans for b*

GGA	Position (cM)	1st scan LOD1	2nd scan LOD
24	130	47.53	-
7	94	2.54	2.43
4	1	2.12	0.61
8	37	1.06	3.10

LOD, logarithm of odds; QTL, quantitative trait locus.

All QTL for b* with LOD scores more than 1.0 in the initial scan are shown. The 2nd linkage scan was conditional on the major QTL on GGA24.

The peak LOD score for a* and b* on GGA 24 occurred at the BCDO2 locus, an obvious candidate gene for shank skin color variations (Eriksson et al., 2008). The BCDO2 gene encodes an enzyme named beta-carotene dioxygenase 2, which cleaves carotenoids (colorful substance) to apo-carotenoids (colorless substance) via an asymmetrical cleavage reaction (Kiefer et al., 2001), and thereby can affect the shank skin color variation. No obvious positional candidate gene other than BCDO2 was identified in this investigation for skin pigmentation. However, each QTL contains genes coding ribosomal proteins (QTLs on GGA2, 4, 8, and 24) and microRNAs (QTLs on GGA2, 4, 7, 8, and 24), which are implicated in skin color variation (McGowan et al., 2008; Tian et al., 2012).

Given the strong relevance of the BCDO2 gene to the shank skin color variation together with indisputable evidence for linkage, we conducted the measured-genotype association analysis and family-based tests of association (i.e., QTDT) to interrogate linkage disequilibrium between the QTL and BCDO2 g.9367 A>C, which was previously used as a marker for quantitative trait linkage analysis, in the presence of the linkage (Table 3). First, the measured-genotype association analysis revealed highly significant differences among the three genotypes in their shank skin color values a* (p = 1.69×10−28) and b* (p = 3.56×10−66). Interestingly, negative dominance implying a decreased genotypic value in the heterozygotes was found for both a* and b* at BCDO2 g.9367 A>C. Second, the QTDT also showed highly significance for associations between the BCDO2 g.9367 A>C variant and the two shank skin color traits; The p-values for a* and b* generated by QTDT were 2.40×10−25 and 1.68×10−65 respectively.

Table 3

Quantitative association between BCDO2 g.9367 A>C with color shank skin color values in Korean native chicken

Trait	BCDO2 g.9367 A>C			MGA, p	QTDT, p
	AA (28 )	AG (146 )	GG (401)
L*	54.49±0.73	55.09±0.46	52.90±0.39	1.94×10−10	5.15×10−8
a*	−0.76±0.15	−0.87±0.09	0.48±0.07	1.69×10−28	2.40×10−25
b*	0.79±0.65	0.43±0.40	11.15±0.33	3.56×10−66	1.68×10−65

BCDO2, beta-carotene dioxygenase 2; MGA, measured genotype association; QTDT, quantitative transmission disequilibrium test.

Although there was lack of evidence for linkage between BCDO2 locus and L*, significant associations of BCDO2 g.9367 A>C with L* were detected by measured-genotype association analysis (p = 1.94×10−10) and QTDT (p = 5.15×10−8). Similar case also occurred for the MC1R locus on GGA11. We previously reported association of the MC1R locus with quantitatively measured shank skin color traits in KNC (Jin et al., 2014). However, we could not detect any significant or suggestive linkage signals between the MC1R locus and the shank skin color traits. The possible reason could be that association between the loci and trait causing loci is not maintained by linkage but by other factors such as natural selection for particular combinations of alleles, population admixture or recent mutation (Chung et al., 2007).

Previous studies reported that Id locus on Z chromosome is another major locus affecting shank color traits (Siwek et al., 2013). However, we could not detect any suggestive or significant QTL on Z chromosome. In this study, the genetic linkage map did not include Id locus. Hence, it is necessary to conduct further studies to investigate effect of Z chromosome on shank color traits in KNC. We also need to take into account that the genetic linkage map in this study only represents 26 autosomal linkage groups and Z chromosome; we could not cover 13 chromosomes in the chicken genome. Hence, the QTLs identified in the current study do not account for the whole QTL effects since our genome scan could not cover several micro-chromosomes. Additionally, there is a lack of coverage of genetic markers in some macro-chromosomes. Therefore, we cannot exclude the possibility that these un-covered genome regions may harbor additional significant QTL for traits of interests. On the other hand, the two classes of melanin (eumelanin/pheomelanin) can be quantified by a combination of chemical reaction and HPLC if we can collect skin samples (Hennessy et al., 2005). Hence, we may identify other loci involved in shank pigmentation with these novel phenotypes provided a sufficient number of samples.

To determine if the BCDO2 g.9367 A>C could be the unique functional variant underlying the identified linkage signal for a* and b*, we carried out a conditional linkage analyses that jointly accounted for association with the BCDO2 g.9367 A>C variant. Conditional on the SNP in the BCDO2, evidence for linkage to BCDO2 locus still remained (i.e., LOD score for a* [from 14.2 to 2.68]; LOD score for b* [from 47.5 to 10.4]) (Figure 2). This indicated that there should be one or more additional variants in the region that affect phenotypic variations of a* and b*.

In conclusion, this study represents the first direct evidence of the presence of a major QTL on GGA24 that affects quantitatively measured shank skin color traits. Moreover, our result based on the additional analyses including family-based test of association also showed that the BCDO2 g.9367 A>C variant is likely to account for a major QTL effect. The results presented in this study should be helpful to understand the genetic structure that underlies genetic variation of skin color-related traits.