Morphological and genetic characterization of the Graciosa donkey breed

ABSTRACT

The Graciosa donkey has been associated with husbandry and social activities in Graciosa Island for decades. Although donkeys played an essential role in the development of the island, their numbers are severely reduced. The Graciosa donkey is characterised by its small stature, strong bony constitution and strong hooves. They are extremely gentle, patient and submissive fully integrated into the island's ecosystem. As so they must be considered as elements of an environmentally sustainable ecosystem. This study aimed to define the morphotype of the Graciosa donkey and to access its genetic diversity using microsatellite markers. Biometric data consisted of 15 morphological traits from 39 individuals. The Graciosa donkey is dolichocephalic, hypermetric according to the dactilothoracic index, hypometric regarding the approximate live weight and far from ground. Twenty-four individuals, which showed no kinship, based on results of 13 markers were considered as founders and were genotyped for two additional markers and further analysed in this study. All microsatellites used were polymorphic, but low levels of genetic diversity were identified. As the Graciosa donkey gains popularity, the information obtained from this study will be very useful for conservation and management purposes by safeguarding its genetic diversity.

KEYWORDS:

• Equus asinus
• Azores
• morphotype
• genetic diversity
• conservation

Introduction

Donkeys were crucial in the development of civilizations worldwide, but it was the livestock species most affected by the industrial revolution as they became, in most countries, obsolete as pack animals and in agriculture (Kugler et al. 2008).

In the Azores archipelago, donkeys played an essential role in the settlement of the islands, starting in the fifteenth century. The donkey was the domestic species best adapted to Graciosa’s features, characterized by small plots of land of difficult access where the donkey is easier to handle. The small land area, small altitudes with sprayed properties and the scarce resources of the inhabitants at that time provided all the conditions for the use of this population of donkeys across the island. Besides helping in agricultural practices, they were also used to carry food for livestock, agricultural products, wood, water and milk, in addition to transporting men, women and children, and were often companions in leisure activities (Mendonça 2005). Due to an accurate knowledge of the animals and of the local conditions, breeding focused on the more rustic, resistant, proportional and balanced individuals, resulting in harmonious animals apt for work. The breeding was so successful that these donkeys were supplied to the surrounding islands, mainly Faial and São Jorge, where they were much appreciated for the transport of milk, as they walk regularly without sudden movements (Bettencourt 1950). Although 1103 donkeys have been reported as existing on the island in 1926 (Bettencourt 1950), numbers have gradually declined with the introduction of bicycles and motor bikes and agricultural mechanisation (Ferreira 1987) to an astonishing 26 individuals reported in 2001 (Kugler et al. 2008). The increasingly limited utility of the Graciosa donkey as farm animals has compromised its future sustainability, which now depends on a shift towards new market demands, such as recreation, tourism, therapies and milk production. It is noteworthy that nowadays most owners are farmers around 70 years old that still possess donkeys for their usefulness as draught and pack animals. However, the attitude towards these animals is changing and some are being maintained as companion animals (Mendonça 2005).

The Graciosa donkey shows typical phenotypic characteristics of the extant wild African ass subspecies the Equus asinus africanus (convex profile, dolichocephaly, less than 1.20 m at the withers and clear coat with a shoulder stripe on the spine) (Supplementary Figure 1), although some present zebra stripes on the front limbs typical of the E. asinus somaliensis (Kugler et al. 2008). The Graciosa donkey is an extremely gentle, patient and submissive animal. They have been known for decades as the Dwarf Graciosa donkey, a confusing nomination. The origin of this name comes from the fact that inhabitants from the Azorean islands used to call donkeys to horses and dwarf donkeys to donkeys, but they are not true dwarfs, because, according to the American metrics of donkeys, they are in the range of standard donkeys (lovelongears.com).

The population of Graciosa donkeys was never subjected to a systematic breeding or selection program. Presently they are just kept by private enthusiasts, which gathered as an association in 2013 to preserve the existing individuals and increase their numbers – the Association of Breeders and Friends of the Graciosa Dwarf Donkey (ACABAIG). According to data from February 2021 provided by the ACABAIG, a total of 35 individuals were registered in the Studbook of the Graciosa donkey and 4 were in the process of being registered. Despite being an endangered breed, these animals are not covered by any governmental support as a farm animal and no conservation actions have been implemented, meaning that there is no state assistance to establish a studbook or register the individuals, and no support in the case of infectious diseases or epidemic outbreaks. Nowadays maintaining the Graciosa donkey is dependent upon enthusiasts, but hobby keeping leads to different breeding purposes in character, type and biological attributes than breeding for agricultural purposes (Kugler et al. 2008), as breeding will depend on the owners’ preferences.

This study aimed to define the morphotype of the Graciosa donkey and to study its genetic diversity, as the still-existing individuals should be protected not only as genetic stock but also as a cultural heritage. Also the knowledge of morphometric measurements in donkeys together with genetic data information is of great importance for the preservation of genetic diversity and can assist in implementing conservation and breeding plans, while preventing genetic erosion.

Materials and methods

Ethics statement

In this study, all sample collection procedures (measurements for the biometric study and blood for the genetic study) were carried out in strict accordance with the Portuguese decree-law N113/2013 from August 7th, Chapter I, Article 2, paragraph 7, subparagraphs a, b, d, e, f. It was also in accordance with Directive 2010/63/EU of the European Parliament and of the Council of 22 September 2010 on the protection of animals used for scientific purposes, Chapter I, Article 1, point 5 (b and f), and therefore an ethical review by the Statement Animal Experiment Committee was not required.

Biometric study

It has been suggested that the capacity of performance of donkeys can be assessed by the description of their morphological characteristics and that draft power is directly proportional to size parameters (Nengomasha et al. 1999).

To access the morphotype of the Graciosa donkey 39 adult Graciosa donkeys (>4 years old), 19 jacks and 20 jennies were measured, in farm, for the biometric study. A total of 13 measurements (Figure 1) were taken, with zoometric stick and tape, from the left side of the animal while standing on a horizontal surface in standard position, except for measurements of the head. Biometric data were collected by a veterinarian following good veterinary practices. Additionally, chest depth was obtained by subtracting the sub-sternal flank height from the height at withers. Approximated live weight (W) was estimated from body measurement data as defined by Pearson and Ouassat (1996), as W = TC^2.65/2188. These authors determined that the most appropriate variable for estimating the live weight of donkeys is the thoracic perimeter (TC).

Figure 1. Measurements taken to the Graciosa donkey for the biometric study.

To determine if samples were generated from a population with a normal distribution, the Shapiro–Wilk’s W test was performed using the Normality Test Calculator (Georgiev, online) followed by Levene’s test for the determination of the homogeneity of variances. Variables that fulfilled the conditions for parametric were compared between genders through the t-test implemented in MS Excel software. For the remaining variables the Mann–Whitney U-test was used. Levene’s test and Mann–Whitney U-test were calculated using the calculators available at www.socscistatistics.com. For each morphological variable mean, standard deviation, median, maximum and minimum values (Table 1), and number of individuals that fall in each quartile (Supplementary Table 1) were assessed with Microsoft Excel software.

Table 1. Maximum, minimum, mean and standard deviation of body measurements and live weight of jennies and jacks of the Graciosa Donkey.

Measurements		Jennies					Jacks
		Mean ± SD	Median	Max	Min	CV (%)	Mean ± SD	Median	Max	Min	CV (%)
Heights (cm)	Withers (WH)^t	106.55 ± 3.88	107	115	99	3.64	106.53 ± 4.27	106	116	99	4.01
	Croup (CH)^t	111.35 ± 4.23	110.50	120	104	3.80	110.05 ± 4.86	111	120	102	4.42
	Sub-sternal flank (SSF)^t	58.74 ± 3.88	58	66	50	6.61	59.28 ± 2.88	58.5	65	54	4.86
Depths (cm)	Chest (CD) (calculated)^U	47.53 ± 3.47	47	58	42	7.30	47.22 ± 2.92	47	52	43	6.18
Lengths (cm)	Head (HL)^t	48.55 ± 2.36	48.50	53	43	4.86	48.21 ± 1.96	48	53	45	4.07
	Ears (EL)^U	25.85 ± 1.74	25	30	24	6.73	25.00 ± 1.59	25	29	23	6.36
	Back (BL)^t^,*	110.00 ± 5.18	108.50	123	102	4.71	105.26 ± 6.96	104	122	94	6.61
	Croup (CL)^U^,*	36.03 ± 2.04	36	39	31.50	5.66	33.95 ± 3.58	35	39	24	10.54
Widths (cm)	Head (HW)^U	20.30 ± 0.78	20	22	19	3.84	20.21 ± 1.06	20	22	18	5.24
	Bi-iliac (BIIW)^t^,*	33.15 ± 2.24	33.50	37	28	6.76	31.05 ± 3.33	30	38	26	10.72
	Bi-ischial (BISW)^t^,**	26.25 ± 3.42	27	33	20	13.03	21.95 ± 3.27	22	28	15	14.90
Circumferences (cm)	Forelimb cannon bone (FCBC)^U	14.45 ± 0.59	14	16	14	4.08	14.79 ± 1.00	15	17	13	6.76
	Carpus (CaC)^U^,*	23.40 ± 1.39	23	28	21	5.94	24.05 ± 1.47	24	27	21	6.11
	Thorax (TC)^t^,*	127.00 ± 6.44	127	140	115	5.07	122.42 ± 9.26	122	142	106	7.56
Approximate live weight (W) (kg) (calculated)^U^,*		172.77 ± 23.22	171.80	222.43	132.07	13.79	157.84 ± 32.32	154.46	230.95	106.42	20.48

^t t-test.

^U Mann–Whitney U-test.

*Significance level (p < 0.05).

**Significance level (p < 0.01).

Based on these 15 morphological traits, body indexes and functionality indexes (Table 2) in which the relations between some elements of width, height, length and weight were evidenced, were calculated to evaluate the proportions of the individuals and deduce proportions that indicate functional aptitudes. These indices also allow comparison between breeds and studies in terms of proportions of body segments, regardless of length differences (Komosa and Purzyc 2009). Values were statistically significant when P-values were < 0.05. Pearson correlation coefficient (R) analysis was used to identify possible associations between the morphological characters for genders separately with StatPlus:mac, AnalystSoft Inc. (statistical analysis program for macOS. Version v8.). For this study, a strong positive correlation was considered when r > 0.50, a strong negative correlation was considered when r < −0.50. Additionally, to evaluate the proportion of variance of a variable in another variable, the coefficient of determination (R²) was calculated.

Table 2. Mean values for Conformation and functional indices determined based on the measurements taken for jennies and jacks of the Graciosa donkey.

Conformation and functional indices	Jennies	Jacks
	Mean ± SD	Mean ± SD
Chest index (CI) = (CD-SSF)^a^,^U	<0	<0
Cefalic index (CfI) = (HW/HL)^b^,^c^,^U	0.42 ± 0.01	0.42 ± 0.02
Pelvic index (PI) = (BIIW/CL)^b^,^U	0.92 ± 0.05	0.92 ± 0.11
Dactilothoracic index (DTI) = (FCBC/TC)^a-c^,^U^,**	0.11 ± 0.01	0.12 ± 0.01
Weight in cannon index (WCI) = (FCBC/W)*100^d^,^U^,*	8.49 ± 1.02	9.76 ± 2.02
Conformation index or Index 1 (CoI or I1) = (TC/WH)^c^,^e^,^f^,^t^,**	1.19 ± 0.04	1.15 ± 0.06
Body ratio or Index 2 (BR or I2) = (WH/CH)^c^,^f^,^t^,**	0.96 ± 0.01	0.97 ± 0.02
Index 3 (I3) = (SSF/WH)^c^,^U	0.55 ± 0.03	0.56 ± 0.02
Index 4 (I4) = (FCBC/SSF)^c^,^U	0.25 ± 0.01	0.25 ± 0.02
Chest height index (CHI) = (SSF/TC)^c^,^U^,**	0.46 ± 0.02	0.49 ± 0.03
Compactness index (CpI) = W/WH^e^,^t	1.62 ± 0.18	1.48 ± 0.26
Back load index at walk (BLW) = (TC^2*98)*(WH/100)^e^,^t^,*	148.53 ± 11.81	138.27 ± 17.49
Back load index at canter (BLC) = (TC^2*56)*(WH/100)^e^,^t^,*	84.87 ± 6.75	79.01 ± 9.99

^a Following the description by Martin-Rosset (1983).

^b Following the description by Martin (2006).

^c Following the description by Fonseca Jiménez et al. (2016).

^d Following the description by Cabral et al. (2004).

^e Following the description by Mariz et al. (2014).

^f Following the description by MacManus et al. (2008).

^t t-test.

^U Mann–Whitney U-test.

*Significance level (p < 0.05).

**Significance level (p < 0.01).

Genetic study

Fresh tissues are ideal sources for DNA extraction by common methods, but blood is one of the least invasive ones. Hairs are another source for DNA extraction, but the yield of DNA can be very low and quickly degrade after extraction (Xiangyu and Zhang 2006). As the ACABAIG intends to create a DNA genebank for future analysis and with different technologies, blood was chosen as the DNA source. Blood samples were collected by a veterinarian by jugular venipuncture to refrigerated sterile tubes with EDTA, following good veterinary practices and is a non-experimental, agricultural practice.

DNA was extracted from fresh blood following standard protocols (Miller et al. 1988).

In a preliminary study for the recognition of the Graciosa donkey as an autochthonous breed the same 39 individuals were genotyped using the 13 microsatellite markers of the core panel recommended for E. asinus genotyping by the International Society for Animals Genetics (ISAG 2017): AHT4 (Binns et al. 1995); ASB23 (Breen et al. 1997); HMS2, HMS3, HMS6, HMS7 (Guérin et al. 1994); HMS18 (Godard et al. 1997); HTG7, HTG10 (Marklund et al. 1994); TKY297, TKY312, TKY337, TKY343 (Tozaki et al. 2001).

After a preliminary analysis, 24 (14 jacks and 10 jennies) of these 39 individuals showed no kinship based on the results of the 13 microsatellite makers instead of declared parentage-offspring relationships and were therefore considered as founders and further analysed in this study. Due to low levels of genetic diversity detected with the ISAG core panel recommended for E. asinus, these 24 individuals were further genotyped for two additional microsatellites from the ISAG core panel for horse genotyping, namely AHT5 (Binns et al. 1995) and HTG6 (Ellegren et al. 1992).

Polymerase chain reaction amplification was performed in 20 µl volumes containing 20 ng of total DNA, 200 µM of each dNTP, 2–5 pmol of each primer and 1 U Taq DNA polymerase (Thermofisher) in reaction buffer containing MgCl₂. Forward primers were fluorescence labelled at the 5’ end with FAM, VIC, NED or PET, allowing multiplexing and simultaneous separation of amplified products. Reactions were conducted in a UNO II Biometra thermocycler as described by Lopes et al. (2015). The PCR products were size fractionated by capillary electrophoresis using an automated sequencer (SeqStudio Genetic Analyzer, Thermofisher) and fragment lengths were determined with the help of the internal size standard GeneScan 500 LIZ (Applied Biosystems). Alleles were scored according to ISAG nomenclature.

Genetic variability was measured by estimating total number of alleles (TNA), effective number of alleles (Ne), observed (Ho), expected (He) and unbiased expected heterozygosities (uHe), fixation index (F) and deviations from Hardy–Weinberg equilibrium with GenAlex 6.5 (Peakall and Smouse 2012). Polymorphism information content (PIC) and probability of occurrence of null alleles (r) were determined with Cervus 3.0.3 (Marshal et al. 1998). Samples in which a single allele per locus was detected were considered as homozygous genotypes, instead of heterozygous with a null allele, to compute genetic diversity parameters. The software Identity (Wagner and Sefc 1999) was used to calculate the probability of paternity exclusion (PE). A principal component analysis (PCA), to condense the genetic variation revealed with the 15 microsatellite loci by individual genotypes, was performed with GenAlex 6.5. Five samples of Terceira Pony genotyped for the same markers were used as outgroups.

Results

Biometric study

Mean, standard deviation, median, maximum, minimum and coefficient of variation of morphological traits from the Graciosa donkey jennies and jacks are presented in Table 1. All morphological traits analysed were homogeneous with CV ≤ 20.48%. The distribution of samples per quartile is presented in Supplementary Table 1. Table 2 shows the conformation indexes determined based on the previously described measurements. A matrix with Pearson’s correlation coefficients, coefficient of determination and p-values for each pair of variables are presented in Supplementary Table 2. Correlation results of jacks are presented above the diagonal and of jennies below the diagonal. Results indicated that W has the strongest positive and significant correlation with the other biometric traits for both genders. In general, jennies have a greater number of strong and significant correlations compared to jacks. For jennies, WH, CH, CL and CaC presented highly and significant correlations with most variables. For jacks CH, WH and BL presented highly and significant correlations with most variables. On the other side, CD for jennies showed no strong correlations with the remaining variables. Phenotypic correlation results of jennies for CD did not present any significant or strong correlation. Also, BIISW only showed strong correlation with TC (r = 0.62), W (r = 0.608) and SSF (r = 0.556). Contrarily to jennies, for jacks FCBC and BIISW presented strong correlations only with CaC (r = 065 and r = 0.571, respectively) and CL presented strong correlations with CaC (r = 0.552) and BIISW (r = 0.521). No negative high correlations were detected.

As expected, the correlation coefficient of W (r = 0.99) with TC corresponds to a coefficient of determination of 1, as the approximate live weight of donkeys is estimated based on the thoracic perimeter.

Genetic study

For the 13 microsatellite markers of the core panel recommended for Equus asinus by ISAG a total of 48 alleles were identified, ranging from 2 for ASB23 and HMS7 to 6 for HTG7, and with a mean of 3.69 alleles per locus. Nevertheless, the least informative locus for this population is HMS6 with the lowest values obtained for Ne, I, Ho, He, uHe and PIC, and the most informative microsatellite marker is HTG7 with the highest TNA, Ne, Ho, He, uHe and PIC values (Table 3). The mean observed heterozygosity was slightly higher than the expected heterozygosity (0.471 against 0.437, respectively). Three out of the 13 markers revealed observed heterozygosity values lower than the expected ones (AHT4, HMS6 and HMS7), with HMS6 showing significant (p < 0.05) deviation from Hardy–Weinberg equilibrium with the corresponding F_is value of 0.464. By using the donkey core marker set the obtained PE was low (97.60%) indicating that for parentage testing for this breed the number of markers has to be increased.

Table 3. Descriptive statistics of the genetic diversity of the Graciosa donkey. Total number of alleles (TNA), number of effective alleles (Ne), information index (I), observed heterozygosity (Ho), expected heterozygosity (He), unbiased expected heterozygosity (uHe), fixation index (F), polymorphism information content (PIC) and estimated frequency of occurrence of null alleles (r).

Locus	TNA	Ne	I	Ho	He	uHe	Fis	PIC	r
AHT4	4	2.756	1.165	0.583	0.637	0.651	0.084	0.583	0.033
ASB23	2	1.332	0.415	0.292	0.249	0.254	−0.171	0.218	−0.034
HMS2	3	1.567	0.608	0.458	0.362	0.370	−0.266	0.310	−0.071
HMS3	5	2.648	1.158	0.625	0.622	0.636	−0.004	0.552	−0.002
HMS6	3	1.184	0.334	0.083	0.155	0.159	0.464	0.148	0.062
HMS7	2	1.229	0.334	0.125	0.187	0.191	0.330	0.169	0.052
HMS18	3	2.254	0.883	0.583	0.556	0.568	−0.048	0.457	−0.017
HTG7	7	3.728	1.480	0.875	0.732	0.747	−0.196	0.686	−0.083
HTG10	5	1.488	0.710	0.375	0.328	0.335	−0.143	0.313	−0.036
TKY297	4	1.296	0.503	0.250	0.228	0.233	−0.095	0.219	−0.018
TKY312	3	2.612	1.024	0.708	0.617	0.630	−0.148	0.544	−0.056
TKY337	4	1.602	0.712	0.458	0.376	0.384	−0.219	0.340	−0.060
TKY343	3	2.698	1.039	0.708	0.629	0.643	−0.126	0.553	−0.048
Mean^+	3.692	2.030	0.797	0.471	0.437	0.446	−0.041	0.392	−0.021
AHT5	8	4.939	1.780	0.864	0.798	0.816	−0.083	0.772	−0.037
HTG6	5	1.760	0.897	0.478	0.432	0.442	−0.107	0.408	−0.032
Mean^++	4.067	2.206	0.870	0.498	0.461	0.471	−0.048	0.418	−0.023

⁺ Mean values for the 13 microsatellite markers of the core panel recommended for Equus asinus by ISAG.

⁺⁺ Mean values for all 15 loci analysed.

To increase the accuracy of parentage testing and to detect higher levels of genetic diversity within the Graciosa donkey breed, two additional markers from the core panel recommended by ISAG for E. caballus genotyping, namely AHT5 and HTG6, were used to genotype the same animals. With the two additional markers, a total of 61 alleles were identified. The most informative marker detected was now AHT5 with the highest values for TNA, Ne, He uHe and PIC. The highest observed heterozygosity still was detected at HTG7. The addition of these two markers increases the PE from 97.60% to 99.30%, providing more confidence when performing parentage testing.

The PCA showed two distinct groups, one comprising all Graciosa donkey individuals with no clear evidence of sub-grouping at the left side of the plot, and another clustering all Terceira Ponies at the right side of the plot (Figure 2). The first three principal axes of the PCA accounted for 46.17% of the total variation of the two-dimensional projection. The first principal component accounted for 37.87% of the underlying variation and the second condenses 8.30% of the variation.

Figure 2. Two-dimensional projection of the principal component analysis of the 24 Graciosa donkey (GD) and five Terceira Pony (TP) individuals based on 15 microsatellite markers along the first two principal axes.

Discussion

The results obtained from the 15 morphological traits and the calculated indexes are important for the characterization of a breed and to define its taxonomic affiliation (Komosa and Purzyc 2009; Ayad et al. 2019). The present study indicated that from a morphological and zoometric point of view, the Graciosa donkey is: dolichocephalic [CfI < 1; (Martin 2006; Fonseca Jiménez et al. 2016)]; hypermetric [DTI = 0.11–0.12; (Mariz et al. 2014)] as a result of a strong bony constitution and therefore suitable for traction; hypometric regarding the approximate live weight [W < 350 kg; (MacManus et al. 2008)] and ‘far from ground’ [CI < 0; (Ribeiro 1989; MacManus et al. 2008)]. They are well proportional as I2 values are close to 1 (I2 = 0.96 for jennies and 0.97 for jacks), and I4 values are 0.25 for both jennies and jacks. An animal is considered proportional when the height at withers is equal, or similar to the height of the croup (I2) and when the height to the sternum is four times the forelimb cannon bone circumference (I4) (MacManus et al. 2008).

The Graciosa donkey is suitable for work supported by indexes I3 (at the upper limit of range) and WCI (8.49 for jennies and 9.76 for jacks). Index 3 is an indicator of the proportionality of individuals and their relationship with the aptitude for work (Martin 2006). The WCI relates the perimeter of the cannon bone with the weight and indicates the capacity of the limbs to displace the body mass (Cabral et al. 2004). The value reported for the Graciosa donkey is higher than the 8.08 reported for the Nordestina donkey breed (Almeida 2009) and much higher than the ones reported for the Mangalarga Marchador horse (4.16 for males and 4.96 for females) (Cabral et al. 2004) and for the Alter horse (3.89 for males and 3.63 for females) (Oom and da Costa Ferreira 1987).

An average of 143 kg for backload index at walk and 79 kg for backload index at canter has been described for the Nordestina donkey breed, a breed mostly used for racing and traction (Pimentel et al. 2014). For the Graciosa donkey average values of 143.4 and 81.94 kg respectively were detected supporting their function as pack animals.

Sexual dimorphism is observed by the higher values for jennies for CH, TC, W, BL, CL, BIIW and BISW and for SSF and CaC for jacks (Table 1). From these TW, W, BL, CL, BIIW and BISW show statistically significant difference between genders. Consequently, several body indices also showed significant differences between genders: DTI, WCI, BR and CHI towards the males and I1, BLW and BLC in favour of the females (Table 2). This is a common finding, as these body regions are associated with the reproductive function of the females (Di Rosa et al. 2007; Sargentini et al. 2009; Liotta et al. 2014; Pimentel et al. 2014; Fonseca Jiménez et al. 2016). Also, CaC showed statistically significant difference between genders.

When compared with other breeds or populations, the Graciosa donkey showed an intermediate stature between donkeys with African origin (Abyssinian, Afar, Hararghe, Ogaden and Omo - Kefena et al. 2011; Algeria - Ayad et al. 2019; Nigeria - Khaleel et al. 2020; Macedónia - Bunevski et al. 2018) and donkeys from European countries (Amiata - Sargentini et al. 2009; Catalã - Folch and Jordana 1997; Banat - Stanisic et al. 2020; Romagnolo - Beretti et al. 2005; Calabrese - Liotta et al. 2014), as revealed by the values of WH, CH and TC. The Graciosa donkey is also lower at withers than the donkey from Miranda, the other Portuguese breed (Quaresma et al. 2005).

Using the core panel recommended for E. asinus, the total number of alleles (48) obtained for the Graciosa donkey was low. This number increased to 61 when two markers recommended for equine genotyping by ISAG (AHT5 and HTG6) were added to our analysis. Also, the other statistic parameters, namely Ne, I, Ho, He, uHe, PIC and PE showed an increase in values with the two additional markers. Nevertheless, our values were always lower than those reported by other studies (Almeida 2009; Jordana et al. 2011; Bordonaro et al. 2012; Matassino et al. 2014; Zhang et al. 2016; Stanisic et al. 2017a, 2017b; Dang et al. 2019; Zeng et al. 2019; Yatkin et al. 2020). However, as the different authors used different marker sets and a different number of individuals were genotyped, these values can be biased. To the best of our knowledge, apart from this study, only Dang et al. (2019) used the complete core panel recommended by ISAG for genotyping 150 donkeys from a local population in Liaoning Province, China. Their study reported a TNA of 101, a mean number of 7.77 alleles per locus, a mean Ho of 0.646, a mean He of 0.608 and a mean PIC value of 0.608. When compared with the Graciosa donkey these values are more than 2 times higher for TNA and mean number of alleles per locus, more than 1.3 times higher for Ho and He and more than 1.5 times higher for PIC. These differences are not surprising as the Graciosa Donkey suffered a major bottleneck in the recent past. Dang et al. (2019) reported PIC values higher than 0.5 for all but three markers (HMS6, HMS7 and TKY337). On the contrary for the Graciosa donkey only five loci (AHT4, HMS3, HTG7, TKY312 and TKY343) of the donkey core panel plus the AHT5 were highly informative with PIC values higher than 0.5. In general, other studies on donkey microsatellite genotyping also reported HMS6 and HMS7 as the least informative markers (Jordana et al. 2011; Bordonaro et al. 2012; Matassino et al. 2014; Zhang et al. 2016; Dang et al. 2019) and HTG7 and HTG10 as the most informative ones (Almeida 2009; Jordana et al. 2011; Bordonaro et al. 2012; Matassino et al. 2014; Stanisic et al. 2017a, 2017b; Dang et al. 2019; Yatkin et al. 2020).

Our data strengthen the need for more studies to define a set of markers to perform paternity testing in a broader range of breeds or populations, as alleles and their frequencies vary between breeds and between populations of the same breed.

As expected low values of genetic diversity were detected. These data are not surprising as the original population of animals brought to the island is most probably still represented in the actual population and until now, phenotypic selection has been based on only a few isolated jennies and jacks.

Conclusion

The Graciosa donkey has been wisely used over centuries by locals as pack animals, for work and traction, as corroborated by the calculated indexes, and as so the interest in preserving such an important genetic resource.

The core panel recommended by ISAG for genotyping E. asinus did not provide enough confidence for performing parentage testing for the Graciosa donkey, thus two additional markers were analysed. This is not surprising as the Graciosa donkey suffered a major bottleneck in the recent past and has its current foundation on a small number of animals. The low levels of genetic diversity detected are of concern if we want to ensure, in the next generations, its sustainability and availability as cultural and genetic heritage. As so breeding strategies, for conservation purposes and in the management of the studbook should include the follow-up of the breed, both morphologically and genetically, to avoid genetic erosion and the search for donkeys, that correspond to the defined morphological type, in Graciosa and in other islands, as these may become an important source of genetic diversity.

Author contributions

ACM conceived the experiment. FC and CCB contacted the breeders and organized the visits. HR, VC and BM collected measurements and blood. MSL and ARA performed the genetic analysis. MSL, ARA and DM analysed the data. MSL, ARA, DM and ACM wrote the manuscript. All authors read the manuscript and agreed to submit it.

Acknowledgments

The authors are grateful to all donkey breeders and owners for participating in this study.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available from the corresponding author [A da Câmara Machado], upon reasonable request.

Additional information

Funding

CBA-UAc is supported by national funds through FCT, IP (Portuguese Foundation for Science and Technology) in the frame of the projects UIDB/05292/2020 and UIDP/05292/2020.