A genetic variation map for chicken with 2.8 million single nucleotide polymorphisms

Summary

Introduction

The generation of a high quality draft sequence for the genome of chicken (Gallus gallus) is an important advance¹. Chickens are good models for studying the genetic basis of phenotypic traits, because of the extensive diversity among domestic chickens selected for different purposes. Monogenic traits are well-studied^2–4, but many interesting traits are complex and determined by an unknown number of genes. Quantitative trait loci (QTLs) have been mapped for a range of traits, including ones for growth, body composition, egg production, antibody response, disease resistance, and behaviour⁵. Determining causative genes is difficult, since each locus controls only a fraction of the phenotypic variance. We will describe a survey of the genetic variation between 3 domestic chickens and their wild ancestor. The 2.8 million single nucleotide polymorphisms (SNPs) that we identified will facilitate mapping of complex traits in many ways. First, improved marker density allows researchers to take advantage of the higher recombination rates in chicken¹, which are 2.5 to 21 cM/Mb depending on the chromosome, as compared to 1 cM/Mb for human and 0.5 cM/Mb for mouse. The previous linkage map used 2000 markers^6,7, but only 800 of these were microsatellites or SNPs, which are the most useful⁸. More importantly, our new data allow researchers to construct detailed haplotypes that segregate in different QTL crosses. Because any mutation underlying a QTL must once have originated from a single founder animal, haplotype comparisons will facilitate the fine mapping of QTLs⁹. To this end, we conduct a genomewide search for evidence of selection due to domestication, and provide an initial characterization of the expected magnitude of these effects.

Genetic variation and utility

Our experiment is outlined in Figure 1. SNPs are generated by partial sequencing at ¼ coverage for each of 3 domestic breeds (a male broiler, a female layer, and a female Silkie), and comparison of the resultant reads to the 6.6x genome for the wild ancestor of domestic chickens, Red Jungle Fowl (RJF). We expect marked heterozygosity within the 3 domestic lines, but not within RJF because the sequenced bird for the genome project is from a highly inbred line that is essentially homozygous.

Figure 1

SNP discovery experiment. We sampled 3 domestic chickens at 1/4 coverage each and compared the resultant sequence to the 6.6x draft genome of Red Jungle Fowl (RJF). Chicken photographs shown here are provided by Bill Payne (RJF), Paul Hocking (broiler), Leif Andersson (layer), and Ning Yang (Silkie).

Comparing the sequence reads for broiler, layer, and Silkie to the genome of RJF, we identified nearly a million SNPs in each instance, at mean rates of about 5 SNP/kb, as shown in Table 1. Notice that all of the “SNP rates” quoted in this paper are computed as nucleotide diversities (π), and given in units of π×10³. After correcting for SNPs detected in more than one line, there are 2,833,578 variant sites, or one potential marker every 374 bp along the 1.06 Gb genome. To assess the reliability of these data, we resequenced 295 SNPs in the same bird in which it was detected (Table S1). As many as 94% of the SNPs were confirmed. However, confirmation rates are sensitive to the functional context (e.g., coding versus non-coding) and SNPs in rare categories are less likely to be confirmed. In fact, only 83% of the non-synonymous SNPs were confirmed. Small indels of a few base pairs in length (mean of 2.3 and median of 1) are detected at rates that are well correlated with the corresponding SNP rates, but smaller by about a factor of 10.

Chicken autosomes are sorted by size into 5 large macrochromosomes (GGA1-5), 5 intermediate chromosomes (GGA6-10), and 28 microchromosomes (GGA11-38). SNP and indel rates are independent of chromosome size, as shown in Figure 2. GGA16 is the sole exception, because it contains the highly variable MHC¹⁰. This result is surprising, as recombination rates on microchromosomes are much higher than on macrochromosomes¹ and studies in other organisms exhibit a positive correlation between recombination rates and polymorphism rates^11–12. We expect that higher gene densities on microchromosomes likely counteract the effect of higher recombination rates.

Figure 2

SNP and indel rates versus chromosome number. We excluded all sequences with “random” chromosome positions. Because of the assembly problems on W, it is not shown. The rates are computed as an average of all 3 domestic lines.

SNP rates between and within chicken lines can be determined from the overlaps between reads. Table 1 demonstrates that almost every pairwise combination gives a SNP rate of just over 5 SNP/kb, except for broiler-broiler and layer-layer, which show about 4 SNP/kb, as expected since the sequenced broiler and layer are from closed breeding lines. To ensure that there are no confounding factors from the single read nature of our data, or the complexities of the overlap analysis, we used comparisons to 3.8 Mb of finished BAC sequence of a different White Leghorn¹³ from the same breed but not the same line as the layer sequenced herein. 15 chromosomes were sampled, and the results confirm our rates of 5 SNP/kb. In another study of 15 kb of introns in 25 birds from 10 divergent breeds of domestic chickens¹⁴, an autosomal rate of 6.5 SNP/kb was reported.

To quantify SNP and indel rate variation versus functional context, we considered three gene sets representing 3868 confirmed mRNA transcripts, 995 chicken orthologs of human disease genes, and 17,709 Ensembl annotations from the RJF analysis¹. Complete details for all 3 lines are tabulated in the supplements (Table S2). An excerpt for broiler is shown in Table 2. Within genes defined by mRNA transcripts, the SNP rates are 3.5, 2.1, 5.7, and 3.4 SNP/kb in 5′-UTR, coding exon, intron, and 3′-UTR regions respectively. In coding regions, indel rates are 43 times smaller than SNP rates. Ka/Ks is 0.098, similar to what is typically seen in vertebrate comparisons. We also studied “conserved non-coding regions” from the RJF analysis¹. SNP rates are similar to those of coding exons, but indel rates are intermediate to those of coding exons and UTRs, which supports the notion that these regions are functional, but may not encode proteins.

Utility of these SNPs depends on their frequency of occurrence in commonly used chicken populations. Hence, we typed 125 SNPs (including coding and non-coding SNPs, randomly distributed across the chicken genome) in 10 unrelated individuals from each of 9 divergent lines representing an assortment of European breeds. This collection includes commercial broiler and layer breeds, standardized breeds selected for their morphological traits, and an unselected breed from Iceland (Table S3). Both alleles segregated in 73% of 1113 successful marker-line combinations (out of 1125 possible combinations). Averaged minor allele frequency is 27%, but it decreases to 20% if marker-line combinations where one of the two alleles is fixed are included. This indicates that a majority of the SNPs are common variants that predate the divergence of modern breeds. Only 12% of the markers had a minor allele frequency of less than 10% in the 90 animals tested.

We demonstrate by example how these data can be used to target specific genome regions. Details of our experiments are in Supplement E (Examples). First, we consider a body weight related QTL on GGA4 that was previously mapped to a 150 cM interval^15,16. After a year of effort, where every known microsatellite (>50) was tested, 26 informative markers were developed. Further progress would have required the laborious sequencing of multiple chickens to find additional polymorphisms in this target region. With the SNP map, we selected 47 random broiler-layer SNPs, and ABI SNPlex assays were developed to type an experimental F₂ cross (n = 466). 28 (60%) of these SNPs were informative, but none had breed specific alleles, confirming that most variations predate domestication. In just one month, we doubled the number of markers, and resolved the initial QTL into two QTLs that affect body weight at 3 and 9 weeks of age.

In addition to providing markers for fine mapping, these SNPs are a rich source of candidate polymorphisms for the causative differences underlying important traits. As an example, candidate genes for disease resistance often include TGF-β^17,18, cytokines¹⁹, and the major histocompatibility complex (MHC). We thus identified 40 SNPs from the SNP map in the coding or promoter regions of 12 cytokine genes. When typed against 8 inbred layer lines, 32 of these SNPs were informative. Cytokine genes on GGA13, including IL4 and IL13, two genes that are expressed in T helper-2 (Th2) cells, drive antibody response. Four of the six SNPs that were polymorphic among lines were in IL4 and IL13, and these SNPs were fixed for different alleles in lines N and 15I, which show differential antibody response to vaccination²⁰. These SNPs therefore allow us to test whether the IL4 and IL13 loci directly determine the observed differential antibody response.

Domestication and selection

Domestic animals are useful models of phenotypic evolution under selection. The challenge is to find not only those loci that determine phenotypic differences, but also the causative alleles. We adopt two approaches, searching for evidence of selective sweeps²¹, and for non-synonymous amino acid substitutions at highly conserved sites. One example of a selective sweep is the IGF2 locus in pigs²². Given the available data, determining the exact haplotype structure is difficult, because blocks of shared alleles can be erroneously disrupted by heterozygosity of the domestic lines and by sequencing errors. However, we can still search for the local reductions in heterozygosity that accompany selective sweep, as long as we are mindful of the sequencing error rate.

We did 3-way comparisons of RJF and all possible combinations of two domestic lines. Given the limited coverage of the latter, we only examined 100 kb segments with at least 10 SNP sites where each qualifying site must have read coverage from every line. In practice, these segments contained an average of 25 to 28 SNPs. Then, we computed how often 80% or more of the SNP sites are identical in the two domestic lines but different in RJF. In Table S4, we show that 0.4 to 1.5% of the segments qualified. However, when we searched for shared alleles between RJF and one of the domestic lines, 1.2 to 2.6% of the segments qualified. Heterozygosity of the domestic lines is more of a confounding factor in searching for blocks of shared alleles between two domestic lines, versus between RJF and one domestic line. This could explain the difference, but if so, then heterozygosity of the domestic lines is the dominant factor in determining such blocks of shared alleles, not selective sweeps. Hence, selective sweeps that occurred before the divergence of modern domestic breeds must have left behind footprints that are much smaller than 100 kb. This would however be entirely consistent with the historically large effective population size of domestic chickens, and the reported high recombination rates.

For a glimpse of the true haplotype patterns, one can compare the aforementioned 3.8 Mb of finished BAC sequence, from the second layer line (L2), to the genome of RJF. These results are overlaid alongside the primary SNP data set in Figure 3. Short RJF-type fragments can be seen in all 4 lines. Shared domestic-type fragments can also be seen, but at sizes of 5 to 15 kb. This is consistent with our inability to detect footprints of selective sweeps at length scales of 100 kb and suggests that a better choice is 10 kb. However, our data are insufficient for such a genomewide analysis.

Figure 3

Detailed haplotype patterns in 3 regions, each covered by 2 overlapping BACs from the second layer line (L2). The primary SNP data are labeled B (broiler), L1 (layer), and S (Silkie). All comparisons are to RJF, and we show only those sites where a SNP is identified in at least one of the 4 lines. Hence, the horizontal scale is linear in the number of SNP sites, but non-linear for size. BLUE colors indicate where a particular line agrees with RJF, while RED colors indicate where it does not. Overlapping BACs on GGA1 and GGA7, but not GGA14, are clearly from different haplotypes.

It has been hypothesized that loss-of-function (LOF) mutations have accumulated in domestic animals, as the result of relaxed purifying selection and selection for adaptive benefits²³. An example of the latter is the deletion in the myostatin gene in cattle selected for muscularity²⁴. Such deletions are rare, and so we looked for non-synonymous SNPs at highly conserved sites using SIFT​²⁵. Every substitution is thus classified as being likely to affect function (intolerant) or not (tolerant). For genes defined by mRNA transcripts, 26% of testable SNPs are intolerant, although only 11% are intolerant if we restrict this to high confidence assessments (Table S5). Usually, it is the domestic allele that is intolerant, but we would emphasize that intolerant SNPs are rare, and only 59% were confirmed by PCR resequencing. Given that the domestic allele is represented by a single read, as opposed to 6.6 for the wild allele, much of this effect is likely due to sequencing errors. However, we noticed the same effect in 424 non-synonymous SNPs that we identified from an analysis of 330,000 ESTs, where every allele was seen in two or more ESTs. We conclude that the LOF hypothesis remains intriguing, but any effect is likely to be small.

Some of the experimentally confirmed SIFT intolerant SNPs could be functionally important. We show one example in Figure 4, from the ornithine transcarbamylase (OTC) gene. It substitutes glycine in RJF to arginine in layer and broiler. This SNP is identical to the G188R substitution associated with hyperammonaemia in humans²⁶. Resequencing of additional domestic birds revealed a high frequency for the intolerant allele in both White Leghorns (p=0.65, n=20) and in broilers (p=0.75, n=6). In mammals, OTC is expressed in the liver and catalyzes the second step of the urea cycle. Chicken OTC is expressed in the kidney and exhibits a low enzymatic activity, with substantial variability among breeds²⁷. Preservation and sequence conservation of OTC, along with all other enzymes in the urea cycle¹, was unexpected because avian species excrete uric acid (not urea) as their primary component of nitrogenous waste, and were believed to be lacking a functional urea cycle. The deleterious nature of human G188R makes this an attractive candidate for phenotypic studies of avian-specific adaptations in the urea cycle.

Figure 4

Multi-species alignments for ornithine transcarbamylase (OTC), indicating non-synonymous substitutions relative to human protein. SIFT intolerant position is indicated by site number and bold-faced lettering. WT=wild type. MUT=mutant.

Discussion

This analysis has provided the first global assessment of nucleotide diversity for a domestic animal in comparison to a representative of its wild ancestor. The small number of birds sequenced is compensated for by the vast number of sites examined. We detected surprisingly little difference in diversity in comparisons between RJF and domestic lines, between different domestic lines, and within domestic lines. The total rates are typically 5 SNP/kb, with the only exception being a slight reduction to 4 SNP/kb in broiler and layer lines that are maintained as closed breeding populations. In comparison, 5 SNP/kb is 6 to 7-fold larger than humans²⁸ and domestic dogs²⁹, 3-fold larger than gorillas³⁰, but similar to the diversity between different mouse subspecies³¹.

Most of the nucleotide diversity observed between and within domestic lines must have originated prior to the domestication of chickens 5,000 to 10,000 years ago. Given a neutral substitution rate of 1.8×10^-9 sites per year for galliform birds³², we estimate that a coalescence time of 1.4 million years would be required to account for the observed rates of 5 SNP/kb. Considering that the rates observed between RJF and domestic lines are not much higher than those between domestic lines, it would seem that domestication has not resulted in a substantial genomewide loss of diversity, as would be expected had a severe population bottleneck occurred. This is important, because it contradicts the assertion that animal domestication began from a small number of individuals in a restricted geographic region³³. That is still a possible scenario for the very earliest phases of domestication, but if so, our data imply that subsequent crossing with the wild ancestor (in the first thousand years until more developed breeds were established) restored this diversity. Nevertheless, extensive diversity is consistent with the ongoing improvements in agricultural traits that have been achieved over the last 80 years, in layer and broiler lines³⁴.

The most important application for this SNP map will be in analysis of QTLs and other genetic traits. Although the density of markers far exceeds what is needed for initial mapping, the principal challenge is not in the detection of linkage but in the identification of genes underlying QTLs⁹. By itself, our SNP map is not adequate. It must be combined with novel strategies and novel resources (like mapping populations specifically designed for fine mapping). The essential problem is the lack of a one-to-one relationship between genotype and phenotype, as the latter is influenced by multiple genetic and environmental factors. This can be overcome, in experimental and domestic animals, by progeny testing and segregation analysis, which permit detailed characterization of haplotypes associated with different QTL alleles, and may eventually lead to the identification of the underlying causative mutations²². This SNP map will facilitate fine mapping.

As an example, the major Growth1 QTL on GGA1 explains about one third of the difference between RJF and White Leghorn in adult body weight and egg weight³⁵. Initial mapping assigned this locus a ∼20 cM confidence interval. Selective back-crossing using sires that have recombinant chromosomes, and QTL analysis using subsequent intercross generations, are currently employed to refine the localization to a few cM, expected to be less than ∼1 Mb. This also establishes a collection of chromosomes of known QTL status. Our SNP map can then be used for haplotype analysis, assuming that the White Leghorns share a chromosomal segment – identical by descent (IBD) – with the causative mutation. The small haplotype blocks detected in this study underscore the need for a larger number of SNPs to identify such IBD segments. Although these small blocks may require greater marker density and more recombinants to identify the causative haplotype, less effort will be required to resolve the actual QTL alleles once the haplotype is found.

Materials and methods

Our broiler and layer lines are from European breeds with dramatic differences in meat and egg production traits. This specialization started only during the first half of the 20^th century³⁶. The sequenced male White Cornish-type broiler is from a closed breeding population commonly used in the production of commercial meat-type hybrids (Aviagen, Newbridge, Scotland); effective population size is about 800. The female White Leghorn layer is from a closed line developed at Swedish University of Agricultural Sciences³⁷; its effective population size has been 60 to 80 birds for the past 30 years. The Chinese Silkie is used in meat/egg production and traditional Chinese medicine³⁸. Selection intensity has been low, and the sequenced female is from a large outbred population.

DNA was extracted from erythrocytes of a single bird, sheared by sonication, and size fractionated via agarose gels. Fragments of 3-kb size were ligated to SmaI-cut blunt-ended pUC18 plasmid vectors. Single colonies were grown overnight, and plasmids were extracted by an alkaline lysis protocol. Sequences were read from both ends of the insert, with vector primers and Amersham MegaBACE 1000 capillary sequencers. Roughly one million reads were generated for each bird. For broiler, layer, and Silkie, we got a total of 841,790, 841,555, and 870,556 successful reads, whose Q20 lengths add to 380,729,199-bp, 372,263,344-bp, and 397,831,117-bp, respectively.

To minimize sequencing errors, we use the Phred quality, Q^39,40. This is related to the single base error rate by the equation: -10×log₁₀(Q). We use more stringent thresholds than normal⁴¹, with Q>25 for the variant site and Q>20 in both flanking 5-bp regions. For an insertion-deletion (indel), the variant site in the shorter allele is given the quality of its two flanking bases. We originally found many artifactual deletions relative to RJF, which upon a closer examination of the sequence reads were due to doublet peaks that got called as singlet peaks. This is an unavoidable flaw of the base caller software. Hence, we raised the indel thresholds to Q30 and Q25. We must still advise caution, and to that end, indels in simple repeats are flagged and none are counted in our summary tables.

Paralog confusion is detected in the course of the genome level BlastN search that determines where the read is supposed to go. Once this is known, the detailed alignments are done within CrossMatch​⁴². Analysis of the RJF genome¹ shows that recent segmental duplications typically agree to 2%. When the best and second best BlastN hits were more than 2% apart, and the best hit was not to a known segmental duplication, the best hit was taken. When either rule was violated, clone-end pairs information was used to resolve the ambiguity. Every alignment had to incorporate 80% of the read. Mapped back to the RJF genome, the amount of usable data for broiler, layer, and Silkie covered 190,513,980-bp, 165,154,746-bp, and 210,214,479-bp respectively.

Polymorphism rates are normalized to the length of the sequence on which we can detect SNPs. To correct for heterozygosity within a line, we compute nucleotide diversity using the approximation⁴³:π=K/∑n−1i=1Li, where K is the number of variant sites found by sequencing n chromosomes in a region of length L. When comparing RJF to one of the 3 domestic lines, n can only be 2 or 3, and it is a stochastic variable, because there is a 50% chance that any two overlapping reads are from the same chromosome. When there are m overlapping reads, the denominator is L2m−1⋅(1+(2m−1−1)⋅(1+12)). We then sum over all possible regions, with different L and m for each region, to get what we call the “effective length”. Similar considerations are used to compute SNP rates within a line, except that n is 1 or 2, and as a result, the denominator becomes L2m−1⋅(2m−1−1).

We compute gene context relative to 5 different data sets. The first 3 are based on experimentally derived genes and the last 2 are based on computer annotations. Riken1 is a set of 1758 full-length cDNAs taken from bursal B-cells of a two week old CB inbred⁴⁴. GenBank refers to 1178 chicken genes with “complete CDS” designation, downloaded as version 2003-12-15. BBSRC is a set of 1184 cDNAs, taken from a larger group of 18,034 cDNAs⁴⁵, which are full-length using a TBlastX mapping to vertebrate Refseq and BlastX mapping to SWALL. Merging all 3 data sets, we have 3868 non-redundant genes. For the detailed gene models, we do a genome level search in BLAT​⁴⁶ and use SIM4​⁴⁷ to compute the exon-intron boundaries. The last two data sets are for 995 chicken orthologs of human disease genes and 17,709 non-redundant Ensembl genes.

Additional details are in Supplement M (Methods).

Supplementary Material

Acknowledgments

International Chicken Polymorphism Map Consortium

Footnotes

References