Samples of hair follicles from native donkeys (Brazilian donkey: 10 samples; Nordestino donkey: 10 samples; Pêga donkey: 10 samples) were studied (Table S1). The 524-bp mitochondrial DNA D-loop sequences of the present study were deposited in GenBank (accession numbers OM416433 - OM416462).
We collected biological samples (hair follicles) from 30 donkeys (Equus asinus) of three Brazilian breeds: Brazilian donkey (10 samples); Nordestino donkey (10 samples), and Pêga donkey (10 samples) (Figure S1). For phylogenetic inferences, we selected the species Equus caballus (NC_001640.1; Xu and Arnason 1994) as outlier group.
Collection of molecular data
Total DNA was extracted from the hair follicle samples using the DNA NucleoSpin® Tissue kit (Macherey-Nagel). The DNA samples of the Brazilian donkey were kindly provided by the Brazilian Agricultural Research Corporation (Embrapa) through a Material Transfer Agreement (23066.043532/2019-91). The PCR assays were conducted in a Veriti 96 Well Thermal Cycler (Applied Biosystems) in a final volume of 25 µl containing 2.0 µl DNA (9 to 15 ng), 0.5 µl of each primer (10 pmol/µl), 6.6 µl of Taq mix (dNTP, buffer, MgCl2 and Taq polymerase), and 15.4 µl Milli-Q water.
The primers described by Lopes et al. (2005) were used for amplification and sequencing of fragments of a mitochondrial locus, control region (displacement loop, D-loop) (F: CTGGTCTTGTAAACC and R: ACAGTTATGTGTGAGCATGG). The sequences were amplified between positions 15351 and 15949 bp. Donkey sequences available in GenBank were added, as well as sequences of Equus caballus as outlier group, Equus africanus of clade I, and Equus africanus somaliensis of clade II (Table S1).
The amplification conditions were one cycle of initial denaturation at 94oC for 5 min, followed by 35 cycles at 94oC for 1 min, 54oC for 1 min and 72oC for 1 min, and a final extension at 72oC for 5 min. The PCR products were purified using precipitation with 20% polyethylene glycol (Sambrook et al. 1989) and sequenced in an ABI PRISM 3500 automatic sequencer (Applied Biosystems). All sequences generated in this study were deposited in GenBank under the following accession numbers: Brazilian donkey (OM416453 - OM416462); Nordestino donkey (OM416433 - OM416442); Pêga donkey (OM416443 - OM416452).
The BioEdit v.184.108.40.206 software (Hall 1999) was used to edit the sequences. The sequences were aligned using the ClustalW Multiple alignment algorithm implemented in BioEdit (Hall 1999; Thompson et al. 2002) and the CLUSTAL algorithm in MEGA 5 (Tamura et al. 2011). The final alignment length was 525 bp.
Phylogenetic reconstruction, population structure and genetic diversity
After editing, the sequences were used for analysis of the number of haplotypes (H), haplotype diversity (Hd), nucleotide diversity (π), molecular analysis of variance (AMOVA), haplotype structure, and phylogenetic reconstruction using Bayesian inference and the maximum likelihood method. The following softwares were used for these analyses: DnaSP (Librado and Rozas 2009), PopART (Population Analysis with Reticulate Trees) (Leigh and Bryant 2015), Arlequin, MrBayes v. 3.2.6 (Ronquist et al. 2012), and RAxML-HPC v.8.2.12 (Stamatakis 2014), respectively. All analyses were performed remotely through the CIPRES Science Gateway 3.3 (http://www.phylo.org/index.php/portal/) (Miller et al. 2011). The evolutionary model HKY + 1 was calculated by Bayesian inference using the Kakusan4 software (Tanab, 2011).
For Bayesian inference, we performed two independent runs of 10 million iterations, with four Markov Monte Carlo chains (MCMC), sampling a tree every 10,000 iterations. Stationarity and convergence of the runs (effective sample size - ESS > 200) were evaluated with Tracer v.1.6 (http://beast.bio.ed.ac.uk/Tracer). The potential scale reduction factor was also used to check convergence of the chain and burn-in (Gelman and Rubin 1992). For maximum likelihood estimation, we used 1000 bootstrap replicates.