Characterization is the first approach to sustain use of a genetic resource (Yakubu et al. 2010). Ideally, it should encompass both the phenotypic descriptions based on a number of indicators and morphological measurements (FAO, 2012), and the genetic characterization through molecular markers (FAO, 2011). Microsatellites Markers have been increasingly used in genetic characterization studies of farmed animals (Sunnucks, 2001) and particularly in the assessment of genetic relationships (FAO, 2008). These DNA markers are still used in such studies (Aljumaah et al. 2015; Seilsuth et al. 2016; Lenstra et al. 2017). The Moroccan goat population is over 6.2 million heads and is ranked in the 13th position worldwide in terms of goat’s numbers (http://faostat3.fao.org/browse/Q/QA/E). Until the beginning of the 21st century, only the Draa goat has been considered as a breed, while other goat populations have traditionally been named as the names of their geographical location, e.g. “goat of the North” or “black goat of the mountains”. In the last decade, some subpopulations of the black goat were described as separate on the basis of slightly different morphologies (Atlas, Barcha and Ghazzalia) (Ibnelbachyr et al. 2015). However, these phenotypic different should be verified with genetic analysis evidences. The highly distributed genetic variability in Moroccan goats may reflect the levels of the ancestral genetic diversity existing in the first individuals that colonized Morocco and/or the existence of an extensive and recurrent inter-population gene flow (Benjelloun et al. 2011). In fact, livestock populations which are geographically, ecologically or culturally isolated may have been affected by different natural and artificial selection (FAO, 2003). Therefore, it can be expected that there are a degree of likeness between Draa breed and other goat populations of Morocco. The aim of this study was to use microsatellite markers to characterize Draa breed and four other Moroccan goat populations.
MATERIALS AND METHODS
The screening was done in the cradle of the Draa goat (South-Eastern Morocco) and its neighbouring areas in south-eastern and southern Morocco (Figure 1). The sampling location was established based on the geo-referenced cells of 2500 km². Within each cell, three herds were randomly selected and 3 to 4 unrelated adult goats were selected per herd, until 30 to 50 samples were collected per population, as recommended by FAO for genetic diversity studies (FAO, 2011). The investigated goat populations were Draa, Atlas, Barcha, Ghazzalia and other populations that were grouped together and named as “undefined goats”. Phenotypic traits were recorded, and then blood samples were taken from the jugular vein using vacutainer tubes. At the end, 150 genotypes were analyzed: 35 from Atlas, 31 from Barcha, 36 from Draa, 32 from Ghazzalia and 16 from ‘’Undefined goats’’. These animals, 55 males and 95 females, were aged approximately from 2 to 6 years.
Microsatellite amplification and analysis
Genomic DNA was extracted using the salting out procedure (Miller et al. 1988). Sixteen microsatellite markers, which are recommended by FAO-ISAG panel were amplified and analyzed (Table 1). The PCR amplifications were performed in a final volume of 28 µL containing 2 μL of genomic DNA (on average 50 ng/μL), 1 μL of each primer (80 nmol/mL) and 22 μL of the following PCR Mastermix solution (SuperMix from Invitrogen, Life Technologies): 22 mM Tris-HCl, 55 mM KCl, 1.65 mM MgCl2, 22 U/mL Taq polymerase and 220 μM of each dNTPs. The following PCR protocol was used: an initial denaturation step (95 ˚C for 5 min) followed by 30 cycles of a denaturation step (94 ˚C for 15 s), an annealing step (annealing temperature of the primers, varied from 55 to 65 ˚C, for 45 s) and an extension step (72˚ for 90 s). The 30 cycles were followed by a final extension (72 ˚C/30 min). The PCR products were resolved on a 3500 DNA fragment analyzer (Applied Biosystem) and the obtained data were analyzed with GeneMapper 5.0 software. Standard DNA samples provided by the ECONOGENE Project (www.econogene.eu) were also genotyped and used to adjust allele length.
The measure of genetic polymorphism such as number of alleles, the observed (Ho) and expected (He) heterozygosity, were calculated using GENETIX software version 4.03 (Belkhir et al. 2004). Allelic richness (Na), which is an unbiased estimator of the number of alleles, was calculated using FSTAT software version 2.9.3 (Goudet, 2001). The GENETIX software was used to estimate the power of these markers to discriminate between populations via the coefficients of genetic differentiation GST (Nei, 1978) and to carry out a factorial correspondence analysis at the population level. Population structure was analyzed through a Bayesian clustering procedure implemented in the version 2.3.4 of STRUCTURE package (Pritchard et al. 2000). For each value of K (K=2 to K=10), five independent runs were carried out under an admixture model with correlated allele frequencies and with no prior information on the population of origin. A burn-in-period of 50000 was followed by 100000 Monte Carlo Markov Chains iterations. Two parameters were calculated to estimate the best fitting number of clusters, i.e. K value: i) the mean values of log likelihood, lnPr(X|K), calculated for each K over 10 runs, and ii) the posterior probabilities of K, Pr(X|K), from the mean lnPr(X|K) values according to Pritchard et al. (2000). Graphical representations of structure analyses were obtained using DISTRUCT software (Rosenberg, 2004). Version 3.5 of ARLEQUIN package (Excoffier et al. 2005) was used to perform an Analysis of Molecular Variance (AMOVA; Excoffier et al. 1992) and to calculate the value of the coefficient of inbreeding (FIS) and FST (Wright, 1968) distances between populations. This FST distance matrix was subsequently used to build a neighbour-network with version 4.13.1 of SPLITSTREE software (Huson and Bryant, 2006).
RESULTS AND DISCUSSION
Polymorphism and heterozygosity of the studied loci
In total, 160 alleles were found and allelic richness values varied between 2.770 and 9.669. Twelve loci were polymorphic, with mean effective number of alleles ranged from 5.853 to 9.669 and heterozygosity from 0.546 to 0.885 (Table 1).
Figure 1 Goat sampling locations are shown on the map of Morocco
Different colors correspond to different breeds / populations
Some sampling points are overlapping since they correspond to individuals from the same or different populations in the same location
Table 1 Allelic richness (Ne), expected heterozygosity (He), observed heterozygosity (Ho), coefficients of differentiation (GST and FST) of microsatellite DNA markers used
NED: fluorescent dye of ABI yellow color; VIC: fluorescent dye of ABI green color; PET: fluorescent dye of ABI red color and FAM: fluorescent dye of ABI blue color.
These loci had GST and FST values ranged from 0.001 to 0.126. The least polymorphic loci were MAF209 (3.185 alleles on average), MCM527 (3.524 alleles), ETH10 (3.488 alleles) and ILSTS05 (2.770 alleles), and were excluded from the subsequent analyses, as recommended by Barker (1994).
Genetic diversity between and within the populations
The average number of alleles was similar in the five populations (Table 2), while the “undefined” goats showed the lowest average followed by the Draa breed. The expected heterozygosity (He) was ranged from 0.635 to 0.705, while the observed heterozygosity (Ho) was 0.579 in Draa, 0.613 in Barcha, 0.614 in Ghazzalia, 0.635 in undefined goats and 0.644 in Atlas. Most of the total variance (88.4%) was distributed within individuals and only 1.85% of variance was due to differences between populations (Table 3). The inbreeding coefficient (FIS) was ranged from 0.161 for Draa breed to 0.037 for the “undefined” group. All values were also statistically significant (P<0.001) with the sole exception of the latter population. For FST coefficient, it varied from 0.005 between Barcha and Ghazzalia to 0.040 between Ghazzalia and Draa. The graphical representation of the factorial correspondence analysis results (Figure 2) showed the distribution of the individual genotypes along the axes corresponding to the three first major components of inertia (total inertia explained 85.3%). The scatter of the points corresponding to the individuals from Draa breed occupied a very large area of the multivariate space of the graph, stretching far from the remaining populations, which, in turn, largely overlapped on the central part of the plot surrounding the origin of the axes.
The distinctiveness of Draa breed and the weak structure of the other populations were also highlighted by structure analyses (Figure 3). At K= 4 – identified as the best fitting resolution according to Pritchard et al. (2000); (Table 4) – a clear split appeared within the group of individuals assigned to the blue cluster at K= 2 and separating half of Draa population from the rest. Also, the appearance of the purple cluster did not create any clear split, but rather highlighted some distinctiveness of Ghazzalia breed in which the red component remained predominant (Figure 3). Conversely, both Barcha breed and the undefined group showed the lowest level of population structure what was confirmed by the average proportions of membership. Within these two populations, the genomic clusters shared by all five populations (red, purple and green) occurred almost at the same percentage, while Draa, Atlas and Ghazzalia were characterized by the prevalence of the blue, green and red components, respectively. The average proportions of membership at each sampling location showed a geographical pattern underlying the distribution of the clusters. In particular, the blue cluster is localized mainly in a specific area in the central part of the country that corresponds to almost the whole distribution area of Draa breed.
Genetic distances and phylogenetic relationships
Nei’s genetic distances between the Draa goat and the other local goats (Table 5) showed that the Draa breed had the highest distance from the other populations. Even if all these genetic distance values of Draa are almost the same, the Neighbour-network graph (Figure 4) confirmed the evidence already revealed by Structure software results: Atlas and Draa branches stem from the same basal edge, thus suggesting an ancestral genetic relationship between these breeds. In addition, the genetic distinctiveness and relative isolation of Draa breed were highlighted by the major length of its branch. Among the studied populations, Barcha breed seems to be the less differentiated, being positioned in the middle of the graph with a branch length close to zero value. The distribution of the major variance component (about 98.2%) among and within individuals is in agreement with the findings of Benjelloun et al. (2011). Low levels of differentiation within sheep and goat populations have already been highlighted by several studies (Cañόn et al. 2006; Peter et al. 2007). The highly distributed genetic variability in Moroccan goats may reflect the levels of the ancestral genetic diversity existing in the first individuals that colonized Morocco and/or the existence of an extensive and recurrent inter-population gene flow (Benjelloun et al. 2011). The observed heterozygosity, which varied from 0.579 in Draa to 0.644 in Atlas goats, are close to those found in the Canary breeds (Martínez et al. 2006), in some local Egyptian goats (Agha et al. 2008) and others African ones (Missohou et al. 2011), but less than Gaddi goat (Singh et al. 2015). In all studied breeds, the observed heterozygosity (He) was lower than the expected heterozygosity (Ho) and the Draa breed showed the most severe departure from Hardy-Weinberg equilibrium due to a deficit of heterozygoty as testified by the FIS value of 0.161. This value is close to that reported in Baladi breed (Agha et al. 2008), in some Chinese populations (Li and Valentini, 2004) and in some breeds of India (Dixit et al. 2012), but less than Berari breed (Kharkar et al. 2015). Likewise, Bruno-de-Sousa et al. (2011) found a high value of inbreeding coefficient (0.124) in Preta breed in Portugal associated with reduced population size.
Table 2 Number of effective alleles (Na), expected heterozygosity (He), observed heterozygosity (Ho) and inbreeding coefficient (FIS) with significance level
NS: non significant.
Table 3 Analysis of molecular variance (AMOVA)
Figure 2 Factorial correspondence analysis showing multivariate relationships among individual genotypes of Moroccan goats (Atlas in yellow, Barcha in red, Draa in white, Ghazzalia in bleu and undefined goats in rose color)
A deviation from equilibrium indicates a departure from random mating and suggests that the higher than expected number of homozygotes probably due to inbreeding. The reduced numbers of Draa goats and the selection for prolificacy and milk under which the breed has been submitted for a longtime might be the cause of its deviation from equilibrium. Thus, the possibility that Draa subpopulations may have suffered from genetic bottlenecks, at least in part of the breed’s distribution area can’t be ignored. Indeed, this event has been reported in several studies in sheep (Fatima et al. 2008) and goats (Korkmaz Ağaoğlu and Ertuğrul, 2012). In the case of Draa breed, its closed geographical area, its reduced population size, the small size of herds and the selection performed by breeders are factors that may have caused its high level of inbreeding. Moreover, as revealed by Structure results plot, a strong population substructuring is present within this breed, with almost half of its individuals forming a clearly defined cluster, which may represent a separate gene pool existing within this breed and may have caused a Wahlund effect that means a reduction of heterozygosity in our population due to subpopulation structure.
Table 4 Mean values of ln Pr(X|K) and posterior probabilities of K (Pr(X|K)) calculated according to Pritchard et al. (2000)
Figure 3 Bayesian clustering performed with STRUCTURE software on Moroccan goat microsatellite data
Table 5 Genetic distances matrix between identified populations (Nei, 1978)
Figure 4 Neighbour-network based on FST distance matrix between the studied populations (Atlas, Barcha, Draa, Ghazzalia and Undefined goats)
Figure 5 Proportion of membership of each sampling location on the results of Structure software run at K= 4
Each individual is represented by a pie chart that represents their genetic clusters
Pie charts with black / white outline correspond to locations where single / multiple individuals have been sampled, respectively
The genetic distances confirm the differentiation level between Draa breed and the other populations. The FST values obtained between our breeds are closed to those obtained by Al-Atiyat et al. (2015) between Jordan goats. Our results indicate that similarities between the Draa and other populations are low even if Draa and Atlas branches stem from the same basal edge. An evidence which was confirmed by the findings of Tadlaoui Ouafi et al. (2002) who reported low values for the distance between the Draa breed and the black goat called ‘’R'halia’’, a term identifying the black goat of the Middle Atlas in Morocco and which corresponds in our case to Atlas and Barcha populations. Even though the genetic makeup of the Moroccan goat populations may have been influenced by gene flow and introgression coming from different areas of the Mediterranean basin during the past (Benjelloun et al. 2011; Pereira et al. 2009), nevertheless the vast majority of Moroccan goat populations likely derive from the first founders that reached Morocco in the late Neolithic. In addition, the results showed a lack of population sub-structuring for most breeds and the presence of two genomic components (i.e. red and purple in Figures 3 and 5) spread over most of the sampling area, fit well with the hypothesis of an ancestral wide genetic basis still shared today due to extensive gene flow, with the sole exception of a sub-population of Draa breed. On the other hand, the presence of a geographical pattern in the distribution of the major genetic components could be explained either by an uneven and gradual introgression of genomic contributions from gene pools of non-native origin, or by the effect of increasing geographical distance which progressively prevents panmixia. Indeed, a geographical gradient in genetic diversity might also reflect the adaptation of animals to the environment and their response to natural selection (Ciani et al. 2013). But to investigate both such effects and the real extent of gene flow, isolation by distance and introgression, a more dense set of molecular markers (e.g. medium or high SNP panels or whole genome sequences) is needed.
The studied Moroccan goat populations were characterized by a substantial genetic diversity weakly structured between populations, with the exception of Draa breed which showed a certain degree of differentiation due either to reproductive isolation or inbreeding, possibly caused by geographical isolation, reduced population size and / or man mediated selection processes.
The authors acknowledge the technical contribution of Mrs. Sekkour M., Lberji A., Dadouch A., Errouidi C., Bouaziz B. and Elmaata M. in carrying out this study. They acknowledge also the technical support of Mr. Bechchari A. in elaborating sampling maps.