Skin biopsies were taken by remote system delivery of biopsy darts. We used a CO2 powered DanInject (Børkop, Denmark) Model JM with both in house and manufactured (Palmer Capshur, Atlanta, GA, USA) 6 mm biopsy darts. Where observed, we attempted to sample distinct groups within each subspecies. The presence of so many mtDNA haplotypes within many subspecies suggests multiple matralines were sampled (see results). Subspecies assignments for each sampled giraffe were based on geographic location and pelage following Dagg and Foster . The skin samples were placed in 0.5 ml room temperature tissue preservative buffer for preservation. Samples were transferred to the same buffer but with 0.2% gluteraldehyde for sterilization before export/import to the USA. The sampling performed by HDZ researchers was performed under Kenyan permit KE911780-1, Ugandan permits UWA/PMR/RES/50 and Ugandan National Council for Science and Technology permit #EC549, Niger Interior Ministry Permit 731 and Namibian Ministry of Environment and Tourism Research/Collection Permit #597/2002. All samples were imported under USDA/APHIS Import Permit #43686. Detailed permit information is available on request to the authors. We extracted genomic DNA from giraffe biopsy samples using a standard phenol chloroform/isoamyl alcohol extraction protocol.
A 654-nucleotide fragment of mtDNA was amplified and sequenced in 266 giraffes and one okapi (Okapia johnstoni). We amplified and sequenced this fragment using the primers L15774 and H16498 . Polymerase chain reaction amplification was performed in a 50 μl reaction using an MWG-Biotech Primus 96 Plus thermal cycler with 35.7 μl sterile double-distilled water, 5 μl 10 × PCR buffer, 5 μl of 25 mM MgCl2, 1 μl of 10 mM dNTP mix, 1 μl of both 25 pM/μl forward and reverse primers, 0.3 μlTaq polymerase (Sigma-Aldrich, St Louis, MO, USA), and approximately 50 ng of genomic DNA. The PCR amplification profile was 94°C for 3 min, followed by 30 cycles of 94°C for 30 s, a primer-specific annealing temperature for 35 s, 72°C for 45 s, ending with a single extension of 72°C for 5 min. All PCRs included a negative control (no DNA). PCR products of expected size were excised from 1% agarose/Tris/acetic acid/EDTA gels and purified using an Ultra Clean Kit (MoBio Laboratories, Solana Beach, CA, USA). The mitochondrial fragment was sequenced in both forward and reverse directions on a Beckman CEQXL2000 capillary sequencer (Beckman Coulter, Fullerton, CA, USA). Sequences were aligned using Sequencher 3.0 (Gene Codes Corp., Ann Arbor, MI, USA).
The mtDNA data matrix (n = 266 sequences) was collapsed to 35 haplotypes using the program Collapse v1.1 . To ensure proper phylogenetic resolution and robust support for relationships among haplotypes, the rest of theCYTb gene was amplified and sequenced from 35 giraffes, representing each of the 35 unique haplotypes, and the okapi. Primers L14724 , L15162, and H15915  were used to amplify and sequence theCYTb gene using the same protocols described above. This sequence was then concatenated with the 654 nt fragment, resulting in an alignment length of 1709 nt (with okapi) or 1707 nt (without okapi). These sequences were deposited in GenBank (accession numbers EU088317–EU88352).
Phylogenetic relationships among the 35 giraffe haplotypes (1709 nt) were estimated using maximum parsimony (equally weighted) (MP), maximum likelihood (ML), and minimum evolution (ME) methods. The HKY85+I+G model of DNA substitution was selected  and used in ML and ME analyses that included only giraffe haplotypes. For analyses that included the okapi, the HKY85 (without accounting for site heterogeneity) model was used. Gaps (insertions/deletions) were coded as a fifth base in MP analyses. Maximum parsimony and ME analyses were executed in PAUP* 4.0b10 . For these analyses, heuristic searches were performed using 100 random sequence additions, with one tree held at each step during stepwise addition, tree-bisection-reconnection branch swapping, steepest descent option not in effect, no upper bound for MaxTrees, and MulTrees option in effect. Maximum likelihood analyses were conducted with TREEFINDER  and parameters of the HKY85+I+G or HKY85 model were estimated along with the tree topology. For each phylogenetic method, robustness of clades was assessed using 1000 bootstrap pseudoreplicates. The okapi sequence was used to root the phylogenies of the giraffe haplotypes. However, due to the large sequence divergence between giraffe and okapi (and therefore the potential for signal saturation), phylogenetic trees including only giraffe haplotypes were also midpoint rooted. Regardless of rooting method or optimality criterion used, clades with a ≥ 80% bootstrap value were maintained across all analyses (Figure 1A and Additional files 2, 3, 4, 5, 6, 7).
Minimum-spanning network between haplotypes
We also constructed a minimum-spanning network of absolute distances between control region haplotypes using the molecular-variance parsimony algorithm as implemented in Arlequin v3.1 .
Genetic structure and diversity
Population structure was deduced with an analysis of molecular variance (AMOVA) using Arlequin V3.1 . In order to identify groups of populations based on genetic differences, we grouped sampling localities to maximise the among-group variance component (Φct). Haplotype and nucleotide diversity indices were calculated with Arlequin V3.1 using mtDNA control region data.
Estimation of divergence times using MDIV 
We generated maximum likelihood estimates of θ, twice the effective female population size (N
) times the mutation rate (u);T, the divergence time between two populations scaled by population size; andM, the gene migration rate between the two populations, also scaled by population size. We assumed uniform prior distributions and applied an HKY model of mutation  to allow for the possibility of multiple mutations per site. We ran Markov chains of 4000000 cycles preceded by a "burn-in" period of 500000 cycles for each pairwise population comparison, set maximum values forT andM of 10 and 30, respectively, and ran the analysis three times for each population comparison using different random seeds. We calculated divergence time (t) using the formulat =T *θ/(2u)*g, where T and θ are generated by the program,u is the mutation rate, andg is generation time in years. We calculatedu as 2*μ*k, whereμ is the mutation rate per nucleotide andk is the length of the sequence. Given the higher mutation rate found in the control region relative to the cytochromeb gene, we used a range of estimated mutation rates for control region sequence that span values found previously in other large mammal species. These included 0.05, 0.10 and 0.15 substitutions/site/MY, and a generation time of 4 years, which is the approximate age of first breeding for giraffes .
We amplified 13 published  and one novel (see Additional file 18) giraffe-specific microsatellite loci to generate multilocus genotypes for the 381 individuals. We performed the PCR amplification in 25 μl reaction volumes using an ABI 480 thermocycler (Perkin-Elmer; Foster City, CA, USA) with approximately 50 ng of genomic DNA as template. Final amplification conditions consisted of 12.5 pmol unlabelled reverse primer, 12.5 pmol fluorescently labeled forward primer, 1.5 mM MgCl2, 200 μM each dNTP, and 0.5 units ofTaq DNA polymerase (Promega; Madison, WI, USA). The thermal profile for PCR amplification was 95°C for 5 min, followed by 35 cycles of denaturing at 95°C for 30 s, a annealing at primer-specific temperature for 30 s, and elongated at 72°C for 30 s, ending with a single extension of 72°C for 10 min. We separated the PCR products on either a 7% polyacrylamide gel electrophoresed on an ABI 377 or through POP4 capillary buffer electrophoresed on an ABI 3100 DNA Analyzer (Applied Biosystems, Inc; Foster City, CA, USA). We assigned allele fragment lengths relative to the GeneScan-500 (TAMRA; Applied Biosystems, Inc; Foster City, CA, USA) size standard using the ABI GeneScan software program. We checked and corrected the data set for errors using MICRO-CHECKER 2.2.3  and MSA 4.00 .
Allelic diversity, and expected (He) and observed (Ho) heterozygosity were calculated using the program GENALEX . Each locus was tested for deviation from Hardy-Weinberg equilibrium and linkage disequilibrium with other loci (p < 0.05) using the program Genepop . Bonferroni corrections to significance values  were applied to account for multiple tests (see Additional files 12A and 12B).
Allele-sharing neighbour-joining network
We generated the neighbor-joining network tree using 14-locus genotypes of 381 individuals. The network was created using the allele-sharing distance Ds  and the program POPULATIONS v1.2.28 .
Bayesian clustering analysis
We used the program STRUCTURE  to infer genetic population structure using genotypes from 14 microsatellite loci of 381 individuals. All individuals were combined into one dataset for analysis, without anya priori population assignments and admixture was allowed. We evaluated K values, the number of assumed populations, from 1–16 using a burn-in of 50000 iterations followed by 500000 iterations for each value of K. Each value of K was run a minimum of three times to evaluate stability (see Additional file 19). We then calculated the posterior probability of population assignment to one of the six subspecies using initial assignments based on thea priori K = 13 cluster proportion results, with the migration parameter set to γ = 0.1 (see Additional file 13A). We also utilized the program GENECLASS2  for a comparative estimate of population assignments using the same K = 13 cluster proportion results for initial population assignments. We used the Rannala and Mountain  Bayesian assignment method with the simulation method of Paetkau  and an assignment threshold level of 0.05 (see Additional file 13B).
Population differentiation and inbreeding coefficients
We calculated pairwise Fst values using microsatellite results for population comparisons at the subspecies, population and sampling site levels using the program Fstat  (see Additional file 16A). Significant values of Fst were determined using the G test in Fstat (α = 0.05). We also calculated Nei's genetic distance  for population comparisons at the subspecies level (GENALEX ) (see Additional file 16B). Population inbreeding coefficients (Fis) were also calculated using Fstat and significant values determined using α = 0.05 (see Additional file 20).
Recent migration rates between the six giraffe subspecies were estimated using a Bayesian MCMC analysis of microsatellite genotypes (BayesAss 1.3 ). Individuals were preassigned to the six subspecies based on sampling location. We used 3000000 iterations, a sampling frequency of 2000, a burn-in length of 999999 iterations, and delta values for allele frequency, migration rate and level of inbreeding of 0.15 (see Additional file 17).
Isolation by distance
We tested for isolation by distance between subspecies, populations and sample locations using a comparison of genetic distance (Fst/(1-Fst)) with geographic distance, applying the Mantel test in GENALEX  (999 permutations, significance level p < 0.01) (see Additional file 21).
Molecular analysis of variance – microsatellites
We calculated molecular analysis of variance (AMOVA) for microsatellite data at the subspecies (Q = 6) and population (Q = 10) levels (999 permutations) using the program GENALEX  (see Additional file 22).