A synergistic combination of two next-generation sequencing platforms with a detailed comparative BAC physical contig map provided a cost-effective assembly of the genome sequence of the domestic turkey (Meleagris gallopavo). Heterozygosity of the sequenced source genome allowed discovery of more than 600,000 high quality single nucleotide variants. Despite this heterozygosity, the current genome assembly (∼1. 1 Gb) includes 917 Mb of sequence assigned to specific turkey chromosomes. Annotation identified nearly 16,000 genes, with 15,093 recognized as protein coding and 611 as non-coding RNA genes. Comparative analysis of the turkey, chicken, and zebra finch genomes, and comparing avian to mammalian species, supports the characteristic stability of avian genomes and identifies genes unique to the avian lineage. Clear differences are seen in number and variety of genes of the avian immune system where expansions and novel genes are less frequent than examples of gene loss. The turkey genome sequence provides resources to further understand the evolution of vertebrate genomes and genetic variation underlying economically important quantitative traits in poultry. This integrated approach may be a model for providing both gene and chromosome level assemblies of other species with agricultural, ecological, and evolutionary interest.
The rapid and continuing development of next-generation sequencing (NGS) technologies has made it feasible to contemplate sequencing the genomes of hundreds—if not thousands—of species of agronomic, evolutionary, and ecological importance, as well as biomedical interest [1]. Recently, a draft genome of the giant panda was described, based solely on Illumina short read sequences [2]. Below, we describe the genome sequence of the turkey (Meleagris gallopavo) determined using primarily NGS platforms. In this case, however, a combination of Roche 454 and Illumina GAII sequencing was employed. While this approach presented unique assembly challenges, the turkey sequence benefits from the particular advantages of both platforms. In addition, unlike the case for the panda, this novel approach allowed us to use a BAC contig-based physical and comparative map, along with the turkey genetic map [3] and the chicken genome sequence [4], to align the turkey sequence contigs and scaffolds to most of the turkey chromosomes. Such an alignment is essential for making long range evolutionary comparisons and employing the sequence to improve breeding practices using, for example, genome-based selection approaches, where chromosome locations are critical. The high throughput and low cost of NGS technologies allowed sequencing the turkey genome at a fraction of the cost of other recently reported genomes of agricultural interest (bovine and equine) [5], [6]. The draft turkey genome sequence represents the second domestic avian genome to be sequenced, and this permits a genome-level comparison of the two most economically important poultry species. When added to the recently published zebra finch genome [7], analysis of the three avian genomes reveals new insights into the evolutionary relationships among avian species and their relationships to mammals. Turkeys, like chickens, are members of the Phasianidae within the order Galliformes. One estimate [8] is that the last common ancestor of turkeys and chickens lived about 40 million (M) years ago; however, other estimates are more recent [9], [10]. Comparison of the turkey genome to that of the chicken provides the opportunity for high resolution analysis of genome evolution within the Galliformes. The turkey has 2n = 80 chromosomes (chicken has 2n = 78) and, as for most avian species, the majority of these are small “microchromosomes” that cannot be distinguished by size alone. Although most turkey chromosomes are syntenic to their chicken orthologues, the chicken chromosome GGA2 is orthologous to two turkey chromosomes, MGA3 (GGA2q) and MGA6 (GGA2p), due to fission at or near the centromere, while GGA4 is orthologous to MGA4 (GGA4q) and MGA9 (GGA4p) [10], [11].
Generally, DNA from a single inbred animal is preferential for sequencing to minimize polymorphism. For the turkey, however, such an option is not available, and thus we sequenced DNA from “Nici” (Nicholas Inbred), a female turkey, which is also the source DNA for the two BAC libraries that have been characterized [12]. Nici is from a subline (sib-mating for nine generations) originally derived from a commercially significant breeding line, but her genome is still extensively heterozygous. A side benefit of this approach was the concomitant identification of extensive and commercially relevant single nucleotide polymorphism (SNP) data, as discussed below. With the exception of the BAC end sequences (BES) used only for chromosome alignment, the sequence data used for this assembly came solely from two sequencing platforms: the Roche/454 GS-FLX Titanium platform (454 Life Sciences/Roche Diagnostics, Branford, CT) and the Illumina Genome Analyzer II (GAII; Illumina, Inc. , San Diego, CA). The 454 data were generated using the latest “Titanium” protocol at Roche and the Virginia Bioinformatics Institute (Virginia Tech) and included both unpaired shotgun reads and paired-end reads produced from two libraries with estimated 3 kilobase pair (Kbp) and 20 Kbp fragment sizes. The 454 runs yielded approximately 3 M read pairs from the 3 Kbp library (average usable read length 180 bases), 1 M read pairs from the 20 Kbp library (average length 195 bases), and 13 M shotgun reads (average length 366 bases). The Illumina sequencing data were generated at the USDA Beltsville Agricultural Research Center and the NIH National Institute on Aging from both single and paired-end read libraries with a 180 bp fragment size for the paired reads. Details on the sequence data are presented in Table 1. These data represent approximate 5× genome coverage in 454 reads and 25× coverage in GAII reads, assuming a genome size similar to that of the chicken at 1. 1 billion bases [4]. In addition, BACs used to generate the 40,000 BES alignment markers by traditional Sanger sequencing spanned ∼6× clone coverage of the genome. Since female DNA was used, coverage of the Z and W sex chromosomes was half that of autosomes; therefore the assembly of both these chromosomes was poor. A modified version of the Celera Assembler 5. 3 [13], [14] was used to produce the contigs and scaffolds in the assembly (see Methods for details). The initial assembly contained 931 Mbp of sequence in 27,007 scaffolds with N50 size of 1. 5 Mbp. The span of the scaffolds was 1. 038 Gbp. The scaffolds contained 145,663 contigs with N50 size of 12. 6 Kbp. The assembled scaffolds were then ordered and oriented on turkey chromosomes using a combination of two linkage maps and a comparative BAC contig physical map. The first turkey linkage map [3] had 405 chicken and turkey microsatellite sequences that mapped to the assembled scaffolds. The second linkage map, based on segregation of SNPs in a different population [15], had 442 SNP markers mapped to the scaffolds. The comparative chicken-turkey physical map [16] provided turkey chromosome positions for 30,922 BES found in scaffolds. Comparison of scaffolds to the marker map resulted in splitting only 39 scaffolds due to inconsistencies between the assembled scaffolds and marker positions on the chromosomes. A total of 28,261 scaffolds containing 917 Mb of sequence were assigned to chromosomes (Table 2). Included in this number were 7,080 single-contig scaffolds that represented repetitive sequences but that could be linked to non-repetitive scaffolds. The remaining 5,858 scaffolds were pooled to form ChrUn (unassigned) which contains 19 Mb of sequence in comparison to about 64 Mb on the current chicken chr_Un. Analysis of the assembled contigs showed that 4. 6% of the sequence was covered only by reads from a single sequencing platform, with 2. 3% covered exclusively by each. If the reads covered the genome uniformly, one would expect to have missed only 0. 67% of the genome with Roche/454 and 0. 0006% with Illumina. The distribution of regions of exclusive coverage for both platforms (Figure S1) shows there was a large number of short (<20 bp) gaps in coverage by Illumina sequencing, whereas the Roche/454 coverage gaps tended to be larger. Mean sequencing gaps were 46 bases for Illumina reads and 72 for the Roche/454 coverage. Coverage biases previously have been shown for both platforms [17], but fortunately, the biases are relatively orthogonal. Therefore, it is definitely beneficial to use data from both platforms in de novo assemblies. The draft turkey assembly was compared to the chicken genome assembly (2. 1), which was sequenced and assembled using traditional Sanger sequencing [4]. Table 3 illustrates that assembly of NGS sequence data, although feasible, does not produce contigs and scaffolds as large as those expected from an assembly based on Sanger sequencing. However, the relatively low cost of NGS sequencing (<$250,000 for the turkey) makes such projects feasible for species with more focused interest groups and facilitates for resources to be directed toward genome analysis and interpretation as opposed to generating raw sequence data. However, chromosome assemblies currently still require the integration of multiple data types including shotgun reads and contigs, genetic linkage maps, BAC maps and BES, and cytogenetic assignments. The challenge was to develop databases and software to achieve this goal. Integrity of the assembly was validated by mapping the assembled turkey scaffolds to 197 Kbp of finished BAC sequence containing part of the MHC B-locus, GenBank accession DQ993255. 2. The average sequence similarity was over 99. 5% and no inconsistencies in the 21 scaffolds that mapped to that region were observed. The extent of the genome coverage could be estimated both from the total span of the assembled scaffolds and from portions of the chicken genome with syntenic matches to the turkey scaffolds. Both methods produced consistent estimates of the size of the euchromatic portion of the turkey genome at about 1. 05 Gbp. With 936 Mbp of sequence in the final chromosomes, including ChrUn, the assembly encompasses an estimated 89% of the total sequence of the genome. One of the striking observations in the chicken genome sequencing project was the difficulty obtaining sequences for specific regions, including the 10 smallest microchromosomes [4]. For example, the chicken genome lacks sequence orthologous to human chromosome 19q. Remarkably, these sequences appeared to be absent not only from the shotgun clone libraries used to generate the whole genome shotgun (WGS) reads but also from all available BAC libraries [18]. Although these regions have high GC content, it is unclear why this region of the genome is resistant to cloning in E. coli. In general, BAC coverage of microchromosomes is less than macrochromosomes in both chicken and turkey BAC libraries, although the HSA19q orthologues are an extreme example of a missing syntenic region. Since the turkey genome was sequenced without any cloning step, the assembly was tested for representation of HSA19q orthologous sequence. Presence of sequences was verified by performing a BLAT analysis of the complete HSA19q sequence against the turkey and chicken genomes (Table S1). Surprisingly, regions orthologous to HSA19q were not represented at a higher frequency in the turkey assembly versus the chicken assembly. As was observed in the chicken, regions orthologous to HSA19p and a small syntenic region from HSA19q are covered well in the turkey assembly (MGA30 and 13, respectively). These results suggest that absence of HSA19q orthologous sequences is not due to the high GC content, in that Illumina sequences show a bias towards higher coverage of GC rich regions [19], [20]. The identification of a single BAC clone that hybridizes across the entire length of a single microchromosome in chicken [21] suggests that the occurrence of microchromosome-specific repeats might be a more likely explanation for the absence of these sequences using both traditional Sanger sequencing as well as NGS technologies. Heterozygous alleles, including both SNPs and single nucleotide insertions and deletions (indels), were detected by scanning the assembled contigs for positions where the underlying reads significantly disagreed with the consensus base [22]. A previous study cataloging heterozygous alleles from assembled shotgun reads within an individual human genome used a similar approach, augmented with a set of quality criteria used to distinguish genuine biological variations from sequencing error [23]. Following this approach, a set of quality criteria was developed and implemented within the assembly forensics toolkit [24]. Two classes of SNVs were catalogued: (1) those with abundant evidence, called strong SNVs (601,490 SNVs), and (2) a more inclusive set called weak SNVs (920,126 SNVs total). In the turkey genome, transitions were roughly 2. 4× more common than transversions: 295,055: 122,731 for strong SNVs and 466,629: 200,743 for all SNVs. Many single base indel positions were detected: 183,215 of 601,490 strong SNVs, and 249,512 out of all 920,126 SNVs. A very small number of SNVs (489 strong, and 3,242 all) were detected with more than two well-supported variants, suggestive of unfiltered sequencing errors or collapsed repeats. The depth of coverage for strong SNVs ranged between 6 and 30 with mean and standard deviation of 15. 3±5. 3, while the depth of coverage for all SNVs ranged between 4 and 5,319 with mean and standard deviation of 41. 4±134. 6. The very high coverage regions are highly likely to be due to collapsed near-identical repeats. Annotation of the turkey genome sequence identified a total of 15,704 genes (Table S2) of which 15,093 were distinct protein coding loci and 611 non-coding RNA genes. In addition, multiple distinct proteins produced by alternative splicing were identified for some loci, giving a total of 16,217 distinct protein sequences. Orthologs between turkey, chicken, and human proteins were defined using sequence homology, phylogenetic trees, and conservation of synteny. All gene annotations are available from the Ensembl genome browser version 57 (http: //e57. ensembl. org). The draft turkey genome assembly was used to test the distribution of nucleotide diversity across the turkey genome by aligning SNPs covering ∼3. 97% of the genome identified through resequencing a reduced representation library from commercial turkeys [15]. Substantial deviations were observed between regions in the genome. Chromosome Z showed the lowest nucleotide diversity, about half (θ = 0. 000273) that of the autosomes, which is likely the result of a lower effective population size of this chromosome and lower recombination rate (Figure S2) [25]. The five largest chromosomes had similar nucleotide diversities as the microchromosomes. Given the higher recombination rate on the microchromosomes, the ensuing higher mutation rate [26], and lower susceptibility to hitchhiking effects, equal rates of nucleotide diversity between micro- and macrochromosomes may seem unexpected. However, these findings are in line with observations in the chicken [27] and may be explained by higher gene density and higher purifying selection on the microchromosomes. Within chromosomes, extended regions of low nucleotide diversity were detected, many of which coincided with centromeres. Comparisons of gene family assignment statistics for the turkey and chicken genome assemblies are shown in Table S4. Although the draft turkey sequence has fewer genes than the current chicken genome build (2. 1), part of the difference may be due to cutoff values used by annotation groups resulting in variation in gene number. Even with this caveat, more than half of the gene families show no change in copy number between them (Table S5a–d). Overall, most families exhibiting variation have general regulatory functions related to transcription, metabolism, cation transport, cell-cell signaling, and cell development or differentiation (Figure 3). Distinct keratin families, encoding major structural proteins of chicken feathers, claws, and scales, have undergone uneven expansion or contraction with considerable variation in number among species. More than half of the innovation families (found in turkey but not chicken) have unknown functions, are singletons, and were annotated by mapping to the zebra finch protein prediction. Species-specific gene families in birds and mammals are summarized in Tables S6, S7. Of these, 881 are specific to turkeys and chickens and 271 specific to birds. The inference for bird-specific functions is of relatively high quality since the likelihood that a bird gene is not simultaneously found in all 13 non-bird species is low. Most of the rapidly evolving gene families in birds have unknown functions. Approximately 83% of the turkey/chicken-specific families and 71% of the bird-specific families have unknown functions. For the remaining families, most have well-defined roles (Table 5). Families related to egg formation (such as avidin, ovocalyxin, and vitellogenin) and scavenger receptors were identified as avian specific in the present and previous analyses [4]. Examination of gene family sizes between the avian species and the platypus, an egg laying mammal, found two egg-related gene families [egg envelop protein (ENSFM00500000271806) and vitellogenin, an egg yolk precursor protein (ENSFM00250000000813) ] to be conserved among the four egg-laying species. Both of these gene families are absent from eutherians. Other gene families specific to egg-laying species (birds and platypus) are mainly related to protein metabolism, cell-cell communication, and regulatory functions. Several other proteins related to egg formation, such as avidin and ovocalyxin, are found in birds but not in platypus. In contrast to unique gene families, only 70 families were completely absent in both the turkey and chicken (33 in all birds) compared to the non-avian species. These include the gene family associated with enamel formation (ENSFM00250000008876, an enamelin precursor related to teeth), which is completely lost in the three avian species. Genes encoding the vomeronasal receptors and several casein related families are also completely absent in the avian species. Several olfactory receptor families specific to mammals are either absent or dramatically reduced in birds. Interestingly, the olfactory receptor 5U1 and 5BF1 gene families, reported to be dramatically expanded in chicken as compared to humans and flies [4], is contracted in turkey. Lineage events in the turkey, chicken, and zebra finch genomes reveal significantly higher synonymous substitution rates on microchromosomes than macrochromosomes (Figure 4a), with a clear inverse relationship with chromosome size. This suggests that genes on the microchromosomes are exposed to more germ-line mutations than those on other chromosomes [38]. However non-synonymous mutation rates do not seem to vary so widely and when combined show the dN/dS ratio (a measure of selection) to increase with chromosome size. These results are consistent with the prediction that the higher synonymous substitution rates of microchromosomes combined with the “Hill-Robertson” effect [39] of higher recombination rates on these smaller chromosomes increases purifying selection [40] on the microchromosomes (Figure 4b and 4c). Theory predicts natural selection to be more efficient in the fixation of beneficial mutations in mammalian X-linked genes than in autosomal genes, where hemizygous exposure of beneficial non-dominant mutations increases the rate of fixation. This “fast-X effect” should be evident by an increased ratio of non-synonymous to synonymous substitutions (dN/dS) for sex-linked genes. As shown in Figure 5, there is solid confirmation of the predicted rapid evolution in the sex-linked genes based on turkey, chicken, and zebra finch genome-wide data. These results confirm that evolution proceeds more quickly on the Z chromosome [41], where hemizygous exposure of beneficial non-dominant mutations increases the rate of fixation. Based on the analysis of differentially evolved genes, 428 and 257 genes were identified as being under accelerated evolution in the turkey and chicken lineages, respectively. Most of the accelerated genes in the turkey lineage have gene ontology (GO) terms related to DNA packaging and regulation of transcription (Figure 6a). In contrast, a large proportion of the accelerated genes in the chicken lineage have GO terms related to negative regulation of cellular component organization and biogenesis, proteolysis, interphase, and cell cycle arrest (Figure 6b). The enrichment of KEGG pathways using DAVID supports the GO term analysis (Table S8). These results suggest that genes with a role in transcriptional regulation are key in the evolution of the turkey, whereas genes involved in protein turnover and cell proliferation have been more important in the evolution of the chicken. For genes classified as innate immune loci by InnateDB (www. innatedb. ca), dN/dS ratios were calculated for each pair of species (turkey-chicken, turkey-zebra finch, etc.) and then compared with ratios for non-immune genes. Innate immune genes showed lower dN/dS ratios than other genes in all species-pairs of mammals and birds, except between turkey and chicken where the values are essentially equal (Figure 7). Using Wilcoxon rank sum test, it is obvious from the comparisons that the innate immune-related genes have been under more purifying selection than non-immune-related genes except between turkey and chicken (Table S9). Evolution of genes of the innate immunity system is thought to be continuous and under balancing selection [42]. However, purifying selection under the same conditions may be the dominant force acting on the vast majority of genes that function within the innate immune system [43]. Although only innate immune genes are under purifying selection by functional constraints, they are also more constrained than other genes. This relationship supports the view that the ancient innate immune system has a highly specialized function, critical for the recognition of pathogens and thus should be under purifying selection. However, unlike other species, the dN/dS ratios for innate immune genes between turkey and chicken are similar to other genes. Perhaps the adaptation of turkey and chicken to different ecological niches has exposed them to new pathogenic environments with potentially lethal pathogens having exerted selective pressures on their genomes. This thesis would suggest that there was a period of accelerated evolution of the innate immunity system after the divergence of these species 30–40 M years ago. The availability of the turkey genome for comparison to the chicken [4] and zebra finch [7] allows for interrogation of the immune gene repertoire. In general, homologs for all the innate immune gene families were found (Table 6), with smaller gene families present in birds. This finding is consistent with earlier comparisons of mammalian with the chicken genome [44] and provides greater evidence of an avian-wide phenomenon. Examples include the chemokines, TNF superfamily, and pattern recognition receptors. Inflammatory CCL chemokines, which occur in all avian and mammalian species, fall into two multigene families (MIP and MCP; Figure S4). There are four MIP family members in the chicken and the zebra finch (CCLi1–4), yet only three family members in the turkey genome build (CCLi2–4). For the MCP family, there are six (CCLi5–10), three (CCLi5–7), and five (CCLi5–7 and 9–10) members in the chicken, zebra finch, and turkey genomes, respectively. The chicken genome sequence lacks TNFSF-family members TNFSF1 and TNFSF3 [44]. Presence of these lymphotoxins controls lymph node formation in mammals [45]; however, lymph nodes are absent in birds [46]. Therefore, it was not surprising that these genes were not found in any of the three avian genomes. In contrast, lack of TNFSF2 (TNFA) was unexpected, since it is found in many fish species [47], and there are several reports of TNF-alpha-like activity in chickens [48]. A sequence homology search in the three avian species only detected TNFSF15, a close relative of TNFSF2. Loss of TNFSF1,2, and 3 (as well as TNFSF14) in the avian lineage could explain these observations (Figure S5). Absence of specific genes from the three avian genomes further implies that particular genomic regions are intrinsically difficult to clone and/or sequence with the traditional Sanger and NGS methods. Finally, clear differences between birds and mammals exist in the size of the pattern recognition receptor families. For example, there are only six NODLR family members in each of the three avian species, in contrast to 22 and 32 in human and mouse, respectively (Table 6 and Figure S6). These are cytoplasmic receptors that recognize a range of ligands that activate caspases, and elicit an inflammatory response. A recent analysis revealed hundreds of NODLR genes in fish [49] with homologs of all mammalian genes. It is therefore clear that NODLR genes were lost during the evolution of the avian genomes. In contrast, while similar numbers of TLRs are found in birds and mammals, evolutionary histories of gene gain, loss, and conversion are complex (Figure S7) [50]–[52]. The avian TLR1A/B and TLR2A/B genes are orthologs of mammalian TLR1/6/10 and TLR2, respectively. All three birds have lost TLR8 and 9 but retained TLR7. The avian TLR21 is the ortholog of mouse TLR13, which was lost in the human lineage, and TLR15 appears to be unique to the avian lineage. Approximately 6. 94% of the turkey genome consists of interspersed repeats, most of which belong to three groups of TEs, the CR1-type non-LTR retrotransposons, the LTR retrotransposons, and the mariner-type DNA transposons (Table 7 and Dataset S1). The CR1 group of TEs is the most abundant, occupying 4. 81% of the genome, which is likely an underestimate because a number of highly degenerate and low copy number CR1-type elements remain to be characterized. Overall, the turkey and chicken genomes are very similar with respect to repeat content and the types of predominant TEs [4], [53] with high sequence similarities between major TEs. For example, CR1_B in turkey and chicken share ∼91% nucleotide identity over a 2 Kbp region, the Birddawg_I LTR retrotransposons share ∼89% identity over a 3. 6 Kbp region, and the mariner transposon Galluhop shares ∼91% identity over the entire 1. 2 Kbp of the full-length element. Similar to the chicken, the Galluhop repeat in turkey is associated with a deletion derivative of ∼550 bp. Repetitive sequences are among the fastest evolving sequences in the genome. Therefore, the conservation of the repeat elements and sequences between the turkey and chicken is indicative of very stable genomes given a divergence time of 30–40 M years. Genome projects enable the collection of large supermatrices of alignable nucleotide sequences for phylogenetic analysis. Galliform phylogeny was re-examined by collecting sequences from the turkey and chicken genomes for 42 loci. These sequences were assembled into the largest supermatrix available for the order, containing 83 galliform species representing 73 genera, with three anseriform outgroup species. With several whole mitochondrial sequences, two genomes, and repeated use in multiple studies, 37 taxa were represented by 11 or more loci, and 12 taxa by more than 20 loci, providing data-rich anchor points that bridged locus sets throughout the tree. For the turkey, the main finding was its close relationship with the Central American ocellated turkey Agriocharis (Meleagris) ocellata (94% bootstrap support) and the relation to the grouses within the phasianids (Figure S10). The turkey-grouse clade has been recovered in several [59]–[61] but not all previous multi-locus studies. The average bootstrap support for the nodes was high and the topology reproduced many features of previous studies, with monophyly of the megapodes, cracids, numidids and odontophorids, and polyphyly of the Perdicinae and Phasianinae within the phasianids. Grouping of an African bird (Ptilopachus petrosus) traditionally classified as a phasianid with the New World quails as recently observed [59] is supported, with the three loci independently reproducing this clustering. The same was true when P. nahanii was used instead of P. petrosus. Polyphyly of Francolinus was expected [62]; however, the implied polyphyly of Lophura was not. Increased throughput and decreased costs of NGS technologies facilitate cost- and time-effective sequencing of genomes. The turkey genome sequence described herein represents the first eukaryotic genome completely sequenced and assembled de novo from data produced by a combination of two NGS platforms, Roche-454 and Illumina-GAII. This genome project is a first where the majority of the production cost was invested in analysis and interpretation rather than generating sequence, and that the assembly is comparable in genome coverage to the predominantly Sanger-based sequences of the chicken and zebra finch. The sequence assigned to the chromosomes covers approximately 93% of the turkey genome. The quality of this sequence makes it a valuable resource for comparative genomics including identification of thousands of SNVs amenable to whole genome analyses. The turkey sequence confirms and extends the previously known high synteny between the turkey and chicken genomes [3]. These two avian species are remarkably similar with only 30 predicted rearrangements (mainly small inversions) distinguishing their genomes, despite last sharing a common ancestor about twice as long ago as the common ancestor of mice and rats or humans and gibbons. Chromosome rearrangements that occurred show a trend towards more acrocentric chromosomes in the turkey than in the chicken. The stability of galliform genomes is further confirmed by the overall conservation of gene sequences and repeat families. At less than a third the size of mammalian genomes, a greater proportion of the turkey genome (∼10%) is under selective constraint versus mammals where the fraction of conserved nucleotides is approximately 5%. This also reflects the reduced percentage of the turkey genome comprised of interspersed repeats (7%). Whereas genomes of close relatives allow for analysis of rapidly changing sequence, those of distant species help elucidate regions conserved during vertebrate evolution. Gene families present only in birds provide a broad perspective on lineage-specific evolution. For example, variation in gene content between birds and an egg-laying mammal (platypus) shows functions shared by egg-laying animals in general as well as those unique to egg-laying birds. Likewise, genes specific to mammalian characteristics such as tooth formation have been lost in avian species. Some gene families such as TLRs of the innate immune system show complex evolutionary histories of gene gain, loss, and gene conversion between mammalian and avian species. The adaptive immune system is a relatively recent innovation peculiar to the vertebrates and provides a valuable framework for genome comparisons [63]. Genes involved in the control and regulation of the immune response towards invading pathogens are subject to strong selective pressures: the so-called “arms race” between pathogen and host. The result has been exceptional sequence divergence between the immune genes of vertebrate species, in particular those between birds and mammals [64]. Additionally, many immune genes belong to gene families that have been subject to lineage specific expansions and contractions, facilitating the evolution of new functions to combat pathogenic challenges. There are many fundamental differences between the immune systems of birds and mammals, including the major histocompatibility complex (MHC) structure [65], absence of lymph nodes in birds [46], and different mechanisms of somatic recombination in the generation of antibody diversity [66]. From an evolutionary perspective, the turkey and chicken provide an interesting case for comparative study. These two genomes have undergone intense artificial selection in recent decades for similar production traits, yet their differentially evolved genes included more functioning in transcriptional regulation in turkey, and more functioning in protein turnover and cell proliferation in chicken. Comparative genomics can provide additional insights into the response of the galliform genomes to this recent period, as well as to their longer histories of domestication. The turkey genome sequence can enhance the discovery of genetic variations underlying economically important quantitative traits, further maximizing the genetic potential of the species as a major protein source.
Vertebrate whole genome sequence assembly is aided by decreased variability in the target genome. To this end, a female turkey “Nici” (donated by Nicholas Turkey Breeding Farms) identified as NT-WF06-2002-E0010 was chosen for sequencing; Nici is also the source DNA for the two BAC libraries that have been characterized [12]. Nici is from an inbred sub-line (i. e. , sib-mating for nine generations) originally derived from a commercially significant breeding line. From her pedigree, Nici has an increased inbreeding coefficient of 0. 624 relative to the founder breeding line. As a prelude to initial genome sequencing, heterozygosity of Nici was compared with that of individuals from several breeder lines by genotyping 147 randomly distributed microsatellites [12]. Mean heterozygosity for Nici was determined to be 0. 31 compared to 0. 33 for other commercial birds. Further SNP genotyping results found Nici was homozygous at 293 of 333 SNPs (87. 99%) compared to an average of 275 (81. 73%) for birds from a Beltsville Small White flock closed for 30 years. Of note, all sequence data accumulated to date suggest that Nici is monomorphic at the MHC, typically the most polymorphic region of the genome [67]. It is noteworthy that NGS depth of coverage allowed for the use of a genome that was only partially inbred. Celera Assembler release 5. 3 was used to produce the assembly. The assembly process can be summarized to the following major stages: Stage 1 (gatekeeper): input of reads and quality control Stage 2 (overlapper): computation of read overlaps and trimming of poor quality sequence based on the overlaps Stage 3 (unitigger): initial assembly of uniquely-assemblable contiguous chunks of sequence based on the overlaps Stage 4 (cgw): scaffolding of unitigs based on mate pair data, followed by merging overlapping unitigs into contigs Stage 5 (consensus): computation of consensus sequences for the contigs There are multiple choices of the assembler modules available for overlapping and unitigging. The traditional OVL overlapper was originally designed for Sanger reads. The more advanced MER overlapper was designed to account for the homopolymer errors that are common in 454 read data. The MER overlapper is more accurate, although several times slower, on pure 454 assemblies. Surprisingly, with the combined Illumina and 454 Titanium data, the MER overlapper had no advantage over OVL, which suggests that the homopolymer errors are less pronounced in the latest Titanium data. Because BOG (best overlap graph) is more tolerant of highly variable read sizes (74 bp to 366 bp), the newer BOG unitigger was used instead of the original unitig module. Three maps were used to produce a Combined Map (CMap) for alignment of assembled sequence to chromosomes. The CMap had 31,769 markers that mapped both to the assembly and to the turkey chromosomes MGA1 through MGA30. Maps for the sex chromosomes W and Z were not used due to fragmentary marker coverage. Instead, scaffolds that aligned only to chicken W and Z chromosomes were identified and then ordered and oriented according to the chicken coordinates. A detailed comparative BAC contig physical map for turkey [16] was generated based on over 43,000 BES, over 80,000 BAC fingerprints, and over 34,600 BAC locations assigned by hybridization to overgo probes corresponding to 2,832 loci [70]. Two different BAC libraries were used: CHORI-260 generated by the Children' s Hospital of Oakland Research Institute and 78TKNMI generated at Texas A&M University. Comparative BAC contigs were assembled based on: (1) consistent (correct strandedness and separation distance) alignments of mate-paired BES to the chicken genome sequence (Build 2, May 2006, http: //genome. ucsc. edu), (2) hybridization to unique overgo sequence probes aligned with the chicken genome, and (3) BAC fingerprint-based contigs [16]. The BAC contig physical map, along with the BES, provides a tool for aligning scaffolds from the turkey sequence to turkey chromosome regions as well as for identification of rearrangements between the chicken and turkey genomes. (Regularly updated versions of this map are available at http: //poultry. mph. msu. edu/resources/resources. htm#TurkeyBACChicken, and it can also be accessed in graphical form at http: //birdbase. net/cgi-bin/gbrowse/turkey09/, see Text S1.) The current number of contigs, end sequence matches to the chicken genome and lengths are provided in Table S12. Most gaps between contigs are due to regions of low BAC density (particularly on microchromosomes MGA18 and 24–30 and on the sex chromosomes that are underrepresented in the BAC libraries and, in some cases, poorly assembled in the chicken sequence). However, some gaps are due to repetitive regions and likely sites of CNV [10]. The average size of comparative map BAC contigs on the autosomes is over 10. 5 Mb with the N50 average autosomal contig size being about 31 Mb. Twelve chromosomes (MGA6,9, 11,12,13,16–21, and 26) are spanned by only a single BAC contig and another five are spanned by two contigs (MGA2,5, 15,23, and 30). In addition to the previously known centric split of GGA2 to MGA3 and 6 and fusion of acrocentric MGA4 and 9 to create the metacentric GGA4, the comparative map suggests movement of more interstially positioned chicken centromeres to positions at or near the telomeres on MGA2,7, 10,11,12, and 13. Although MGA3 and 7 both contain short p arms visible in the turkey karyotype [10], no evidence of centromeric breaks internal to sequences orthologous to that of the chicken were found on these chromosomes, although there are a couple of short terminal contigs on MGA7 that could comprise a very small p arm. Of course, there may also be repetitive sequences or other sequences that were refractory to assembly in the chicken sequence that may be located on p arms of MGA3 and 7 (MGA8 and 14 are difficult to resolve near the likely p end telomere due to multiple rearrangements, and microchromosomes MGA18 and 25–30 tend to be fragmented in the BAC map and less well assembled in the chicken sequence due to poor BAC coverage). A total of 768 SNPs were genotyped on a F2 population of two genetically different commercial turkey lines that consisted of 18 full sib families with a total of 948 offspring. SNPs were genotyped using the Illumina Golden Gate assay. Of the 768 SNPs, 458 were eventually used to build linkage maps for 27 chromosomes (MGA1–17,19–26,28,30) (Table S13). The linkage map was constructed with a modified version of CRI-MAP software kindly provided by Drs. Liu and Grosz of Monsanto. All markers were checked for non-Mendelian inheritance errors using the option “prepare. ” Linkage maps for the individual chromosomes were constructed in a number of iterative rounds using the “build” option within CRI-MAP starting with a threshold of LOD = 5 with subsequent stepwise lowering the LOD threshold until LOD = 0. 1. Closely linked markers not separated by recombination events were ordered according to their location on the chicken sequence map (build WASHUC2). The order of markers in the final map was verified using the “flips” option. Strong SNVs are positions at which: (1) at least three reads support each nucleotide variant, (2) the sum of the top three quality values for each variant is at least 60, and (3) the overall depth of coverage is at most 30. For this analysis, gaps in the multiple alignments were assigned a quality value equal to the minimum quality of the flanking bases. In addition, if the SNV was an indel within a homopolymeric run, at least one Illumina read was required supporting each variant. The support and quality value thresholds should reduce chance sequencing errors to 1/1,000,000, and the 30-fold depth of coverage threshold should filter out apparent variations caused by near-identical repeats. The requirement for Illumina reads verifying homopolymer indels was used to filter out well-known 454 sequencing biases [71]. Weak SNVs are similar to strong SNVs but with relaxed thresholds. For weak SNVs, at least two reads had to support both variants, and the sum of the top two quality values for each variant had to be at least 45. The restriction on the depth of coverage was removed. As with strong SNVs, if the variant was an indel within a homopolymer, at least one Illumina read supporting each variant was required. Draft annotation was generated using two independent methods. First, a draft annotation of 12,206 putative protein coding loci was generated by combining evidence from multiple sources using JIGSAW [72]. Evidence for genes included spliced alignments of known proteins and mRNAs from multiple vertebrates, and expressed sequence tags (ESTs) from chicken and turkey. Ab initio gene predictions by Twinscan [73] and GlimmerHMM [74] (trained on chicken genes) were also considered by JIGSAW but with a very low relative weight, such that no gene models were based solely on ab initio predictions. The protein mappings were given the highest weight, followed by full-length mRNA alignments and then EST alignments. Proteins and mRNAs were taken from the most recent Ensembl gene builds for chicken, zebra finch, and green Anole lizard, and from the GenBank RefSeq database (the “other vertebrates” section plus mouse and zebrafish genes). Many thousands of genes and gene fragments were eliminated from the initial predictions if the computational evidence was insufficient. Second, the gene-finding pipeline at Ensembl [75], which also uses a combination of known proteins, ESTs, and cDNAs to annotate genes, was used to generate an independent set of protein-coding genes and noncoding RNAs. After combining the two gene lists, the total number of distinct protein coding loci was 15,093 plus 611 noncoding RNA genes, for a total of 15,704 genes. Some loci were identified as producing multiple distinct proteins due to alternative splicing, giving 16,217 distinct protein sequences. Orthologs between turkey, chicken, and human proteins were defined using sequence homology, phylogenetic trees, and conservation of synteny. Homologous pairs and orthology type are available from the version 57 Ensembl Compara database (http: //e57. ensembl. org). It was assumed that all 1∶1 orthologs were correct and were used to define conserved syntenic regions. Further orthologs were then defined from the one-to-many and many-to-many relationships, if the homologs mapped to a conserved syntenic region. This allowed for a 7%–8% increase in the number of defined orthologs for all species (Table S14). Multiple (three-way) alignments were built on the turkey, chicken [4], and zebra finch [7] genomes using Pecan [33]. Pecan is a global multiple sequence aligner that assumes no major rearrangements in the input sequences. Thus, sets of collinear segments were defined before aligning the sequences. Searches were based on the turkey-chicken and chicken-zebra finch pairwise BLASTZ-net alignments [76]. The genomes were repeat-masked using RepeatMasker (www. repeatmasker. org), followed by BLASTZ analysis, to find all highly similar regions, which were grouped in chains using the axtChain software and refined using the netChain software [77]. GERP [34] was used to get both per-base conservation scores and conserved elements. Searches were combined based on similarities and de novo repeat analysis approaches. For de novo repeat analysis, RepeatScout [80] and LTR_Struc [81] were used under default conditions. Overall, 944 repeat elements were identified, many of which were redundant. LTR_Struc uncovered three LTR retrotransposons that had both LTRs and the reverse transcriptase domain. During similarity searches, known chicken repeats (Repbase Update, http: //www. girinst. org) as well as representatives of different types of TEs were used as query. Repeats identified by RepeatScout were classified by comparison with known chicken repeats as well as with representative LTR retrotransposon protein sequences, non-LTR retrotransposon (or LINE) protein sequences, and DNA transposases. Whenever possible, the chicken homolog was used to assign names for the turkey TEs. Sequences other than those derived from the turkey genome project were collected from NCBI (www. ncbi. nlm. nih. gov) using BLAST [97], retrieved through searches in Ensembl (www. ensembl. org), and identified by BLAT searches on the UCSC Genome Bioinformatics database (genome. ucsc. edu). ESTScan v. 2 [98], [99] was used to correct predicted coding regions for frameshift and other sequencing errors. The alignments of amino acid and coding sequences used for the analysis of gene conversion and the construction of phylogenetic trees were generated with MUSCLE v. 3. 7 [100]. Gene trees reconciled with species trees were calculated using TreeBeST v. 1. 9. 2 [101] and trees were visualized using Archaeopteryx [102]. Perl scripts and modules from Bioperl [103] were used to manipulate sequence and phylogenetic data. DNA sequences for 42 high-coverage loci (11 mitochondrial) were collected from GenBank, for all Galliformes and for three outgroup Anseriformes species. One representative species of each genus was selected based on frequency of coverage, and additional representatives were chosen for the known polyphylous genus Francolinus and in cases of complementary locus coverage. The 86 final species were represented by 10. 6 kb on average. For Coturnix coturnix, loci from the complete mitochondrial sequence of C. japonica (often considered a subspecies of the former) were used. Northura maculosa was classified as a cracid at NCBI but did not join the other cracids in preliminary trees; this GenBank entry was apparently a misspelling and misclassification of the tinamou N. maculosa and was removed from the study; GenBank was notified of the problem. For each locus, sequences were aligned and the unmasked alignments were concatenated, partitioned by gene (except that the two mitochondrial rRNA genes were fused as were the eight mitochondrial coding sequences). A maximum likelihood tree was constructed using RA×ML with the GTR-GAMMA model and 100 full bootstraps were taken.
In contrast to the compact sequence of viruses and bacteria, determining the complete genome sequence of complex vertebrate genomes can be a daunting task. With the advent of “next-generation” sequencing platforms, it is now possible to rapidly sequence and assemble a vertebrate genome, especially for species for which genomic resources—genetic maps and markers—are currently available. We used a combination of two next-generation sequencing platforms, Roche 454 and Illumina GAII, and unique assembly tools to sequence the genome of the agriculturally important turkey, Meleagris gallopavo. Our draft assembly comprises approximately 1. 1 gigabases of which 917 megabytes are assigned to specific chromosomes. Comparisons of the turkey genome sequence with those of the chicken, Gallus gallus, and the zebra finch, Taeniopygia guttata, provide insights into the evolution of the avian lineage. This genome sequence will facilitate discovery of agriculturally important genetic variants.
genetics and genomics/comparative genomics
genetics and genomics/genetics of the immune system
computational biology/comparative sequence analysis
evolutionary biology/genomics
genetics and genomics/genome projects
computational biology/genomics
genetics and genomics/bioinformatics
Multi-Platform Next-Generation Sequencing of the Domestic Turkey (Meleagris gallopavo): Genome Assembly and Analysis