Genetic structure and variability within and among populations of the fat-tailed Barbarine sheep breed using microsatellites markers

This study investigates the genetic diversity and the structure of the most dominant native fat-tailed Tunisian sheep breed (Barbarine, BAR) using microsatellite markers. Blood samples from 183 BAR animals, belonging to 4 subpopulations according to phenotypic traits, were collected across all regions in Tunisia. BAR animals and 31 Appenninica Italian sheep breed (APP) used as an out-group were genotyped at 17 microsatellites loci. A total of 270 alleles were identified with average gene diversity equal to 0.812. The mean observed heterozygosity (0.745) and allelic richness (8.09) estimates were high within BAR breed highlighting notable levels of genetic diversity. The low FIS (0.078) and FIT (0.084) values indicate low level of inbreeding within this breed while a low FST estimate (0.007) shows that the subpopulations are not genetically differentiated. The clustering analysis performed with ‘structure’ detected the absence of substructures and the clear uniqueness of the BAR. Tomiuk and Loeschcke’s DTL genetic distance values confirmed the distinction between APP and BAR breeds. Results arising from our microsatellites analysis represent a starting point for the valorization of this indigenous Tunisian sheep breed. A suggestion was made to monitor its genetic variability and for the preservation of this breed for the next generations.


INTRODUCTION
Sheep genetic resources are playing a very important role in developing countries, especially throughout the Near East and North Africa (NENA) regions where most of the local breeds have not been sufficiently characterized while their share expressed in animal units represents 15% of the world small ruminants (Galal, 2010).In Tunisia, sheep farming is the most important domestic livestock activity.According to 2010 Food Agricultural Organization (FAO) data, among the 109.897 million sheep heads reared in North Africa, 7.234 million are reared in Tunisia showing a growth of 4.44% in the last 10 years (FAOSTAT, 2012).Tunisia holds a numerous and diverse set of native and exotic sheep breeds adapted to a range of dry environment with extreme climate fluctuations during the year.There are four native sheep breeds, three meat breeds (Barbarine, reach up to 15% of the total carcass weight in wellshaped adult animals (Khaldi, 1984).This fat is considered an effective means of resistance to adverse climate conditions.It can undergo enormous changes in weight after food deprivation without significant alteration of the physiological state of the animal.However, the very large fat tail began to be an annoying problem due to the necessity of shepherd assistance at mating time and a change in the eating habits of Tunisian consumers who start to prefer less fatty meat.Indeed, an intense competition between BAR and other thin-tail breed has appeared (Bedhiaf-Romdhani et al., 2008b) and the BAR is beginning to be threatened by the loss of its genetic traits through uncontrolled crossbreeding with other native and exotic thin tail breeds.Thus, the proportion of BAR to the national sheep population has decreased from 87% in 1964to 63% in 2006. Bedhiaf-Romdhani et al. (2008a) considered that the breed is composed by ten different populations defined as ecotypes.However, it should be noted that, during our sampling campaigns covering all areas of Tunisia, the majority of herds throughout Tunisia were mixed with different ecotypes.In the northern region, two ecotypes or phenotypes were identified: the black face called "Barbarine à tête noire" (NBTN) and the brown face called "Barbarine à tête rousse" (NBTR).Another black face variety with denser and longer wool nominated "Chalfie" (CH), more adapted to cold climates according to the breeders, was found in the North.In the centre, three ecotypes were recognized: the brown face (CBTR), the brown face with a white line in the middle called "Barbarine à Liste Frontale Blanche" (LFB) and the Barbarine with brown eyes' areas, muzzle and legs named "Sagaa" (SG).
In the southern region, three ecotypes were encountered: the brown face (SBTR), the black face (SBTN) and the Barbarine with black eyes' areas, muzzle and legs named "Sardi" (SR).It should be pointed out that SR, SG and LFB individuals have been found in the majority of private herds throughout Tunisia.Only brown face and black face homogeneous herds were encountered among big breeders.This is mainly due to the selection made in the sixties following the opinion of Palian (1966), who advised to differentiate animals with black face from those with brown face.
In this study, we grouped the brown face ecotypes and black face ecotypes present in northern, central and southern Tunisia in two separated groups called BTR and BTN, respectively.Appenninica (APP) is a native Italian sheep breed reared mainly in central Italy.It is a mediumcoarse wool breed with semi-lopped ears kept primarily for meat production (Associazione Nazionale della Pastorizia, 1997).In the present work, it was used as an out-group for comparison purposes.Despite the importance of molecular information in the establishment of genetic improvement programs, in the management and conservation of breeds, the genetic diversity of Barbarine and other indigenous sheep in Tunisia has not been sufficiently studied using molecular tools; only RAPD-PCR markers were performed to investigate Tunisian sheep breeds (Khaldi et al., 2010;El Hentati et al., 2012).The present study is aimed at investigating the genetic variability and the population structure of the Barbarine breed using microsatellites markers and testing the hypothesis that this breed represents genetically distinct ecotypes.

MATERIALS AND METHODS
Research protocols followed the guidelines stated in the Guide for the Care and Use of Agricultural Animals in Agricultural Research and Teaching (FASS, 1999) and followed the protocols and approaches of other published researches on genetic variation and population structure on several other livestock species (Dalvit et al., 2008(Dalvit et al., , 2009;;Zanetti et al., 2010;Maretto et al., 2012).

Animal sampling
The sampling strategy in agreement with the organization of livestock sector in Tunisia accounted for private herds and public herds of the Office of Public Lands (OTD).The public herds belonging to the Tunisian government included only two ecotypes of Barbarine: brown face (OBTR) and black face (OBTN).A total of 183 individual blood samples of BAR were collected from unrelated animals belonging to different populations distributed in the Tunisian territory: BAR brown face (BTR, n = 63) which includes NBTR, CBTR, SBTR, OBTR and LFB; BAR black face (BTN, n = 60) which includes NBTN, CBTN, SBTN, OBTN and CH; SG (n = 30) and SR (n = 30).Because of the absence of herd books, the animals of private herds were chosen as three unrelated animals from each farm or small flock based on the information provided by the farmer to avoid sampling of closely related individuals.Figure 1 shows the geographical areas from which the BAR individuals were sampled.31 individual blood samples were collected from several flocks of APP.The DNA extraction was carried out using the Wizard Genomic DNA Extraction Kit (Promega, USA) starting from 300 µL of whole blood.

Amplification and genotyping of microsatellite markers
A panel of 17 microsatellite markers was established according to ISAG/FAO Standing Committee (2004) recommendations and to previous studies (Baumung et al., 2006;Dalvit et al., 2008Dalvit et al., , 2009)).The following loci were used: OarAE54, OarFCB20, URB58, McM527, INRA23, TGLA53, MAF65, OarCP49, MAF214, HSC, INRA63, OarAE119, OarAE129, ILSTS087, OarFCB304, OarCP34 and CSRD247 (Table 1).Genotypes for all 17 microsatellite markers were determined by means of three multiplex fluorescent PCR reactions.Amplification was performed using standard PCR reactions in a GeneAmp 9700 thermal cycler (Life Technologies, USA) starting from 50 ng of purified DNA.The 17 microsatellites were amplified with the following conditions: initial denaturation step of 5 min at 95°C, 35 cycles of 30 s at 95°C, 1 min 30 s at 61°C and 30 s at 72°C and a final extension of 30 min at 60°C using the Type-IT Microsatellite PCR Kit (Qiagen, Hilden, Germany).Allele size was determined with a CEQ 8000 Genetic Analysis System (Beckman Coulter, Fullerton, CA, USA).

Statistical analysis
The total number of alleles per locus (TNA), allelic frequencies, observed (H o ) and expected (H e ) heterozygosity were calculated using GENETIX version 4.05.2 (Belkhir et al., 1996(Belkhir et al., to 2004)).Exact tests for deviation from Hardy-Weinberg equilibrium (HWE) were applied using a Markov Chain Monte Carlo simulation (100 batches, 5,000 iterations per batch and a dememorization number of 10,000) as implemented in GENEPOP version 3.4 (Raymond and Rousset, 1995).The MSA software (Dieringer and Schlötterer, 2003) was employed in calculations of allelic richness (AR, mean number of alleles per locus corrected by sample size), gene diversity (GD) (Nei, 1987) and Wright's fixation index [F IS (Weir and Cockerham, 1984)].Allelic richness and private alleles per population were calculated using rarefaction method to adjust for different population sizes using ADZE (Szpiech et al., 2008).Molecular coancestry coefficients (f ij ) and Tomiuk and Loeschcke (1995) genetic distances (D TL ) were measured using MOLKIN 3.0 (Gutierrez et al., 2005).D TL distances among populations were represented by a Neighbor-Net tree using SplitsTree4 (Huson and Table 1.Analyzed loci, chromosome (chr), fragment size (bp), total number of detected alleles (TNA), allelic richness (AR), gene diversity (GD) and fixation indices (F IS , F IT and F ST ) according to Weir and Cockerham (1984).
Fragment Bryant , 2006).Moreover, a factorial correspondence analysis (FCA) was performed based on individual genotypes using GENETIX version 4.03.The analysis of molecular variance (AMOVA) was performed by the ARLEQUIN software (Excoffier et al., 2005) using the codominant allelic distance matrix with 1000 permutations.
To study the population structure and to detect the most likely number of clusters (K) in the dataset, the software STRUCTURE version 2.2 (Pritchard et al., 2000) was used.To choose the appropriate number of inferred clusters to model the data, 50 independent runs were performed for each K cluster (2 < K < 10).All analyses used a burn-in period of 30,000 and 150,000 iterations for data collection.The best number of clusters fitting the data was established by plotting the Ln Pr(X|K) over the 50 independent runs for each K, as suggested by Pritchard et al. (2000) and also following Evanno et al. (2005).The output obtained was used directly as input by the cluster visualization program DISTRUCT (Rosenberg, 2004).

Genetic variability at microsatellite loci
The total number of alleles detected in the whole dataset was 270 and all markers were polymorphic in each population.Descriptive statistics on the variability of the investigated loci are reported in Table 1.The largest number of alleles was found at locus OarCP49 (29) and the smallest at loci OarAE129 and OarCP34 (9).The mean number of alleles per locus was 15.88 ± 5.15, and the mean allelic richness was 8.09 ± 0.27.The gene diversity across all loci was 0.812 ± 0.094, ranging from 0.538 (OarAE129) to 0.922 (OarCP49).The distribution of markers along the genome was quite wide involving 14 chromosomes with three of them (14, 19 and 5) interested by two markers (Table 1).Considering all the populations, F-statistics parameters were: F IS = 0.080 ± 0.010, F IT = 0.101 ± 0.010 and the F ST index was equal to 0.023 ± 0.003 (P < 0.001).To test if these parameters were affected by the presence of APP, which is an outgroup, they were computed again after removing this breed; the estimates were even smaller in BAR as F IS was equal to 0.078 ± 0.066, F IT = 0.084 ± 0.072 and F ST = 0.007 ± 0.010.

Breed variability and differentiation
The genetic variability of each subpopulation was initially studied in terms of the number of observed alleles and allelic richness, as shown in Table 2. BTN showed the largest number of alleles per locus (11.19 ± 3.99), followed by SG (11.06 ± 4.12) while APP (7.59 ± 2.43) and LFB (7.94 ± 2.33) showed the lowest values.Due to differences in sample sizes, the rarefaction method (Szpiech et al., 2008) was used in calculations of AR and private allelic richness (P AR ).BTN and SG showed the largest AR values calculated for 13 individuals (8.47 and 8.37, respectively) and APP and LFB the least (6.69 and 7.68, respectively).P AR was uniformly small and less than one in all BAR subpopulations and the highest P AR value was found in APP (0.49).In general, the BAR subpopulations showed considerable H o values: OBTN exhibited the highest H o (0.799 ± 0.161) followed by CH (0.772 ± 0.111) while APP and LFB exhibited the lowest H o (0.663 ± 0.141 and 0.703 ± 0.164, respectively).H o was always lower than the H e (Table 2).Most loci within populations tended to be within HWE (Table 2), the highest number of loci departing from HWE was found in APP (4 loci: MAF65, ILST087, MCM527, TGLA53) while in BAR populations only one locus was departing from HWE in BTR, SR and SG (OarAE129, CSRD247, MAF214, respectively).
Inbreeding (F IS ) was rather high for APP, LFB and SG showing an excess of homozygotes among loci (Table 2).Another way to measure within-populations diversity is the estimation of molecular coancestry, a measure of relatedness among individuals.Molecular coancestry estimates varied from 0.187 ± 0.023 (CH) to 0.198 ± 0.023 (BTN) in the BAR breed and were higher than in the APP breed (0.178 ± 0.024) (Table 2).The between populations F ST and the D TL genetic distances are given in Table 3.The higher F ST distance (always higher than 0.050) was found between APP and the group of BAR subpopulations whereas lower F ST estimates, ranging from 0.000 to 0.020 were found within BAR breed.
Tomiuk and Loeschcke's D TL genetic distance showed low values between BAR ecotypes ranging from 0.043 to 0.148.Between APP and BAR, this distance ranged from 0.186 and 0.296.The neighbor-net obtained from D TL distances (Figure 2) showed a clear distinction between APP and BAR subpopulations with APP placed at the end of the longest branch and BAR ecotypes placed in proximity to one another: OBTR and OBTN, BTR and BTN, SR and SG, CH and LFB.Variation via AMOVA procedure was 0.66% among BAR ecotypes, and 6.63% between BAR and APP.The factorial correspondence analysis was performed including all populations and loci using the corresponding allele frequencies (Figure 3).The first three components explained the 59.94% of the total variation, 36.85% of which explained by Axis 1 that clearly separates the BAR breed from the APP breed; 12.07% explained by Axis 2 that slightly separates the SR subpopulation; and 11.02% explained by Axis 3 that separates BAR individuals.It can be seen from the results of all the analyzed parameters that genetic diversity between BAR sub-populations is not significant.This breed seems to be genetically uniform population showing a high result between individual genetic diversity.

Population structure
Assignment test was performed using the program STRUCTURE with the number of expected  clusters (K) ranging from 1 to 10.The Ln Pr(G|K) increased from K = 2 to K = 4 and then reached a "plateau" at K = 6.According to Evanno et al. (2005), we assumed K = 3 as the most likely number of clusters describing our dataset (Figure 4).At K = 2, APP separates from BAR, while for K = 3, the diagram showed an identity between BAR strains.All individuals of all BAR subpopulations were represented as admixed individuals without showing any notable subdivision among the different groups.This revealed that there was a pattern of miscegenation during the development of BAR breed.At K = 4, a clear subdivision inside the APP population was identified but the BAR subpopulations maintained the previous clusterization (Figure 5).

DISCUSSION
In the present work, we investigated for the first time the genetic variability and the population structure of the most representative Tunisian native sheep breed, the Barbarine.This breed showed a high genetic variability.The high number of alleles for each locus as well as the high GD values confirmed that all microsatellite markers used were appropriate to analyze diversity in this Tunisian native breed.The BAR showed high withinbreed genetic diversity values.The genetic diversity estimates: mean number of alleles (MNA), AR, H e and H o were higher compared to APP.The BAR has not been genetically characterized before; therefore, it was not possible to compare our results with that reported in literature.However, previous analyses evaluating breeds originated from BAR, have shown lower genetic diversity measurements in US Tunis breed (Blackburn et al., 2011a, b) and Italian Barbaresca breed (Tolone et al., 2011).This loss of neutral genetic variability is probably due to a reduced effective population size of the founder populations.Compared to sheep breeds closer to the center of domestication, the BAR showed significantly  (Soma et al., 2012).Previous studies presented the same range of within-breed genetic variability in Italian sheep breed (Bozzi et al., 2009;Lasagna et al., 2011;Tolone et al., 2011) and in Austrian and Spanish sheep breeds (Baumung et al., 2006;Legaz et al., 2008), whereas Arora andBhatia (2004, 2006) and Sodhi et al. (2006) found lower values while investigating Indian sheep breeds.The large genetic diversity of the BAR was also demonstrated by Khaldi et al. (2010) and El Hentati et al. (2012) using RAPD markers.Estimates of F ST confirmed the very low level of differentiation within BAR, only 0.7% of the genetic variability was explained by difference among subpopulations or ecotypes.These findings are also supported by the AMOVA analysis; only 0.66% of the total variation was present among BAR ecotypes.El Hentati et al. ( 2012) showed that 6.06% of total variance was recorded between populations of BAR using RAPD markers.Moreover, results obtained with STRUCTURE evidenced the absence of clear clusters within the BAR breed (Figure 4) and also FCA analysis, f ij values found in the BAR ecotypes and the low values of D TL genetic distances between them further support the hypothesis of breed homogeneity.The F ST value increased to 2.35% when including APP indicating that genetic differentiation between these two breeds which originated from different geographic area, is notable.The D TL genetic distances supported the large diversity of APP from all BAR ecotypes.Moreover, APP showed the highest F IS value, the lowest AR, H o and H e values indicating the limited genetic variability present in this breed compared to BAR ecotypes.
This difference in the within-breed variability level can be partly explained by the differences in the management of flocks between Italy (Europe in general) and Tunisia.In fact, BAR is a large traditional population characterized by a traditional farming system with animals living mainly in extensive conditions and without any major selection pressure.This will explain this high within-breed degree of genetic variability (Lauvergne et al., 2000).While analyzing the APP breed, Lasagna et al. (2009) also found a similar value of F IS (0.118) confirming a high deficit of heterozygotes in this breed.The heterozygosity values within BAR showed differences i) between brown  (Pritchard et al., 2000).Mean ln Pr(G|K) values within each K (among 50 runs) are presented by solid circles.(b) ΔK values calculated following Evanno et al. (2005).
face and black face phenotypes: H o was higher in the black face group consisting of BTN, CH and OBTN compared with brown face group formed by BTR and OBTR and ii) between private animals and animals belonging to public herds (OBTR and OBTN).The Neighbor-Net based on D TL distances further confirms the distinctness of the public herds.Their closeness is probably due to the common management strategy followed by the OTD farms.Even though the samples of OBTN were collected from two geographically separated farms and those of OBTR from three separated farms, OBTR showed a lower genetic variability and a slightly higher inbreeding with respect to OBTN.This difference could be explained by the selection done in the sixties of the previous century to differentiate animals with black face from those with brown face (Palian, 1966).
Moreover, according to the breeders and Rekik and Ben Hammouda (2000), the more favorable areas of the north of Tunisia are suited for the black face strain of the Barbarine breeds.The level of inbreeding in BAR (F IS = 0.058) compared with that of other breeds is considered low: 0.300 in Egyptian fat tailed sheep (El Nahas et al., 2008); 0.469 in Hamdani fat tail Irakian sheep (Al-Barzinji et al., 2011); 0.030 to 0.080 in Italian (Sicilian) breeds (Tolone et al., 2011); 0.078 to 0.118 in native sheep breeds undergoing in situ conservation (Dalvit et al., 2009).This decreased level of inbreeding in BAR could be mainly explained by the mating management of this breed.In fact, the shepherd assistance is needed to lift the fat tail at copulation time.Consequently, this morphological trait of BAR should confer to breeders a paternity controlling convenience and mating between close relatives is generally avoided.However, the rather high level of inbreeding in LFB and SG populations could be due to the farming practices carried out in the towns (Centre Tunisia) where the samples of theses populations were collected: small flocks (rams and ewes) belonging to neighboring breeders are usually reared together in extensive pastures allowing for frequent mating with close relatives.Mean molecular coancestry estimates within BAR were comparable to that of APP but rather low if compared with the results obtained by Dalvit et al. (2008) in Alpine sheep breeds.These results also emphasize the higher variability found in the BAR breed.
The FCA analysis, STRUCTURE and D TL genetic distances suggested a clear genetic separation between APP and BAR.The neighbor-net (Figure 2) obtained from D TL distance estimates grouped the BAR ecotypes together showing the close relationship among them.APP appeared at the end of the longest branch confirming its use as an out-group.Clustering obtained by STRUCTURE revealed the high homogeneity of the BAR breed since no clear subdivision in ecotypes was found.At K = 3, several BAR individuals showed admixed pattern which can be considered as an introgression in different proportions of genetic material from other population.This admixed pattern may be related to the old miscegenation of BAR breed with a very long thintailed breed which replaced the BAR since 300 B.P., as depicted in Phoenician and Roman monuments (Khaldi, 1989).Crossbreeding between BAR and Thin-tail Algerian breed in particular, commonly practiced in central and northern Tunisia in the past few decades, could further explain this introgression.This uncontrolled crossbreeding raises an important threat to the BAR integrity; hence, the necessity to take more appropriate actions.Maintaining population size and genetic variability through planned mating is therefore needed to preserve this genetic resource.
The valorization of BAR breed is also a step towards its sustainability.In fact, this breed continues to be the main source for commercial meat providers and it will significantly contribute to meet the objectives of the national strategy in red meat self-sufficiency (Ben Salem et al., 2011).A survey done by Bedhiaf-Romdhani et al. (2008b) revealed that consumers do still prefer the BAR's meat for its tenderness, flavor and smell.Socio-cultural value of this breed still preferred as sacrificial animal in religious and family celebrations is to be considered for conservation priorities purposes.Furthermore, Ben Salem et al. (2011) asserted that the recent mutations of the production systems induced by climate change and social constraints are in favour of the BAR breed.This molecular investigation can be helpful to provide the information necessary for an appropriate management of this autochthonous breed in order to improve BAR breed productive performance.

Conclusion
This study presents the first results regarding the characterization of genetic variability using microsatellite markers of the most dominant sheep breed in Tunisia.Traditionally, classification of BAR breed was based on visible phenotypic traits thus several ecotypes were identified.However, molecular characterization proved the homogeneity of BAR breed and its high level of genetic diversity.In fact, there is no basis for considering varying ecotypes of BAR sheep; they could be considered as one population.This consolidation has several appealing aspects for conservation, in that flocks across the country can be more easily managed to maintain genetic variability.This breed seems to have a large reservoir of genes able to be oriented in several selection directions.However, the BAR seems to be threatened by losing its gene pool by crossbreeding and it starts to be considered as a less productive breed.This breed is increasingly neglected by farmers who either change it through crossbreeding or replaced altogether with more productive ones.To minimize the potential risk of genetic diversity loss in this breed, it is necessary to give it a priority for conservation in order to preserve its global native genetic diversity.Therefore, urgent steps must be taken in the near future to ensure its maintenance and appropriate genetic management.

Figure 3 .
Figure 3. Spatial representation of the two breeds Barbarine (BAR) and Appenninica (APP), in green, as defined by the factorial correspondence analysis.

Figure 4 .
Figure 4.Estimated posterior probabilities of ln Pr(G|K).(a) ln Pr(G|K) values are presented as a function of the number of clusters(Pritchard et al., 2000).Mean ln Pr(G|K) values within each K (among 50 runs) are presented by solid circles.(b) ΔK values calculated followingEvanno et al. (2005).

Table 2 .
Number of analyzed samples, mean number of alleles (MNA), allelic richness obtained with rarefaction method (AR), private allelic richness (P AR ), observed (H o ) and expected (H e ) heterozygosity, within-population heterozygote deficiency (F IS ), number of loci deviated from Hardy-Weinberg equilibrium (HWE) and within-population molecular coancestry (f ij ).

Table 3 .
Tomiuk andLoeschcke's distances (D TL ) below diagonal and F ST (P < 0.05 in italic) above diagonal between the analyzed populations.