Genetic Diversity and Variability in two Italian Autochthonous Donkey Genetic Types Assessed by Microsatellite Markers

Abstract

Since 13^rd century, Italian domestic autochthonous donkey population has been characterised by Mediterranean grey mousy cruciate ancestral phenotype, currently typical of Amiata donkey (AD) genetic type. This phenotype persisted up to the 16^th century when a marked introduction of Hispanic and French big sized and dark bay or darkish coloured sires occurred. In the context of a safeguard programme of Latial Equide resources, the aim of this research was to evaluate the genetic diversity and similarity between the AD breed and an autochthonous donkey population native from Lazio, the Viterbese donkey (VD), using molecular markers. A total of 135 animals (50 AD and 85 VD) were genetically characterised by using 16 short tandem repeat markers. A high genetic differentiation between populations (F_ST=0.158; P<.01) and a low between-breeds genetic similarity (0.233±0.085) were observed. Correspondence analysis, the result of STRUCTURE software analysis and analysis of molecular variance would seem to indicate genetically different entities as well. It would be desirable to increase the number of comparison with other breeds to better understand the origin of VD. Moreover, results obtained in this study suggest that the loss of genetic variation observed in VD could mainly derive from unnoticed sub-population structuring (Wahlund effect), rather than to other factors such as inbreeding, null alleles or selection influence.

Key Words:

• Microsatellite
• Donkey
• Autochthonous population
• Genetic differentiation

Introduction

Nowadays in Europe, about 55 breeds or varieties of donkey are still present, 15 of which are present in Italy (FAO, 2013). Originally, from the phylogenetic point of view, Italian donkeys were divided into four ancestral breeds: Pugliese, Siciliana, Pantelleria and Sarda (Mascheroni, 1929). Besides, within these breeds it was possible to recognise some subdivisions or local breeds; for example, Siciliana breed was divided into two populations, one of the West and the other of the East of the island. The Pugliese breed, the most diffused in Italy in the past, included these populations: Calabrese, Basilicata, Leccese, Martina Franca, Marchigiana and Romagnolo. After a period of abandonment of many donkey populations, during the last few years, due to exploitation of donkey’s products (milk and meat) many local populations are growing in census. Moreover, some small populations are undergoing safeguard programmes (Rizzi et al., 2011). The use of donkeys as riding and pack animals has been always important in Central Italy, because of the landscape characterised by rugged hills. In particular, in this area there are two donkey populations, the Amiata donkey (AD) and the Viterbese donkey (VD). Amiata donkey breed was present in the area of the Amiata mountain (a large inactive volcano) since the end of the 19^th century and it was characterised by uniform grey mantle with typical striping of the limbs, slim shape, and withers height up to 140 cm. Actually the height of its withers was about 129 cm in females and 132 cm in males (Cecchi et al., 2007). After a bottleneck in the years 1945 to 1970 due to the advent of rural motorisation, the population has started to increase again because of specific actions carried out both by private and public organisations. More than 152 farms are active today in Central Italy and they maintain a population of about 600 registered animals (Cecchi et al., 2006), with a relatively low genetic variability (Ciampolini et al., 2007).

Up to the 18^th century, VD was mainly characterised by grey mousy-crusader or dapple gray or sable or dark sable coat, middle size (withers height of 120 cm), rusticity, robustness, and longevity (Mascheroni, 1929). In the following years, crosses with Marchigiani or Martinesi stud (Mascheroni, 1929; Tortorelli, 1973) characterised by dark sable or blackish coat and great size (withers height more than 130 cm) and following selective interventions on the obtained crosses gave rise to some ecotypes denominated Laziale and Ciociaro, which, although characterised by a wide somatic variability, demonstrated to be good sires giving robust pack-mules, particularly suitable for wood hauling (Tortorelli, 1973).

Until 2011, VD was considered extinct because it was not registered in Herd Book. From 2001 a safeguard programme of Latial Livestock Resources has been carried out involving VD too. In the context of this programme, thanks to the molecular characterisation, this population was registered in the Zootechnical Book of Equine and Ass populations at limited diffusion with 0009742 (7.05.2012) decree of the Italian Ministry of Agriculture, Food and Forestry Policies with the denomination of Asino Viterbese or Asino di Allumiere. In 2012, 146 heads were estimated distributed throughout about 15 farms.

Genetic typifying is an important preliminary step in any safeguard biodiversity programme (Matassino et al., 2010). Literature reports studies on the genetic variability within and between breeds using microsatellite markers in some Spanish breeds (Andaluza, Catalana, Mallorquina, Encartaciones, and Zamorano-Leonesa; Aranguren-Mendez et al., 2001, 2002), one French donkey breed (Baudet du Poitou; Bellone et al., 1998), the Catalonian donkey breed (Jordana et al., 1999, 2001) and three Croatian donkey populations (Ivankovich et al., 2002). In Italy, studies were carried on some donkey Sardinian populations (Cosseddu et al., 2001; Colli et al., 2013); Martina Franca and Romagnolo breeds (Blasi et al., 2005; Colli et al., 2013); Ragusano, Pantesco and Grigio Siciliano genetic types (Blasi et al., 2005; Guastella et al., 2007; Bordonaro et al., 2012; Colli et al., 2013); Amiata breed (Ciampolini et al., 2007; Colli et al., 2013) and Asinara (Colli et al., 2013).

Given the breeding area of AD and VD and the similarity in the coat colour, it could be argued that AD can be an ancestor of VD. Hence, the aim of this study, inserted in a wide safeguard programme of Equide Latial biodiversity (with particular reference to the Regional Law n. 15 of 1 March 2000), was to contribute to the evaluation of genetic diversity and population genetic structure of VD compared with AD, using 16 microsatellite markers.

Materials and methods

Animals

A total of 135 animals – 50 ADs (15♂♂, 35♀♀) and 85 (34♂♂, 51♀♀) VDs – were sampled from different flocks trying to avoid closely related individuals. In particular, 15 farms located in the provinces of Viterbo, Roma and Frosinone for VD, and 35 farms located in different provinces of Tuscany and Lazio for AD.

DNA, short tandem repeat and statistical analyses

Genomic DNA was extracted from 5 mL of peripheral blood samples and DNA was isolated using GenElute Blood Mammalian Genomic DNA Minipreps Kit^® (Sigma-Aldrich, St. Louis, MO, USA). Sixteen microsatellites were investigated (ATH04, HTG07, HTG10, HMS02, VHL20, COR71, ASB23, HMS07, VHL209, COT58, HMS03, HTG04, HTG06, HMS06, HMS045, SGC28) and they were amplified in multiplexes in 4 independent PCR reactions. Primer sequences for the microsatellites are available at the Domestic Animal Diversity Information System web site (FAO, 2013). Genotyping was carried out by using an Applied Biosystems 310 DNA sequencer (Applied Biosystems, Foster City, CA, USA). Data analysis was performed using the GENEMAPPER software (ABI, v. 4.0, 2005).

Number of alleles (Na), allelic frequencies, expected and observed heterozigosity were estimated using Genepop v.3.3 (Raymond and Rousset, 1995). Exact tests for deviation from the Hardy-Weinberg (HW) equilibrium were performed using the ARLEQUIN package v. 3.11 (Excoffier et al., 2005).

The following parameters were computed at the population level by using the programme MolKin (Gutièrrez et al., 2005): molecular coancestry coefficients (Caballero and Toro, 2002), kinship distance, and F-statistics. The molecular coancestry between two individuals – i and j (f_ij) – is the probability that two randomly sampled alleles from the same locus in two individuals are identical by state (Caballero and Toro, 2002). The molecular coancestry of an individual i with itself is self-coancestry (s_i), which is related to the coefficient of inbreeding of an individual i (F_i) by the formula F_i=2s_i–1. In turn, the kinship distance (D_k) between two individuals i and j is D_k=[(s_i+s_j)/2]-f_ij (Caballero and Toro, 2002). MolKin computes within- and between-breed molecular coancestry and D_k by simply averaging the corresponding values for all the within or between-population pairs of individuals.

Genetic similarities were investigated by comparing the individual multilocus genotype of each individual with each other (Ciampolini et al., 1995). Genetic similarity is defined as P=A/2L, where P is the proportion of common alleles (A) in relation to the 2L possibilities (L=number of considered loci). The similarities between each pair of individuals were then averaged over the whole population. In addition, the factorial correspondence analysis performed with the option 3D over populations was estimated using Genetix software (Belkhir et al., 2000). Moreover, a likelihood test of breed assignment was performed using the ARLEQUIN package v. 3.11 (Excoffier et al., 2005).

Further assignment test and analysis of molecular variance were carried using GenAlEx software v. 6.5 (Peakall and Smouse, 2012) which offers both frequency-based and distance-based analyses.

Breed differentiation was also investigated using the Bayesian clustering algorithm implemented in the STRUCTURE software v. 2.2 (Falush et al., 2007). The test was run for K from 2 to 6 and for each K value five repetitions were performed. The test was carried: i) using an admixture model; ii) using 100.000 burn in period; iii) using 100.000 Markov Chain Monte Carlo iterations; iv) giving no prior information; v) choosing the best fitting K value by the likelihood score [LnP(D)] and Evanno method (Evanno et al., 2005).

Results and discussion

All 16 microsatellite markers resulted to be polymorphic in the whole sample and in each breed. A total of 100 alleles were found in VD and 95 in AD, with Na ranging from 2 to 10 in VD [mean value 6.25±2.57; coefficient of variation (CV)=41%] and from 3 to 13 in AD (mean value 5.94±2.89; CV=49%). The mean effective number of alleles (Ne) was 2.99 (±1.20; CV=40%) in VD and 3.00 (±1.49; CV=50%) in AD. The most polymorphic loci were: HTG7 in both genetic types (13 alleles in AD and 10 in VD), HTG10 (9 alleles in AD and 10 alleles in VD) and AHT4 (10 alleles in AD and 9 alleles in VD); the less polymorphic loci were: HMS06 in both breeds (2 alleles in VD and 3 in AD); HTG4 (3 alleles in both populations), HMN045 (3 alleles in AD and 4 alleles in VD) and VHL209 in AD (3 alleles) (Table 1).

Only 10 markers are common to those used by other authors (ATH4, HTG07, HTG10, HMS02, VHL20, HMS07, HMS03, HTG04, HTG06, HMS06) (Aranguren-Mendez et al., 2001, 2002; Jordana et al., 2001; Bordonaro et al., 2012; Colli et al., 2013); ASB23 marker has also been used by Bordonaro et al. (2012), while Ivankovic et al. (2002) used a completely different set of markers. Ciampolini et al. (2007) used two markers more than that reported in this research (ASB2 and HMS1). So, although a comparison with other breeds can be biased due to the different marker sets used by different authors, we highlight that in terms of mean Na the genetic variability observed in VD was higher than that reported in three Sicilian donkey breeds (Bordonaro et al., 2012). On the contrary, AD showed the lowest Na. Furthermore, both populations had an Na lower than that reported by Aranguren-Mendez et al. (2001) on five Spanish donkey breeds (7.0-7.5). In both populations AHT4 and HTG7, markers have the highest Na as reported by Aranguren-Mendez et al. (2001) and by Bordonaro et al. (2012), while HTG4 marker has only three alleles as reported by Bordonaro et al. (2012).

The two populations VD and AD were characterised by a different allelic pattern within locus (Figure 1). A total of 52 breed-specific alleles were observed (28 in VD and 24 in AD), most of them at a frequency lower than 10%. On the contrary, 7 breed-specific alleles (5 in AD and 2 in VD) had a frequency higher than 20%: allele 137 of ATH4 (50%), allele 231 of HMS02 (42%), allele 169 of HTG04 (65%), allele 150 and 152 of marker HMS06 (39 and 56%, respectively) in AD, and allele 165 of HMS03 (45%) and allele 159 of ASB23 (23%) in VD (Figure 1).

The significant overall loci F_ST index was 0.158 (P<.01). This value reveals that 15.8% of the genetic variation observed in the sample is explained by population differences, whereas the remaining is due to the differences within population. This value is generally higher than those observed in studies carried out on European donkey breeds (Arangures-Mendez et al., 2002; Ivankovich et al., 2002; Bordonaro et al., 2012; Colli et al., 2013) and in other species such as for example Italian horse breed (Di Stasio et al., 2008), Italian cattle breeds (Cecchi et al., 2012), and Italian sheep breeds (Ciani et al., 2013).

Observed heterozigosity (Table 1) ranges from a minimum value=0.400 for the locus HMS7 to a maximum value=0.785 for locus ASB23, with an average value of 0.578±0.155 in AD and from 0.380 for locus HMS06 to 0.798 for locus ASB23 in VD with an average value of 0.607±0.162. Expected heterozigosity ranges from a minimum value=0.310 for locus HTG4 to a maximum=0.797 for locus ASB23 with an average value of 0.607±0.163 in VD and from 0.335 for locus VHL209 to 0.874 for locus HTG7 in AD with an average value of 0.602±0.165 (not showed data). The results revealed moderate degree of heterozigosity in the two populations. The expected mean heterozigosity was lower in both populations than that reported in a study on Catalonian donkeys (0.712; Jordana et al., 2001), on three donkey populations of the Croatian coastal region (0.66-0.70; Ivankovic et al., 2002), and on the Baudet du Poitou donkey population (0.623) (Bellone et al., 1998), but higher than that reported on three Sicilian donkey breeds (0.581) (Bordonaro et al., 2012).

Of the total 16 analysed microsatellite, 4 markers were in significant HW disequilibrium in VD (AHT4, VHL209, HMS7 and ASB23) and 5 markers in AD (AHT4, HMS2, HTG10, HTG07 and ASB23) (data not shown); the disequilibrium is given by an excess of homozygotes. The polymorphism information content (PIC) per locus showed nine markers with values over 50% in both populations and an average value of 55.52% in AD and 56.81% in VD, respectively (Table 1). The PIC was originally introduced by Botstein et al. (1980). It refers to the value of a marker for detecting polymorphism within a population, depending on the number of detectable alleles and the distribution of their frequency and it has been proved to be a general measure of how informative a marker is; the higher is the PIC value, more informative a marker is. In the present study, ATH4 and HTG07 and HTG10 microsatellites appeared informative, as reported also by Bordonaro et al. (2012) on three Sicilian donkey breeds and by Aranguren-Mendez et al. (2001) on five Spanish donkey breeds.

The highest values of F_IS was observed in VD (0.070) and the lowest was observed in AD (0.056). Molecular coancestry (f_ij), D_k and F_i are similar in the two breeds ranging from 0.392 (VD) to 0.397 (AD) for f_ij, from 0.317 in AD to 0.324 in VD for D_k, and from 0.428 in AD to 0.434 in VD for F_i. These values are so close to each other that the difference could be due to stochastic effects rather than to an actual difference between the two populations.

The genetic similarity within the populations were 0.467 (±0.086) and 0.473 (±0.107) respectively in VD and AD. These values were higher than those reported in other species such as cattle [0.281; D’Angelo et al. (2006); 0.374-0.420; Cecchi et al. (2012)], sheep [0.318-0.370; D’Angelo et al. (2009)], and dog [0.455; Ciampolini et al. (2011); 0.412; Ciampolini et al. (2013)], showing a low genetic variability in the two populations.

Also within-breed mean coancestry and Kinship distances show that both populations are homogeneous; in fact, the values observed in our study for f_ij and for F_i were clearly greater than that reported in literature on other species such as horse (Marletta et al., 2006), sheep (Dalvit et al., 2008, 2009; D’Angelo et al., 2009; Ciani et al., 2013) and cattle (Cecchi et al., 2012), while D_k was smaller than reported in literature. The observed values might suggest that the low level of genetic variability have arisen as a possible consequence of mating among relatives.

The between-breed D_k and the mean f_ij values were 0.323 and 0.336, respectively, while the between-breeds genetic similarity was 0.233 (±0.085), showing a low genetic similarity between these two populations.

By adopting a likelihood test of population assignment, it emerged that all the animals were correctly allocated to the true breed of origin as showed in Figure 2, thus highlighting a clear differentiation between the two donkey populations. This result was supported by FCA (Figure 3) and by the bayesian clustering analysis (Figure 4).

Correspondence analysis highlighted a clear absence of transvariation area between the two genetic types and an apparent greater dispersion of VD population against AD; indeed, the surface leakage concerning the two populations was higher in the VD with a value equal to 1.41 (41%). Analysing the results derived from STRUCTURE (Figure 4) using K=2, it would seem that the two populations are separated in different groups as observed from plot 1 of Figure 4; for K=3 and K=4 it would seem that VD population was constituted by two sub-populations. Moreover, for K=5 and K=6 both AD and VD would seem divided into different sub-populations. From LnP(D) showed in Table 2 and Evanno method reported in the Appendix the best K value is found to be 6; in fact, the K6 plot shows VD population seemingly characterised by more genetic heterogeneity as it was marked by four clusters, while AD population was marked by two clusters.

The results of STRUCTURE software analyses, together with F_IS value make us hypothesise the presence of a Wahlund effect that affects the calculation of the inbreeding coefficient.

Even the analysis of molecular variance allows to define those two populations as different; in fact, the 26% (P= 0.001) of molecular variance is explained by the differences between them (Figure 5).

Conclusions

In the limits of the observation field, the results show that actually a high genotypic diversity between the two populations exists. Indeed, the two populations, although having the same level of heterozygosity, are characterised by quanti-qualitative allelic diversity, consisting into a different allelic pattern within locus. All distance measures (F_ST and genetic similarities) would seem to indicate genetically different entities. It would be desirable to increase the number of comparison with other breeds to better understand the origin of VD. In particular, research will continue comparing VD with Ragusano, Martina Franca and Asino Romagnolo donkey breeds. Concerning the low genetic variability observed in VD, inbreeding can be one cause but we must remember that there are other factors that can also cause a lack of heterozygotes in the population. Firstly, null alleles (i.e. non-amplifying alleles) may be present and lead to a false observation of excess homozygotes. Despite the fact that pedigrees were not available for analysis, it was not possible to demonstrate the presence of null alleles (usually caused by a mutation in the primerbinding site leading to an allele that will not amplify). Secondly, the presence of population substructure within the breed may lead to a Wahlund effect. The selection influence could not be proved because production data were not available. Results obtained in this study suggest that the loss of genetic variation in VD could mainly derive from unnoticed sub-population structuring (Wahlund effect). It would be very interesting to increase the number of VD samples genotypically typed using a larger number of microsatellite markers. In fact, the analysis of a large number of loci increases the power of detecting population substructure because each locus will contain an independent history of the population depending on the amount of random drift, mutation, and migration that has occurred. An important contribution to better understand its division into sub-populations could be provided using mitochondrial DNA analysis.

Figure 1. Alleles and their frequencies in Viterbese and Amiata donkey populations for each marker.

Figure 2. Results of population assignment.

Figure 3. Factorial correspondence analysis between Amiata and Viterbese donkeys.

Figure 4. Plot of genomic components estimated for the total sample using the software STRUCTURE from K=2 to K=6. Different colours indicate different clusters.

Figure 5. Results of the analysis of molecular variance.

Table 1. Number of alleles, observed heterozigosity, polymorphism information content, and effective allele number within each microsatellite locus in AD, VD and the total.

Marker	AD				VD				Total
	Na	Hobs	PIC (%)	Ne	Na	Hobs	PIC (%)	Ne	Na	Hobs	PIC (%)	Ne
HMS06	3	0.540	46.64	2.14	2	0.380	30.65	1.61	4	0.676	62.87	3.09
AHT4	10	0.620	67.21	3.32	9	0.797	77.13	4.92	15	0.854	83.79	6.84
HTG7	13	0.840	85.22	7.45	10	0.734	71.77	3.76	13	0.803	78.86	5.08
VHL20	6	0.420	38.08	1.68	5	0.535	49.60	2.15	9	0.504	48.44	2.02
HMS7	6	0.400	38.99	1.69	6	0.432	41.17	1.76	7	0.452	40.82	1.74
HTG10	9	0.580	71.84	3.94	10	0.784	79.04	4.62	10	0.852	83.60	6.75
HMS3	5	0.700	62.99	3.17	6	0.705	66.27	3.39	6	0.786	75.45	4.68
HTG4	3	0.490	43.62	2.00	3	0.310	28.08	1.45	4	0.604	54.44	2.52
HMS2	8	0.640	58.15	2.83	7	0.640	57.62	2.80	10	0.732	69.26	3.73
HTG6	4	0,680	64.80	3.30	4	0.732	68.76	3.74	7	0.737	68.98	3.8
VHL209	3	0.235	19.18	1.27	7	0.391	36.86	1.64	7	0.329	30.72	1.49
HMS45	3	0.580	47.68	2.25	4	0.535	47.76	2.15	4	0.547	48.13	2.21
COR71	7	0.700	67.35	3.45	8	0.674	62.23	3.06	9	0.609	64.62	3.23
COR58	5	0.620	58.67	2.84	7	0.722	67.18	3.60	7	0.708	65.68	3.43
SGV28	4	0.440	43.56	2.14	3	0.541	48.04	2.18	4	0.547	47.62	2.21
ASB23	6	0.785	74.32	4.51	9	0.798	76.77	4.99	9	0.808	78.13	5.21

Na, number of alleles; Hobs, observed heterozigosity; PIC, polymorphism information content; Ne, effective allele number.

Table 2. STRUCTURE: likelihood scores.

K	LnP(D)
2	-5413.04
3	-5351.32
4	-5289.88
5	-5250.76
6	-5153.92

LnP(D), likelihood score.

APPENDIX

Appendix. Evanno method.

Table

Run	K	LnP(D)	Var[LnP(D)]	SD	L’(K)=L(K)–L(K-1)	IL”(K)=L’(K+1)-L’(K)	AK=m(IL”(K)I)/s[L(K)]
1	2	-5412.7	113.7	10.663
2	2	-5413.3	114.7	10.710
3	2	-5412.8	113.7	10.663
4	2	-5413.1	114.5	10.700
5	2	-5413.3	114.7	10.710
Mean value		-5413.04	114.26	10.689
1	3	-5325.5	190.2	13.791
2	3	-5453.4	440.9	20.998
3	3	-5326.5	192.6	13.878	61.72	-0.28	-0.018351285
4	3	-5326.8	193	13.892
5	3	-5324.4	188.5	13.729
Mean value		-5351.32	241.04	15.258
1	4	-5288.9	284.2	16.858
2	4	-5292	290.5	17.044
3	4	-5294	291.7	17.079	61.44	-22.32	-1.31998642
4	4	-5286.7	280.6	16.751
5	4	-5287.8	282.7	16.814
Mean value		-5289.88	285.94	16.909
1	5	-5267.8	368	19.183
2	5	-5274.2	378.7	19.460
3	5	-5289.3	378.9	19.465	39.12	57.72	3.051219357
4	5	-5266.2	366.5	19.144
5	5	-5156.3	300.4	17.332
Mean value		-5250.76	358.5	18.917
1	6	-5152.8	414.2	20.352
2	6	-5153.5	415.2	20.377
3	6	-5145.6	400.2	20.005	96.84	-96.84	-4.750260842
4	6	-5154.4	416.8	20.416
5	6	-5163.3	431.9	20.782
Mean value		-5153.92	415.66	20.386

LnP(D) likelihood score; SD, standard deviation.