Genetic and Population Structure of Croatian Local Donkey Breeds

: The two native Croatian donkey breeds (Littoral-Dinaric donkey and Istrian donkey) were marginalized in the second half of the 20th century and were on the verge of biological extinction. The aim of this study was to analyze the demographic and genetic status of two donkey breeds, two decades after the start of protection by analyzing their pedigrees and genetic structure. The average generation interval was higher for the Istrian donkey (7.73) than for the Littoral-Dinaric donkey (7.27). The rate of the effective number of founders compared with the effective number of ancestors in the Littoral-Dinaric donkey (1.03; 325/316) and in the Istrian donkey (1.08; 70/65) revealed no evidence of a genetic bottleneck. The inbreeding coefﬁcient ( F ) and the average relatedness coefﬁcient ( AR ) was lower in the Littoral-Dinaric donkey population (0.99%; 0.13%) than in the Istrian donkey population (1.77%; 1.10%). Genetic microsatellite analysis showed relatively high genetic diversity in Littoral-Dinaric donkey and Istrian donkey breeds, expressed by mean allele number (5.92; 5.85) and expected heterozygosity (0.650; 0.653). Genetic differentiation between the Littoral-Dinaric donkey and the Istrian donkey has not signiﬁcantly increased in the last two decades ( F ST = 0.028). Genetic analysis also showed no evidence of high inbreeding or genetic bottleneck in both breeds. A total of 11 haplotypes including 28 polymorphic sites were found in 30 samples. Analysis of mtDNA has shown that the Littoral-Dinaric donkey and Istrian donkey breeds belong to the Equus asinus africanus group. The study conﬁrms the need to use different analytical approaches to get a regular and complete insight into the situation and trends within and between breeds, so that the existing diversity can be fully preserved.


Introduction
Donkey breeds, similarly to other domestic animal species, are an important part of national, regional, and global genetic resources. The total number of donkey breeds in the world is a relatively small proportion (178; 2.8%) of the world's mammalian breeds [1]. In the 19th and 20th centuries, local donkey breeds came under extinction pressure as they lost their function as working animals in the community, especially as pack animals and riding animals. The number of extinct donkey breeds in relation to the total number of donkey breeds is relatively low (4; 2.2%), a large number of breeds are "at risk" (22.5%), and for most donkey breeds, their endangered status is unknown (128; 71.9%) [1]. The donkey is an animal of arid areas, whose importance may be significant if climate change and global warming occurs.

Pedigree Analysis
The Studbook data of the Istrian donkey and Littoral-Dinaric donkey up to 2020 were analyzed. The genealogical information was traced back to the founder animals, the oldest of which was born in 1972. The created database consists of the pedigree information for 1233 Istrian donkeys and 5030 Littoral-Dinaric donkeys. For each animal in the database, information such as the UELN (Universal Equine Life Number), name of the individual, date of birth, sex, and pedigree information were collected. The reference populations consist of those animals that are potentially contributing to the next generation. In this study, animals with two known parents were chosen as the first reference population (REF 01), and active (living) donkeys of both sexes were selected as the second reference population (REF 02).
Pedigree analysis was conducted using the ENDOG v4.8 software [11]. Several analyses were carried out to evaluate the pedigree and population parameters. Pedigree completeness was observed through the mean number of complete generations, the mean number of maximum generations, and the mean number of equivalent generations [32].
Parameters related to inbreeding (mean inbreeding coefficient (F), mean average relatedness coefficient (AR), effective population size (Ne)), probability of gene origin (effective number of founders (f e ), effective number of ancestors (f a ), effective number of founder herds (f h )), genetic contributions of the founders and ancestors (genetic conservation index (GCI), and generation interval (GI)) were assessed. The average relatedness coefficients and individual inbreeding coefficients were calculated based on Gutiérrez and Goyache [11]. Boichard et al. [32] define f e as the number of equally contributing founders that will produce the same genetic diversity as in the assessed population, and f a as the minimal number of ancestors necessary to explain the genetic diversity in the reference population. The effective number of founders (f e ) was computed by the formula f e = 1/∑ (p i 2 ), where p i is the proportion of the genes of the descendant population contributed by founder i [33]. The effective number of ancestors (f a ) was calculated on the basis of q i as the marginal genetic contribution of ancestor i (f a = 1 / ∑ q i 2 ) [32]. The ratio between the effective number of founders (f e ) and effective number ancestors (f a ), as an indicator of a population bottleneck, were calculated according to Boichard et al. [32] The genetic conservation index was calculated on the basis of the proportion of genes of founder i (p i ) in the pedigree of an animal (GCI = 1 / ∑ p i 2 ) [34]. Generation intervals (GI) were calculated as the average age of the parents at the birth of the progeny subsequently used for reproduction, as well as the average age of the parents of all offspring [35,36], and were calculated for the following four paths: father-daughter/son, and mother-daughter/son. The four pathways were compared pairwise using independent sample t-tests.

Microsatellite Analysis
Hair root samples were collected from 60 Littoral-Dinaric donkeys (45 female and 15 male animals in age from 4 to 9 years) and 60 Istrian donkeys (43 female and 17 male animals in age from 4 to 10 years) reared on different breeding farms in the Mediterranean part of Croatia (Littoral-Dinaric donkey in Dalmatia, Istrian donkey in Istria). Sampling was achieved among minimally related individuals using pedigree information and breeder's information. DNA extraction and genotyping were performed in an authorized laboratory. Thirteen microsatellite loci (AHT4, HMS2, HMS7, TKY297, TKY343, ABS23, HMS3, HTG10, TKY312, HMS18, HMS6, HTG7, TKY337) which have been recommended for individual identification and parentage verification of equines, were screened from the set recommended by the International Society for Animal Genetics (ISAG) and previous reports [37,38].
Genetic diversity parameters were estimated for each microsatellite locus and across all loci for each population by total (N A ), effective (N AE ), and private number of alleles (N P ), unbiased estimates for observed (H O ) and expected (H E ) heterozygosity, and migration rate (Nm), using GenAlEx 6.5 [39]. Effective (N AE ) are calculated as follows: where is p i is frequency of the i-th allele in a locus and ∑ p i 2 is the sum of the squared population allele frequencies.
Unbiased estimates for observed (H O ) heterozygosity is calculated as the proportion of N samples that are heterozygous at a given locus. Expected (H E ) heterozygosity is calculated as: where H E is the expected heterozygosity (i.e., the proportion of heterozygosity expected under random mating), and p i is the allele frequency of the i-th allele.
Migration rates (Nm) are calculated as follows: where F ST is Wrights fixation index calculated on the basis of the average expected heterozygosity of the subpopulations (H E ) and the expected heterozygosity of the total population (H T ). Cervus software version 3.0.7. [40] was used to calculate the null allele frequency (F (null) ) and polymorphic information content (PIC) using the following equations: where p ij is the frequency of the j-th allele for marker i and the summation extends over n alleles.
The Wright's F statistics (F ST , F IS and F IT ) as proposed by Weir and Cockerham [41], and the allelic richness (AR), were computed using FSTAT Version 2.9.3 [42]. Fisher's exact test was performed to test possible significant departures from the Hardy-Weinberg (HW) proportions using GenePop 4.3 [43]. p-values of the heterozygote deficit and the excess for each locus were obtained simultaneously.
Genetic relationships among individuals were represented using factorial correspondence analysis (FCA) performed with the option of using 3D over populations with Genetix 4.05 software [44]. Likelihood test of breed assignment and analysis of molecular variance were performed using the GenAlEx 6.5 [39]. Heterozygosity-excess was quantified using the approach of Cornuet and Luikart [45], implemented in the program Bottleneck 1.2.02 [46], and shifted away from an L-shaped distribution of allele frequencies.

Mitochondrial DNA Analysis
MtDNA was extracted from blood samples of 15 unrelated older female animals (older than 14 years, bred in different and geographically dispersed herds) per donkey breed (30 samples in total). Based on the complete donkey mtDNA sequence (GenBank X97337) [21], two pairs of primers were designed (5 -AGTCTCACCATCAACCCCCAAAGC-3 and 5 -CCTGAAGTAGGAACCAGATG-3 ) to amplify a 358 bp fragment of the hypervariable D-loop region comprised between sites 15,479 and 15,837. Amplification was performed on an MJ Research PTC-100 thermal cycler with the following conditions: initial denaturation at 95 • C for 5 min, followed by 32 cycles of 94 • C (1 min) at 52 • C (30 s) and 72 • C (1 min), and a final extension at 72 • C (5 min). Fragments were sequenced using the ABI PRISM BigDye Terminator Cycle Sequencing Ready Reaction Kit and an ABI PRISMTM ® 310 Genetic Analyzer. Multiple alignments of mtDNA sequences were performed using the Clustal-W program (Version 1.82, 2001), and were analysed with MEGA software version X [47].

Pedigree
The demographic data of the Littoral-Dinaric donkey and Istrian donkey populations are shown in Figure 1. The number of animals per year constantly grew in the two last decades. The growing trend in the size of populations is the result of the establishment of programs for their protection, the definition of breeding programs, the establishment of Studbook registers (opened during the inventory of populations), the promotion of breeds (tourism, hobby, milk, meat), the establishment of breeding associations, breeding exhibitions, and other activities. In the period from 2012 to 2014, stagnation is observed, mainly due to changes in the support system for breeders of native breeds (adaptation to the EU legal framework; Croatia joined the EU in 2013) and the Studbook registers closed.
following conditions: initial denaturation at 95 °C for 5 min, followed by 32 cycles of 94 °C (1 min) at 52 °C (30 s) and 72 °C (1 min), and a final extension at 72 °C (5 min). Fragments were sequenced using the ABI PRISM BigDye Terminator Cycle Sequencing Ready Reaction Kit and an ABI PRISMTM ® 310 Genetic Analyzer. Multiple alignments of mtDNA sequences were performed using the Clustal-W program (Version 1.82, 2001), and were analysed with MEGA software version X [47].

Pedigree
The demographic data of the Littoral-Dinaric donkey and Istrian donkey populations are shown in Figure 1. The number of animals per year constantly grew in the two last decades. The growing trend in the size of populations is the result of the establishment of programs for their protection, the definition of breeding programs, the establishment of Studbook registers (opened during the inventory of populations), the promotion of breeds (tourism, hobby, milk, meat), the establishment of breeding associations, breeding exhibitions, and other activities. In the period from 2012 to 2014, stagnation is observed, mainly due to changes in the support system for breeders of native breeds (adaptation to the EU legal framework; Croatia joined the EU in 2013) and the Studbook registers closed. In the 2020 Littoral-Dinaric donkey population, 43.1% of adult males and 68.0% of females are reproductively active, and in the Istrian donkey population, reproductive activity is similar (50.9% males and 65.7% females). About half of the adult males occasionally participate in reproduction to reduce inbreeding (stallion rotation). The age structure of the donkey populations is favorable. In the living populations of Littoral-Dinaric donkeys and Istrian donkeys (REF 02), 32.7% and 33.9% were the three years of age, respectively ( Figure 2). In the populations of the Littoral-Dinaric donkey and Istrian donkey, only 11.0% and 10.5% were older than 20 years, respectively. The number of active Littoral-Dinaric donkey and Istrian donkey breeders in 2020 is 809 and 158, respectively, and the average herd size is 4.1 and 4.8 animals, respectively. In the 2020 Littoral-Dinaric donkey population, 43.1% of adult males and 68.0% of females are reproductively active, and in the Istrian donkey population, reproductive activity is similar (50.9% males and 65.7% females). About half of the adult males occasionally participate in reproduction to reduce inbreeding (stallion rotation). The age structure of the donkey populations is favorable. In the living populations of Littoral-Dinaric donkeys and Istrian donkeys (REF 02), 32.7% and 33.9% were the three years of age, respectively ( Figure 2). In the populations of the Littoral-Dinaric donkey and Istrian donkey, only 11.0% and 10.5% were older than 20 years, respectively. The number of active Littoral-Dinaric donkey and Istrian donkey breeders in 2020 is 809 and 158, respectively, and the average herd size is 4.1 and 4.8 animals, respectively.
An analysis of all pedigree entries of the Littoral-Dinaric donkey population revealed 2724 animals for which both parents were known (REF 01, Table 1 The effective number of founders (f e ) and ancestors (f a ) for the Littoral-Dinaric donkey was 325 and 316, respectively ( Table 1). The ratio between f e and f a is 1.03, indicating the absence of a bottleneck in the populations. A similar ratio between f e and f a is observed in Istrian donkeys (1.08), and does not indicate a bottleneck, which is confirmed by the results of microsatellite analysis.  The effective number of founders (fe) and ancestors (fa) for the Littoral-Dinaric donkey was 325 and 316, respectively ( Table 1). The ratio between fe and fa is 1.03, indicating the absence of a bottleneck in the populations. A similar ratio between fe and fa is observed in  The average inbreeding coefficient (F) in the Littoral-Dinaric donkey population was 0.99%, and ranged from a minimum of 0.02% (2000, the beginning of breeding consolidation and a larger number of founders) to a maximum of 0.99% (2020) (Figure 3). In the last decade, the average inbreeding range was constantly growing, from 0.28 to 0.99, indicating some problems in mating schemes. The average inbreeding in the population of the Istrian donkey was 1.77%, ranging from a minimum of 0.11% (2000, the beginning of breed consolidation) to a maximum of 1.77% (2020). In recent years, the average inbreeding coefficient ranged from 1.79% to 1.77%. The mean AR coefficient in the Istrian donkey population was 1.10%. The effective population size calculated via individual increase in inbreeding is 25.99 and is lower than in the Istrian donkey population (32.78). decade, the average inbreeding range was constantly growing, from 0.28 to 0.99, indicating some problems in mating schemes. The average inbreeding in the population of the Istrian donkey was 1.77%, ranging from a minimum of 0.11% (2000, the beginning of breed consolidation) to a maximum of 1.77% (2020). In recent years, the average inbreeding coefficient ranged from 1.79% to 1.77%. The mean AR coefficient in the Istrian donkey population was 1.10%. The effective population size calculated via individual increase in inbreeding is 25.99 and is lower than in the Istrian donkey population (32.78).    The REF 02 populations of the Istrian donkey had a high completeness value in the first two generations (96.94% and 51.68%, respectively), whereas in the third generation, completeness decreased to 18.59%, and in the population of the Littoral-Dinaric donkey, the completeness value is lower (60.21%, 24.28%, and 5.52%, respectively). Beyond the fifth parental generation, the completeness of both donkeys' REF 02 populations was close to zero.
Estimates of the average age of parents at the birth of their offspring in the Istrian donkey and Littoral-Dinaric donkey populations ranged from 7.41 to 7.79 years, respec- the completeness value is lower (60.21%, 24.28%, and 5.52%, respectively). Beyond the fifth parental generation, the completeness of both donkeys' REF 02 populations was close to zero.
Estimates of the average age of parents at the birth of their offspring in the Istrian donkey and Littoral-Dinaric donkey populations ranged from 7.41 to 7.79 years, respectively (Table 2). Higher values are observed in the Littoral-Dinaric donkey population in the mother-daughter pathway (7.79 ± 4.296), compared with the father-daughter pathway (p < 0.02). Considering all pathways, the average age of parents at the birth of their offspring was 7.50 ± 4.067 years in the Istrian donkey populations, and 7.59 ± 4.057 years in the Littoral-Dinaric donkey populations. The average generation interval in the Istrian donkey population was 7.73 ± 3.94 years and ranged from 7.53 ± 3.84 (mother-daughter) to 7.89 ± 3.91 (father-daughter). The average generation interval in the Littoral-Dinaric donkey population was 0.47 years less (7.27 vs. 7.73 years).

Microsatellite Variability and Genetic Diversity Indices
All 13 microsatellite markers were polymorphic in the whole sample and in each breed. The basic parameters of genetic diversity across microsatellites are presented in Table 3. A total of 83 alleles (N A ) were found, ranging from 3 (HMS6) to 10 (HTG7) per locus, with a mean of 6.31. The average values for E NA and AR were 3.317 and 6.017, respectively.
The abundance of genetic variation in microsatellite loci was indicated by a mean value of PIC index (PIC = 0.624) and 86% of them were highly polymorphic [48]. The estimated null allele frequency (F (null) ) of all 13 microsatellites has values below 20% (varying from −0.030 HMS18 to 0.114 HMS6) and was suitable for genetic analyses. The Ewens-Watterson neutrality showed that none of the tested markers favored any kind of selection, as the F values (sum of square of allele frequency) were within the upper and lower limits of the 95% confidence interval. This shows the suitability and utility of these markers for genetic diversity studies The average observed heterozygosity was lower (H O = 0.648) than the expected heterozygosity (H E = 0.673), resulting in a positive, but not significant, inbreeding coefficient (F IS ) of 1.4%. Of the 13 loci examined, all were within the Hardy-Weinberg equilibrium, with the exception of HMS3, which had a significant excess of heterozygotes of 0.048 (p < 0.01) ( Table 3). The overall inbreeding index F IT averaged 0.058 and no loci made a significant contribution. The overall F ST index revealed that 4.5% of the total genetic variation observed in the sample is explained between breeds. Since most of the total genetic variability comes from differences among individuals (95.5%), and only 4.5% is due to differences between breeds, a low level of genetic differentiation exists between the populations studied. The locus that contributed most to the differentiation of the samples was HMS6, whereas ABS23 proved to be a nondiscriminatory marker (Table 3). A low level of population differentiations contributes to a high migration rate of the samples and loci overall (N m = 12.431).  Summary statistics describing microsatellite marker polymorphisms and genetic diversity per breed are presented in Table 4. Allelic variability between the two breeds was almost identical, with an average N A and E NA of 5.85 and 3.33 for Istrian donkeys and 5.93 and 3.31 for Littoral-Dinaric donkeys, respectively. Of the total 83 alleles, 11 alleles were private (P A ), which were evenly represented between the Istrian donkey (P A = 5) and Littoral-Dinaric donkey (P A = 6) breeds. The highest frequency of P A was observed in the Littoral-Dinaric donkey breed (AHT 4, P A = 0.233), whereas the frequencies of the other P A ranged from 0.008 to 0.05 for both breeds. The mean H O and H E heterozygosity between breeds were very close, 0.640 and 0.653 for Istrian donkeys and 0.656 and 0.650 for Littoral-Dinaric donkeys. A loss of heterozygotes was observed in the Istrian donkey (2.9%), to which two loci gave a relevant significant contribution (HMS3 and HMS6). On the contrary, an excess of heterozygotes (0.2%) in the Littoral-Dinaric donkey population significantly contributed AHT4 and HMS3 markers. The scatter plot of the factorial correspondence analysis (FCA) summarizes the individual relationships in metric space, in which the first axis and second axis contribute to total variation with 5.87% and 5.51%, respectively (data not shown). FCA did not reveal a defined grouping between Istrian donkey and Littoral-Dinaric donkey populations, although there is some indication of separation. The results of the GenAlEx assignment test revealed that 93% of the animals could be assigned to the population from which they were sampled. The highest percentages of individuals assigned to their original population occurred in the Istrian donkey (94.8%) and Littoral-Dinaric donkey populations (91%) (Figure 5). The Wilcoxon test showed no significant results for population bottleneck under the stepwise mutation model (SMM) for Istrian donkeys (p = 0.729) and Littoral-Dinaric donkeys (p = 0.773). The qualitative graphical method based on the allele frequency spectra showed no shift in the allele frequency distribution, and a normal L-shaped curve was observed; thus, these results point to the population resilience of the Istrian donkey and Littoral-Dinaric donkey.

Mitochondrial Variability and Relationship with Other Donkey Breeds
Based on the 358-bp fragment of the D-loop region of the mtDNA, seven different haplotypes were identified in the Littoral-Dinaric donkey population and four haplotypes were identified in the Istrian donkey population. Twenty-seven polymorphic sites were caused by transitions or transversions, although one polymorphic site was the result of an insertion. Nucleotide sequence diversity among all haplotypes ranged from 0.27% to 5.77%. The mean group distance within the Littoral-Dinaric donkey haplotypes was 0.031, and within the Istrian donkey haplotypes it was lower, at 0.005.
To further analyse the phylogenetic relationship of the donkeys studied, they were compared with sixty-nine publicly available mtDNA D-loop sequences of donkey breeds from other countries (https://www.ncbi.nlm.nih.gov/gene/ (accessed on 17 January 2022)). The basic observation is that the Littoral-Dinaric donkey and the Istrian donkey belong to the phylogenetic group Equus asinus africanus (d = 0.030, d = 0.040), and a clear separation from the group Equus asinus somalicus can be observed (d = 0.067, d = 0.049) ( Figure 6). The mtDNA sequences of the Istrian donkey are structured by two phylogenetic groups (haplo-group A, sequence ID 01 and ID 02; haplo-group B, ID 03 and ID 04), but the genetic distance between the groups is small (d = 0.006). The mtDNA sequence of the Littoral-Dinaric donkeys were also divided into two phylogenetic groups (haplo-group C, sequence LDD 01, LDD 02, LDD 03, LDD 04, and LDD 05; haplo-group D, LDD 06, and LDD The Wilcoxon test showed no significant results for population bottleneck under the stepwise mutation model (SMM) for Istrian donkeys (p = 0.729) and Littoral-Dinaric donkeys (p = 0.773). The qualitative graphical method based on the allele frequency spectra showed no shift in the allele frequency distribution, and a normal L-shaped curve was observed; thus, these results point to the population resilience of the Istrian donkey and Littoral-Dinaric donkey.

Mitochondrial Variability and Relationship with Other Donkey Breeds
Based on the 358-bp fragment of the D-loop region of the mtDNA, seven different haplotypes were identified in the Littoral-Dinaric donkey population and four haplotypes were identified in the Istrian donkey population. Twenty-seven polymorphic sites were caused by transitions or transversions, although one polymorphic site was the result of an insertion. Nucleotide sequence diversity among all haplotypes ranged from 0.27% to 5.77%. The mean group distance within the Littoral-Dinaric donkey haplotypes was 0.031, and within the Istrian donkey haplotypes it was lower, at 0.005.
To further analyse the phylogenetic relationship of the donkeys studied, they were compared with sixty-nine publicly available mtDNA D-loop sequences of donkey breeds from other countries (https://www.ncbi.nlm.nih.gov/gene/ (accessed on 17 January 2022)). The basic observation is that the Littoral-Dinaric donkey and the Istrian donkey belong to the phylogenetic group Equus asinus africanus (d = 0.030, d = 0.040), and a clear separation from the group Equus asinus somalicus can be observed (d = 0.067, d = 0.049) ( Figure 6). The mtDNA sequences of the Istrian donkey are structured by two phylogenetic groups (haplo-group A, sequence ID 01 and ID 02; haplo-group B, ID 03 and ID 04), but the genetic distance between the groups is small (d = 0.006). The mtDNA sequence of the Littoral-Dinaric donkeys were also divided into two phylogenetic groups (haplo-group C, sequence LDD 01, LDD 02, LDD 03, LDD 04, and LDD 05; haplo-group D, LDD 06, and LDD 07), and the genetic distance between haplo-group C and D is moderate (d = 0.039). The Littoral-Dinaric donkey haplo-group D vs. haplo-group C are closer to the Istrian donkey haplo-group A and haplo-group B (d = 0.018 vs. d = 0.046).   Two phylogenetic groups of the Istrian donkey (haplo-group A and haplo-group B) are phylogenetically closely related to the Ragusano donkey (d = 0.006, d = 0.009), and the distance from the Romagnolo donkey, Asinara donkey, and Amiata donkey increases (d = 0.014-0.020). The highest distance is between the Istrian donkey and the Sardo donkey /Martina Franca donkey (d = 0.029, d = 0.034). Between two phylogenetic haplo-groups of the Littoral-Dinaric donkey, haplo-group C is more distant from the Italian donkey breeds (d = 0.026-0.044) than haplo-group D (d = 0.020-0.031). The Istrian donkey is also phylogenetically more closely related to the Catalonian donkey than to the Littoral-Dinaric donkey (d = 0.029-0.033 vs. 0.042-0.061).

Discussion
Breeding consolidation activities of the Littoral-Dinaric donkey and Istrian donkey populations were initiated at the end of the 20th century. In the first phase, an inventory of the remaining donkey population in Croatia was conducted, as well as research on phenotype. Once the breed standards were established, animals that met the requirements were entered into the Pre-Book during the inventory process over the next decade. In 2010, a main breed studbook was established, in which the animals from the Pre-Book were entered. In this way, the populations were systematically divided into two breeds after carefully checking the breeding status of each animal.
Positive population trends can be observed in the donkey populations studied (Littoral-Dinaric donkey, Istrian donkey), as the number of males, females, and registered offspring increases (Figure 1). The number of donkey breeders per breed has also increased over the past five years by 22% in the Istrian donkey population (from 129 to 158 herds) and 29% in the Littoral-Dinaric donkey population (from 628 to 810 herds). The age structure of the two breeds is also favorable, and about 3 4 of the population is less than ten years old ( Figure 2). The observed population trends are the result of a polyvalent and careful approach to breed conservation. After the start of the program for the protection of breeds, recognizing the importance of educating the public about the benefits of preserving donkey breeds initiated activities to promote the preservation of donkey breeds at local, regional, and national levels. After two decades of activities for the promotion and economic reaffirmation of donkeys, they are more valued animals, and are involved in different activities (leisure, agrotourism, folklore, exhibitions, etc.), economic programs, through the production of donkey milk and meat (people on the Istrian peninsula traditionally consume equine meat), and ecoservice function (biodiversity conservation). The female/male ratio was 3.1/1 in the Littoral-Dinaric donkey and 4.5/1 in the Istrian donkey populations, respectively. The female/male ratio was similar in Pêga donkeys (4.97/1) [9], Miranda donkeys (3.63/1) [8], and Amiata donkeys (2.90/1) [6], but higher than in Martina Franca donkeys (1.48/1) [7]. The reason for the difference in ratio between females and males is in the herd size and breeding technology.
In the Istrian donkey and Littoral-Dinaric donkey populations, the greatest number of traced generations was six and five, respectively, and the average maximum of complete and equivalent generations were 1.80, 1.24, 0.88, and 1.01, 0.74, 0.56, respectively. In the Catalonian donkey, the number of complete and equivalent generations was 1.23 and 1.96 [5]. In Amiata donkeys, the average maximum of complete and equivalent generations were 1.4, 0.53, and 0.78, respectively [6]; in Asinina de Miranda donkeys, 0.33, 0.22, and 0.28, respectively [8]; and in Andalusian donkeys, 1.09, 0.52, and 0.75, respectively [10]. Some donkey breeds have a much deeper pedigree, such as the Martina Franca donkey, for which the average maximum, complete and equivalent generations were 4.67, 1.97, and 3.01, respectively [7]. In Littoral-Dinaric donkeys and Istrian donkeys, pedigree depth has been relatively modest, but has increased in recent years, and is becoming an increasingly reliable tool for population management. Giontella et al. [49] suggest that knowledge of relationships between animals is essential for the genetic management of breeds, and represents one of the principal tools for optimizing their conservation strategies.
Generation intervals are important factors of population management measures [50]. The GI in the Littoral-Dinaric donkey population was lower than in the Istrian donkey population (7.27 vs. 7.73 years) and was similar in the Andalusian donkey (7.34) [10]. The GI for the four different pathways were similar. In Catalonian donkeys [51] and Amiata donkeys [6], lower GI were observed (6.74 and 6.65 years), whereas higher GI were observed in the Asinina de Miranda donkey (8.18 years) [8], the Martina Franca donkey (9.09 years) [7], and the Pêga donkey (10.70 years) [9] populations. The lower GI in the Littoral-Dinaric donkey population is a consequence of the population increasing more quickly and donkeys entering into reproduction early (three to four years). Folch and Jordana [51], and Cecchi et al. [6] observe higher GI in maternal pathways than in paternal pathways (7.32/6.16; 7.0/5.9) in the Catalonian and Amiata donkey populations, but the GI was longer in the sire-offspring pathways than in the dam-offspring pathways [7]. Prolonging the generation interval may be useful to increase the number of males and females selected for breeding, which would thereby progressively increase the effective population size, which is inversely proportional to the rate of inbreeding [52,53], and therefore, it could be a useful tool to maintain the genetic diversity of the population. The average age of parents at the birth of their offspring was similar in the Littoral-Dinaric donkey and Istrian donkey populations (7.59 and 7.50 years, respectively) and lower than in the Asinina de Miranda population (9.32 years) [8]. Navas et al. [10] observed in the Andalusian donkey populations a slightly higher mean age of parents in the dams-sons pathway (8.23 years) than dams-daughters pathway (7.84 years).
The effective number of founders and ancestors in the Littoral-Dinaric donkey population was relatively high (325, 316), which is typical for populations that do not have long historical pedigree records [54]. The effective number of ancestors and founders in the population for Istrian donkeys was lower (70, 65), and the number of ancestors explaining 50% of the genetic variability was 23. In the Catalonian donkey, the effective number of founders and ancestors was 70.6 and 27, respectively [5], in the Amiata donkey 114.9 and 42 [6], in the Martina Franca donkey 22 and 18 [7], in the Asinina de Miranda 38 and 34 [8], and in the Andalusian donkey 153.2 and 142 [10]. The rate of the effective number of founders compared with the effective number of ancestors can be used to determine the bottleneck in the population. If the ratio is 1, the population is stable. A larger ratio reflects a stronger bottleneck effect [32]. The f e /f a ratio in the Littoral-Dinaric donkey and Istrian donkey populations was 1.03 and 1.08, respectively, which is more favorable compared to the Catalonian donkey (2.61) [5]. In the Amiata, Andalusian, Asinina de Miranda and Martina Franca donkeys, the f e /f a ratio is similar (0.85, 1.08, 1.12, 1.22) [6][7][8]10]. The high f e /f a ratio suggests a disproportionate use of some breeding animals, presumably stallions, resulting in a loss of genetic diversity compared with that expected under random mating conditions [55]. The results of the current study did not indicate a bottleneck, so the Littoral-Dinaric donkey and the Istrian donkey populations were stable.
The trend and level of the inbreeding coefficient (F) in the population are crucial for maintaining the genetic diversity of the breed. The level of the inbreeding coefficient in the Littoral-Dinaric donkey population was lower (0.99%) than in the Istrian donkey population (1.77%), but the trend of the inbreeding coefficient is more favorable in the Istrian donkey population (in the last five years the level of F is stable) than in the Littoral-Dinaric donkey population (the level of inbreeding has increased in the last decade). The inbreeding coefficients of the Catalonian and Martina Franca donkeys were higher (3.36; 3.95) [5,7] and those of the Asinina de Miranda, Amiata, and Andalusian donkeys (0.08, 0.29, 0.70) were lower [6,8,10]. The inbreeding coefficients increased with the knowledge of the pedigree, but the poor pedigree completeness levels result in the underestimation of the inbreeding level [56], and therefore, the values found must be carefully analyzed [8]. The observed F value was followed by a low AR value in the Littoral-Dinaric donkey and Istrian donkey populations (0.13, 1.10), indicating a reduced representation of each individual in the whole population. Donkeys with the lowest AR, used as stallions and mares, can reduce inbreeding, and balance the gene contributions of the founders in the population, and thus genetic variability [57]. In the Catalonian and Martina Franca donkeys, the AR coefficients were higher (3.76, 7.35) [5,7], and in the Amiata, Andalusian and Asinina de Miranda donkeys, they are higher (0.94, 0.81, 0.33) [6,8,10].
The average GCI in the Littoral-Dinaric and Istrian donkey populations were relatively low (1.76, 2.42). In general, the trend in GCI in the Littoral-Dinaric donkey and Istrian donkey populations increased over time. Donkeys with higher GCI values have a greater balance in the number of founders and are thought to contain genes transmitted by the founders. From this point of view, animals with higher GCI values are crucial for breed conservation and the genetic variability they possess. Quaresma et al. [8] observed lower values for a GCI (1.30) in the local Portuguese donkey breed Asinina de Miranda and found that the average GCI was negatively associated with the age of the donkeys.
The  [14], four Chinese breeds (range from 6.0 to 6.80) [58], or donkey breeds from the Balkans (9.3) [19], but they are higher than two Sicilian native breeds (4.05) [59]. Comparing observed and expected heterozygosity, the values of Istrian donkey and Littoral-Dinaric donkey populations are mostly between the reported values of previous studies [20,58], with averages of H O = 0.560 and H E = 0.70, and H O = 0.63-0.71, H E = 0.93-0.70. Although, results of heterozygosity for Croatian donkey breeds were higher than those for Sicilian donkey [59] and Italian donkey [16] breeds, but lower than the heterozygosity in the Balkan donkey breed research [19]. Previous research results by Ivanković et al. [13], which included both breeds, showed slightly higher values for H E (H E(ID) = 0.68, H E(LDD) = 0.70) than the results from this study. Comparing the observed and expected heterozygosity of the Croatian donkey populations with 59 horse breeds (included warmblood, coldblood and pony breeds) [60] found that it was below the average for all breeds (H E = 0.80, H O = 0.70).
The results of the preserved genetic diversity in both native breeds are the large breeding area and the geographical distribution of the breeds, especially in the case of the Littoral-Dinaric donkey. In case of the Littoral-Dinaric donkey, its breeding area is Dalmatia and Dalmatian Zagora (12,200 km 2 ), whereas the breeding area for the Istrian donkey is in the Istrian peninsula (3500 km 2 ). Both breeds have a wide genetic base, which is the result of historical migration and the purchasing of animals between villages that contributed to the exchange of genetic material.
The These results indicate that there are no appreciable differences in the level of genetic variability among the Croatian donkey breeds. A comparison of genetic diversity using microsatellites in animal species can be generally biased due to the different marker sets used, but they still represent the genetic structure of populations very well; therefore, these results are comparable to those of previous studies. Littoral-Dinaric donkey and Istrian donkey populations have a lower than average number of alleles compared with five local Spanish breeds (range from 7.0 to 7.5) [14], four Chinese breeds (range from 6.0 to 6.80) [58], and donkey breeds from the Balkans (9.3) [19], but they are higher than two local Sicilian breeds (4.05) [59]. Microsatellite genotyping results revealed a relatively high degree of heterozygosity in the Croatian donkey populations included in this study. Comparing the observed and expected heterozygosity, the values of the Littoral-Dinaric donkey and Istrian donkey populations are mostly between the reported values of previous studies, such as those reported by Di et al. [58] (average H O = 0.560 and H E = 0.70) and Yatkın et al. [20] (H O = 0.63-0.71, H E = 0.93-0.70). Although, results of heterozygosity for Croatian donkey breeds were higher than those for Sicilian donkeys [59] and Italian donkeys [16], but lower than the heterozygosity in the Balkan donkey breed [19]. Previous research results by Ivanković et al. [13] which included both breeds, showed slightly higher values for H E (H E(LDD) = 0.70, H E(ID) = 0.68) than the results of this study. A very low percentage of homozygous individuals in the Istrian donkey population was confirmed by a low, non-significant inbreeding coefficient (F IS(ID) = 0.029). A small excess of heterozygotes (F IS(LDD) = −0.002) was observed in the Littoral-Dinaric population. For both populations, this outcome is expected, and can be explained with several reasons. A large breeding area, and thus, a large geographical dispersion of breeds, especially in the case of the Littoral-Dinaric population (see File S1), contributes to a lower possibility of mating between relatives, and consequently, to the preservation of greater genetic diversity. For decades, residents followed the practice of purchasing (or exchanging) male and female breeding animals of the same breed from other villages rather than from neighboring households, which helped control inbreeding. Reports of the inbreeding coefficient in donkey populations from previous studies vary considerably, from 0.09 to 0.225 [20,58]. The large ancestral gene pool, which has a very large number of founders, lack of selection, and planned mating (avoiding mating of relatives), have undoubtedly contributed to the preservation of genetic variability in Croatian local donkey breeds.
The high rate of gene flow between the studied populations was confirmed by the high number of migrants (N m = 12.431), which reduces genetic differentiation among populations. These results are supported by the low genetic differentiation between Littoral-Dinaric donkey and Istrian donkey (F ST = 0.028) populations, suggesting that 97.2% of the total genetic variation resulted from genetic differentiation within breeds. This value is almost identical to the F ST = 0.0212 from the previous research of Ivanković et al. [13] Low genetic differentiation for the Spanish donkey population was reported by Aranguren-Méndez et al. [14] and for the Turkish donkey population by Yatkin et al. [20]. The weak differentiation between the two donkey breeds was highlighted by a population assignment where 93% of individuals were correctly assigned to the true breed of origin. It is likely that the retention of ancestral variation contributes to a large gene pool (populations spread from centers of origin possess larger genetic diversity), rather than weak differentiation between breeds.
Although the population size of the Littoral-Dinaric donkey and Istrian donkey breeds have declined significantly throughout history (in 1937 there were~39,000 animals) [61], these populations have not experienced a bottleneck. It is not unusual to reveal the absence of bottleneck events despite the reduction in population size [62].
The results of the mtDNA study show genetic variability within and between the two donkey breeds in Croatia, but no geographic clustering was observed. Cozzi et al. [28] observed a similar situation in six Italian donkey breeds. Genetic variability within and between breeds was not reflected in their geographical clustering. Phylogenetic analysis based on mtDNA classified two Croatian donkey breeds into the group Equus asinus africanus. Moreover, in other studies concerning the mtDNA D-loop, the donkey breeds included in the study are also in the phylogenetic group Equus asinus africanus, which was established by other authors [25,28,29,31]. All identified haplogroups of the Istrian donkey population (A, B), and one of the Littoral-Dinaric donkey population (D), are close to the Italian donkey breeds, which is important for understanding the breed's phylogenetic position and further breeding consolidation. For the future program of preserving the uniqueness of two Croatian local donkey breeds, this observation is important, especially the knowledge of phylogenetic closeness with some Italian donkey breeds; however, when reconstructing phylogenetic relationships within and between breeds using mtDNA, we must consider that mtDNA is transmitted from mother to offspring, and therefore, it can only give insights in maternal lineages. This might be discordant with nuclear DNA, and in some cases provide a biased perspective on the evolutionary history of the breeds (species).

Conclusions
The preservation of genetic diversity is one of the primary tasks of the program for the protection of local donkey breeds, especially when they are threatened with extinction. After the biological survival of the Littoral-Dinaric donkey, and when the Istrian donkey was threatened, the launch of a program for inventorying, and a breeding consolidation of the remaining populations two decades ago, led to a biological recovery of the populations. Breed indicators observed in a complex study that the pedigree and genetic structure of donkey breeds point to certain population trends, but also to risks that need constant monitoring. Weaknesses in pedigree indicators due to lower pedigree depth were partially verified using genetic structure indicators. In the Littoral-Dinaric donkey population, the level of inbreeding and average relatedness is lower than in the Istrian donkey population, but the trends are less favorable, which requires constant monitoring. If it is necessary to slow down the growth of inbreeding, the equal contribution of animals to reproduction should be ensured, especially those that are less represented in the pedigree, so that the relatedness between mated animals is minimized, the generation interval is extended, and so on. Analysis of pedigrees and microsatellites has not revealed any bottleneck in the populations. Regular monitoring should be continued, precisely to optimize further development of donkey breeds. It is certainly necessary to promote the breeds and develop new sustainable strategies for their conservation (e.g., through tourism, asinotherapy, milk production (and meat), ecoservice function, and so on). Only a complete conservation program for local breeds is sustainable in the long run.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available to preserve privacy of the data.