Methods

Whole Genome Amplification of Ancient DNA.

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.

Sample Selection.

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.

SNP Selection.

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.

Genotype Calls in Outgroup Species.

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.

Phylogenetic Analysis.

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.