Phylogenetic Relationships of Turkish Indigenous Donkey Populations Determined by Mitochondrial DNA D-loop Region

Simple Summary This paper represents the first fundamental report of mtDNA diversity in Turkish indigenous donkey breeds and presents findings for the origin and genetic characterization of donkey populations dispersed in seven geographical regions in Turkey, and thus reveals insights into their genetic history. The median-joining network and phylogenetic tree exhibit two different maternal lineages of the 16 Turkish indigenous donkey populations. Abstract In this study, to analyze the mtDNA D-loop region and the origin of the maternal lineages of 16 different donkey populations, and to assess the domestication of Turkish indigenous donkeys in seven geographical regions, we investigated the DNA sequences of the D-loop region of 315 indigenous donkeys from Turkey. A total of 54 haplotypes, resulting from 35 polymorphic regions (27 parsimoniously informative and 6 singleton sites), were defined. Twenty-eight of these haplotypes are unique (51.85%), and 26 are shared among different Turkish indigenous donkey populations. The most frequent haplotype was Hap 1 (45.71%), followed by two haplotypes (Hap 4, 15.55% and Hap 7, 5.39%). The breed genetic diversity, evaluated by the haplotype diversity (HD) and nucleotide diversity (πD), for the Turkish donkey populations ranged from 0.533 ± 0.180 (Tekirdağ–Malkara, MAL) to 0.933 ± 0.122 (Aydin, AYD), and from 0.01196 ± 0.0026 (Antalya, ANT) to 0.02101 ± 0.0041 (Aydin, AYD), respectively. We observed moderate-to-high levels of haplotype diversity and moderate nucleotide diversity, indicating plentiful genetic diversity in all of the Turkish indigenous donkey populations. Phylogenetic analysis (NJT) and median-joining network analysis established that all haplotypes were distinctly grouped into two major haplogroups. The results of AMOVA analyses, based on geographic structuring of Turkish native donkey populations, highlighted that the majority of the observed variance is due to differences among samples within populations. The observed differences between groups were found to be statistically significant. Comparison among Turkish indigenous donkey mtDNA D-loop regions and haplotypes, and different countries’ donkey breeds and wild asses, identified two clades and which is named Somali (Clade IV) and Nubian (Clade V) lineages. The results can be used to understand the origin of Turkish donkey populations clearly, and to resolve the phylogenetic relationship among all of the different regions.

Over the last 30 years, as farm mechanization and traffic automation have become popular, the number of donkeys has reduced strikingly in the world [14].
The objective of this study was: (i) to determine the genetic variability between seven different geographical regional donkey populations in Turkey, (ii) to investigate the genetic structure of Turkish native donkey populations and to explore the maternal origin, with mtDNA D-loop sequences, and explore the domestic history of Turkish indigenous donkey populations, (iii) to evaluate genetic breed sub-structuring, and (iv) the comparison of mitochondrial DNA sequences of Turkish populations with other countries' domestic donkey breeds and wild asses. The results will be used to better understand the Turkish donkey population's origin, and to resolve the phylogenetic relationship among all the different regions.

Population Determination, Sample Collection, and DNA Isolation
According to the FAO database (Domestic Animal Diversity Information System, accessed 15 September 2019); the estimated number of donkeys in Turkey is about 150,000. Three types of donkey are recognized in Turkey (FAO/DADIS 2012). The Anatolian is found all over the country, not just on the central plateau, and is usually grey or black in color. The Merzifon or Marsovan, from the town and district of the same name in Amasya province in the central Black Sea region, is at risk. The Karakaçan type, presumably named after the Karakaçan nomads of the Balkans and Thrace, is also considered at risk (FAO-DADIS 2012) [20].
In this study a total of 315 blood samples were collected from Turkish indigenous donkeys raised in seven different geographical regions, from 16 different provinces of Turkey ( Figure 1). These locations were selected to represent the expected native Turkish donkey breeds: Anatolian, Merzifon, Karakaçan, etc. The provinces and the potential geographical dispersion of some Turkish indigenous breeds, with the number of animals, are given in Table 1. These provinces of Turkey were known for their relatively large donkey populations.
The objective of this study was: (i) to determine the genetic variability between seven different geographical regional donkey populations in Turkey, (ii) to investigate the genetic structure of Turkish native donkey populations and to explore the maternal origin, with mtDNA D-loop sequences, and explore the domestic history of Turkish indigenous donkey populations, (iii) to evaluate genetic breed sub-structuring, and (iv) the comparison of mitochondrial DNA sequences of Turkish populations with other countries' domestic donkey breeds and wild asses. The results will be used to better understand the Turkish donkey population's origin, and to resolve the phylogenetic relationship among all the different regions.

Population Determination, Sample Collection, and DNA Isolation
According to the FAO database (Domestic Animal Diversity Information System, accessed 15 September 2019); the estimated number of donkeys in Turkey is about 150,000. Three types of donkey are recognized in Turkey (FAO/DADIS 2012). The Anatolian is found all over the country, not just on the central plateau, and is usually grey or black in color. The Merzifon or Marsovan, from the town and district of the same name in Amasya province in the central Black Sea region, is at risk. The Karakaçan type, presumably named after the Karakaçan nomads of the Balkans and Thrace, is also considered at risk (FAO-DADIS 2012) [20].
In this study a total of 315 blood samples were collected from Turkish indigenous donkeys raised in seven different geographical regions, from 16 different provinces of Turkey ( Figure 1). These locations were selected to represent the expected native Turkish donkey breeds: Anatolian, Merzifon, Karakaçan, etc. The provinces and the potential geographical dispersion of some Turkish indigenous breeds, with the number of animals, are given in Table 1. These provinces of Turkey were known for their relatively large donkey populations. From those donkey samples, 10 mL of whole blood was taken from the jugular vein with EDTAcoated vacutainer tubes. Genomic DNA was isolated by using a standard phenol-chloroform-isoamyl alcohol extraction method [29]. The concentration of DNA was analyzed to compare with the  From those donkey samples, 10 mL of whole blood was taken from the jugular vein with EDTA-coated vacutainer tubes. Genomic DNA was isolated by using a standard phenolchloroform-isoamyl alcohol extraction method [29]. The concentration of DNA was analyzed to compare with the standard DNA marker concentration on agarose gels. The quality and quantity of DNA was checked on 0.6% agarose gels, prepared with a Tris-Boric acid-EDTA buffer. The isolated DNA was diluted with TE buffer (10:1) and stored at +4 • C till polymerase chain reaction (PCR) analysis.

Amplification of mtDNA Control Regions (D-Loop) and Sequencing
The 383 bp segment, from the mitochondrial DNA control region (D-loop), from the donkey genome, was amplified by PCR, as described by Aranguren-Méndez et al. [4]. The primer sequences were amplified between the sites 15,387 to 15,769 bp. The PCR conditions were optimized for all the primer pairs selected for the study. mtDNA D-loop was amplified by PCR, carried out in 20 µL reaction volume, containing 50 ng of genomic DNA, 2.0 mM MgCl 2 , 0.2 mM each of dNTP, 0.5 µM of each primer, 1 X PCR buffer, and 1 U of Taq DNA polymerase. Amplification was performed using a BIORAD MyCycler Thermal Cycler (Bio-Rad Laboratories, Inc, Berkeley, CA, USA) with the following conditions: initial denaturation at 94 • C for 5 min followed by 35 cycles of denaturation at 94 • C for 30 s, annealing at 58 • C for 30 s, extension at 72 • C for 2 min, and final extension at 72 • C for 15 min. The amplified products were electrophoresed, visualized under UV radiations, and purified to be sequenced using a ABI prism 3500 genetic analyzer. PCR products were sequenced for both directions. The raw sequence trace files were checked for the presence of ambiguous bases using ChromasPro v. 1.7.4 software.

Data Analyses
The obtained sequences were controlled for the presence of ambiguous bases by using the software, ChromasPro v. 1.7.4. For identifying the haplotypes (Hap), all sequences were aligned with the complete donkey mtDNA sequences, MK896308 [30] and X97337 [11], commonly used as reference sequences, using the ClustalW algorithm implemented in BioEdit v7.0.9 [31] and MEGA v.7.0.26 software [32]. The number of haplotypes (N H ), number of private haplotypes (P H ), haplotype diversity (H D ), total polymorphic sites (S), parsimony informative sites (SPI), singleton site (S S ), shared haplotype (S H ), the nucleotide diversity (π D ), and average number of nucleotide differences (k) were estimated using DNASP 6.12.03 ×64 software [33]. Pairwise population genetic differentiation (F ST ) and the average number of pairwise differences within and between populations were calculated using Arlequin v.3.5.2.2 [34]. In order to define the genetic relationships between Turkish indigenous donkey populations mtDNA sequences and those available in Genbank for other donkey populations originated from European, Asian, and African origin domestic donkeys, as well as Ethiopian donkeys, and wild asses, a total of 429 sequences were considered. All the sequences were trimmed at 344 bp to maximize the number of sequences included. The mtDNA sequences from the GenBank database used for haplogroup definition and relationship investigation are reported in Table S1. To evaluate the genetic variance among donkeys within and between populations, we performed the analysis of molecular variance (AMOVA) using Arlequin v.3.5.2.2 [34]. MEGA v.7.0.26 software was used to analyze the relationships among the haplotypes identified in Turkish native donkeys [32].
The T-REX web server was used to analyze the relationships among the haplotypes identified in Turkish native donkey populations by geographical regions. Genetic relationships among populations were reconstructed using median-joining networks (MJN) using Network v.10.1.0.0 software [35]. Whole mtDNA D-loop sequences identified in this research were submitted to the GenBank with the accession numbers MH683672-MH683725.

Genetic Diversity of the mtDNA D-Loop Analyses of Turkish Indigenous Donkeys
The mtDNA D-loop sequences were taken for all analyzed 315 samples from seven different regions, and from 16 different provinces of Turkey. In Table 2 Table S1. The haplotype distribution, by clusters corresponding to the Turkish indigenous donkeys, is shown in Figure 2b. Construction of the phylogenetic tree by the maximum likelihood method showed two different clusters. The maximum composite likelihood estimates of the nucleotide variation pattern, determined for the 315 sequences, were 33.00% (A), 27.3% (T), 28.6% (C), and 11.1% (G). Positions containing gaps and missing data were eliminated.   In Table 3, the observed haplotypes, their annotation, and frequencies in studied populations are presented.  (21), Marmara (19), Mediterranean (18), and Black Sea (17) regions. Southeast Anatolia populations have higher mtDNA diversity with 21 haplotypes. Eleven of these haplotypes were observed only in Southeast Anatolian populations. It was observed that the lowest haplotype diversity is in Central Anatolia and Eastern Anatolia. It is thought that the reason for this decrease was due to the low number of samples, and the sampling of only one province in these regions. The haplotype number in each population varies from 4 to 15. The highest haplotype number was determined in Mardin (MRD) donkeys followed by Amasya-Merzifon (MER), Kırklareli (KIR), and Antalya (ANT) donkey populations. The lowest haplotype number was determined in (: Tekirdag-Malkara (MAL), Kütahya (KUT), Kahramanmaraş (KRM)) donkeys followed by Tokat (TOK) and Aydın (AYD) donkey populations.   The breed genetic diversity (Table 2), evaluated by the haplotype diversity (H D ) and nucleotide diversity (π D ), for the Turkish indigenous donkey populations ranged from 0.533 ± 0.180 (MAL)-0.933 ± 0.122 (AYD), and 0.01196 ± 0.0026 (ANT)-0.02101 ± 0.0041 (AYD), respectively. The overall values for haplotype and nucleotide diversities were 0.763 ± 0.023 and 0.01656 ± 0.00046, respectively ( Table 1). Analysis of mtDNA D-loop region for each population revealed that the AYD donkey population had the highest haplotype diversity followed by the KAS donkey population. In MAL donkeys, the lowest haplotype diversity was obtained, followed by KAR and SAN donkey populations. AYD donkeys had the highest estimates of nucleotide diversity followed by the ISP donkey population. In general, nucleotide diversities ranged from 0.01196 in ANT donkeys to 0.02101 in AYD donkeys. The highest number of private haplotypes (n = 8) was observed in the KRM population. Additionally, the average number of nucleotide variation (k) ranged from 4.172 (ANT) to 7.333 (AYD) with an average of 5.745.
The 54 haplotypes were defined by 29 parsimony informative sites. Parsimony informative sites [32] are defined as mutations that have a minimum of two nucleotides that are present at least twice in the sampled population, whereas noninformative sites are singleton sites. Of the polymorphic sites, there were 35 transitions and two transversions, and two sites had both a transition and transversion, suggesting a strong bias toward transitions. No insertions or deletions were observed. At a population level, all the breeds showed one or more singleton sites. In Table 2, the CAT population counted the highest number of non-informative singleton sites (S S ), while the lowest number of polymorphic sites (S) was identified in MAL and KRM populations. In addition, the mean number of nucleotide differences (k) was the highest for AYD donkeys (7.333) followed by ISP donkeys (6.689), while ANT donkeys had the lowest mean number of nucleotide differences (4.172). The AYD breed observed five haplotypes, the highest π value (0.02101), the highest number of haplotype diversity (0.933), and the highest nucleotide diversity (k = 7.333). The high k value for the AYD breed is due to the high difference among the five haplotypes identified in the population. Despite the low number of samples, a high H D and (π) value, with five haplotypes (Hap 1, Hap 4, Hap 7, Hap 14, and Hap 25) for the AYD was found (Table 1). However, ANT have the lowest π D value (0.01196) and k value (k = 4.172), but haplotypes have a high H D value (0.717).
The matrix of pairwise F ST (fixation index) value and D A distances value between Turkish indigenous donkey populations from seven geographical regions are shown in Table 4. The F ST values range from 0.0302 (Aegean vs. Central Anatolia) to 0.2523 (Marmara vs. Central Anatolia). The highest differences in F ST values were found between the Marmara and Central Anatolia region populations (F ST : 0.2523, p < 0.01), followed by between the Mediterranean and Central Anatolia region donkey populations (F ST : 0.2130, p < 0.01). In contrast, there is little genetic differentiation between Southeastern Anatolia and Black Sea donkeys (F ST :0.0004), followed by between Mediterranean and Black Sea (F ST :0.0068). In the study, some F ST values were observed to be negative (Table 4). Negative F ST values should be effectively seen as zero values. A zero value for F ST means that there is no genetic subdivision between the populations considered. In addition, the matrix of pairwise F ST values between Turkish indigenous donkey populations from 16 different provinces is revealed in Figure 3.  Table S2 Table S3) for 315 tested sequences from different regions was constructed using the T-REX web server. The NJT clearly demonstrated three distinct clades in Turkish donkey populations from seven geographical regions ( Figure 4). Three geographical regions, including CAR, AER, and SAR, clustered in clade I, while a different three geographical region populations, including MDR, MRM, and EAR, clustered in clade II. Clade III included only the Black Sea region populations (BSR). The topology of this tree reflected the patterns of Reynolds' linearized pairwise genetic distances, where MRM and MDR, MRM and EAR, BSR and EAR, CAR and AER, MDR and EAR, as well as BSR and SAR region donkeys were genetically much closer to each other than to the rest of the pairwise genetic distances (F ST : Supplementary Table S2).  The median-joining network of mtDNA D-loop region was described between Turkish donkey populations ( Figure 5). The size of the provided node is proportional to the number of samples represented in a haplotype, with the smallest node representing a single individual. Branch length is proportional to the mutational distance; only mutational distances greater than one are indicated. In   The median-joining network of mtDNA D-loop region was described between Turkish donkey populations ( Figure 5). The size of the provided node is proportional to the number of samples represented in a haplotype, with the smallest node representing a single individual. Branch length is proportional to the mutational distance; only mutational distances greater than one are indicated. In  The median-joining network of mtDNA D-loop region was described between Turkish donkey populations ( Figure 5). The size of the provided node is proportional to the number of samples represented in a haplotype, with the smallest node representing a single individual. Branch length is proportional to the mutational distance; only mutational distances greater than one are indicated. In Figure 5, two major clades were obtained from the reconstructed median joining network. The clade  Figure 6 show that Clade 1 and Clade 2 distribution more clearly.

Turkish Donkey Populations Comparing with Other Donkey Populations
Comparison of the 427 reference (the GenBank database samples, Supplementary Table S1) Table 5.

Turkish Donkey Populations Comparing with Other Donkey Populations
Comparison of the 427 reference (the GenBank database samples, Supplementary Table S1)   The phylogenetic relationship among the 145 identified haplotypes was analyzed through the median-joining network. The MJN showed that there are five distinct lineages, as indicated in the phylogenetic tree and the star-like phylogeny (Figure 7a). The haplotypes of each population are detected by the color code, the abundance by the relative size of the symbol, and the diffusion among the breeds with the pie division of the different colors (Figure 7b,c).   [16], Green: [26], Yellow: [15], Orange: [30], Purple: [17], White: [11], Turquoise: [28], Grey: [38] Four haplotypes were determined in E. asinus somalicus, and were grouped in the same (III) lineage. Lineages IV and V consist of domestic donkey populations. In Figure 6, Lineages IV and V, which have two different clades, are obviously defined in domestic donkey populations: clade 1 contains 50.68% of sequences, whereas clade 2 contains 49.31% of sequences in domestic donkey populations. Table 6 shows that genetic distances and nucleotide diversity of the two clades. The AMOVA analysis is a useful tool to check how maternal genetic diversity is distributed within, and among, populations whose structure is quantified by FST values. In this study, we checked possible structures by creating and comparing different groups of populations. The analysis was  [16], Green: [26], Yellow: [15], Orange: [30], Purple: [17], White: [11], Turquoise: [28], Grey: [38]).
The MJN analysis detected 145 haplotypes from 742 sequences (427 deposited GenBank database sequences (Supplementary Table S1) and 315 Turkish donkeys) showing a high variability (H D = 0.890 ± 0.009). The median joining network showed that there are five distinct lineages, as shown in the phylogenetic tree and the star-like phylogeny. Three samples, Hap 126-127 and Hap 63, were horses (E. cabullus) that grouped together as an outgroup to the donkeys in the network. Equus cabullus, Equus kiang, and Asiatic wild asses (E. hemionus luteus, E. hemionus kulan, E. Hemionus onager) shared the (I) and (II) lineage, and had their respective lineages. Lineage (I) (E. cabullus) has three haplotypes, lineage (II) (E. hemionus luteus, E. hemionus kulan, E. Hemionus onager, E. kiang) has eight haplotypes. Four haplotypes were determined in E. asinus somalicus, and were grouped in the same (III) lineage. Lineages IV and V consist of domestic donkey populations. In Figure 6, Lineages IV and V, which have two different clades, are obviously defined in domestic donkey populations: clade 1 contains 50.68% of sequences, whereas clade 2 contains 49.31% of sequences in domestic donkey populations. Table 6 shows that genetic distances and nucleotide diversity of the two clades. Table 6. Genetic distance and nucleotide diversity of the two clades.

Population
Pairwise Distance ± SD π D ± SD The AMOVA analysis is a useful tool to check how maternal genetic diversity is distributed within, and among, populations whose structure is quantified by F ST values. In this study, we checked possible structures by creating and comparing different groups of populations. The analysis was made with two hypotheses: Hypothesis 1; the geographical distribution of the countries was considered. The domestic donkey populations were clustered in five groups: (1) Turkey's indigenous donkey populations (n = 315), (2) European countries data: Balkan donkeys [26], Italian donkeys [28], Serbian donkeys [15] (n = 112), (3) Asian countries data: China, Tadzhikistan, Kyrgyzstan, Iran, Mongolia donkeys (n = 219) [17,30], (4) African origin domestic donkeys, and Kenyan donkeys, (n = 85) [16], (5) wild asses and others (n = 11) [11,[36][37][38][39][40][41]. In hypothesis 2, the F ST values and genetic distances (Reynolds' linearized genetic distances in Turkish donkey populations from different regions) were considered between Turkey's domestic donkey populations. (1) CAR, AER, and SAR region donkeys; (2) MDR, MRM, and EAR region donkeys, and (3) BSR region donkeys. Table 7 shows the results of the AMOVA analysis according to these hypotheses. Φ ST : 0.07705 0.000 *** a Φ CT : variation among groups divided by total variation; Φ SC : variation among sub-groups divided by the sum of variation among sub-groups within groups and variation within sub-groups; Φ ST : the sum of variation groups divided by total variation. b ns = p > 0.05; * p < 0.05; *** p< 0.001. Hypothesis 1; the AMOVA analyses results revealed that the variation among groups, among populations within groups, and within populations were 7.95%, 1.72%, and 90.34%, respectively. The findings showed that most of the obtained variance is due to differences within populations. Variance components among groups were nonsignificant, demonstrating non-significant geographical distribution in the analyzed donkey populations. Moreover, the variance component among populations, within groups (p < 0.05), and within populations (p < 0.001) were significant. Hypothesis 2; the AMOVA analysis findings indicated that the differentiation among groups, among populations within groups, and within populations were 7.31%, 0.39%, and 92.29%, respectively. The findings revealed that most of the obtained variance is based on differences within populations. Variance components among groups were significant (p < 0.05), which showed important geographical distribution in the analyzed Turkish donkey populations by using Reynolds' linearized genetic distances. In addition, variance components among populations within groups were non-significant statistically, but were statistically significant within populations (p < 0.001).

Discussion
As exposed in some of genetic diversity studies in donkeys [13,19,24,25,30,42,43], the D-loop region of maternally ancestral mtDNA supplies adequate evidence to assess population genetic variation, evolutionary relationships, and matrilineal genetic origin of a species under consideration. This study presents the substantial analysis of mtDNA diversity in Turkish indigenous donkeys and supplies knowledge about the maternal lineage origins and genetic diversity of donkey populations, and so understanding into their genetic history.

Phylogenetic Relationships of Turkish Indigenous Donkey Populations
In the present study, mtDNA D-loop region from 315 Turkish indigenous donkeys were identified to clarify their phylogenetic relationship and haplogroups, and to define mtDNA haplotypes and their maternal ancestry. The results showed quite high genetic diversity within and between 16 Turkish donkey populations from the seven different geographical regions. By sequence alignment, 54 distinct haplotypes were defined among Turkish donkey populations, 28 of which are unique and 26 are shared among different donkey populations. No congruence of haplogroup to a population's geographic ancestry was observed. Similarly to other research, the content of A+T was the richest in the mtDNA D-loop in Turkish donkey populations [28,44]. The haplotype distribution by clusters, corresponding to the Turkish native donkeys, showed two different clusters by using the maximum likelihood method (Figure 2b). These results were similar to Çınar Kul et al. [27] and Yalçın's [45] studies.
The MJN analysis retrieved two major clades among the 16 different Turkish donkey populations analyzed in the present study. Two clades, separated from five nucleotide substitutions each, were found. These clades comprise haplotypes for which common origins are assumed since they share a characteristic pattern of mutations. As revealed in the network, it clearly showed two lineages and indicated a star like phylogenetic pattern, in which two large haplotypes, Hap 1 and Hap 4, are in the center of clade 1 and clade 2 lineages, respectively. The clade 2 lineage predominated slightly (60.635%, 191/315). D-loop gene lineages in Turkish donkeys appear mixed with those of other donkeys.
The haplotype and nucleotide diversity within Turkish native donkeys were (0.763 ± 0.023) and (0.01680 ± 0.00046), respectively. The haplotype and nucleotide diversity values of other countries' donkeys ranged from (0.890 ± 0.013), (0.01759 ± 0.00053) in the Chinese donkeys to (0.983 ± 0.007), and (0.02843 ± 0.00195) in the African origin domestic donkeys, respectively. The haplotype and nucleotide diversity values of the Turkish donkey populations was found to be lower than that of the other donkey populations (Table 5; Balkan, Chinese, Serbian, Italian etc.), indicating a relatively low diversity. This finding is consistent with other genetic diversity studies in different countries' donkey populations [15,26,28]. The haplotype diversity found here was higher than that found in a previous study in three different Turkish donkey populations [27], and the nucleotide diversity found here was lower than that found in that previous study. These results indicate a relatively higher level of genetic diversity in the 16 Turkish donkey populations compared with other countries' donkey populations. For example, the haplotype diversity and nucleotide diversity values of 10 Balkan donkey populations were 0.982 ± 0.002 and 0.017 ± 0.009 [26]. However, according to Walsh's work on the required sample size for the diagnosis of conservation units a sample of 59 individuals is necessary to reject the hypothesis that individuals with unstamped ("hidden") character states exist in the population size. Thus, the sample size necessary to reject a hidden state frequency of 0.05 is 56, when sampling from a finite population of 500 individuals [48]. Our genetic diversity estimation is therefore a precise reflection of Turkish indigenous donkeys due to the large sample size used in this study.
The MJN analysis identified 145 haplotypes from 742 mtDNA D-loops (427 Genbank sequences and our 315 Turkish donkeys) based on 344 bp of the control region sequence, and all these haplotypes grouped into five lineages.  (11,(15)(16)(17)26,28,30,37,38). The MJ network analysis showed that both Somali and Nubian lineages had made a genetic contribution to Turkish indigenous donkey evolution, as 19 haplotypes of lineage Somali (n = 191 samples) and 30 haplotypes of lineage Nubian (n = 124) were identified in Turkish donkey populations. These findings showed that 16 different Turkish native donkey populations and other domestic donkey populations possessed abundant mtDNA diversity, and indicated two different maternal origins. Clade 1 includes 49.59% of sequences, whereas clade 2 includes 48.52% of sequences. Out of 145 haplotypes, 53 were referred to the Somali lineage (clade I), whereas 77 haplotypes belonged to Nubian lineage (clade II). As shown in the network, two lineages were clearly revealed, and showed a star like phylogenetic pattern, in which two large haplotypes, H2 and H5, are located in the center of the clade 1 and clade 2 lineages, respectively. The clade 1 lineage appeared to predominate slightly (49.59%, 368/742) in the domestic donkeys. The results showed that the Nubian (E. africanus africanus) and Somali (E. asinus somali) wild asses were the probable progenitor for ancient Turkish donkeys and domestic donkeys. These results support previous studies on the origins of the domestic donkey [15][16][17]26,28,37].
In our study, to assess a hierarchical structure among Turkish native donkey populations and other countries' donkeys (Italian, Serbian, Chinese, etc.) and wild asses, AMOVA analysis was performed under two different hypotheses, grouping populations in different clusters. In both hypotheses, the AMOVA analysis results indicated that the majority of the observed variance was due to differences among individuals within populations. The greater part of the variation was observed within individuals (90.34% hypothesis 1 and 92.29% hypothesis 2) whereas the differences among groups represented only the 7.95% (hypothesis 1) and 7.31% (hypothesis 2) of the variation, respectively. Our results are similar to other studies on the mtDNA D-loop region [15,26,28,42]. In hypothesis 1; among populations within groups and among individuals within populations there was statistically significant variation (p < 0.001, p < 0.05). In hypothesis 2; among groups and among individuals within populations there was also statistically significant variation (p < 0.001, p < 0.05). Our AMOVA analysis and phylogenetic analysis revealed low genetic variation between the groups that were defined geographically among Turkish donkey populations. It was observed that there is a small difference in the donkey population of the BSR compared to the other region's populations. In order to explain the reason for this difference, it is considered appropriate to conduct new studies with more samples in BSR.

Conclusions
The present study demonstrated the abundant mtDNA diversity existing in Turkish indigenous donkey populations. The detection of 54 haplotypes in 315 donkeys suggests that abundant genetic variety exists in Turkish donkey populations. We confirmed two different maternal lineages (Somali and Nubian) of domestic donkeys reported by other researchers. No obvious geographical distribution was found among Turkish donkey populations. The study found high nucleotide and haplotype diversity values, no haplotypes clearly distinct from other populations, and no clear clustering on the median joining tree in Turkish donkeys. In summary, our results imply the moderately high genetic diversity of the 16 Turkish donkey populations in mtDNA D-loop region, and give an insight into the origin of the analyzed populations. Our results obviously exclude the Asian wild donkey as ancestors of Turkish indigenous donkey populations, and the two African wild donkeys (Somali and Nubian) are the likely ancestors of Turkish indigenous donkeys. The results provided here could be regarded as a genetic structure of the maternal ancestors of Turkish donkey populations.