Genetic diversity and structure of 2 indigenous sheep breeds (Kotel and Teteven) in Bulgaria using microsatellite markers

ABSTRACT

This study analyzed the genetic diversity and structure of two Bulgarian indigenous sheep (Ovis aries) breeds – Kotel and Teteven, each presented by 4 flocks, using 15 microsatellite markers. A total of 195 alleles were identified. Among the markers used, the most informative were INRA5 and OarFCB20 with PIC values of 0.87. Both breeds showed a high level of genetic diversity with average values of He and Ho, respectively 0.74 and 0.76. The observed high level of genetic diversity is combined with the absence of heterozygous deficiency reported as Fis values of 0.04 and 0.07, respectively for Kotel and Teteven breeds. Low level of genetic differentiation between the two breeds was observed, as evidenced by the low value of Fst of 0.09, and analysis of variance (ANOVA) data showing that only 9% of the total variation is due to variation between breeds, whereas 91% is due to variation within breeds. Furthermore, we revealed a breeding practice that comprises an exchange of animals between the breeds which reflects in the clustering obtained by STRUCTURE analysis. According to it, Kotel breed is a more homogenous population and forms cluster 1, whereas Teteven breed is more heterogenous and consists of two subpopulations that form two distinct clusters. The results obtained in the present study confirm the high potential of microsatellite markers in studies related to both breed affiliation and purity and can serve as a basis for the implementation of future management and conservation programs of these indigenous breeds.

KEYWORDS:

• Sheep breeds (O. aries)
• microsatellite markers (SSR)
• genetic diversity
• genetic structure

Introduction

Recent studies in livestock genetic resources in Europe unambiguously emphasize the importance of local breeds for maintaining genetic diversity in order to respect different production, cultural and region needs [1]. Local livestock genetic resources are important sources of genetic variation for providing genetic improvement, preventing diseases and loss of genetic diversity [2]. Diversity of indigenous breeds contributes significantly to both national and regional food quality and variety. In this line, the conservation and sustainable utilization of farm animal genetic resources are of vital importance. Knowledge about the genetic diversity of a population is a crucial parameter for the implementation of successful genomic selection and conservation of livestock genetic resources. The ability of the genetic assignment of individual members to a target breed is essential to the success of the conservation initiatives. It should be based on in-depth knowledge of the genetic resources of the specific breed [3]. Therefore, it is important to genetically characterize indigenous breeds using the recent advantages of molecular marker techniques.

Microsatellite markers (SSR, simple sequence repeats) have been extensively used for characterization of genetic diversity of indigenous breeds because of their high degree of polymorphism, co-dominant inheritance, locus specificity and extensive genome coverage. Recently, microsatellite markers were successfully applied for characterization of sheep (Ovis aries) genetic diversity, population structure and genetic differentiation, aiming development and application of genetic conservation strategies in different European countries like Austria, Hungary, Italy, Spain, Switzerland, [4–9] as well as in the countries from the Balkan Peninsula like Greece [10–15].

In Bulgaria, sheep breeding is the most ancient branch of animal husbandry. Despite its relatively small territory, the country has a rich diversity of autochthonous domestic livestock breeds, of which seventeen are local sheep breeds [16–17]. These breeds are well adapted to the various climatic conditions of the country and they represent a reservoir of important genes in breeding programs under the changeable climatic conditions in recent years. They have been phenotypically described and characterized, but these data do not give information about genetic variability within and between studied breeds which is important for their preservation. In Bulgaria, microsatellite markers were used for the first time in sheep with the aim to determine genetic relationship and genetic structure among 7 indigenous Bulgarian sheep breeds in order to facilitate and plan their sustainable development, utilization and conservation [1, 18]. However, there is still lack of information about the genetic diversity and structure of the remaining 12 indigenous sheep breeds which is important for their effective preservation and conservation.

The objective of this study was to assess the genetic diversity of two indigenous sheep breeds in Bulgaria, namely Teteven and Kotel and to determine the genetic structure of two populations. In addition, the purity and assigning of animals of several flocks to a specific breed is discussed.

Materials and methods

Sampling and DNA extraction

Ninety-six blood samples were collected from two indigenous sheep breeds, of which 48 were from Teteven and 48 were from Kotel sheep breeds. The breeds have different habitats in the region of the Balkan Mountains in Bulgaria. Teteven and Kotel breeds belong to long, thin-tailed group of sheep breeds. The main area of Teteven breed distribution includes the Central Stara Planina (Balkan Mountains) near the town of Teteven and the villages of Ribaritza, Cherni Vit, Golyam Izvor and Galata. Kotel breed is a representative of the sheep of the Tsakel type. The area of distribution is in the mountainous and semi-mountainous regions of Kotel Mountain and the adjacent territories.

Each sheep breed was presented by 12 unrelated individuals of 4 independent flocks. Genomic DNA was extracted using DNeasy Blood & Tissue Kit (Qiagen) according to the manufacturing instructions. The isolated DNA samples were assessed for their quantity and quality by spectrophotometric measurements (NanoDrop Thermo Scientific, Wilmington, DE) and were checked through electrophoresis in 0.8% agarose gel.

Microsatellite amplification

A total of 15 microsatellite markers recommended by the Food and Agriculture Organization (FAO) (http://www.fao.org/dad-is/) and International Society for Animal Genetics (ISAG) (http://www.isag.org.uk/Docs/2005_Panels Markers Sheep Goats. Pdf) based on their level of allelic diversity were selected as recommended to be used in genotyping and paternity testing analyses (Table 1).

Table 1. SSR locus information.

Marker name	Chromosome location	Origin	Primer sequences F / R
CSRD247	14	Ovine	F: GGACTTGCCAGAACTCTGCAAT
			R: CACTGTGGTTTGTATTAGTCAGG
HSC	20		F: CTGCCAATGCAGAGACACAAGA
			R: GTCTGTCTCCTGTCTTGTCATC
INRA0063	14	Bovine	F: ATTTGCACAAGCTAAATCTAACC
			R: AAACCACAGAAATGCTTGGAAG
MAF0214	16	Ovine	F: GGGTGATCTTAGGGAGGTTTTGGAGG
			R: AATGCAGGAGATCTGAGGCAGGGACG
OarAE0129	5	Ovine	F: AATCCAGTGTGTGAAAGACTAATCCAG
			R: GTAGATCAAGATATAGAATATTTTTCACACC
OarCP0049	17	Ovine	F: CAGACACGGCTTAGCAACTAAACGC
			R: GTGGGGATGAATATTCCTTCATAAGG
OarFCB0011	2	Ovine	F: GCAAGCAGGTTCTTTACCACTAGCACC
			R: GCCTGAACTCACAAGTTGATATATCTATCAC
OarFCB0304	19	Ovine	F: CCCTAGGAGCTTTCAATAAAGAATCGG
			R: CGCTGCTGTCAACTGGGTCAGGG
D5S2	5		F: TACTCGTAGGGCAGGCTGCCTG
			R: GAGACCTCAGGGTTGGTGATCAG
INRA0005	12	Bovine	F: TTCAGGCATACCCTACACCACATG
			R: AAATATTAGCCAACTGAAAACTGGG
INRA0023	3	Bovine	F: GAGTAGAGCTACAAGATAAACTTC
			R: TAACTACAGGGTGTTAGATGAACTC
MAF0065	15	Ovine	F: AAAGGCCAGAGTATGCAATTAGGAG
			R: CCACTCCTCCTGAGAATATAACATG
McM0527	5		F: GTCCATTGCCTCAAATCAATTC
			R: AAACCACTTGACTACTCCCCAA
OarFCB0020	2	Ovine	F: AAATGTGTTTAAGATTCCATACATGTG
			R: GGAAAACCCCCATATATACCTATAC
SPS0113	10		F: AAAGTGACACAACAGCTTCTCCAG
			R: AACGAGTGTCCTAGTTTGGCTGTG

The primers for amplification of the corresponding microsatellite loci were 5’-labeled with either VIC, PET or FAM. For amplification of 14 of the SSR loci 4 mulitiplex reactions were performed as follows:

Multiplex А: D5S2(FAM), INRA5(VIC), MAF65(PET);

Multiplex B: OarAE129(PET), OarFCB304(VIC), OarFCB11(FAM);

Multiplex C: INRA23(PET), OarFCB20(FAM), SPS113(FAM), McM527(VIC);

Multiplex D: CSRD247(VIC), HSC(FAM), MAF214(PET), OarCP49(VIC)

Primer pair INRA63 (FAM) was used separately in polymerase chain reaction (PCR).

PCR was performed in a total volume of 20 µL containing 20 ng DNA template, 10 µL 2XMyTaq™ HS Mix (BIOLINE), 10 pmol of each primer (F and R) and miliQ H₂0 up to the final volume.

PCR amplification was performed on an Axygen™, MaxyGene II Thermal Cycler using the following PCR conditions:

For multiplexes B and D: 1 cycle of denaturation at 95 °C for 12 min, followed by 31 cycles each consisting of 20 s at 95 °C, 1 min at 63 °C, 1 min at 72 °C, and 5 min extension at 72 °C.

For multiplexes A and C: 1 cycle of denaturation at 95 °C за 10 min, followed by 31 cycles each consisting of 30 s at 95 °C, 30 s at 55 °C, 1 min at 72 °C, and 5 min extension at 72 °C.

For amplification of the INRA63 locus, the PCR program includes 1 cycle of denaturation 95 °C for 12 min, followed by 31 cycles each consisting of 20 s at 95 °C, 1 min at 58 °C, 1 min at 72 °C, and a final step – extension at 72 °C for 5 min.

Fragment analysis

Genotyping of the microsatellite markers was carried out on an automated ABI PRISM 3130XL Genetic Analyzer (Applied Biosystem, USA) using LIZ 500 (Applied Biosystem, USA) as the internal size standard. Electropherograms were genotyped automatically with the GeneMapperTM software v5.0.

Statistical analysis

GeneAlex v6.5 [19] was used for calculation of Nei’s genetic distance, expected (He) and observed heterozygosity (Ho), fixation indices (Fis, Fit and Fst) as well as for performing AMOVA and Principal Coordinate Analysis. The polymorphic information content (PIC) per each SSR marker was calculated using Powermarker 3.25 [20]. Structure v 2.3.4 was used for analysis of the genetic structure [21] where Admixture was used as ancestry model, Length of Burnin Period was set to 100000 and number of MCMC Reps after Burnin was set to 200000. The number of presumed clusters (K) was set from 1 to 16. Ten iterations were performed for each K value. The method by Evanno et al. [22] was used to estimate the most probable number of clusters (K) using Structure harvester v 0.6.94 (http://taylor0.biology.ucla.edu/structureHarvester/). Clumpak (http://clumpak.tau.ac.il/) was used to obtain an average result over the 10 simulation runs at the identified most probable K.

Results

Population genetic variability

A total of 195 alleles were identified in the entire sheep population consisting of 96 animals genotyped at 15 microsatellite loci. The number of alleles (Na) varied from 5 in locus OarAE129 to 25 in locus OarCP49 with an average of 13 alleles/locus. The average Ne is 4.33. Parameters revealing the genetic variability in the studied population for each locus are given in Table 2.

Table 2. Allele range, number of identified alleles (Na), number of effective alleles (Ne), mean number of alleles per locus (Mn), heterozygosities: observed (Ho), expected (He), Polymorphism Information Content (PIC) and F-values in each locus.

	Allele range	Na	Mn	Ne	Ho	He	PIC	Fis	Fit	Fst
D5S2	188 − 202	7	4.63	3.24	0.65	0.67	0.71	0.04	0.13	0.10
INRA5	116 − 152	13	7.88	5.30	0.80	0.80	0.87	0.00	0.09	0.09
MAF65	125 − 143	10	5.63	3.86	0.75	0.73	0.75	−0.03	0.05	0.07
OarAE129	139 − 169	5	3.50	2.47	0.52	0.58	0.56	0.11	0.18	0.08
OarFCB304	156 − 186	13	6.38	3.72	0.77	0.70	0.74	−0.11	−0.01	0.09
OarFCB11	120 − 144	12	5.75	4.23	0.85	0.76	0.81	−0.12	−0.03	0.08
INRA23	200 − 224	13	7.13	5.23	0.81	0.78	0.84	−0.04	0.05	0.09
OarFCB20	90 − 116	13	7.25	5.51	0.88	0.81	0.87	−0.08	0.01	0.08
SPS113	236 − 251	11	6.25	4.27	0.73	0.75	0.79	0.02	0.10	0.07
McM527	161 − 248	12	6.38	4.86	0.82	0.77	0.8	−0.07	0.00	0.07
CSRD247	209 − 265	13	6.38	3.82	0.65	0.73	0.8	0.11	0.21	0.11
HSC	267 − 297	14	7.88	5.27	0.74	0.80	0.86	0.07	0.15	0.09
MAF214	173 − 265	16	6.25	3.37	0.70	0.65	0.71	−0.07	0.04	0.11
OarCP49	80 − 134	25	9.00	4.88	0.85	0.77	0.84	−0.10	0.00	0.09
INRA63	157 − 209	18	8.38	4.93	0.85	0.78	0.84	−0.10	0.00	0.09
Mean		13	6.58	4.33	0.76	0.74	0.79	-0.02	0.06	0.09
Total		195

The expected heterozygosity (He), which is considered the best estimator of genetic diversity in the population, ranged from 0.81 in locus OarFCB20 to 0.58 in locus OarAE129 with a population mean of 0.74 for the 15 loci analyzed in the population. High level of He was also observed in locus INRA5 (0.80), HSC (0.80) and loci INRA23 (0.78) and INRA63 (0.78). On the other hand, the observed heterozygosity (Ho) ranged from 0.52 in locus OarAE129 to 0.88 in locus OarFCB20 with a population mean of 0.76, indicating high level of genetic variability in the breeds. The polymorphic information content (PIC) ranged from 0.87 for OarFCB20 and INRA 5 to 0.56 for the marker OarAE129. The average PIC for the 15 microsatellite markers was 0.79 and no marker showed PIC less than 0.5, indicating that all loci were highly polymorphic.

The average Fis value was −0.02 (p < 0.05). Two markers, namely OarAE129 and CSRD247, showed Fis values higher than 0.1 (Fis = 0.11), indicating an excess of homozygous individuals in the populations for these markers. The Fit fixation index, which measures the heterozygosity loss of the individuals with respect to the overall population, was 0.06 (p < 0.05), indicating a 6% general deficit of heterozygous individuals in the sheep population. The mean Fst index, which measures the degree of genetic variability explained by differences between both breeds, was 0.09. Therefore, it is evident that the total genetic variation mostly corresponds to differences among individuals (91%) and only 9% is the result of differences between breeds. This result indicates that genetic variation measured by the microsatellite markers showed low level of genetic differentiation between the studied breeds.

Genetic variability among Kotel and Teteven breeds

The parameters used to evaluate the genetic diversity for each breed are presented in Table 3. The number of identified alleles for Kotel and Teteven breed were 170 and 139 respectively. Kotel breed had a higher mean number of alleles per locus [11.33] in comparison to Teteven breed (9.27). Higher values for Ho and He (0.78 and 0.81, respectively) were observed in Kotel breed, whereas these values in Teteven breed were slightly lower (0.72 and 0.77, respectively). There is no large difference between Ho and He values which means that both breeds showed high genetic diversity in general. The inbreeding parameter (Fis) is close to zero in both sheep breeds which means low level of heterozygote deficit. Fis is slightly higher in Teteven (0.07) than in Kotel breed [0.04].

Table 3. Number of alleles, mean number of alleles, mean observed (Ho) and expected (He) heterozygosities, coefficient of inbreeding (Fis) in Kotel breed (KOT) and Teteven breed (TET).

		Mean
Breed	Number of alleles	Na	Ne	Ho	He	Fis
KOT	170	11.33	5.77	0.78	0.81	0.04
TET	139	9.27	4.89	0.72	0.77	0.07
Mean	154.50	10.3	5.33	0.75	0.79	0.05

In general, both breeds showed high genetic diversity for all loci analyzed in this study. Hardy Weinberg equilibrium tests for the two breeds from 15 microsatellite markers showed low deviations (p < 0.001) in only 3 of 30 HW tests performed (Table 4), which is due to the absence of heterozygote deficit in each breed.

Table 4. Hardy Weinberg (HW) equilibrium test by breed.

Breed	Locus
	D5S2	INRA5	MAF65	OarAE129	OarFCB304	OarFCB11	INRA23	OarFCB20	SPS113	McM527	CSRD247	HSC	MAF214	OarCP49	INRA63
KOT	0.432	0.456	0.998	0.847	0.991	0.876	0.865	0.437	0.654	1.000	0.000***	0.427	0.834	0.999	0.393
TET	0.053	0.558	0.829	0.220	0.712	0.638	0.771	0.727	0.227	0.001***	0.043	0.000***	1.000	1.000	0.166

*** significant at p < 0.001.

Genetic variability within and among the flocks of Kotel and Teteven breeds

Both breeds are presented by 4 flocks respectively (Table 5). The mean number of identified alleles (Na) varied from 5.27 in Flock 8 (Teteven breed) to 7.47 in flock 2 (Kotel breed) (Table 5). Na was higher in all studied flocks (FL1-FL4) of Kotel breed in comparison to the FL5-FL8 of Teteven breed, as is shown in Table 3. Moreover, our study showed the absence of 6 alleles in 5 out of 15 microsatellite loci as well as the presence of new alleles with the frequency of 0.17-0.21 in 3 out of the tested here microsatellite loci in flock 8 which were not found in the flocks of Kotel breed and the remaining 3 flocks of Teteven breed.

Table 5. Number of identified alleles (Na), heterozygosities: observed (Ho), expected (He), and mean coefficient of inbreeding (Fis) in each flock.

Flocks/breeds*	Na	Ne	Ho	He	Fis
FL1 - Kot	7.33	4.80	0.76	0.77	0.02
FL2 - Kot	7.47	5.33	0.79	0.78	−0.01
FL3 - Kot	6.87	4.60	0.77	0.77	−0.01
FL4 - Kot	7.07	4.94	0.80	0.78	−0.03
FL5 - Tet	5.67	3.70	0.70	0.71	0.01
FL6 - Tet	6.07	3.54	0.67	0.69	0.03
FL7 - Tet	6.87	4.38	0.78	0.75	−0.03
FL8 - Tet	5.27	3.36	0.79	0.68	−0.16

* FL1-FL4 indicates the studied 4 flocks of Kotel (Kot) breed; FL5-FL8 indicates the studied 4 flocks of Teteven (Tet) breed.

The expected heterozygosity (He) ranged from 0.68 in Flock 8 (Teteven breed) to 0.78 in flock 2 and flock 4 (Kotel breed). The lowest level of observed heterozygosity (Ho) (0.67) was found in flock 6 (Teteven breed) and the highest (0.80) in Flock 4 (Kotel breed). Generally, Ho and He were higher in all four flocks of Kotel breed. Mean Fis varied from −0.16 in flock 8 to 0.03 in flock 5 of Teteven breed. Positive but lower than <0.1 Fis values were observed for FL1 (0.02, Kotel breed) as well as for FL5 (0.01) and FL6 (0.03) of Teteven breed, which does not indicate a risk of inbreeding in any of the flocks.

Genetic differentiation and admixture analysis

To study the genetic differentiation and admixture among the flocks of both breeds, genetic structure analysis using the Structure software and Principal Coordinate Analysis (PCoA) were performed. The results from applying the method by Evanno et al. [22] showed that the most probable number of genetic clusters is three (K = 3) (Figure 1a). Figure 1b illustrates the genetic structure at K = 3 for all 96 samples of the 8 analyzed flocks. As is seen from Figure 1b the flocks from Kotel breed (FL1-FL4) showed highly homogenous structure with almost all samples belonging predominantly to a single cluster. On the contrary high heterogeneity of the population was observed in the flocks from Teteven breed. Although the majority of animals from flocks 5, 6 and 7 belonged to the second major cluster, which is separated from the Kotel breed cluster, five of the animals from these three flocks (one from flock 6 and four from flock 7) belonged to the same genetic cluster as the samples from the Kotel breed. Moreover, all animals of flock 8 from Teteven breed belonged to a third major cluster with a genetic structure which was different from the rest of the flocks of the two analyzed breeds. Similar results were obtained applying PCoA as shown in Figure 2 where again three major clusters were observed with flock 8 of Teteven breed separated from the other two identified clusters.

Figure 1. Genetic structure of 8 sheep flocks representing two Bulgarian sheep breeds based on SSR marker data. a) Estimation of the most probable K using the delta K method by Evanno et al. [22], b) Bar plot representing the genetic structure at K = 3. FL1-FL4 represent data for flocks from Kotel breed, FL5-FL8 represent data for flocks from Teteven breed.

Figure 2. Principal coordinate analysis of 8 sheep flocks representing two Bulgarian sheep breeds based on SSR marker data. FL1-FL4 represent data for flocks from Kotel breed, FL5-FL8 represent data for flocks from Teteven breed.

We further compared the Fst values between the studied flocks as a measure for their genetic differentiation (Table 6). In general, the flocks from Kotel breed showed low genetic differentiation with Fst values ranging between 0.024 and 0.035 suggesting homogenous structure within the breed. The flocks from Teteven breed showed higher Fst values ranging between 0.046 and 0.092. As seen from Table 6 the highest Fst values were scored between flock 8 and the rest of the flocks from the Teteven breed (0.067-0.092). Surprisingly, the Fst values between flock 8 and the flocks from Kotel breed (0.053-0.066) were lower. As a whole the results from the pairwise Fst comparison between the flocks were in complete consistence with the analysis of the genetic structure (Figures 1 and 2).

Table 6. Pairwise comparison of Fst values among 8 flocks from Kotel and Teteven breed.

FL1	FL2	FL3	FL4	FL5	FL6	FL7	FL8
0.000								FL1
0.024	0.000							FL2
0.032	0.035	0.000						FL3
0.033	0.034	0.030	0.000					FL4
0.062	0.057	0.056	0.048	0.000				FL5
0.055	0.060	0.056	0.059	0.054	0.000			FL6
0.050	0.043	0.040	0.033	0.046	0.050	0.000		FL7
0.065	0.066	0.062	0.053	0.089	0.092	0.067	0.000	FL8

Discussion

Population genetic diversity

The study of genetic variation plays an important role in developing rational breeding strategies for economical animal species [23]. In the present study the genetic diversity of two indigenous sheep breeds in Bulgaria was assessed by using 15 microsatellite markers from the ISAG/FAO recommendation list for genotyping and parentage testing in sheep.

The level of variation depicted by the number of alleles at each locus serves as a measure of the impact of each of the studied breeds on differentiation within livestock populations. The average Na found in this study is higher than that observed by Nanekarani et al. [24] with 15 microsatellite markers in a population of 360 sheep from four breeds in Iran and by Mastranestasis et al. [11] with 11 SSR markers in the Levsos breeds population in Greece. However, this value is lower than the observed average Na value of 13.4 in 362 sheep that originated from three Colombian sheep breeds distributed over 43 farms in 11 Colombian provinces [25] using 10 microsatellite markers. The comparison of the Ne with the number of observed alleles per locus provides evidence for the predominance of certain alleles in each breed.

The polymorphic information content (PIC) is an indicative parameter for informativeness of microsatellite markers and their usefulness in diversity analysis of a breed. In our study, the high values of PIC (> 0.5) [26] and the high average number of alleles per locus indicate that the panel of 15 microsatellite markers used is suitable for studies of genetic diversity in Bulgarian sheep. In this sense, INRA5, OarFCB20 and HSC could be considered the most informative markers of our test panel. Previous studies on the genetic diversity and structure in 7 indigenous sheep breeds [1] using 6 SSR markers also showed high level of polymorphism in locus OarFCB20. An average value of 0.79 for the PIC once again indicated abundant genetic diversity in the Bulgarian sheep population as well as the utility of these markers for population assignment [27] as well as genome mapping [28] studies in addition to genetic diversity analysis.

The expected heterozygosity (He), which is considered the best estimator of genetic diversity in a population [29] showed a value of 0.74, indicating that the indigenous sheep population in Bulgaria has high genetic variability. This is very promising for the conservation of the indigenous sheep breeds as a genetic resource since high levels of genetic variability allows a better adaptation to environmental changes which have become more evident in recent years and reduces the negative effects of inbreeding that are common in livestock farming. The observed heterozygosity (Ho) ranged between 0.52 (OarAE129) and 0.88 (OarFCB20), with an average of 0.76 for the 15 loci analyzed in the population. It is observed that in 8 out of 15 microsatellite loci studied here Ho was higher than He, which means an absence of heterozygous deficit. Overall, the observed heterozygosity is higher than in other sheep breeds from Columbia [25], China [30], Sudan [31], India [27], similar to that reported for 10 native Greek breeds [10] but lower than that observed in the Levsos sheep breed from Greece [11].

The average Fst over the loci and samples was 0.9 showing that most genetic differentiation is explained within breeds. The low Fst value suggests that local breeds are not well differentiated, which could be explained by the common history and breeding practice. Additionally, the admixture could be result of the high flow of genes between farms due to the practice to borrow or lease rams and lack of effective control in breeding programs [25] in Bulgaria. Similar results have been reported by Stahlberger-Saibetkova et al. [7] and Al-Atiyat et al. [32]. However, lower estimates of genetic differentiation between sheep breeds were found by Ocampo et al. [25], Bozzi et al. [5], Zhong et al. [30], Álvarez et al. [33], Ben Sassi et al. [34], Hoda et al. [14] and Abdelkader et al. [35], who reported Fst ranging between 0.05 and 0.3. The Fis value for the overall population (-0.02) is lower than the values reported in Turkish sheep breeds [36], Tunisian sheep breeds [37], Moroccan sheep breeds [38] and seven indigenous Bulgarian sheep breeds not included in this study [1]. The Fis values obtained according to the breeds studied here indicated absence of a loss of heterozygosity in contradiction to many other studies in sheep that reported the presence of heterozygous deficit that could be generally explained by subdivision among flocks (Wahlund effect) and uncontrolled cross breeding without respect to the distribution area of each breed by the breeder.

Within-breed genetic differentiation

Both breeds showed a high mean number of alleles. The Kotel breed had higher MNA (11.33), which in turn is higher than the values reported by Hristova et al. [1] for Breznishka (10 alleles/locus), Copper-Red Shumenska (10.5 alleles/locus), Karakachanska (9.5 alleles/locus), Local Karnobatska (9.17 alleles/locus), Blackhead Plevenska (10.33 alleles/locus) and Starozagorska (9.5 alleles/locus0 sheep breeds in Bulgaria.

The absence of deviation found in HW tests and the low values of Fis found in each breed indicates a very low deficit of heterozygote individuals in contradiction to the study of Hristova et al. [1] in seven Bulgarian indigenous sheep breeds and of Kusza et al. [39] on the genetic relationship among five Bulgarian sheep breeds using 16 SSR markers. The absence of heterozygote deficit is due to the negative Fis values in 3 out of 4 tested flocks of Kotel breed and in 2 out of 4 tested flocks of Teteven breed. In these flocks, the Ho was higher or equal to the He. The high values of the observed heterozygosity and the negative values of Fis could be explained by the high gene flow between the flocks of the breeds. Very low Fis values were also observed in the study of Mastranestasis et al. [11], on genetic diversity and structure of the Levsos sheep breed in Greece including 10 flocks on the island of Levsos. The authors reported significantly negative Fis for 5 out of the 10 flocks and consistently higher observed heterozygosity than that expected for all flocks.

Furthermore, our results revealed a high rate of animal exchange not only between neighbouring and more distant farms in the frame of the breed but also between breeds, which suggests a lack of appropriate management of the conservation programs of local sheep breeds. Our study showed that the genetic differentiation of the studied flocks of both breeds is lower in Kotel breed than in Teteven breed. Among the studied four flocks of Teteven breed, the highest Fst was observed for flock 8, which means that this flock is isolated and the exchange of animals with other flocks of Teteven sheep bred is restricted. Moreover, this flock is characterized by a high proportion of alleles not observed in any of the flocks of Kotel and Teteven breeds. The presence of alleles specific only for this flock was also observed. This is also reflected in the clustering obtained by the PCoA and STRUCTURE analyses. The STRUCTURE analysis showed a specific grouping of the eight flocks in three clusters in which all of them had different proportions of membership. The first cluster corresponds to the flocks of the Kotel breed, which is the most homogenous, consistent with the observed low genetic differentiation between flocks. The second cluster includes three flocks of Teteven breed. It is the most heterogeneous one because of the higher observed gene flow from the flocks of the Kotel breed and from an unknown sheep breed previously referred to as flock 8 of the Teteven breed. This flock comprises animals that are more genetically distant from the Teteven breed and forms an additional cluster, which is almost homogenous. This indicates that this flock consists of animals not belonging to both of the studied breeds. Our study showed that this flock has overlapping habitats with other local sheep breeds like Karakachanska and Central Stara Planina sheep breeds. Additional studies including the breeds having overlapping habitats will reveal the possible assigning of the animals of that flock to any of the breeds in the region.

The study showed that the accurate management of the purity of the breeds is an important tool for the implementation of the management practice for the proper conservation of indigenous sheep breeds in Bulgaria. In this case, the use of molecular markers is an appropriate tool for both the study of genetic diversity and admixture analysis and for the determination of the purity of the breeds. They allow the assigning of individual animals to specific breeds, which sometimes is impossible with the use of phenotypic tools.

Conclusions

The study showed that the genetic variability in Bulgarian sheep breeds is ascribed by the within-breed variability. The low Fst value of 0.09 in the overall population in this study, suggests that local breeds are not differentiated enough. Lower differentiation was observed in Teteven than in Kotel breed with the exception of flock 8 of Teteven breed. The number of unique alleles, as well as the specific grouping of this flock in a separate cluster not belonging to the clusters of both breeds studied here, suggests its different origin. This shows that molecular markers are an appropriate tool for assigning animals/flocks to specific breeds and determining the admixture processes in breeding practice where there is no strong control on the proper management of the local sheep breed conservation. Further analysis of the neighbour breeds from the region will allow the correct assigning of the animals of flock 8 of Teteven breed.

Disclosure statement

The author declares no conflicts of interest.

Data availability

All data that support the findings reported in this study are available from the corresponding author upon reasonable request.