G C A T T A C G G C A T genes Article Maternal Phylogenetic Relationships and Genetic Variation among Rare, Phenotypically Similar Donkey Breeds Andrea Mazzatenta 1,2, * , Massimo Vignoli 1 , Maurizio Caputo 1 , Giorgio Vignola 1, * , Roberto Tamburro 1 , Francesco De Sanctis 3 , Jordi Mirò Roig 4 , Roberta Bucci 1 , Domenico Robbe 1 and Augusto Carluccio 1 1 Faculty of Veterinary Medicine, University of Teramo, SP 18, 64100 Teramo, Italy; mvignoli@unite.it (M.V.); mcaputo@unite.it (M.C.); rtamburro@unite.it (R.T.); rbucci@unite.it (R.B.); drobbe@unite.it (D.R.); acarluccio@unite.it (A.C.) 2 Physiology and Physiopathology Section, Neuroscience, Imaging and Clinical Sciences Department, University G. D’Annunzio of Chieti-Pescara, Via Dei Vestini, 31, 66100 Chieti, Italy 3 Section of Immunology, Department of Medicine, University of Verona, Piazzale L.A. Scuro 10, 37134 Verona, Italy; Francesco.desanctis@univr.it 4 Servei de Reproducció Equina, Departament de Medicina i Cirurgia Animals, Autonomous University of Barcelona, Edifici V, 08193 Bellaterra, Spain; jordi.miro@uab.es * Correspondence: amazzatenta@unich.it (A.M.); gvignola@unite.it (G.V.) Abstract: The mitochondrial DNA (mtDNA) D-loop of endangered and critically endangered breeds has been studied to identify maternal lineages, characterize genetic inheritance, reconstruct phyloge- netic relations among breeds, and develop biodiversity conservation and breeding programs. The aim of the study was to determine the variability remaining and the phylogenetic relationship of Martina Franca (MF, with total population of 160 females and 36 males), Ragusano (RG, 344 females and  30 males), Pantesco (PT, 47 females and 15 males), and Catalonian (CT) donkeys by collecting genetic  data from maternal lineages. Genetic material was collected from saliva, and a 350 bp fragment of Citation: Mazzatenta, A.; Vignoli, M.; D-loop mtDNA was amplified and sequenced. Sequences were aligned and evaluated using standard Caputo, M.; Vignola, G.; Tamburro, bioinformatics software. A total of 56 haplotypes including 33 polymorphic sites were found in R.; De Sanctis, F.; Roig, J.M.; Bucci, R.; 77 samples (27 MF, 22 RG, 8 PT, 19 CT, 1 crossbred). The breed nucleotide diversity value (π) for Robbe, D.; Carluccio, A. Maternal all the breeds was 0.128 (MF: 0.162, RG: 0.132, PT: 0.025, CT: 0.038). Principal components analysis Phylogenetic Relationships and grouped most of the haplogroups into two different clusters, I (including all haplotypes from PT and Genetic Variation among Rare, CT, together with haplotypes from MF and RG) and II (including haplotypes from MF and RG only). Phenotypically Similar Donkey Breeds. Genes 2021, 12, 1109. https:// In conclusion, we found that the primeval haplotypes, haplogroup variability, and a large number of doi.org/10.3390/genes12081109 maternal lineages were preserved in MF and RG; thus, these breeds play putative pivotal roles in the phyletic relationships of donkey breeds. Maternal inheritance is indispensable genetic information Received: 24 June 2021 required to evaluate inheritance, variability, and breeding programs. Accepted: 13 July 2021 Published: 22 July 2021 Keywords: genetics; mitochondrial DNA; donkey; Martina Franca; Ragusano; Pantesco; Catalonian Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affil- 1. Introduction iations. The genus Equus is the only remaining member of the family Equidae, which includes both extant and fossil species [1]. The noncaballine forms include the African wild ass Equus africanus; zebras Equus quagga (formerly Equus burchellii), Equus grevyi, and Equus zebra (with two subspecies Equus zebra zebra of South Africa and Equus zebra hartmannae Copyright: © 2021 by the authors. of Namibia and Angola); and the Asian wild asses Equus kiang and Equus hemionus (with Licensee MDPI, Basel, Switzerland. subspecies Equus hemionus kulan and Equus hemionus onager) [2,3]. The domestic donkey This article is an open access article E. africanus is accepted as a subspecies of the African wild ass [4,5]. distributed under the terms and Geographically isolated donkey populations, referred to as breeds here, including conditions of the Creative Commons Martina Franca (MF) of the continental Puglia region, Ragusano (RG) of the island province Attribution (CC BY) license (https:// of Sicily, Pantesco (PT) of the small island of Pantelleria (85 km2 ), and Catalonian (CT) of the creativecommons.org/licenses/by/ Spanish Catalonian region, share similar phenotypes due to genetic inheritance/ecological 4.0/). Genes 2021, 12, 1109. https://doi.org/10.3390/genes12081109 https://www.mdpi.com/journal/genes Genes 2021, 12, 1109 2 of 14 constraints and the outcomes of biodiversity conservation programs. The natural history of the donkey is of interest because the desire to safeguard genetic biodiversity is growing [6] and because the three autochthonous Italian breeds are considered endangered, including MF with 160 females and 36 males (Associazione Nazionale Asino Martina Franca) and RG with 344 females and 30 males (Domestic Animals Diversity Information System), which are critically endangered, as is PT with 47 females and 15 males (Associazione Italiana Allevatori). The phylogeny of the domestic donkey is not yet well understood. Historical files show that the Roman Empire dominated the Iberian Peninsula for ~613 years from 218 BC to 395 [7]. At that time, the MF donkey was commonly used for transport and military operations; thus, it is presumed that it was widely distributed across Spain [8]. Later, Spain dominated central Southern Italy for ~148 years from 1559–1707 BC, and CT was likely introduced during this time in Italy. However, phenotypic traits and anthropological documents are often insufficient to ascertain breed history, origin, and the occurrence of genetic exchange [9]. Instead, mitochondrial DNA (mtDNA) sequencing can determine intra- and inter-species historical, biogeographic, and phylogenetic rela- tionships [10]. The extrachromosomal mitochondrial genome, unlike the nuclear one, is inherited only through the maternal lineage, is haploid, and does not undergo genetic recombination [11,12]. The application of clonal polymorphisms to study the genetics of domestic animals is valuable [8,10]. Variation in the D-loop region of mtDNA and the lack of recombination in mtDNA make it a highly informative tool for matrilineal studies, for determining intraspecies phylogenetic relationships, and for characterizing intrabreed variation [12–16]. mtDNA studies of dog breeds, which have greater phenotypic and working variability compared to the donkey, which is relatively uniform, have revealed genetic information on their domestication, evolution, and hereditary diseases [17,18]. mtDNA studies of equine breeds were used to investigate their origin [19–26] and to track breed migration and distribution by comparing the maternal lines in different populations [27,28]. The complete donkey mitochondrial genome sequence was essential to date the divergence from the horse between 8 and 10 MYA [29,30], which is earlier than paleontological data [24,31] and data from restriction endonuclease analysis [32]. Interestingly, two lineages of the domestic donkey were identified using mtDNA: Clade 1 for the Nubian lineage (E. a. africanus) and Clade 2 for the Somali lineage (E. a. somaliensis). These lineages resulted from two separate domestication events among two wild ancestral populations located in (1) the Atbara region and Red Sea Hills (NW Sudan) and in (2) southern Eritrea, Ethiopia, and Somalia [5,33–36]. However, the existence of another ancestor of the domestic donkey belonging to an unrecognized extinct African wild population has been suggested [6,36,37]. Genetic studies on the biodiversity of the Italian donkey are limited and have primarily focused on variability among protein markers and microsatellites [38–40]. Recently, whole genome sequencing [41] and mtDNA [42] were used to study the evolution and genetic diversity of Italian donkey populations. In this study, we evaluated the mtDNA D-loops of endangered and critically endan- gered Italian donkey breeds. mtDNA sequences, single-nucleotide polymorphisms (SNPs), and haplotypes were identified and analyzed to investigate the matrilineal assortment within and between asinine breeds with similar phenotypes and to investigate the ori- gin and phylogenetic relationships between asinine breeds to better manage rare donkey breeds by establishing proper breeding and conservation programs. 2. Materials and Methods One-hundred and twenty-three salivary samples were collected from eight official breeding stations (Supplementary Table S1 additional material; Istituto Incremento Ip- pico associated with Facoltà di Medicina Veterinaria Università degli Studi di Teramo (Fondo Rustico Chiareto) and Centro di Conservazione del Patrimonio Genetico dell’Asino 2. Materials and Methods One-hundred and twenty-three salivary samples were collected from eight official Genes 2021, 12, 1109 3 of 14 breeding stations (Supplementary Table S1 additional material; Istituto Incremento Ippico associated with Facoltà di Medicina Veterinaria Università degli Studi di Teramo (Fondo Rustico Chiareto) and Centro di Conservazione del Patrimonio Genetico dell’Asino della razza Martina della razza FrancaFranca Martina (Azienda agricola (Azienda RussoliRussoli agricola Crispiano); Istituto Istituto Crispiano); Incremento Ippico Re- Incremento Ip- pico Regione gione Campania Campania Santa MariaSanta MariaVetere Capua Capua andVetere Aziendaand Agricola AziendaCiro Agricola SchiròCiro Schirò Corleone- Corleone-Monreale Monreale unfortunately unfortunately have no successful have no successful sequencedsequenced samples) samples) in accordance in accordance with the with the standards standards for careforandcareprotection and protection of animals of animals used used forfor scientificpurposes scientific purposes Directive Directive 2010/63/EU. 2010/63/EU. This This study study waswas approved approved bybythethe Ethics Ethics Committee Committee (Protocol (Protocol No.No. 62128 62128 of of 27 27 April April 2018). 2018). TheThe samples samples collected collected werefrom were fromfree-range free-rangeanimals animalswithwith certificates certificates of origin, which origin, whichwerewereused usedtotoexclude excludeanimals animals ofof thethe same same maternal maternal descent, descent, in order in order to to in- increase crease genetic genetic variability variability in the in the sample sample set.set. Seventy-seven Seventy-seven samples, samples, including including MF MF= 27,=RG 27, =RG22,= PT 22, =PT8,=CT 8, CT = 19, = 19, andand 1 Italian 1 Italian crossbreed,were crossbreed, weresequenced sequencedsuccessfully successfully(Figure (Figure 1), 1), whereas the remaining 46 samples were whereas the remaining 46 samples were corrupted. corrupted. Figure 1. Breed Figure 1. Breeddistribution. distribution.The Thered red perimeters perimeters pinpoint pinpoint the the original original Mediterranean Mediterranean distribution distribution of theof thedonkey four four donkey breeds breeds investigated: Martina Franca (MF), Ragusano (RG), Pantesco (PT), and Catalonian investigated: Martina Franca (MF), Ragusano (RG), Pantesco (PT), and Catalonian (CT). (CT). Genetic Genetic material material was was collected collected from from saliva saliva using using aa sterile sterile oral swab, transferred to FTA FTA mini-cards, and stored in multibarrier multibarrier pouches (Whatman Labware pouches (Whatman Labware Products, Products, U.K.). U.K.). The The reference material is is available available at at the theO.V.U.D. O.V.U.D.(University (UniversityVeterinary VeterinaryHospital) Hospital)Centre Cen- tre for for thethe breeding breeding of of large large animalsatatthe animals theFaculty FacultyofofVeterinary VeterinaryMedicine, Medicine, University of Teramo, Teramo, Italy. Based on the complete donkey mtDNA mtDNA sequence sequence (GenBank (GenBank X97337) X97337) [29], [29], two pairs two pairs of ofprimers primerswere were designed designed to amplify to amplifythe hypervariable the hypervariableregion between region sites 15390 between sites and 15750 15390 [43], which and 15750 is a fragment [43], which of the D-loop is a fragment mtDNA of the D-loop (http://bioinfo.ut.ee/primer3, mtDNA (http://bioinfo.ut.ee/pri- accessed on 4 April 2019). After extraction from the FTA mer3, accessed on 4 April 2019). After extraction from the FTA mini-card, mini-card, DNA DNAwas amplified was am- by PCRby plified inPCR a 25 inµLareaction containing 25 μL reaction 50 ng of containing 50DNA, ng of 2.5 DNA,mM2.5 MgCl2, 0.2 mM0.2 mM MgCl2, each mMdNTP,each 0.5 µM 0.5 dNTP, PERμM 50 - PER CC AAG5′- CC GAC AAG TATGACCAA TATGGACAAAG-3 0 and FOR 50 -TTG GAG GGA TTG CTG GGA AG-3′ and FOR 5′-TTG GAG GGA ATT CTG TTG TC-30ATT primers, × PCR buffer, TC-3′1primers, 1× PCR and 1 U of buffer, Taq and 1UDNA polymerase of Taq (Fermentas, DNA polymerase Thermo (Fermentas, Thermo Fischer Scientific). The amplification was performed using the Mastercyclercycler Fischer Scientific). The amplification was performed using the Mastercycler thermal ther- (Eppendorf, mal USA) withUSA) cycler (Eppendorf, the following with the conditions: initial denaturation following conditions: at 94 ◦ C for initial denaturation at 594min °C followed for by 35 cycles 5 min followed by 35 94 ◦ C of of cycles for9430°Cs,for ◦ C for 30 s, 72 ◦ C for 30 s, and then, a final 58 30 s, 58 °C for 30 s, 72 °C for 30 s, and then, a extension at 72 at◦ C for 5 min. final extension 72 °C for 5 min. The raw sequence trace files were checked for the presence of ambiguous bases using the software Chromas v.2.5.1 (http://www.technelysium.com.au/, accessed on 4 April 2019). Sequences were aligned with Muscle; in Supplementary Table S2 of the Supple- mentary Materials is the alignment with other Italian donkey breeds [42]. The number of polymorphic sites (parsimony informative and singleton sites), the number of haplotypes Genes 2021, 12, 1109 4 of 14 (private and shared haplotypes), nucleotide diversity, and the average number of nucleotide differences were calculated according to Tajima (1983) and Nei (1987) using MEGA7, Fu’s neutrality statistic test, and Tajima’s D test with DnaSP 6.12.01, as well as using a maximum parsimony analysis and the maximum composite likelihood method. The median-joining network (for sequences, see Supplementary Table S3 in the Supplementary Materials) and principal coordinates analysis (PCoA) were performed with DARwin software [44–46]; Supplementary Figures S1 and S2 were created by using itol.embl.de/tree [47]. The other statistical analyses were performed with Statistica 7.0 StatSoft. 3. Results 3.1. Breed Haplotype Analysis The successfully analyzed samples from the eight certified breeding centers (Table S1) included MF, RG, PT, and CT donkeys. The mtDNA D-loop hypervariable region between sites 15390 and 15750 (GenBank ID # 2466755, Table S3) was fully sequenced for 77 sam- ples, and 56 haplotypes including 33 polymorphic sites were found. Of these, 14.6% of haplotypes have frequencies greater than 2.7%: Hap 1, 22 and 36; Hap 30, 37, 40 and 53 F = 4.1%; Hap 51 F = 12.2% the most common; while 85.4% of haplotypes are rare and have frequencies of 1.4%. Table 1 shows molecular diversity indices for each breed. The breed genetic diversity for all the breeds, evaluated by the nucleotide diversity value (π), was 0.128. Within subpopulations, π was 0.098, and the mean interpopulation evolutionary π was 0.03. The haplotypes identified in the analyzed breeds included 6 for PT, 14 for CT, and 22 each for MF and RG. Table 1. D-loop nucleotide polymorphisms and molecular diversity indices per breed tested in the study. Breed n NHap SNPs π Martina Franca (MF) 27 22(5 s) 13 0.162 ± 0.022 Ragusano (RG) 22 22(3 s) 15 0.132 ± 0.028 Pantesco (PT) 8 6(3 s) 1 0.025 ± 0.001 Catalonian (CT) 19 13(4 s) 2 0.038 ± 0.009 crossbreed 1 1 1 ALL 77 56 33 0.128 s, shared haplotypes. In this dataset, there is a crossbreed with its own haplotype. n, number of individuals sampled per breed; NHap, the number of haplotypes in each breed with the number of shared haplotypes in parentheses; SNPs, the number of polymorphic sites; π, nucleotide diversity with standard deviation. All the sequenced samples were aligned with a reference sequence (GenBank X97337) to highlight the presence of SNPs. The absolute number of mutations identified for each sequenced sample is shown (Figure 2A, Table 2 and Table S4). Moreover, all identified SNPs were classified according to their positions in the reference sequence and were charac- terized as transitions (purine–purine and pyrimidine–pyrimidine substitutions) or deletion (Figure 2B, left panel; every column represents a sequenced sample). The percentages of samples in which each SNP was identified are shown as histograms (Figure 2B, right panels, Table 2 and Table S4). The multivariate test of significance for nucleotide frequencies (Table S5) showed no differences in nucleotide composition among breeds (p = 0.98). Maximum composite likelihood estimates of the nucleotide substitution pattern per breed, positions containing gaps, and missing data were eliminated as recommended by the literature (Supplementary Table S6) [45]. Genes2021, Genes 12,x1109 2021,12, FOR PEER REVIEW 5 5ofof15 14 Figure 2. SNP analysis performed on mtDNA isolated from Martina Franca (MF), Ragusano (RG), Pantesco (PT), and Catalonian Figure 2. SNP(CT) donkeyperformed analysis breeds. (A)on Absolute mtDNAnumber of from isolated mutations found Martina in each Franca sample (MF), by Sanger Ragusano (RG),sequencing. (B) Every Pantesco (PT), and SNP was characterized Catalonian (CT) donkey according breeds.(A)to the position Absolute numberon the reference found of mutations sequence, the sample in each kind of by mutation, and the percentage Sanger sequencing. of (B) Every samples SNP was containing the according characterized mutation (brown histograms) to the position using on the MEGA7. reference R = reference sequence, the kindsequence (X97337). of mutation, and the percentage of samples containing the mutation (brown histograms) using MEGA7. R = reference sequence (X97337). Genes 2021, 12, 1109 6 of 14 Table 2. Absolute and relative haplotype frequencies in the four analyzed donkey breeds. SNP Martina Franca (MF) Catalano (CT) Pantesco (PT) Ragusano (RG) Position Reference Mutation Samples (n) % Samples (n) % Samples (n) % Samples (n) % 15719 C T 0 0 0 0 0 0 1 5 15699 C T 6 21 0 0 0 0 4 18 15663 A G 6 21 0 0 0 0 4 18 15653 C T 6 21 0 0 0 0 4 18 15646 A G 2 7 0 0 0 0 1 5 15645 G A 6 21 0 0 0 0 4 18 15622 A G 4 14 0 0 0 0 2 9 15600 A G 29 100 19 100 8 100 17 77 15599 C T 6 21 0 0 0 0 4 18 15581 A G 6 21 11 58 0 0 4 18 15570 A G 6 21 0 0 0 0 4 18 15524 C DEL 0 0 0 0 0 0 2 9 15504 T C 6 21 0 0 0 0 4 18 15491 C T 6 21 0 0 0 0 4 18 15485 G A 6 21 0 0 0 0 4 18 15474 C T 0 0 0 0 0 0 2 9 3.1.1. Population Analyses Higher base composition differences within and between breed sequences were found for MF and RG (Table 3). Table 3. Estimates of base composition differences within and between breed sequences using MEGA7. Breed Mean within Groups Martina Franca (MF) 0.289 Ragusano (RG) 0.521 Pantesco (PT) 0.044 Catalonian (CT) 0.020 Mean between Groups MF × CT 0.162 MF × RG 0.387 MF × PT 0.158 RG × CT 0.276 RG × PT 0.278 CT × PT 0.031 Despite the limited number of samples, the following haplotypes were found: 17 in MF (Hap 1, 3, 5, 9, 10, 11, 15, 20, 24, 25, 34, 39, 41, 46, 47, 50, 54); 19 in RG (Hap 4, 6, 8, 12, 13, 16, 17, 18, 19, 21, 23, 29, 38, 42, 43, 44, 48, 49, 55); 9 in CT (Hap 2, 14, 26, 27, 28, 31, 33, 45, 52); and 3 in PT (Hap 7, 32, 35). The most represented seven haplotypes were Hap 51, common to all breeds; Hap 22, 36, and 37 found in MF and shared with PT and/or RG; Hap 30, distinctive of MF and also found in CT; and Hap 40 and 53, characteristic of CT. The sequences were aligned with the reference sequence GenBank X97337 and other similar GenBank sequences representative of the Somali and Nubian African donkey lines, Chinese and other Asiatic E. kiang, E. hemionus, and E. h. kulan lines (Figure S1), and other Italian breeds (Figure S2, Supplementary Materials as in Ref. [42]: Romagnolo donkey (ROD), Amiata donkey (AMD), Sardinian donkey (SAD), Asinara donkey (ASD), Ragusano donkey (RAD)). 3.1.2. Origin and Phylogenetic Relationships The phylogenetic relationships among the 55 haplotypes were calculated using a median-joining network. The haplotypes of each breed are color coded, the abundance of the haplotype indicated by the relative size of the symbol, and diffusion among the breeds is shown in color-coded pie charts (Figure 3). Looking at the distribution of haplotypes by breed, most of them are represented alone. In particular, MF and RG diverge from Genes 2021, 12, 1109 7 of 14 the reference sequence X97337 hap and show the highest variability. Conversely, PT is closely related and even shares haps with MF. Most CTs are grouped into two clusters: CT hap outgroups are closely related mainly to MF and completely share two haps with MF. Genes 2021, 12, x FOR PEER REVIEW 8 of 15 Hap 51 is the most common and is shared among all races. The Fu neutrality statistic test and Tajima’s D test were performed to address the hypothesis of population expansion, and the results were Fs = −2.523, p < 0.05, and D = 2.004, p < 0.05, respectively. Median-joining Figure3. 3.Median-joining Figure network network analysis. analysis. This This analysis analysis was was based based onon 350 350 bpbp ofof the the mtDNA mtDNA D-D-loop loop sequences of 77 donkeys (Martina Franca (MF), Ragusano (RG), Pantesco (PT), Catalonian (CT)), which consists of 55 haplotypes (DARwin 6.0). Each breed is color coded, and for each hap- lotype, the proportions of the different breeds are shown, in black the crossbreed haplotype. The reference sequence X97337 is indicated as x. Genes 2021, 12, 1109 8 of 14 sequences of 77 donkeys (Martina Franca (MF), Ragusano (RG), Pantesco (PT), Catalonian (CT)), Genes 2021, 12, x FOR PEER REVIEW which consists of 55 haplotypes (DARwin 6.0). Each breed is color coded, and for each haplotype, 9 of 15 the proportions of the different breeds are shown, in black the crossbreed haplotype. The reference sequence X97337 is indicated as x. PCoA PCoAanalysis analysis based on on based the the dissimilarity matrix dissimilarity returned matrix two different returned clusters, two different clus- clusters, ters I andIII. clusters andInterestingly, in cluster II. Interestingly, II, there in cluster are only II, there MF (6, are only MF 17,(6,27, 17,35) 27,and 35) RG and(9, RG21,(9,35, 21, 53), while the rest of the haplotypes are grouped into cluster I. However, six 35, 53), while the rest of the haplotypes are grouped into cluster I. However, six haplotypes haplotypes are arenot notincluded includedininclusters clustersI Iand andII:II:they theyare areMFMFHap Hap22, 22,37, 37,and and40; 40;and andRGRGHap Hap13,13,16, 16, and 38 (Figure 4). and 38 (Figure 4). Figure Figure4.4.PCoA. PCoA.The Thetwo twoaxes axesofofthe thePCoA PCoAplot plotwere werebased basedon onthe thedissimilarity dissimilaritymatrix matrixaccording accordingtotoKimura Kimura(1980) (1980)and and generated generatedusing usingDARwin DARwin6.0.6.0.Two Two different different clusters, cluster I and cluster cluster II, II, included includedmost mosthaplotypes haplotypesofofthethebreeds breedsstudied stud- ied (Martina (Martina Franca Franca (MF), (MF), Ragusano Ragusano (RG), (RG), Pantesco Pantesco (PT), (PT), Catalonian Catalonian (CT)). (CT)). Haplotypes Haplotypes thatthat did did not into not fall fall into either either cluster cluster were were MF Hap MF Hap 22,and 22, 37, 37, 40 andand 40 RG andHapRG Hap 13,and 13, 16, 16, 38. and 38. 4.4.Discussion Discussion 4.1. BreedMolecular 4.1. Breed MolecularAnalysis Analysis Donkeybreeds Donkey breedsrepresent representaafascinating fascinatingmodel modelofofdomesticated domesticatedbiodiversity; biodiversity;thus, thus,aa number of studies have analyzed donkey pedigree and genetics. Pedigree number of studies have analyzed donkey pedigree and genetics. Pedigree and reconstruc- and reconstruc- tionstudies tion studiesusually usuallylack lackthe thecorresponding correspondinggenetics genetics[1,48,49], [1,48,49],andandgenetic geneticstudies studieshave have frequentlylacked frequently lackedlineage lineageinformation information[39,42,50]. [39,42,50]. Pedigreestudies Pedigree studies onon Italian Italian (MF, (MF,Amiata) Amiata)andandSpanish Spanish(CT, Andalusian, (CT, Miranda) Andalusian, don- Miranda) keys found dramatic losses in genetic variation due to high rates of donkeys found dramatic losses in genetic variation due to high rates of inbreeding inbreeding [1,48,49,51]. However, pedigree [1,48,49,51]. However,incompleteness and the occurrence pedigree incompleteness of a bottleneck and the occurrence event (e.g., of a bottleneck MF event in 1980) may have led to over- or under-estimations of genetic variation, (e.g., MF in 1980) may have led to over- or under-estimations of genetic variation, which which could affectaffect could breeding strategies breeding [1,51].[1,51]. strategies To overcome this bias, To overcome this in thisinstudy, bias, genetic this study, analysis genetic was analy- performed on subjects with certificates of origin from official authorized breeding sis was performed on subjects with certificates of origin from official authorized breeding centers. This approach centers. allowedallowed This approach us to link usmtDNA data to pedigree to link mtDNA record, population, data to pedigree and breed record, population, and breed to identify the same maternal descendants and to select individuals with pre- sumed higher genetic variability to preserve biodiversity. A further bias in breed studies comes from unbalanced sample size comparisons for populations (e.g., difference of five- fold [39,42]), genetic structure, genetic variability, genetic robustness, average relatedness, Genes 2021, 12, 1109 9 of 14 to identify the same maternal descendants and to select individuals with presumed higher genetic variability to preserve biodiversity. A further bias in breed studies comes from unbalanced sample size comparisons for populations (e.g., difference of five-fold [39,42]), genetic structure, genetic variability, genetic robustness, average relatedness, inbreeding, co-ancestry, the degree of nonrandom mating, and origin [50,52]. In this study, we collected a similar number of samples per breed, except PT because of its limited number of breed lines; consequently, we analyzed a balanced sample. We found higher genetic variability in MF and RG, which disagrees with studies based on pedigree [1]. RG and MF are widely used on farms where natural random mating still occurs, and we believe that is the source of the variability. We found lower variability in PT, as expected, because of the limited number of individuals (n = 62) and because of its isolation. In the PT certificates of origin, eight distinct maternal lineages are attested, which is in agreement with the statement that the recovery of the PT breed started from a small nucleus of founders [39]. Furthermore, based on the number of haplotypes, the genetic robustness of PT is dramatically lower and the pedigree certificate alone cannot predict it. The low variability found in CT was not predicted because of the number of individuals in the population and the common distribution in the large region of Catalonia. This phenomenon could be the result of a breeding program with an unbalanced number of males vs. females (517 vs. 310 respec- tively; Asociacio del Foment dela Raca Asina Catala), rather than the result of the number of subjects and area occupied. Consequently, a potential bottleneck was produced by human artificial selection, which led to a loss in variability [42,49,53]. Because variability, robustness, and the degree of nonrandom mating decrease as average relatedness increases due to inbreeding and co-ancestry, a new reproduction program with multiple approaches is needed. The molecular characteristics analyzed showed distinct nucleotide frequencies among breeds, which is in line with the literature [29,30]. The transition rate was greater for purines, and the transition rate per breed decreased from MF to RG to CT, which is partially discordant with other studies on the Italian donkey [42]. The nucleotide diversity (π) values are in line with a previous study on CT [54]. Cozzi et al. 2017 [42] found high nucleotide diversity in the MF donkey, but found low nucleotide diversity in the RG donkey. We found the highest diversity in MF followed by RG and the lowest variability in CT followed by PT. The differences in our results and a previous study on the Italian donkey could be due to the following biases: (i) a lack of certificate of origin information leading to uncertain breed origin, (ii) unbalanced sampling between breeds with 74% of the breeds being Asinara and Sardo donkeys, both from the island of Sardinia [42]. Overall, molecular indices show greater genetic variability among the Italian than the Spanish donkeys according to the previous literature [42,49,53]. 4.2. Population Analyses The population analyses showed more differences within RG, followed by MF, CT, and PT. The RG genetic population structure is analogous to the highly heterogeneous large maternal Balkan donkey population with a more complex genetic structure than previously thought [53]. The Balkan donkey population is highly genetically diverse despite their severe population decline, probably due to introgression of other related breeds [53]. Introgression probably occurred for RG, but not MF. Our analysis of PT is consistent with a microsatellite genetic variability study that showed lower variability in PT compared to RG and another Sicilian breed, the Grigio Siciliano (GS) [39]. The genetic variability observed in PT, RG, and GS [39] is lower than that reported in five Spanish breeds [54] and three Croatian breeds [55], but higher than that observed in the Amiata donkey from Italy [31] and in Chinese breeds [36]. The comparisons between breeds showed more differences between MF and RG and fewer differences between PT and CT, between RG and PT, and between RG and CT. The fewest differences are between PT and CT. Genes 2021, 12, 1109 10 of 14 Two different studies on the Balkan donkey revealed different interpretations. In the first study, no correspondence between geographical areas and maternal genetic structure was found. Because the difference between the Balkan donkey and the African Burkina Faso donkey outgroup was also low, the authors could not trace routes of expansion in the donkey; consequently, they suggested that the species spread very quickly after domestication [56]. The second study assessed three Balkan donkey populations: Istrian (IS), north Adriatic (NA), and Littoral-Dinaric (LD), and their results suggested that IS is a unique breed, which mixed with LD during sporadic migration events, and that NA and LD are genetically similar [55]. Our study suggests similar effects of migration by MF on CT and PT, which are in accordance with historical reconstructions. 4.3. Origin and Phylogenetic Relationships The well-established identification of two main lineages and the probable existence of another unrecognized extinct wild ancestor in domestic donkeys are believed to be the result of separate domestication events [5,6,33–37]. However, the genetic structure of the Chinese donkey indicates another possible line [36]. Therefore, in the donkey, such as in the dog [17], genetic data support multicentric breed origins. In agreement with this emerging theory is the identification of a potential new clade unique from the MF and RG Italian donkeys. Furthermore, another recent study also suggested multiple breed origins [42]. In Croatian and Serbian donkeys, three haplotype groups were found [55] with distinct nuclear gene pools [53]. A heterogeneous genetic structure of the Balkan donkey was hypothesized because there was no geographical structure; thus, it was difficult to trace the routes of expansion in the donkey [56]. However, other hypotheses for the complex genetic relationships among Italian donkey breeds and breeds living in the Mediterranean and Balkan areas [42] include ancestry and the genetic makeup of modern donkey populations. Our analysis suggests a multicentric domestication phenomena coupled with multiple waves of colonization and counter colonization, such as what occurred when the Roman Empire brought the MF donkey to Spain and when Spain brought the MF donkey to Italy; this is in accordance with the hypothesis suggested by Stanisic and colleagues [53]. There is a rising interest in maintaining genetic diversity in animal populations to safeguard the widest possible genetic resources through conservation programs [57,58]. For conservation of domestic breeds, the preservation of genetic capital is crucial because they are already zootechnical forms with reduced original natural variability. The biodiversity of the domestic form is the result of a genetic pool derived from interaction with semi-artificial environments regulated by human needs and human migratory movements. As a result, genetic and phenotypic changes with respect to the wild species of origin have been addressed slowly based on human necessity. Donkey conservation represents a biological problem with regard to analogous phe- notypes. The similar donkey breed phenotypes may be a result of genetic exchange among breeds, identical origin from an African or Asian population, equivalent climate conditions, and/or similar types of work carried out by the donkeys. Our molecular analyses showed a dramatic loss in variability across all breeds tested, but especially in PT and CT. MF and RG have higher numbers of haplotypes and SNPs than PT and CT. In practice, the MF donkey is an important reservoir of biodiversity that must be preserved with the widest possible range of its genetic heritage, and errors of consanguinity for aesthetic purposes should be avoided with the help of conservation programs. Future conservation programs should include certificates of origin and genetic analyses of the matrilines at least. 5. Conclusions In conclusion, this study: (i) identified D-loop mtDNA characteristics for the MF donkey and three phenotypically similar breeds, (ii) identified different matrilines in the MF donkey and similar breeds, (iii) identified the biodiversity of each breed, and (iv) de- termined the phylogenetic relationships among the breeds. This extensive study on the Genes 2021, 12, 1109 11 of 14 biodiversity and phylogenetic relationships of MF, RG, PT, and CT donkey breeds is useful for future domestication studies. In this study, we analyzed the genetics of a limited number of endangered and highly endangered donkey breeds. The data showed significant loss in variability among all the breeds evaluated, which is in agreement with previous studies that used different methods. The primeval haplotypes, haplogroup variability, and large number of maternal lineages are preserved in the MF and RG breeds; thus, they play putative pivotal roles in the phyletic relationships of the studied donkey breeds. Given the level of endangerment undergone by these breeds, actions are necessary to ensure their long-term survival and conservation. Improving the reproduction and management of existing populations, clarifying their historic interactions by studying their population genetics, and extending and improved monitoring of maternal lineages are valid options. Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/genes12081109/s1: Figure S1: All sequences were aligned with reference sequence GenBank X97337 and other similar sequences present in GenBank representative of African donkey lines Somalicus and Nubianus, Chinese ones, and other Asiatic E. kiang, E. hemionus, and E. h. kulan; Figure S2: The sequences aligned with the reference GenBank X97337 and other similar GenBank sequences representative of other Italian breeds (as in Ref. [42]: Romagnolo donkey (ROD), Amiata donkey (AMD), Sardinian donkey (SAD), Asinara donkey (ASD), Ragusano donkey (RAD)); Table S1: Microchip number of the samples successfully sequenced by breed and official breeding station; Table S2: CLUSTAL multiple sequence alignment by MUSCLE (3.8) of Italian donkey MF, RG, PT, and CT with other sequences of Italian donkey in GenBank (Legend: RAD = Ragusano donkey; AMD = Amiata d.; ASN = Asinara d.; SAD = Sardinian d.; ROD = Romagnolo d.); Table S3: Sequences deposited in GenBank, Submission No. 2466755; Table S4: SNP alignments performed on mtDNA isolated from MF, RG, PT, and CT donkey breeds; Table S5: Nucleotide frequencies per breed. A multivariate test of significance returned no differences in nucleotide composition among breeds (p = 0.98); Table S6: Maximum composite likelihood estimate of the pattern of nucleotide substitution. Each entry shows the probability of substitution (r) from one base (row) to another base (column). Rates of different transitional (A × G and C × T) substitutions are shown in bold and those of transversional (G × T and A × C) substitutions are shown in italics. Author Contributions: The authors contributions are as follows: conceptualization, A.M. and A.C.; methodology, A.M.; software, F.D.S.; validation, A.M. and A.C.; formal analysis, A.M. and F.D.S.; investigation, M.C., J.M.R., D.R. and R.B.; resources, A.C.; data curation, A.M., M.V., R.T., G.V. and A.C.; writing, original draft preparation, A.M., D.R., F.D.S. and A.C.; writing, review and editing, A.C., A.M., M.V., R.T. and G.V.; supervision, A.C.; funding acquisition, A.C. All authors have read and agreed to the published version of the manuscript. Funding: The present study was carried out in the framework of the Project Demetra (Dipartimenti di Eccellenza 2018–2022, CUP_C46C18000530001) funded by the Italian Ministry for Education, University, and Research. Institutional Review Board Statement: The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the by the Ethical Committee (Protocol No. 62128 of 27 April 2018) of the University of Bologna. Informed Consent Statement: Not applicable. Data Availability Statement: Sequences deposited in GenBank Submission No. 2466755 and showed in Table S3. Conflicts of Interest: The authors declare no conflict of interest. References 1. Macfadden, B.J. Evolution. Fossil Horses—Evidence for Evolution. Science 2005, 307, 1728–1730. [CrossRef] 2. Krüger, K.; Gaillard, C.; Stranzinger, G.; Rieder, S. Phylogenetic Analysis and Species Allocation of Individual Equids Using Microsatellite Data. J. Anim. Breed. Genet. 2005, 122 (Suppl. 1), 78–86. [CrossRef] Genes 2021, 12, 1109 12 of 14 3. Groves, C.; Grubb, P. Ungulate Taxonomy; The Johns Hopkins University Press: Baltimore, MD, USA, 2011; p. 336. 4. Beja-Pereira, A.; England, P.R.; Ferrand, N.; Jordan, S.; Bakhiet, A.O.; Abdalla, M.A.; Mashkour, M.; Jordana, J.; Taberlet, P.; Luikart, G. African Origins of the Domestic Donkey. Science 2004, 304, 1781. [CrossRef] 5. Kimura, B.; Marshall, F.B.; Chen, S.; Rosenbom, S.; Moehlman, P.D.; Tuross, N.; Sabin, R.C.; Peters, J.; Barich, B.; Yohannes, H.; et al. Ancient DNA from Nubian and Somali Wild Ass Provides Insights into Donkey Ancestry and Domestication. Proc. R. Soc. B Biol. Sci. 2011, 278, 50–57. [CrossRef] 6. Rizzi, R.; Tullo, E.; Cito, A.M.; Caroli, A.; Pieragostini, E. Monitoring of Genetic Diversity in the Endangered Martina Franca Donkey Population. J. Anim. Sci. 2011, 89, 1304–1311. [CrossRef] 7. Piganoil, A.F.G. La Conquête Romaine; Edité par Félix Alcan: Paris, France, 1927. 8. Rand, D.M. The Units of Selection on Mitochondrial DNA. Annu. Rev. Ecol. Syst. 2001, 32, 415–448. [CrossRef] 9. Clutton-Brock, J. Horse Power: A History of the Horse and the Donkey in Human Societies; Harvard University Press: Cambridge, MA, USA, 1992. 10. Bowling, A.T.; Del Valle, A.; Bowling, M. A Pedigree-Based Study of Mitochondrial D-Loop DNA Sequence Variation among Arabian Horses. Anim. Genet. 2000, 31, 1–7. [CrossRef] [PubMed] 11. Hutchison, C.A.; Newbold, J.E.; Potter, S.S.; Edgell, M.H. Maternal Inheritance of Mammalian Mitochondrial DNA. Nature 1974, 251, 536–538. [CrossRef] 12. Brown, W.M.; George, M.; Wilson, A.C. Rapid Evolution of Animal Mitochondrial DNA. Proc. Natl. Acad. Sci. USA 1979, 76, 1967–1971. [CrossRef] 13. Wallace, D.C. Why Do We Still Have a Maternally Inherited Mitochondrial DNA? Insights from Evolutionary Medicine. Annu. Rev. Biochem. 2007, 76, 781–821. [CrossRef] 14. Mirol, P.M.; Peral García, P.; Vega-Pla, J.L.; Dulout, F.N. Phylogenetic Relationships of Argentinean Creole Horses and Other South American and Spanish Breeds Inferred from Mitochondrial DNA Sequences. Anim. Genet. 2002, 33, 356–363. [CrossRef] 15. Stoneking, M.; Sherry, S.T.; Redd, A.J.; Vigilant, L. New Approaches to Dating Suggest a Recent Age for the Human MtDNA Ancestor. Philos. Trans. R. Soc. Lond. B Biol. Sci. 1992, 337, 167–175. [CrossRef] 16. Vigilant, L.; Pennington, R.; Harpending, H.; Kocher, T.D.; Wilson, A.C. Mitochondrial DNA Sequences in Single Hairs from a Southern African Population. Proc. Natl. Acad. Sci. USA 1989, 86, 9350–9354. [CrossRef] [PubMed] 17. Verginelli, F.; Capelli, C.; Coia, V.; Musiani, M.; Falchetti, M.; Ottini, L.; Palmirotta, R.; Tagliacozzo, A.; De Grossi Mazzorin, I.; Mariani-Costantini, R. Mitochondrial DNA from Prehistoric Canids Highlights Relationships between Dogs and South-East European Wolves. Mol. Biol. Evol. 2005, 22, 2541–2551. [CrossRef] [PubMed] 18. Mazzatenta, A.; Carluccio, A.; Robbe, D.; Giulio, C.D.; Cellerino, A. The Companion Dog as a Unique Translational Model for Aging. Semin. Cell. Dev. Biol. 2017, 70, 141–153. [CrossRef] [PubMed] 19. Achilli, A.; Olivieri, A.; Soares, P.; Lancioni, H.; Hooshiar Kashani, B.; Perego, U.A.; Nergadze, S.G.; Carossa, V.; Santagostino, M.; Capomaccio, S.; et al. Mitochondrial Genomes from Modern Horses Reveal the Major Haplogroups That Underwent Domestication. Proc. Natl. Acad. Sci. USA 2012, 109, 2449–2454. [CrossRef] 20. Cothran, E.G.; Juras, R.; Macijauskiene, V. Mitochondrial DNA D-Loop Sequence Variation among 5 Maternal Lines of the Zemaitukai Horse Breed. Genet. Mol. Biol. 2005, 28, 677–681. [CrossRef] 21. Kavar, T.; Dovč, P. Domestication of the Horse: Genetic Relationships between Domestic and Wild Horses. Livest. Sci. 2008, 116, 1–14. [CrossRef] 22. Lippold, S.; Matzke, N.J.; Reissmann, M.; Hofreiter, M. Whole Mitochondrial Genome Sequencing of Domestic Horses Reveals Incorporation of Extensive Wild Horse Diversity during Domestication. BMC Evol. Biol. 2011, 11, 328. [CrossRef] 23. Cieslak, M.; Pruvost, M.; Benecke, N.; Hofreiter, M.; Morales, A.; Reissmann, M.; Ludwig, A. Origin and History of Mitochondrial DNA Lineages in Domestic Horses. PLoS ONE 2010, 5, e15311. [CrossRef] [PubMed] 24. Guastella, A.M.; Zuccaro, A.; Criscione, A.; Marletta, D.; Bordonaro, S. Genetic Analysis of Sicilian Autochthonous Horse Breeds Using Nuclear and Mitochondrial DNA Markers. J. Hered. 2011, 102, 753–758. [CrossRef] [PubMed] 25. Jansen, T.; Forster, P.; Levine, M.A.; Oelke, H.; Hurles, M.; Renfrew, C.; Weber, J.; Olek, K. Mitochondrial DNA and the Origins of the Domestic Horse. Proc. Natl. Acad. Sci. USA 2002, 99, 10905–10910. [CrossRef] [PubMed] 26. Vilà, C.; Leonard, J.A.; Gotherstrom, A.; Marklund, S.; Sandberg, K.; Liden, K.; Wayne, R.K.; Ellegren, H. Widespread Origins of Domestic Horse Lineages. Science 2001, 291, 474–477. [CrossRef] [PubMed] 27. Kivisild, T.; Reidla, M.; Metspalu, E.; Rosa, A.; Brehm, A.; Pennarun, E.; Parik, J.; Geberhiwot, T.; Usanga, E.; Villems, R. Ethiopian Mitochondrial DNA Heritage: Tracking Gene Flow across and around the Gate of Tears. Am. J. Hum. Genet. 2004, 75, 752–770. [CrossRef] 28. Matisoo-Smith, E.; Robins, J. Mitochondrial DNA Evidence for the Spread of Pacific Rats through Oceania. Biol. Invasions 2009, 11, 1521–1527. [CrossRef] 29. Xu, X.; Gullberg, A.; Arnason, U. The Complete Mitochondrial DNA (MtDNA) of the Donkey and MtDNA Comparisons among Four Closely Related Mammalian Species-Pairs. J. Mol. Evol. 1996, 43, 438–446. [CrossRef] 30. Luo, Y.; Chen, Y.; Liu, F.; Jiang, C.; Gao, Y. Mitochondrial Genome Sequence of the Tibetan Wild Ass (Equus Kiang). Mitochondrial DNA 2011, 22, 6–8. [CrossRef] Genes 2021, 12, 1109 13 of 14 31. Ciampolini, R.; Cecchi, F.; Mazzanti, E.; Ciani, E.; Tancredi, M.; De Sanctis, B. The Genetic Variability Analysis of the Amiata Donkey Breed by Molecular Data. Ital. J. Anim. Sci. 2007, 6, 78–80. [CrossRef] 32. George, M., Jr.; Ryder, O.A. Mitochondrial DNA Evolution in the Genus Equus. Mol. Biol. Evol. 1986, 3, 535–546. [CrossRef] [PubMed] 33. Meadow, R.H. Equus africanus in Arabia. In Equids in the Ancient World; Uerpmann, H.P., Ed.; Ludwig Reichert Verlag: Wiesbaden, Germany, 1991; pp. 12–33. 34. Vila, E. Data on equids from late fourth and third millennium sites in Northern Syria. In Equids in Time and Space; Mashkour, M., Ed.; Oxbow Books Ltd.: Oxford, UK, 2002; pp. 101–123. 35. Rossel, S.; Marshall, F.; Peters, J.; Pilgram, T.; Adams, M.D.; O’Connor, D. Domestication of the Donkey: Timing, Processes, and Indicators. Proc. Natl. Acad. Sci. USA 2008, 105, 3715–3720. [CrossRef] 36. Chen, S.Y.; Zhou, F.; Xiao, H.; Sha, T.; Wu, S.F.; Zhang, Y.P. Mitochondrial DNA Diversity and Population Structure of Four Chinese Donkey Breeds. Anim. Genet. 2006, 37, 427–429. [CrossRef] 37. Kefena, E.; Dessie, T.; Tegegne, A.; Beja-Pereira, A.; Yusuf Kurtu, M.; Rosenbom, S.; Han, J.L. Genetic Diversity and Matrilineal Genetic Signature of Native Ethiopian Donkeys (Equus asinus) Inferred from Mitochondrial DNA Sequence Polymorphism. Livest. Sci. 2014, 167, 73–79. [CrossRef] 38. Matassino, D.; Cecchi, F.; Ciani, F.; Incoronato, C.; Occidente, M.; Santoro, L.; Ciampolini, R. Genetic Diversity and Variability in Two Italian Autochthonous Donkey Genetic Types Assessed by Microsatellite Markers. Ital. J. Anim. Sci. 2014, 13, 3028. [CrossRef] 39. Bordonaro, S.; Guastella, A.M.; Criscione, A.; Zuccaro, A.; Marletta, D. Genetic Diversity and Variability in Endangered Pantesco and Two Other Sicilian Donkey Breeds Assessed by Microsatellite Markers. Sci. World J. 2012, 2012, 648427. [CrossRef] 40. Colli, L.; Perrotta, G.; Negrini, R.; Bomba, L.; Bigi, D.; Zambonelli, P.; Verini Supplizi, A.; Liotta, L.; Ajmone-Marsan, P. Detecting Population Structure and Recent Demographic History in Endangered Livestock Breeds: The Case of the Italian Autochthonous Donkeys. Anim. Genet. 2013, 44, 69–78. [CrossRef] [PubMed] 41. Bertolini, F.; Scimone, C.; Geraci, C.; Schiavo, G.; Utzeri, V.J.; Chiofalo, V.; Fontanesi, L. Next Generation Semiconductor Based Sequencing of the Donkey (Equus asinus) Genome Provided Comparative Sequence Data against the Horse Genome and a Few Millions of Single Nucleotide Polymorphisms. PLoS ONE 2015, 10, e0131925. [CrossRef] 42. Cozzi, M.C.; Valiati, P.; Cherchi, R.; Gorla, E.; Prinsen, R.T.M.M.; Longeri, M.; Bagnato, A.; Strillacci, M.G. Mitochondrial DNA Genetic Diversity in Six Italian Donkey Breeds (Equus Asinus). Mitochondrial DNA A DNA Mapp. Seq. Anal. 2018, 29, 409–418. [CrossRef] 43. Aranguren-Mendez, J.; Beja-Pereira, A.; Avellanet, R.; Dzama, K.; Jordana, J. Mitochondrial DNA variation and genetic relation- ships in Spanish donkey breeds (Equus asinus). J. Anim. Breed. Genet. 2004, 121, 319–330. [CrossRef] 44. Tamura, K.; Nei, M.; Kumar, S. Prospects for Inferring Very Large Phylogenies by Using the Neighbor-Joining Method. Proc. Natl. Acad. Sci. USA 2004, 101, 11030–11035. [CrossRef] 45. Tamura, K.; Battistuzzi, F.U.; Billing-Ross, P.; Murillo, O.; Filipski, A.; Kumar, S. Estimating Divergence Times in Large Molecular Phylogenies. Proc. Natl. Acad. Sci. USA 2012, 109, 19333–19338. [CrossRef] 46. Kumar, S.; Stecher, G.; Tamura, K. MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets. Mol. Biol. Evol. 2016, 33, 1870–1874. [CrossRef] 47. Letunic, I.; Bork, P. Interactive Tree Of Life (iTOL) v5: An online tool for phylogenetic tree display and annotation. Nucleic Acids Res. 2021, 49, 1–4. [CrossRef] 48. Cecchi, F.; Ciampolini, R.; Ciani, E.; Matteoli, B.; Mazzanti, E.; Tancredi, M.; Presciuttini, S. Demographic Genetics of the Endangered Amiata Donkey Breed. Ital. J. Anim. Sci. 2006, 5, 387–391. [CrossRef] 49. Navas, F.J.; Jordana, J.; León, J.M.; Barba, C.; Delgado, J.V. A Model to Infer the Demographic Structure Evolution of Endangered Donkey Populations. Animal 2017, 11, 2129–2138. [CrossRef] 50. Cinar Kul, B.; Bilgen, N.; Akyuz, B.; Ertugrul, O. Molecular Phylogeny of Anatolian and Cypriot Donkey Populations Based on Mitochondrial DNA and Y-Chromosomal STRs. Ank. Univ. Vet. Fak. Derg. 2016, 63, 143–149. 51. Gutiérrez, J.P.; Marmi, J.; Goyache, F.; Jordana, J. Pedigree Information Reveals Moderate to High Levels of Inbreeding and a Weak Population Structure in the Endangered Catalonian Donkey Breed. J. Anim. Breed. Genet. 2005, 122, 378–386. [CrossRef] 52. Ivanković, A.; Ramljak, J.; Konjačić, M.; Kelava, N.; Dovč, P.; Mijić, P. Mitochondrial D-Loop Sequence Variation among Autochthonous Horse Breeds in Croatia. Czech J. Anim. Sci. 2009, 54, 101–111. [CrossRef] 53. Stanisic, L.J.; Aleksic, J.M.; Dimitrijevic, V.; Simeunovic, P.; Glavinic, U.; Stevanovic, J.; Stanimirovic, Z. New Insights into the Origin and the Genetic Status of the Balkan Donkey from Serbia. Anim. Genet. 2017, 48, 580–590. [CrossRef] 54. Aranguren-Méndez, J.; Jordana, J.; Gomez, M. Genetic Diversity in Spanish Donkey Breeds Using Microsatellite DNA Markers. Genet. Sel. Evol. 2001, 33, 433. [CrossRef] 55. Ivankovic, A.; Kavar, T.; Caput, P.; Mioc, B.; Pavic, V.; Dovc, P. Genetic Diversity of Three Donkey Populations in the Croatian Coastal Region. Anim. Genet. 2002, 33, 169–177. [CrossRef] 56. Pérez-Pardal, L.; Grizelj, J.; Traoré, A.; Cubric-Curik, V.; Arsenos, G.; Dovenski, T.; Marković, B.; Fernández, I.; Cuervo, M.; Álvarez, I.; et al. Lack of Mitochondrial DNA Structure in Balkan Donkey Is Consistent with a Quick Spread of the Species after Domestication. Anim. Genet. 2014, 45, 144–147. [CrossRef] Genes 2021, 12, 1109 14 of 14 57. Bjørnstad, G.; Røed, K.H. Evaluation of Factors Affecting Individual Assignment Precision Using Microsatellite Data from Horse Breeds and Simulated Breed Crosses. Anim. Genet. 2002, 33, 264–270. [CrossRef] 58. Toro, M.A.; Barragán, C.; Ovilo, C. Estimation of Genetic Variability of the Founder Population in a Conservation Scheme Using Microsatellites. Anim. Genet. 2003, 34, 226–228. [CrossRef]