Banat donkey, a neglected donkey breed from the central Balkans (Serbia)

Banat donkey population in Serbia, sample size and data collection

According to the Domestic Animal Diversity Information System (DAD-IS) of the FAO (accessed May 25, 2019), the estimated number of donkeys in Serbia is 500–1,000. A small number of adult donkeys and foals is kept in individual households usually in rural regions, while the number of sexually mature individuals registered in 2015 at three localities where large donkey herds are currently maintained, Special Nature Reserve “Zasavica” (ZA), Stara Planina mountain (SP), and in vicinity of the Kovilj village (KO), is 281. These animals are classified to the endangered Balkan donkey breed found in the Balkan Peninsula (Kugler, Grunenfelder & Broxham, 2008; Stanisic et al., 2015; Stanisic et al., 2017). In spring 2017, we visited these localities and recorded all individuals that differ morphologically from the Balkan donkey and, according to local breeders, belong to the Banat donkey.

Since the controlled mating between Banat donkeys over the past 25 years was more or less common only in ZA, we used animals from this locality in our study. These individuals have already been studied genetically by Stanisic et al. (2017), and were characterized by the unique nuclear gene pool and prevalence of the mtDNA Clade 2 haplotypes. However, after re-capturing them, we observed that only eight out of 16 animals were indeed morphologically different from the Balkan donkey (Figs. 1A and 1B). The remaining ones, although assigned to the Banat donkey based on their nuclear genetic profiles, were visually indistinguishable from the Balkan donkey because they were mid-sized animals, with grey to chocolate coat color, dark cross on shoulder and silver bright coat pigmentation on the chest and abdomen, along the medial side of the extremities and around the nose and the eyes (Fig. 1C). Thus, these animals were delineated as potential hybrids (HY) and were separated from eight individuals representing Banat donkeys (BanD). All but one BanD individual were adult females, three to eight years old (average age 5.25 years). The average age of HY, which were all females, was 5.75 years.

Eighteen morphological traits described by Stanisic et al. (2015) were measured in animals used in the present study (Fig. S1) approved by the Ethical Committee of the Faculty of Veterinary Medicine, University of Belgrade, Serbia (approval reference number 01-19/7). All re-captured individuals were carefully handled by trained laborers, positioned on flat and hard grounds with parallel legs, and calmed before the non-aggressive measurement was carried out. A single person recorded all the measurements taken from the right side of the body following the procedure used previously by Stanisic et al. (2015) for assessing body measurements in the Balkan donkey. That is, six body variables (carpal circumference, chest circumference, ear length, head length, tarsal circumference and tibia circumference) were measured with a measuring tape, and 11 (back height, body length, carpal height, chest depth, chest width, croup length, croup width, head width, hip height, tarsal height and wither height) with the Lydtin stick. The body weight (BW) of each animal was calculated following the formula given by Pejić (1996): BW (kg) = [(chest circumference)² × (body length)]/11877.

We analysed morphological traits of BanD and HY together with those available for 53 female Balkan donkeys from Serbia whose morphological features have been studied previously by Stanisic et al. (2015). These animals were studied genetically as well (Stanisic et al., 2017), and were selected among 77 individuals based on the assignment probability (q ≥ 0.60) to one of the two nuclear gene pools obtained in the Bayesian clustering analysis performed with STRUCTURE v. 2.3 (Pritchard, Stephens & Donnelly, 2000). Two gene pools in question, so-called blue gene pool (BGP) and red gene pool (RGP), represent two sub-populations of the heterogeneous Balkan donkey breed (Stanisic et al., 2017). We therefore assembled two groups using selected individuals, BalkD-BGP (average age 5.27 years) and BalkD-RGP (average age 5.03 years), regardless on the origin of individuals (ZA, SP and KO). In that way, we used 26 females belonging to the BGP (14 from ZA, nine from SP, and three from KO), and 27 females belonging to the RGP (five from ZA, five from SP and 17 from KO).

Statistical analyses

For each of the four assembled donkey groups (BanD, HY, BalkD-BGP and BalkD-RGP), mean, median and quartile range of each morphological variable (expressed in cm or kg) as well as coefficient of variation (CV, expressed in %) were assessed with STATISTICA 13 (TIBCO, Palo Alto, CA, USA).

Multivariate analysis of variance (MANOVA) was used to test the impact of differences between studied groups on variation of morphological parameters. The significance of tests was assessed with both Hotelling-Lawley and Wilks approximate F-tests. Hotelling’s T2 test was used for testing the significance of differences between donkey groups based on all measured parameters. Then, one-way Analysis of Variance (ANOVA) was applied to test the null hypotheses about influence of differences between groups on variation of each character. Statistical significance of differences between groups based on ANOVA was assessed using the Tukey’s test (HSD-test). These analyses were performed with R programming language (R Core Team, 2017), and ‘Hotelling’ package (Curran, 2018) implemented in R for Hotelling’s T2 test.

To further describe the importance of the influence of differences between donkey groups on variability on studied parameters, we calculated contribution of variances to the total expected variances (intraclass correlation coefficient). Expected variances were calculated using STATISTICA 13. The same software was used for Principal component analysis (PCA) performed in order to group measured parameters according to their loadings (i.e., correlation coefficients) with selected and rotated principal components (PC). The most informative PCs were selected according to the Kaiser’s role (λ > 1, where λ stands for eigenvalue of particular PC), and then rotated by the Varimax method to maximize variance of loadings between measured parameters and PCs within PCs. Since PCs are in orthogonal position against each other, the correlation between them is 0, and thus, the measured parameters that have their highest loadings with the same PC were considered to be correlated and to belong to the same PCA-group, with weak correlations with parameters from other PCA-groups. Relationships between studied groups were assessed in cluster analysis by calculating a matrix of squared Mahalanobis distances between groups used to build a dendrogram with unweighted pair-group method with arithmetic mean (UPGMA), and by Canonical Discriminant Analysis (CDA). Also, Stepwise Forward Discriminant Analysis was used to evaluate the importance of studied morphological parameters (predictors) for discrimination of examined donkey groups. For each successive discriminative model, the percentage of correct allocation is calculated. Importance of parameters was assessed according to order by which parameters entered the discriminative model and the difference in the percentage of correct allocation between the model with and the model without included parameter. These analyses were performed with STATISTICA 13.

Genetic profiles of individuals belonging to four groups (BanD, HY, BalkD-BGP and BalkD-RGP), assessed previously by Stanisic et al. (2017) with 11 microsatellites (AHT4, AHT5, ASB23, CA425, HMS2, HMS3, HMS6, HTG6, HTG7, HTG10 and VHL20), were used for the assessment of the standard genetic diversity parameters, i.e., number of different alleles (A), effective number of alleles [Ae = 1∕(1 − H_E)], number of private alleles (PA), observed heterozygosity (H_O), expected heterozygosity (H_E), and inbreeding coefficients (F_IS) in each group, with GenAlEx 6.5 (Peakall & Smouse, 2012). The same software was employed for Principal Coordinates Analysis (PCoA) used for assessing genetic affinities of groups, and for estimating pairwise group F_ST values and the number of migrants per generation (Nm) calculated as [Nm = 1∕4(1∕F_ST − 1)]. The matrix of pairwise group F_ST values was summarized by two-dimensional scaling (MDS) using Paleontological Statistics (PAST) ver. 3.0 (Hammer, Harper & Ryan, 2001). The analysis of molecular variance (AMOVA, Excoffier, Smouse & Quattro, 1992), aimed at partitioning the overall molecular variation to within- and among-group variation, was performed with Arlequin ver. 3.0 (Excoffier, Laval & Schneider, 2005). Bayesian clustering analysis was carried out with STRUCTURE in order to assess the optimal number of independent genetic groups (K). Burn-in length and run length were 500,000 iterations each, and ten independent runs for each of the assumed K = 1–6 were performed under the admixture and correlated allele frequencies models. The optimal K value was determined by the highest mean estimated log probability of data (Pritchard, Stephens & Donnelly, 2000).

Based on the variability of the mtDNA D-loop, Stanisic et al. (2017) reported 19 mtDNA haplotypes in Banat and Balkan donkeys (h1–h9 belonging to the Clade 1, h10–h19 to Clade 2, GenBank accession numbers KR081377 –KR081395) of which h12 was predominant in the Banat donkey that harboured also a rare h18. These haplotypes were used together with those found in six Spanish breeds and two populations from Africa (Morocco and Zimbabwe) (GenBank accession numbers AF416593, AF416594, AF416595, AF416596, AF416597, AF416598, AF416599, Aranguren-Mendez et al., 2004), six Italian donkey breeds (KX622700 –KX622727, Cozzi et al., 2017), donkeys from Turkey and Cyprus (Cinar-Kul et al., 2016), ancestral Nubian wild donkeys (HM622634, HM622635, HM622636, Kimura et al., 2011), Somali wild donkeys (AY569545, AY569546, AY569547, Beja-Pereira et al., 2004) and a donkey reference sequence (X97337, Xu, Gullberg & Arnason, 1996), for assessing their genealogical relationship by constructing median-joining (MJ) network with NETWORK 4.6.1.2 (Bandelt, Forster & Rohl, 1999). The matrix comprising 71 mtDNA D-loop sequence was aligned using MUSCLE (Edgar, 2004) in MEGA4 (Tamura et al., 2007), and was truncated to match the length of the shortest available sequence (307 base pairs-bp). The parameter ε was kept at 0, and resulting network was manually assembled for visualization.

Banat donkey population in Serbia

The largest number of individuals that differ morphologically from the Balkan donkey and, according to local breeders, represent the Banat donkey, is currently present at KO: 12 adult individuals (ten females, two males) and seven foals <3 years old, while only three adult females and one foal were recorded at SP. The Banat donkey is maintained together with the Balkan donkey in both KO and SP, and the mating between breeds is allowed at both localities. Local breeder from ZA, however, aimed at preventing interbreed mating in their freely roaming herd by separating females during the mating season, and thus, we used individuals from this locality in our study. At present, seven adult females, one adult male and five foals are found at ZA. This herd was initially established with 10 individuals some 25 years ago, and, over time, old animals have successively been replaced with younger ones, some of which were acquired from nearby local farmers. Altogether, the number of extant adult Banat donkey individuals in ZA, KO and SP is less than 100, with apparent domination of females (20 females vs. three males).

Morphological characterization of the Banat donkey

All 18 morphological traits measured in individuals from four groups were homogeneous (CV ≤ 21%, Table 1). Mean values of all body measurements were larger in Banat donkeys than in hybrids and individuals from both sub-populations of the Balkan donkeys (Table 1).

A significant effect of donkey groups on variation of studied parameters was supported by both Hotelling-Lawley and Wilks approximate F-tests in MANOVA analysis. Hotelling’s T2 test revealed statistically significant difference between BanD and both sub-populations of the Balkan donkey (Table 2). Furthermore, Tukey’s test in one-way ANOVA analysis clearly separated BanD into a distinct homogeneous group by all studied morphological parameters except ear length (Table 3), and revealed that statistically significant differences between two sub-populations of the Balkan donkey were obtained for four parameters: back height, chest depth, hip height and wither height, which are highly correlated (see below). Contribution of expected variances of groups’ main effect to the total variance (i.e., the intraclass correlation coefficient reflecting the ratio between expected group variance and sum of expected group variances and errors) was highest for body weight and hip height (>60%), and lowest for chest width and ear length (<20%) (Fig. S2). Thus, the difference between groups was best demonstrated by two parameters, body weight and hip height.

In the PCA analysis, the first three PCs, which explain 75% of the total variance, were selected for further rotation (Table S1). The first rotated PC explains almost half of the total variance (34.4%), and nine out of 18 studied morphological parameters have their highest loadings with this PC. Since parameters body weight and hip height, which have the highest contribution of factor group to the total expected variance (Fig. S2), belong to the same PCA-group, they were considered to be strongly correlated. It is worth mentioning that this analysis revealed that morphological traits that describe body length and weight, and the size of the head and extremities are highly correlated (first PCA-group), as well as those related to the height of animals (PCA-group 2) and the size of the chest and croup (PCA-group 3).

Morphological distinctiveness of Banat donkeys was evident also in the UPGMA dendrogram (Fig. 2) and in the Canonical Discriminant Analysis in which first two canonical variables jointly explain 93.0% of the total variance (Fig. 3). BanD was separated from other groups along the first canonical variable which explains most of the total variance (74.4%), while HY was separated from BalkD-BGP/BalkD-RGP along the second canonical variable which explains 18.6% of the total variance. The outcomes of this analysis are in accordance with the outcomes of the Hotelling’s T2 test which did not support significance of overall difference between two sub-populations of the Balkan donkey by morphological parameters.

Discriminant model based on all 18 measured parameters achieved 76.5% of classification accuracy, and Stepwise Forward Discriminant Analysis procedure revealed that this accuracy could be also achieved by the model that includes only four parameters: hip height, croup width, body length and chest depth (Fig. 4). Taking into account the loadings of these parameters with first three rotated PCs in the PCA analysis (Table S1), it is evident that they belong to different PCA-groups, i.e., hip height and body length to the first, croup width to the second, and chest depth to the third PCA-group.

Genetic characterization of the Banat donkey

We take advantage of the dataset of Stanisic et al. (2017) (genetic profiles of BanD, HY, BalkD-BGP and BalkD-RGP assessed with 11 microsatellites from the Equine Genotype Panel 1.1 recommended by the International Society for Animal Genetics (ISAG) for equine parentage and identification testing) to characterize genetically Banat donkey and to compare four donkey groups.

The parameters of genetic diversity of Banat donkey, represented by seven females, were as follows: A = 58, Ae = 3.749 ± 0.255, H_O = 0.799 ± 0.075, and H_E = 0.777 ± 0.024 (Table 4). Two private alleles were found in this group characterized also by the significant excess of heterozygotes (F_IS =  − 0.111 ± 0.095).

All pairwise-group F_ST values were statistically significant at 95% level (Table 5). F_ST value between Banat donkeys and potential hybrids was in the range of that observed between the two sub-populations of the Balkan donkey (0.044 vs. 0.037), while those between Banat donkeys and each of the two sub-populations of the Balkan donkey were almost doubled (Table 5). Concordantly, the highest number of migrants per generation, reflecting past gene flow, was observed between Banat donkeys and potential hybrids (Nm = 5.38) and between the two sub-populations of the Balkan donkey (Nm = 6.44). In MDS graph (Fig. S3), Banat donkeys and potential hybrids were separated from the two sub-populations of the Balkan donkey along the first dimension, and the same grouping was observed in the PCoA graph defined with axis 1 (explaining 68.70% of variability) and axis 2 (explaining 29.60% of variability) (Fig. 5). A small but statistically significant portion of the total genetic variability was allocated to the among-population component (6.63%, P = 0.001, Table S2). The outcomes of the Bayesian clustering analysis (Fig. 6) were concordant with those reported previously by Stanisic et al. (2017), with optimal number of three genetic groups (Fig. S4). Banat donkeys and potential hybrids were assigned to the unique green gene pool, and they kept their genetic integrity in analyses at K = 3–5. The same pattern was observed for Balkan donkey individuals assigned to the red gene pool (RGP), while those assigned to the blue gene pool (BGP) displayed admixed nuclear genetic profiles at K = 4–5.

Genealogical relationship of the mtDNA haplotypes of the Banat and other studied donkeys

Majority of studied mtDNA haplotypes were grouped into two lineages corresponding to Clades 1 and 2 defined previously by Beja-Pereira et al. (2004) and Kimura et al. (2011), while those found in Somali wild donkeys were allocated to the third lineage which was genealogically linked to Clade 2 (Fig. 7). Since we used somewhat shorter sequences than those used for the delineation of distinct mtDNA haplotypes in corresponding studies, haplotypes distinguished from each other in those studies were sometimes presented jointly within the same node in our MJ network. For instance, haplotypes RAD-H4 and ROD-H5 found in two Italian donkey breeds occupy the same node because the mutation that distinguishes them was not present in our aligned matrix that was truncated. The node positioned in the centre of Clade 2 comprised two ancient Nubian wild donkey haplotypes, haplotypes h12, h17 and h18 found in Banat and Balkan donkeys, ATI-1 typical for Catalana, Mallorqina and Zamorano-Leonessa breeds from Spain and for two donkey populations from Africa (one from Morocco and the second one from Zimbabwe), and haplotypes found in individuals belonging to five donkey breeds from Italy (AMD-H6 present in Amiata donkeys, ASD-H1 and ASD-H2 found in Asinara donkeys, RAD-H1 and RAD-H2 found in Ragusano donkeys, ROD-H4 from Romagnolo donkey, and SAD-H1 from Sardo donkeys).

Banat donkey population in Serbia

The donkey breed that is distributed today throughout the Balkan Peninsula and represented by a relatively small number of individuals is the heterogeneous Balkan donkey (Stanisic et al., 2017). However, additional donkey breeds may also be present in this large region, such as those found in Croatia (Istrian and North-Adriatic donkeys, Ivankovic et al., 2002) and in the north-eastern part of Serbia, the Banat region. The latter, the Banat donkey, is at present an uncharacterized and highly neglected donkey breed in Serbia. For instance, the common practice of local breeders at three localities in Serbia where large donkey herds are currently found, KO, SP and ZA, was to maintain a couple of individuals that differ from the prevalent Balkan donkey in order to demonstrate the morphological variability of donkeys to the people visiting their farms. Furthermore, over the past 25 years mating between the breeds was mainly allowed in KO and SP, and in a somewhat more limited fashion in ZA. Since in ZA we found potential hybrid individuals that are morphologically indistinguishable from the Balkan donkey and are characterized by a distinct nuclear gene pool (Stanisic et al., 2017; this study), it is evident that well-planned actions aimed at the conservation of this valuable genetic resource are urgently needed.

We found that, out of 281 sexually mature individuals registered in 2015 at KO, SP and ZA, 23 animals may be delineated as the Banat donkey based on their morphological features. Due to such a low number of Banat donkeys in Serbia, this breed may be characterized as critically endangered (FAO, 2015). However, according to local breeders, one or two Banat donkeys per household are still present in the rural parts of Vojvodina, the northern part of Serbia comprising Banat, Bačka and Srem districts. Therefore, future efforts of breed conservation should aim not only to preserve the original traits of the Banat donkey, through well-planned mating strategies in maintained herds (Kugler, Grunenfelder & Broxham, 2008), but also to acquire such individuals from local farmers. These individuals should be used for mating in order to prevent loss of genetic diversity and inbreeding through mating among relatives in small-sized populations (e.g., Gutiérrez et al., 2005; Quaresma et al., 2014). Furthermore, the search for Banat donkeys should encompass a larger area, preferably the entire Balkan region, because of the common transfers of donkeys throughout it. Given the general lack of controlled mating of donkeys in the past in the Balkans (Kugler, Grunenfelder & Broxham, 2008; Stanisic et al., 2015; Stanisic et al., 2017), all individuals should be characterised both morphologically and genetically in order to be assigned unambiguously to the Banat donkey breed.

Morphological characterization of the Banat donkey

Comparative analysis of the 18 morphological traits of Banat and Balkan donkeys revealed that Banat donkeys are >70 kg heavier, >10 cm taller and >15 cm longer than the potential hybrids and individuals from both sub-populations of the Balkan donkey. Their chest circumference is >20 cm wider, and their hips and withers are 8–15 cm higher than those measured in other studied groups. Body measurements of a single studied 8-years old Banat donkey male (Table 1) generally correspond to those recorded for younger females (average age 5.25 years), with a somewhat larger difference observed in body length, carpal circumference, chest depth, hip height and tibia circumference, which are the most common morphological traits used for donkey breed characterization (e.g., Folch & Jordana, 1997; Sargentini et al., 2009; Sargentini et al., 2018; Kefena et al., 2011; Kosťuková et al., 2015; Labbaci et al., 2018; Ayad et al., 2019). This would indicate a possible sexual dimorphism in the Banat donkey which is indeed present in certain donkey breeds (e.g., Catalonian donkey from Spain, Folch & Jordana, 1997). For instance, Folch & Jordana (1997) found that the greatest difference between Catalonian donkey sexes is at the cephalic level, and also in height from the posterior third of the body and extremities perimeter, indicating that females are more elevated and more slender than males. However, lack of sexual dimorphism or explicit anatomical dissimilarities between the sexes is also common in donkeys (e.g., six Ethiopian donkey populations, Kefena et al., 2011), and thus, inferences on the possible Banat donkey’s sexual dimorphism require analyses of additional males.

A profound morphological difference between the Banat and Balkan donkey was demonstrated in all the performed statistical analyses, suggesting that the Banat donkey may be acknowledged as a distinct breed based on its morphological features (Alderson, 2003). The most important parameters for the separation of the Banat donkey from the Balkan donkey and potential hybrids are body weight and hip height, which however are highly correlated. We also found that the classification accuracy in the Stepwise Forward Discriminant Analysis obtained with the model that included all studied morphological traits may also be achieved by applying the model that includes only four traits: hip height, croup width, body length and chest depth. These results highlight the redundancy of the majority of morphological traits used in our study, which are, more or less, similar to those commonly used for morphological characterization of donkey breeds (e.g., Folch & Jordana, 1997; Sargentini et al., 2009; Sargentini et al., 2018; Kefena et al., 2011; Kosťuková et al., 2015; Labbaci et al., 2018; Ayad et al., 2019). Nonetheless, an important finding is that hip height/chest depth, croup width and body length are weakly correlated because these parameters (or parameters with which they are highly correlated) belong to different PCA-groups (Table S1). This would suggest that all these parameters should be used for assessing morphological distinctiveness of Banat and Balkan donkey breeds. It is worth mentioning that we found that two sub-populations of the Balkan donkey, which are distinct at the genetic level (Stanisic et al., 2017; this study), may be distinguished at the morphological level as well, based on one of the four highly correlated parameters, namely back height, chest depth, hip height and wither height.

The fact that extant Banat donkeys are taller, longer, somewhat wider and much heavier than Balkan donkey is in line with speculations of local breeders from Serbia that they are descendants of the heard of more robust, larger and stronger Spanish donkeys transferred from Spain to the Banat region ∼300 years ago for a specific purpose, namely to help farmers work the local vineyards. Their body constitution was apparently maintained by more or less controlled mating over time despite the general lack of such management practice in Serbia (Stanisic et al., 2015) and worldwide (Kugler, Grunenfelder & Broxham, 2008). However, the discovery of a potential hybrid population in ZA, where Banat donkeys are kept together with Balkan donkeys (Stanisic et al., 2017), highlights the need for urgent preservation of the Banat donkey and maintenance of its original traits through well-planned mating strategies (Kugler, Grunenfelder & Broxham, 2008). It is noteworthy that Banat donkey females produce greater quantities of quality milk than Balkan donkey females (∼30 l, data acquired from local breeders), implying that maintenance of this breed by local breeders may be economically beneficial.

Although further genetic studies are required to decipher whether the black-coated Catalonian donkey from the north-eastern Spain (Folch & Jordana, 1997) was indeed introduced to the Banat region, this breed has been commonly used in the past for the improvement of numerous European and American donkey breeds (Jordana & Folch, 1996). The mating between large-sized Catalonian donkey (Folch & Jordana, 1997) and mid-sized and lighter coloured Balkan donkeys (Stanisic et al., 2015) may account for the successive decrease in size and change in the coat colour in extant Banat donkeys. A unique morphological feature of the Banat donkey, black stripes on legs resembling those typical for Somali wild donkeys (Groves, 1986; Moehlman, 2002) present today in Somalia, Ethiopia and Eritrea (Moehlman, 2002), may represent a relict and primitive marking (Johnson & Johnson, 2008) indicative of their origin. According to Sargentini et al. (2018), black stripes on legs and some other body markings characterize also the Italian Amiata donkey breed and their ancestors.

Genetic characterization of the Banat donkey

Although in this study we used genetic profiles at the nuclear DNA level of only seven adult Banat donkey females, all the parameters of genetic diversity calculated for this group were similar or even higher compared to those observed in the two larger sub-populations of Balkan donkey females (Table 4), and also in populations of other donkey breeds distributed worldwide (Jordana, Folch & Sanchezm, 1999; Jordana et al., 2016; Aranguren-Mendez, Jordana & Gomez, 2001; Ivankovic et al., 2002; Blasi et al., 2005; Ciampolini et al., 2007; Zhu et al., 2013; Matassino et al., 2014; Rosenbom et al., 2015). Furthermore, two private alleles were observed in Banat donkeys characterized also by the significant excess of heterozygotes. Therefore, although small, the female Banat donkey population currently kept in ZA is not severely affected by the loss of genetic diversity or by inbreeding. The hybrid group of eight individuals found in ZA is characterized by slightly lower levels of genetic diversity, three private alleles, and statistically insignificant excess of homozygotes. Since they are morphologically indistinguishable from Balkan donkeys but mainly harbour Clade 2 mtDNA haplotypes which are dominant in Banat donkeys, they most likely emerged as a result of mating between Banat donkey females with Balkan donkey males which are abundant in ZA. However, the scores of these individuals in the PCoA plot (Fig. 5) were not intermediate between the Banat and the Balkan donkey, and thus we cannot exclude the possibility of genetic input from sources other than the Balkan donkey in these donkeys. This is supported by the presence of three private alleles in these eight individuals. Given the fact that Banat donkey individuals currently present in ZA were acquired from local farmers which did not keep records on individuals which were used for mating with Banat donkey, the genetic source accounting for the positioning of the score of the hybrid group in PCoA graph remains unknown. Nonetheless, it is very likely that these individuals that display genetic affinity towards Banat donkey and morphological similarity with the Balkan donkey are hybrids with rather complex and unknown ancestry.

Genetic distinctiveness of the Banat donkey was supported by all the analyses of genetic differentiation that were performed in the course of our study. Thus, we support the previous findings of Stanisic et al. (2017) that the Banat donkey differs genetically from the Balkan donkey, and that it may be acknowledged as a distinct donkey breed in Serbia, the Balkans and worldwide (Alderson, 2003). A group of seven Banat donkey females represents a good core population which can be enlarged by introduction of unrelated individuals, both females and males, required for maintaining high levels of genetic diversity essential for successful preservation of the breed (e.g., Gutiérrez et al., 2005; Quaresma et al., 2014). Management practice in ZA should be improved in order to limit the mating between different breeds and enable preservation of original Banat donkey traits (Kugler, Grunenfelder & Broxham, 2008).

Genealogical relationship among the mtDNA haplotypes of the Banat and other studied donkeys

Although our mtDNA data do not have the power to fully resolve the genealogical relationship among the mtDNA haplotypes used in our study, the finding that haplotypes h12 and h18 which are prevalent in the Banat donkey occupy the same node with haplotype ATI-1 typical for Catalana, Mallorqina and Zamorano-Leonessa breeds from Spain and for two donkey populations from Africa (one from Morocco and the other from Zimbabwe), would indeed indicate their close genealogical relationship. On the other hand, the presence of two ancient Nubian wild donkeys’ haplotypes in the same node in MJ network may imply a relict nature of mtDNA haplotypes found in extant Banat donkeys. Thus, the Banat donkey may indeed represent a breed that was established over time in the new environment, the Banat region in Serbia, from a herd that may have been introduced from Spain some 300 years ago. According to Aranguren-Mendez et al. (2004), mtDNA haplotypes ATI-1 and ATI-3, shared by three Spanish donkey breeds and two African donkey populations in more or less similar proportions, represent the most ancestral mtDNA types.

Recent molecular evidence questions the hypothesis of Beja-Pereira et al. (2004) and Kimura et al. (2011) that extinct Somali wild donkeys were the ancestors of the Clade 2 donkeys (e.g., Kefena et al., 2014; Stanisic et al., 2017; Xia et al., 2019) and, in addition, puts forward alternative centres of origin of donkeys belonging to this lineage (Rosenbom et al., 2015). Our data, however, would imply that the origin of Clade 2 donkeys may be much more complex than previously thought. For instance, we found also that a number of haplotypes found in individuals belonging to several donkey breeds from Italy (including the Amiata breed, Cozzi et al., 2017) are positioned within the same central node of Clade 2 as Spanish ATI-1 haplotype, indicating a close genealogical relationship of donkeys from Iberian Peninsula and Italian donkeys reported previously by Cozzi et al. (2017). Also, Spanish ATI-1 haplotype is found in Anatolian and donkeys from the Cyprus (Cinar-Kul et al., 2016). Therefore, our data support a close genealogical relationship of donkeys from all three south European Peninsulas, Turkey/Cyprus and Africa (i.e., the Mediterranean region), and furthermore, reveal that haplotypes found in several Italian donkeys occupy an intermediate position at the backbone of the MJ network between Clade 1 and Clade 2 haplotypes. These haplotypes of Italian donkeys were not available and were not included into the study of Stanisic et al. (2017) in which an intermediate positions between two Clades were occupied mainly by haplotypes found in donkeys from Albania, Bulgaria and Croatia.

Regarding the Banat donkey, which is characterized by the high proportion of Clade 2 mtDNA haplotypes, this breed may indeed represent a lineage evolving from Spanish matrilineal lineages, while a unique morphological feature shared between the Banat donkey, the Amiata donkey (Sargentini et al., 2009; Sargentini et al., 2018) and the Somali wild donkeys (Groves, 1986; Moehlman, 2002) may imply their common origin.