Sample collection and DNA isolation
Fecal, blood and tissue samples were collected from plains zebra (E. quagga/E. burchelli) in two parks of southern Africa: Etosha National Park, Namibia (n = 38) and Kruger National Park, South Africa (n = 33). For the purposes of consistency with historical ELA allele nomenclature, we hereafter refer to the species by its former scientific name, E. burchelli. With fecal samples, three to five pellets were collected from each individual and allowed to dry. Epithelial cells from the outermost mucosal layer were scraped from the desiccated pellets using a sterile razor blade. Tissue samples were preserved in DMSO/EDTA/Tris/salt solution and blood samples in ethylenediaminetetraacetic acid (EDTA). All samples were stored at-20°C until DNA extraction. Sample collection was approved by the Animal Care and Use Committee (Protocol #R217-0510B) at UC Berkeley.
Whole genomic DNA was extracted from blood and tissue using Qiagen kits (Valencia, CA). Non-invasive samples, collected from feces, are subject to contamination, enzyme degradation (e.g. ), and hydrolytic and oxidative damage that may result in lower DNA yield and increased error rates (most commonly allele dropout ). Thus, we used the AquaGenomics protocol (MultiTarget Pharmaceuticals, Inc.) optimized for fecal DNA extraction. A few fecal samples suffered degradation which resulted in failed PCR-amplifications. These degraded samples were re-extracted using the QIAmp fecal extraction kit (Qiagen), also designed specifically for fecal DNA extraction.
PCR-amplification and sequencing
We targeted two MHC loci of the Equine Lymphocyte Antigen (ELA) system, ELA-DRA and DQA, by polymerase chain reaction (PCR)  and genotyped these loci through direct sequence-based typing. We amplified 246 bp of the DRA using equid-specific primers, Be3 and Be4 , and 205 bp of the DQA using the primers DQA-2e and DQA-2f . These primers targeted the functionally significant exon 2 of both genes, a region consisting of antigen binding sites (ABS) as predicted by their human lymphocyte antigen (HLA) equivalent . PCR mixes (total reaction volume of 15 μL) for both genes contained approximately 25-50 ng DNA, 2 uL GeneAmp 10 × PCR buffer (100 mM Tris-Cl, pH 8.3, 500 mM KCl, 15 mM MgCl2, 0.01% (w/v) gelatin), 1 U AmpliTaq Gold DNA polymerase (Applied Biosystems), 0.4 mM dNTPs, 15 μg bovine serum albumin (New England BioLabs) and 0.50 μM of each primer.
Amplification of the DRA locus used the following "touch-down" thermocycling profile: an initial denaturation at 95°C for 10 min; 2 cycles of 94°C for 1 min, 60°C for 1 min, and 70°C for 35 s; 18 cycles of 93°C for 45 s, 59°C for 45 s, and 70°C for 45 s, with the annealing temperature decreasing by 0.5°C with each cycle; 35 cycles of 92°C for 30 s, 50°C for 30 s, and 70°C for 1 min; final extension at 72°C for 10 min to allow for complete amplification of the targeted gene. PCR-amplification of the DQA locus used the following thermocycling profile: an initial denaturation at 95°C for 6 min; 40 cycles of 94°C for 45 s, 56°C for 45 s, and 72°C for 1 min; final extension at 72°C for 5 min.
DRA amplicons were purified prior to sequencing by incubating with Exonuclease I and Shrimp Alkaline Phosphatase at 37°C for 30 minutes. Purified products were cycle-sequenced in both forward and reverse directions using the Big Dye® Terminator v.3.1 kit and run on an ABI 3730 automated sequencer (Applied Biosystems).
Identification of MHC alleles
Sequence chromatograms were edited and aligned using the software Geneious 4.7 . Allelic phase for DRA heterozygous sequences was determined by computational inference with the haplotype reconstruction program PHASE v.2.1 . This program has been found to be accurate in determining allelic phase even in extremely variable loci, such as the MHC  and, therefore, is considered to be a reliable method for allele identification. We conducted five runs, using different initial random seed values, and compared phase results across runs. A threshold posterior probability of 0.9, a value considered significantly higher than the standard (see ), was used to assess the accuracy of the allelic phase determination. Individuals not meeting this threshold were dropped from use in further analyses.
Given the large number of heterozygous sites in the DQA locus and previous evidence for multiple loci , all PCR-amplicons were cloned and sequenced to identify novel haplotypes. PCR products were extracted and purified with the QIAquick Gel Extraction Kit, (Valencia, CA) and cloning was performed using a TOPO-TA® cloning kit with Mach 1™-T1R competent cells (Invitrogen). Amplicons were ligated into pCR®4 TOPO vectors and transformed into E. coli competent cells. Sixteen to twenty-three positive clones per individual were picked with a sterile toothpick and screened by sequencing (protocol described above). The high number of PCR-amplified clones was sufficient to avoid errors, such as recombinant sequences generated during PCR . Each allele was confirmed with at least two observations, meaning that it had to be found in at least one homozygous individual or two heterozygous individuals to be included in the following analyses.
Sequence data and alignments
Novel MHC alleles identified in E. burchelli were compiled with a reference panel of Equidae sequences (GenBank, NCBI), including horse (E. callabus), ass (E. asinus), onager (E. hemionus), kiang (E. kiang), plains zebra (E. burchelli), mountain zebra (E. zebra), Grevy's zebra (E. grevyi) and Przewalski's horse (E. przewalski). A list of ELA-DRA and DQA sequences from each equid species and their respective GenBank accession numbers are listed in Additional file 1. As the ELA-MHC nomenclature is currently in revision, names for previously discovered alleles follow designations given in Janova et al. (2009) and novel sequences discovered here were named based on the recommendations outlined by the MHC allele nomenclature committee . The new nomenclature is expected to be established soon on the IPD-MHC Database (http://www.ebi.ac.uk/ipd/mhc). Identical alleles shared between species were given species-specific numbering. Reference and novel nucleotide, and corresponding amino acid sequences were aligned using the Geneious 4.7 sequence alignment tool and editor .
Statistical analyses of diversity and evolution
Standard descriptive diversity indices for each locus within the genus Equidae were calculated using MEGA4 . These indices included the number of alleles (A), variable nucleotide positions (VNP), parsimony informative positions (PIP), transition/transversion bias ratio (R), Kimura 2-parameter gamma (K2P+Γ) evolutionary distance (d) and Poisson-corrected amino acid distance. The K2P+Γ model accounts for multiple hits, differences in transitional and transversional substitution rates and variation in substitution rates among sites following a gamma-shaped distribution. Estimates of the gamma shape parameter (α) were determined in PAUP*v4.0b0  to be α = 0.9872 for the DRA data and α = 0.4181 for the DQA data. Standard error of distance estimates were obtained by using a bootstrap procedure with 10,000 pseudoreplicates.
Four different methods, implemented in RDP v.3.44 beta package , were used to test for recombination and detect potential recombinant events: (1) RDP, (2) GENECONV, (3) Maximum Chi, and (4) BootScan. The highest acceptable p-value for all methods was set at a conservative value of 0.10, with a Bonferroni correction for multiple comparisons and a window size of 30 variable nucleotides for all approaches except BootScan. For analyses in BootScan, 1,000 bootstrap replicates were conducted under the Kimura model (transition/transversion ratio = 1.341), with a window size of 100 bp, step size of 20 nucleotides and cut-off value of 0.70.
Selection, averaged across the gene, was estimated using MEGA4  in terms of the relative rates of nonsynonymous (d
N) and synonymous (d
S) base pair substitutions, according to Nei and Gojobori (1986) with the Jukes and Cantor correction for multiple hits . Z-tests of selection were performed over all sites, and separately at ABS and non-ABS, under the null hypothesis of neutrality (d
N = d
S) and the alternative hypotheses of non-neutrality (d
N ≠ d
S), positive selection (d
N > d
S), and purifying selection (d
N < d
Site-specific selection analyses
As selection will realistically act on only a small subset of amino acids in a protein, averaging substitution rates over entire gene regions is considered to be a conservative indicator of positive selection . Therefore, we used a more powerful maximum-likelihood based method, implemented in the CodeML subroutine of the software PAML  which allows the rates of ω = d
S to vary among codons [31, 64]. This method has been suggested to be more sensitive than other methods for detection of molecular evidence of selection . The models employed here, called 'random-sites' models, do not require a priori information on the functional significance of each site and estimate the nonsynonymous to synonymous rate ratio (ω) to indicate selective pressure at the protein level (ω < 1: purifying selection, ω = 1: neutral evolution, ω > 1: positive selection). In this analysis, we used the Equidae alignments to assess heterogeneity in ω across the two MHC genes (DRA and DQA) and to identify codons under positive selection. We fit the alignment to the following codon 'random-sites' models, in PAML: M0 (one ratio: best average ω across all sites), M1a (nearly neutral: estimates the proportion of sites that best-fit ω = 0 versus those best-fit by ω = 1), M2a (positive selection: adds a third set of sites to M1a that have ω > 1 and estimates the best-fit for this added ω value and associated proportion of sites), M3 (discrete: fits proportions and ω values assuming three classes of sites labeled 0, 1, and 2 such that ω
0 < ω
1 ≤ ω
2), M7 (beta: ω is beta-distributed on [0, ]) and M8 (beta and omega: a proportion of sites are beta-distributed on [0, ] and the remaining proportion have an average ω
2 > 1 ). M0 is the only model that does not allow for variation in ω across codon sites. Whereas M1a and M7 allow only for neutral evolution and purifying selection at some proportion of sites, M2a, M3, and M8 also allow for the possibility of positive selection at a proportion of sites.
Likelihood ratio tests (LRT) were used to compare nested models based on their log-likelihood . We compared M0 and M3 to test for the significance of heterogeneity in ω across sites, whereas M1a was compared with M2a, and M7 with M8 to test for positive selection. Significant adaptive evolution was inferred if twice the difference in log-likelihood values was greater than the chi-square critical value for the given degrees of freedom. We used the Bayes empirical Bayes (BEB) approach  to estimate mean ω and standard errors across codon positions. Specific sites under positive selection were indicated by estimates of ω > 1 and posterior probabilities > 0.95. This approach accounts for sampling errors in the maximum likelihood estimates of the parameters and has a low false positive rate. Tree files used in PAML analyses were generated using a maximum likelihood approach in PhyML , under the Kimura 3-parameter and the Kimura 2-paramter model of nucleotide substitution for the DRA and DQA locus, respectively. Models of nucleotide substitution and the distribution of rate variation across nucleotide sites (gamma) were estimated in PAUP*v4.0b0 .
Phylogenetic relationships among Equidae DRA and DQA sequences were reconstructed using a Bayesian approach implemented in MrBayes 3.1 . The data set was partitioned and the best-fit models were determined for each codon position using the Akaike Information Criterion (AIC) in MODELTEST v.3.7 . Bayesian inference involved running six Metropolis-coupled MCMC chains (1 cold and 5 heated) simultaneously at n incremental temperature of 0.1, and chains were run for seven and sixteen million generations for the DRA and DQA data, respectively. Trees were sampled every 100 generations and the first 25% of trees found were discarded, leaving the remaining trees to be used for estimating the consensus tree. Two independent analyses were conducted and results were compared to check for convergence by confirming that the average deviation of split frequencies approached 0 (with values less than 0.01). We also checked that the potential scale reduction factor (PSRF) approached 1 and that chains mixed sufficiently (with chain mixing values greater than 0.2 between chain pairs). Finally, we used the program Tracer v1.4  to ensure whether sampling from the posterior distribution of each parameter was sufficient and had reached a large enough effective sample size (ESS > 200) for accurate parameter estimation. Posterior probabilities, representing the probability that a specific node is observed, were recorded. This analysis was run on both non-partitioned and partitioned data, and the optimal model was determined using Bayes Factors.
DRA sequences from Bos taurus (DQ821713), Ovis aries (Z11600) and Sus scrofa (AY754888) obtained from GenBank (NCBI) were used as outgroups. For DQA trees, available sequences from B. taurus (AB548942), O. aries (M33304) and S. scrofa (EU195146) were used as outgroups.