Assessing patterns of genetic admixture between sheep breeds: Case study in Algeria

Abstract In developing countries, cross‐breeding between local breeds and indigene or exotic breeds represents one of the main threats to the livestock diversity, leading to genetic dilution and loss of unique allelic combination underlying essential local adaptive traits. In this study, two Algerian sheep breeds, known to be highly admixed, were considered as a case study, to demonstrate how combination of different methodologies coupled with the use of specific softwares can be efficient to assess the spatial structuration of a hybrid zone, even in a case of extreme admixture. A fine sampling covering distribution areas of both breeds was implemented in order to study the admixture area and adjacent zones from a phenotypic (i.e., 19 quantitative traits were considered) and a genetic point of view (i.e., 21 microsatellites markers were used). Both approaches gave concordant patterns, highlighting areas with sheep most differentiated (or less admixed) for each breed. In detail, the region of Biskra appeared as the most preserved for the Ouled‐Djellal breed and the northwest of Laghouat was identified as the most preserved area for the Rembi breed. The approach proposed in the study offers a low‐cost solution to identify the most representative flocks of a breed, allowing the implementation of efficient conservation plans.


| 6405
HARHAT AT Hal. considering only domesticated populations in spite of the importance of this issue in the livestock conservation field.
According to FAO, across the world, one-third of farm animal breeds face extinction, and cross-breeding with exotic breeds usually exported from developed to developing countries or with indigenous breeds, represents one of the greatest threat for livestock diversity, particularly in developing countries (FAO 2014(FAO , 2015. It is precisely in developing countries that the greatest genetic diversity can be recorded (Da Silva & Benjelloun, 2016). Breeds of these areas are still largely managed in a traditional way (i.e., with management practices characterized by "soft" artificial selection and important place for natural selection (Taberlet et al., 2008)) and then have remained closely connected to their original habitat (Hoffmann, 2013) unlike "industrial" breeds that are submitted to harsh artificial selection. Cross-breeding leads to the homogenization of the livestock genetic pool and induces a genetic dilution of the specific allelic combinations, selected over millennia by natural and artificial selection that confer to the valuable local breeds strong adaptation to their environments.
In Algeria, Ouled-Djellal, one native sheep breed considered more productive, is favored by most of the breeders. Algerian farmers under increasing economic pressure realize unsupervised cross-breeding (i.e., not carried out in the framework of selection plans) between Ouled-Djellal and local Algerian breeds, with the hope to improve their productivity (Madani, Yakhlef, & Abbache, 2003). The practice is so extensive that three breeds have been found highly genetically admixed with the Ouled-Djellal: Rembi, Taâdmit, Berber and to a lesser extent Barbarine (Gaouar, Da Silva, et al., 2015;Gaouar, Kdidi, & Ouragh, 2016). The situation between Ouled-Djellal and Rembi appeared as the most critical: Gaouar, Da Silva et al. (2015), by use of 30 microsatellite markers, recorded a pairwise F ST of 0.009 [0.002-0.017] 95% and by genome-wide single-nucleotide polymorphism genotyping, obtained a value for the pairwise F ST not significantly different from zero (p-value = .001) between the breeds.
In this study, we considered the extreme case of admixture between Ouled-Djellal and Rembi as a case study, in order to analyze at a fine geographic scale, the hybrid zone, from a phenotypic and genetic point of view. The study was conducted in the Algerian region where the most part of Ouled-Djellal and Rembi are encountered, that is, including three zones, the cradles of each breed and the contact zone intermediate between them. Microsatellites markers were used to infer genetic diversity. Our objectives were to assess spatial distribution of genetic admixture in the different zones and to determine in which genetic and phenotypic patterns overlapped. This study demonstrated how combination of different methodologies (Wright's F statistics (1951), discriminant analyses of principal components, spatial analyses of principal components) coupled with the use of softwares specially designed to identify and estimate individual and/or population admixture (FLOCK (Duchesne & Turgeon, 2012), NEWHYBRIDS (Anderson & Thompson, 2002)) can be efficient to dissect the spatial structuration of the hybrid zone, whereas in case of extreme admixture classical approaches are ineffective; giving consequently strong foundation to the development of conservative plans.

| Sampling design
Sheep were sampled in their respective distribution areas (i.e., areas where they are mostly encountered): Ouled-Djellal was sampled in Biskra, M'Sila, Djelfa and in the south of Laghouat; Rembi was sampled in Tiaret, Djelfa and in the north of Laghouat, covering a surface area of around 105,854 km 2 . Tiaret region (considered as the cradle of the Rembi breed) is dominated by Rembi, whereas Biskra (considered as the cradle of the Ouled-Djellal breed) and M'Sila areas are dominated by Ouled-Djellal; the middle area, including Djelfa and Laghouat, represents the overlapping zone supposed to favor strongest admixture (see Figures 2 and 3). The sampling area was largely located between the Saharan Atlas and the high Tellian plateaus. Overall, 54 farms were visited during spring 2013 in order to phenotype sheep of both breeds, and 78 farms were visited during spring 2014 in order to collect blood samples.  Table S1). For each individual, 19 quantitative traits were assessed: body weight, head length, head width, ear length, ear width, neck length, body length, chest girth, cannon perimeter, tail length, trunk length, chest width, internal width of the chest, external width of the chest, chest depth, rump length, hip width, ischium width, withers height (for details considering the process of phenotyping see Harkat et al., 2015 andLaoun et al., 2015).  Table S1).

| DNA extraction, polymerase chain reaction (PCR), and fragment analysis
DNA was extracted from the blood samples according to the salting out procedure (Miller, Dykes, & Polesky, 1988). Twenty-three microsatellites were amplified; all the microsatellites except CSSM66 and

| Data analysis
The mean number of alleles per breed, the average observed (H o ) and expected (H e ) heterozygosity over loci per breed were estimated using ARLEQUIN 3.5 (Excoffier & Lischer, 2010). To calculate allelic richness and the richness of private alleles, we used the rarefaction method (Kalinowski, 2004) implemented in HP-RARE (Kalinowski, 2005) adopting a sample of 44 individuals. Polymorphic information content (PIC) and effective number of alleles (Na e ) were estimated for all markers using the Molkin software (version 2.0) (Gutiérrez, Royo, Alvarez, & Goyache, 2005).
Frequencies of null alleles at each locus and for each breed were estimated with INEST (http://genetyka.ukw.edu.pl/INEst10_setup. exe), in order to take into account simultaneously null allele frequencies at each locus and the average level of the intrapopulation inbreeding as a multilocus parameter (Chybicki & Burczyk, 2009).
The unbiased estimator of Wright inbreeding coefficient, F IS , was calculated following Weir and Cockerham (1984) (f estimator). Its significance was assessed using a permutation method (10,000 permutations) implemented in the GDA (Lewis & Zaykin, 1999).
The extent of population subdivision was examined by calculating the global multilocus F ST value. The index of pairwise F ST of Weir and Cockerham (1984) and their associated 95% confidence intervals were determined using GDA (Lewis & Zaykin, 1999).
To assess the degree to which groups of individuals differ from each other, from a genetic and a phenotypic point of view, we performed discriminant analyses of principal components (DAPC), using the approach implemented in the ADEGENET package (Jombart, 2008) within the statistical package R version 3.2.2 (R Core Team 2015). For the "phenotypic" DAPC analysis, quantitative variables were sizecorrected using the residuals after regression on overall body size.
We conducted spatial analyses: (1) spatial analysis of principal components (sPCA), a multivariate method optimizing the identification of spatial genetic patterns, was performed with the package ADEGENET (Jombart, 2008) using as connection network, the Delaunay triangulation (Upton & Fingleton, 1985); (2) a similar approach was conducted considering phenotypic traits by the use of MULTISPATI-PCA (Dray, Said, & Debias, 2008), that takes into account the spatial location of samples. Computations were conducted with the "ade4" (Chessel, Dufour, & Thioulouse, 2004) and "spdep" packages (Bivand, Pebesma, & Gomez-Rubio, 2008) for the R statistical software package. For both analyses, Monte Carlo tests were used to check the statistical significance of the spatial structures (global and/or local spatial structure) for 10,000 iterations. The sPCA results were visualized by plotting the samples according to their geographic coordinates and coloring them according to their respective scores along the first sPCA components.
The approach of Duchesne and Turgeon (2012) implemented in the software FLOCK was also used to estimate the level of genetic differentiation between groups of sheep. FLOCK is a non-Bayesian method, specially designed to provide spatial and/or temporal admixture maps. This method provides high resolution power even when most genotypes are admixed.

| RESULTS AND DISCUSSION
The current study investigated, for the first time, an area of admixture between two sheep breeds at fine scale, with the principal objective of identifying most preserved flocks of each breed. Ouled-Djellal and Rembi, local sheep breeds of Algeria, were considered as a case study of critical admixture analysis, and the flocks were evaluated from phenotypic and genetic point of view.

| Genetic diversity
The twenty-three microsatellites loci were highly polymorphic (see details in Table S2); a total of 293 different alleles (mean = 12.70 per locus ± 3.51) were found in the two breeds, effective number of alleles ranged from 3.02 for MAF209 to 8.27 for MAF70 (mean = 5.15 per locus ± 1.73,) with a mean PIC of 75.78 ± 8.40. On average, 95.80% of individuals were successfully typed for each microsatellite (the worst score was found for DYMS1 with 86.18% of successful genotypes obtained).
The software INEST suggested possible null alleles for more than 20% of individuals, for two microsatellites (HUJ616 and MAF65); these observations were made in both breeds. Hence, we decided to exclude these microsatellites from the following analyses.
Considering both breeds, the average observed heterozygosity (H o ) was 0.77 ± 0.10 and average expected heterozygosity H e was 0.78 ± 0.08. All the genetic diversity indices considered were found really close between the two breeds ( We failed to detect significant linkage disequilibrium between pairs of loci in each considered breed after False Discovery Rate correction.  Monte Carlo test (10,000 iterations) indicating significant global spatial structure (p-value = .001) (Figure 4).

| Spatial structuration
Given the critical level of admixture recorded in the study, classical Bayesian methods (e.g., implemented in the STRUCTURE software, Pritchard, Stephens, & Donnelly, 2000) used to infer the optimal number of clusters (K) were inappropriate. We then used FLOCK, with the aim to test the significance of the pattern that emerged from the previous analyses; in this perspective, samples were partitioned into three clusters (K = 3) using the "samples mode" to specify that initial partition should include Rembi from Laghouat and Ouled-Djellal from Biskra in single and different clusters. For each sample, the number of allocations to each cluster was computed (Table 3). A majority of Ouled-Djellal from Biskra was assigned to K1 (66.7%); Ouled-Djellal from Biskra were mainly assigned to K2 (69%), with to a lesser extent Ouled-Djellal from M'Sila. Rembi from Laghouat were mainly assigned to K3. The global p-value was highly significant (p-value = .0016) pointing out the existence of the genetic structure tested.
Finally, we used the NEWHYBRIDS program to classify individuals further into genealogical classes. This analysis is a probability-based model, which computes through MCMC, the posterior probability, qi, of individuals belonging to distinct genealogical classes. According to Anderson (2008), the software performs better when "pure" representatives of the parent breeds are specified a priori. Hence, based on previous results, we specified that Ouled-Djellal from Biskra and Rembi from Laghouat had highest likelihood to be composed of "pure" individuals of each breed. The Figure 5 displays posterior probability for each individual to be assigned to one of the classes "pure Ouled-Djellal," "pure Rembi," or "F2 hybrids." Histograms (on the right of the  (3) it was precisely in the overlap zone (i.e., Djelfa and Laghouat) where admixture was supposed to be stronger, that the most preserved area for the Rembi breed was identified (northwest of Laghouat).
Rembi is subjected to intense practices of cross-breeding with Ouled-Djellal (MATET, 2009), considered by the breeders as the most productive breed. The situation reaches such a point that for specialists, initial specimens (characterized by uni-colored bay-fawn fleece, short neck and massive horns, as described by Chellig, 1992)

| CONCLUSION
A fine scale analysis of this type implies considerable sampling effort (in order to cover the distribution ranges of the breeds); F I G U R E 5 Results of NEWHYBRIDS analysis considering 153 sheep genotyped with 21 microsatellites. On the main figure, the posterior probability, qi, of individuals belonging to distinct genealogical classes is represented for each individual; the three classes considered are "pure Ouled-Djellal" in red, "pure Rembi" in green and "F2 hybrids" in blue. On the right, histograms show the proportions for each sample, of individuals assigned to one of the class: in pink proportion of individuals assigned to a class with 0.5 < qi < 0.8, and in red individuals assigned to a class with qi > 0.8 microsatellites, as cost effective-markers, constitute an affordable solution and hence represent valuable option to infer genetic diversity in such case. This study showed that the use of microsatellite markers combined with spatial analyses (using different methodologies and specialized softwares) constitutes relevant solution to assess genetic subtle substructures, even in situation of strong admixture. Cross-breeding is one of the main threats affecting the genetic diversity of livestock (FAO 2014(FAO , 2015. For example, only in North Africa, high levels of admixture were recently reported among Egyptian goats (Elbeltagy et al. 2016), Turkish sheep (Yilmaz, Cemal, & Karaca, 2014), Tunisian sheep , and among Moroccan sheep (Gaouar, Kdidi, et al. 2016). Locally adapted breeds of livestock represent true reservoirs of biodiversity; they are the result of unique evolutionary phenomena, showing strong adaptation to harsh environments. These highly resilient breeds, offering key solution to deal with the global warming issue, are endangered, and in the same way, the worldwide food security is threatened. The ability to analyze hybridization zones and to identify most preserved areas from genetic dilution is essential in such a context, allowing implementing conservation plan and preserving this unique genetic heritage.

ACKNOWLEDGMENTS
We thank Belkacem Mezrou, Hadibi Mohamed, Tarik Tahir, for their assistance in the field. We are also grateful to Algerian breeders for their involvement.

CONFLICTS OF INTEREST
None declared.