Resolving the evolution of extant and extinct ruminants with high-throughput phylogenomics

We have genotyped 16,353 animals representing 61 cattle breeds and 70 species, as divergent from Bos taurus as the Savannah elephant (Table S1), with the Illumina BovineSNP50 BeadChip (11, 12) according to Illumina protocols (13). To examine the quality of genotype calls in these outgroup species, we first sequenced the SNP site and flanking regions for rs17871403 in 14 species, with pronghorn the most divergent of the sequenced species (Table S2). This SNP was chosen because it has been well characterized in cattle and is a member of a SNP panel that is widely used for parentage analysis (14). Of the genotypes produced by the BovineSNP50 assay (Illumina) for this SNP in these species, 99.13% were concordant with the sequence when we allowed for genotype ambiguity (i.e., WW and SS) (see Methods). One of the six genotyped North American mountain goats and one of the eight genotyped caribou had discordant BovineSNP50 and sequence-based genotype calls (Table S2). This analysis of a single SNP across multiple species suggests a genotyping error rate for BovineSNP50 loci of only 0.87%.

We next aligned all 40,843 SNP probe sequences, which are 50 bases in length, to the international sheep genomics consortium (www.sheephapmap.org) genome assembly (available at https://isgcdata.agresearch.co.nz/ and in an annotated form at http://www.livestockgenomics.csiro.au/sheep/oar1.0.php) and found that only 26,098 (63.9%) could be uniquely aligned, primarily due to the incomplete status of the assembly. Of these SNP, 829 had an unknown base (N) identified at the position of the SNP, and for the remaining 25,269 SNPs, there were 308,518 genotypes called in 17 sheep. Genotype calls were in agreement with the genotype predicted from the respective sequence base for 298,311 genotypes (96.7%). There were 1,834 heterozygous genotypes and 8,373 genotypes that were homozygous for an allele not predicted by the sequence assembly. This suggests a BovineSNP50 genotyping error rate of between 2.7 and 3.3% in the outgroup species.

Finally, when minor allele frequencies (MAF) averaged over 40,843 SNPs were plotted against average genotype call rates, samples from outgroup species with the lowest call rates had higher than expected MAF (Fig. S1). This appears to be indicative of DNA quality issues since, for example, DNA for the Capra ibex samples was extracted from irradiated blood samples that had been stored under refrigeration for several years. On removing these samples, there was almost no correlation between MAF and call rate (Fig. S1). This indicates that as genetic distance from cattle increases and call rate decreases, spurious heterozygote and alternate homozygote genotype calls rarely arise, indicating support for the quality of these data.

Using genotypes for 40,843 SNPs scored with the BovineSNP50 BeadChip (see Methods), we produced a completely bifurcating tree with highly supported nodes for 61 Pecoran species, that contains species that diverged up to 29 million years ago (Fig. 1) (15). There were 39,695 parsimony-informative characters using all 678 animals and, remarkably, 21,019 with cattle excluded. Within the Bovidae, only nine nodes had support <100%. We propose 17 relationships and increase the support for 16 previously proposed nodes within the infraorder, when compared to the supermatrix phylogeny of Marcot (3). A striking observation from the phylogeny is that taxonomic classifications of families and subfamilies mirror the topology of the cladogram, since higher taxa form monophyletic groups. This is an improvement over earlier phylogenies, as previously questionable groupings are now shown to be monophyletic.

Currently, PCR-based and non-PCR-based multiple strand displacement amplification (MDA) approaches are used to perform whole genome amplification (16, 17). MDA requires high-quality DNA over 2 kb in length and was found to be inefficient for the ancient bison DNA. Consequently, we used a universal linker-based PCR amplification performed with the GenomePlex Whole Genome Amplification kit (Sigma-Aldrich) to amplify the minute amounts of damaged DNA preserved in bone samples from two ancient Russian Bison priscus specimens and test whether the Illumina iSelect platform could be used to analyze samples derived from extinct species. The first, sample BS662, was collected from permafrost deposits at Alyoshkina Zaimka, Siberia, and is approximately 20,000 years old (18). The second, ACAD012, was collected from Sur'ya 5 cave in the Ural Mountains and has been accelerator mass spectrometry radiocarbon dated to 34,460 ± 290 years BP. Due to the low amounts of DNA from the ancient specimens and the short DNA fragment lengths produced in the whole genome amplification of degraded ancient samples, the genotype call rates for these samples were much lower than for modern bison (Table S1). However, when these ancient samples were included in the Bovini phylogeny (Fig. 1), BS662 was basal to the modern Bison bison clade as expected, but ACAD012 fell within the modern Hereford cattle clade. When we sequenced several overlapping fragments that had been individually amplified from the hypervariable mitochondrial control region of sample ACAD012, we identified variability within the overlapping regions. This is consistent with the sample having been contaminated with modern DNA or being extremely degraded, as also suggested by our genotype data and consequently the sample was removed from the study. A replicate whole genome amplification (library identification KCMU02) was produced from the B. priscus sample used to generate BS662, and when this sample was included in the data set, it was sister to BS662, and both remained sister to modern bison within the phylogeny. However, in the preparation of this library, we avoided the initial DNA fragmentation step within the amplification protocol that appeared to greatly improve the quality and quantity of produced genotypes, as KCMU02 produced a higher genotype call rate (54.9 vs. 45.8%) and far lower heterozygosity (11.5 vs. 39.6%) than did BS662 (Table S3). While only 76.1% of the 12,279 genotypes that were called in both samples were identical, 99.7% of the homozygous genotypes, the only genotype class that has the potential to be phylogenetically informative (see Methods), were identical between the replicates.

Phylogenetic relationships were also inferred for 48 cattle breeds (n = 372 animals) (Table S1) using parsimony, with most nodes being highly supported (bootstrap values >70%). To accommodate heterozygotes, data were first coded with heterozygotes as polymorphic (noninformative) and then as an independent character state (see Methods). When coded as polymorphic, the topology of the cladogram corresponded to the known geographic origins of breeds (Fig. 2A). Interestingly, however, when heterozygotes were coded as distinct characters, the topology changed and no longer clearly reflected the biogeography of breed origins (Fig. 2B).

To further resolve the issue of breed origins, we constructed phylogenetic networks which can reveal conflicting signals in the data (Fig. 3 and Fig. S2). In Fig. S2, Bos taurus indicus and Bos taurus taurus are distinct groups with long edges between the subspecies. Within B. t. taurus, using the Reynolds et al. (19) distance metric and parsimony cladograms (Fig. 2), African taurine cattle were inferred to be more divergent from European cattle than are the Asian B. t. taurus breeds, with 100% bootstrap support in cladograms (Fig. 2 and Figs. S2 and S3). Because SNP were almost exclusively discovered from European B. t. taurus samples (12), there is a strong ascertainment bias toward SNP common within European B. t. taurus on the BovineSNP50 BeadChip, leading to severe biases in estimates of genetic distance that have prevented us from accurately dating the nodes separating European, African, and Asian cattle (Figs. S3 and S4). Furthermore, the data were recalcitrant to correction for ascertainment (see Methods). The network with individuals at node tips (Fig. 3) appears to accurately depict the admixed nature of many populations, for example, the relationship of Belgian Blue to Holsteins and Shorthorns, and Jersey to Iberian and British breeds. The network also reveals pedigree relationships, with sire HO020740 being an interior node to son HO020879.

Ancient DNA was extracted from fossil bison bone specimens using the standard phenol/chloroform/Amicon Ultra-4 method (17). DNA extractions, omniplex library preparations, and PCRs were set-up and performed in a geographically isolated, dedicated ancient DNA facility at the University of Adelaide, Australia. To generate a library of genomic fragments from limited ancient DNA extract, DNA was amplified using the PCR-based GenomePlex Whole Genome Amplification kit (WGA2; Sigma-Aldrich) according to the following protocol: 10 μL DNA were thoroughly mixed with 2 μL library preparation buffer and 1 μL library stabilization solution, and denatured at 95 °C for 2 min. After denaturation, 1 μL library preparation enzyme was added to generate omniplex libraries, followed by a series of incubations at 16 °C for 20 min, 24 °C for 20 min, 37 °C for 20 min, and 75 °C for 5 min in a thermal cycler (Corbett Life Science). The omniplex libraries were next amplified using a limited number of genomic amplification cycles. PCR amplification was conducted in a 75-μL reaction volume containing 14 μL omniplex library, 7.5 μL amplification master mix, 48.5 μL nuclease-free water, and 5 μL WGA DNA polymerase. The PCR amplification conditions were initial denaturation at 95 °C for 3 min, followed by 15 cycles of 94 °C for 15 s and 65 °C for 5 min. GenomePlex-amplified ancient DNA products were finally purified using the GenElute PCR Clean-Up kit (Sigma-Aldrich). Ancient DNA libraries were verified by PCR amplification and sequencing of the hypervariable mtDNA control region before analysis with the BovineSNP50 BeadChip (Illumina). A second amplification, labeled KCMU02, of the sample that produced BS662 was constructed using the same protocol as above, except the genomic fragmentation step within the WGA2 protocol was omitted.

Table S1 shows the numbers of animals genotyped from each species or cattle breed. In taxa or breeds where <10 animals were genotyped, all animals were sampled. If >10 animals were genotyped, animals with the highest genotype call rates and earliest birth dates were selected. When pedigree information was available, closely related animals were avoided, except in Angus and Holstein where 10 old animals (born in the 1950s, 1960s, and 1970s) and 10 recently born animals (born in the late 1990s and 2000s) were selected. When more than 50 animals within a breed had call rates of at least 98% and no pedigree information was available, 10 animals were sampled at random. Samples belonging to recently formed crossbred breeds were removed from the analysis, as these samples distort parsimony phylogenies. Genotypes for the two ancient Bison samples were included despite their much lower genotype call rates, which were expected due to DNA degradation and fragmentation, and the use of whole genome amplification, which affect the fidelity of the Infinium assay. The provenance of all samples included in the analyses is provided in Table S4.

The BovineSNP50 BeadChip (Illumina) consists of SNP primarily discovered by the sequencing of reduced representation libraries (11), the alignment of random shotgun reads from six cattle breeds to the Hereford assembly, or from the draft assembly of the bovine genome (12). To improve genotype quality for B. t. indicus and the outgroup species, we manually adjusted genotype call clusters in Illumina BeadStudio to improve genotype calls. Where pedigree information was available, such as in O. aries and B. bison, the rate of misinheritances was minimized. A set of 40,843 SNP was selected from the 54,693 loci queried by the assay. Loci selected for analysis were all located on autosomes, had a call rate of at least 80% in 36 (75%) B. t. taurus breeds, and were not monomorphic in all breeds. This strategy was effective in selecting informative SNP with fnew genotype errors (Table S5). Data are available at http://animalsciences.missouri.edu/animalgenomics/publications/php.

Almost 96% of the beads on the BovineSNP50 BeadChip query Infinium II SNP, in which adenine and thymine share a fluorescent probe and guanine and cytosine share a different fluorescent probe. For samples in which all four bases are present at a single locus, AA, AT, and TT genotypes produce indistinguishable fluorescence intensities, as do GG, GC, and CC. Thus, A/T or C/G SNP discovered in B. t. taurus were limited in the assay design (1.8 and 2.2%, respectively, and use Infinium I chemistry). However, in species diverged from B. t. taurus where all four bases could be present, genotypes are WW (W is the IUPAC code for A or T bases) for one homozygote class, SS (S is the IUPAC code for G or C bases) for the alternate homozygote, and NN (ambiguous) for the heterozygote class. This ambiguity is evident when sequences and genotypes for outgroup species were compared (Table S2). The WW and SS genotypes were identified in BeadStudio as AA and BB genotype calls.

Most parsimonious trees were inferred from the genotypes using TNT version 1.1 (29). In the analyses involving the outgroup species, phylogenetic signal was obtained only from the homozygous genotypes, and AA homozygotes were coded as “0,” BB homozygotes were coded as “1,” heterozygotes were coded as a polymorphic character state (i.e., “[0,1]”), and missing genotypes were coded as “?.” However, in the analyses of the cattle breeds, an additional data set was created in which heterozygotes were identified by a unique character state (i.e., AA = 0, AB = 1, BB = 2). A heuristic search was conducting using the search technology in TNT, and the search level was initially set to 20. Specifically, we used the SPR-TBR algorithm followed by random sectorial searches, constrained sectorial searches, exclusive sectorial searches, and 10 rounds of tree-drifting. The complete search was replicated 20 times, with 10 rounds of tree fusing at the conclusion of these 20 replicates. A subset of the samples from the tribe Bovini was independently analyzed along with the ancient bison samples to validate the quality of the data generated from these ancient samples. A data set with 714 samples from all taxon groups was first used to construct the most parsimonious trees. After excluding samples with low quality DNA, low bootstrap support, and/or nonsensical placement in the cladogram (i.e., elephant and horse as sister to B. taurus), a final data set with 678 samples was used to construct most parsimonious trees. The cladogram was rooted with Antilocapra americana. Using these 678 samples, bootstrap support was calculated using 1,000 pseudoreplicates, and for expediency, the SPR-TBR heuristic search was used.

Allele frequencies were estimated for 40,843 SNP in 22 breeds (Table S6), and these frequencies were used to estimate pairwise Reynolds distances (19) among the breeds (Fig. S3). Several attempts were made to correct estimates of genetic distance for SNP ascertainment bias. First, distances were calculated from haplotype frequencies. Haplotypes were inferred for the autosomes of all genotyped animals in our collection within each breed group (Table S6) using fastPhase (30). From these haplotyped samples, haplotypes were extracted for the study animals for 885 nonoverlapping loci, each comprising six SNP for which the intermarker distance was <50 kb for contiguous SNP. Haplotype frequencies were estimated for each of the 885 loci within each breed group and were used to estimate Reynolds distances between breeds. Next, we formed weighted distances by averaging individual SNP distances weighted according to the frequency of unascertained SNP (31) possessing the MAF observed in each of the two populations. Finally, we also subsampled approximately 3,000 or approximately 8,000 SNP such that the resulting MAF distribution conformed to the unascertained distribution of bovine SNP (31) in Angus or Holstein, respectively. The subsample size was determined by the severity of underrepresentation of SNP within the MAF range 0.005–0.015 and indicates that ascertainment bias was more severe for Angus than for Holstein. Reynolds and Nei genetic distances corrected for sample size (Table S6) were estimated for each subsample and were averaged across 1,000 bootstrap replicates. Distances were used to construct neighbor-joining and UPGMA trees with Phylip (32). None of the approaches taken to correct for ascertainment bias were able to establish a tree in which branch lengths were clock-like. Biases in the allele frequency spectrum differ within B. t. taurus breeds (Fig. S4) causing the distances between breeds to not be clock-like.

Figures of phylogenies and cladograms were produced in MrEnt3 (33), and phylogenetic networks were constructed using SplitsTree version 4.10 (34). Distances based upon allele frequencies at 40,843 SNP were used to construct a network of 22 breeds. Due to memory limitations in SplitsTree, genotypes at 14,023 SNP were used to construct a network of 372 individuals belonging to 48 breeds. Default settings in SplitsTree were used to construct the networks.