Methods for Genotyping of the Honey Bee (Apis mellifera L.: Hymenoptera: Apidae) in Bulgaria

Abstract

Honey bees are insects of great biological, ecological, and economic importance. They are the subject of a variety of scientific studies. As social insects, they are a preferred and widely used model for clarification of the evolution of social behavior. Because of the haplodiploidy as a mechanism for determining gender, differentiation in the functions of individuals in the bee family, and for its economic importance, Apis mellifera is a significant subject of research in the fields of ontogeny, population genetics, and selective breeding. The biological significance of bees is rooted in the fact that they are main pollinators in the natural environment. About 80% of the pollination of entomophilous plants is carried out by Apis mellifera. In all crops, active pollination significantly increases their yields. Honey bees are a valuable economic asset due to the ensemble of their products, which include honey, bee pollen, propolis, royal jelly, and bee venom, used by humans for food and treatment. The main objective of this chapter was to describe the basic methods used for genotyping of Apis mellifera in Bulgaria. These techniques have been useful to produce a system of population criteria, and taxonomically important molecular markers are applicable in future activities related to the preservation and selection of the Bulgarian honey bee.

1. Introduction

The geographic distribution of honey bee Apis mellifera L. includes Africa, Europe, and the Near East regions. At present, 29 subspecies of A. mellifera are recognized based on morphometric characters [1, 2, 3, 4]. These subspecies are also described as “geographic races” because their distributions correspond to distinct geographic areas.

Five evolutionary lineages have been characterized based on morphometric, molecular, ecological, ethological, and physiological traits [5]. Four of those occur naturally in the Mediterranean Basin: African lineage (A), West and North European lineage (M), Southeast Europe lineage (C), and Near and Middle East lineage (O) [6, 7, 8, 9, 10, 11].

On the basis of morphometric analysis, Ruttner [2] showed that A. m. macedonica existed as a native honey bee in Bulgaria. According to Petrov [12], in Bulgaria, there exists a local honey bee, which is called A. m. rodopica. This local bee inhabits the highest and the central parts of the Rhodope Mountains. Moreover, A. m. rodopica cannot be found in any other regions. Because of this natural isolation, these honey bees are preserved from genetic admixture and influence of other bee races. The Bulgarian honey bee has been adapted to the territory of the country’s landscape and conditions, in which bee families are developing normally and do not show inclinable swarming. Due to the proven biological and productive qualities of the populations and their adaptation to the conditions specific to Bulgaria [19], a gene pool of honey bees needs to be thoroughly studied and preserved.

At present, the subjects of professional and unprofessional beekeeping in Bulgaria are over 800,000 bee colonies. About 13,000 of them are under the control of the National Beekeeping Association. The structure of the sector shows that beekeeping in the country is still extensive and scattered. About 71% of the beekeepers have up to 50 bee families, approximately 24% care for apiaries with 50 to 149 families, and only 5% have apiaries with over 150 families [20].

The genetic diversity of honey bee populations in Europe is threatened by (1) uncontrolled introgression with other subspecies into adapted local populations, (2) stresses from the changing environment due to global climate change, (3) environmental pollution, and (4) the emergence of new pathogens. To respond to these threats and be in full compliance with the requirements of the Convention on Biological Diversity, specific tasks and activities should be organized. This is to preserve the subspecies of honey bees and their ecotypes as a genetic reserve for future selective breeding activities. A lack of efforts and targeted activities for that purpose will have a devastating impact on honey bees and living nature, including the pollination of wild flora and agricultural plants. On a global scale, this would lead to a widespread decline in biodiversity and agriculture yield and would reflect on the planet’s food sources, as 1/3 of human food depends on bee pollination [21].

The Bulgarian local honey bee was studied, in the past, mainly for morphoethological traits [13, 14, 15]. These investigations indicated that in almost all low-elevation parts of the country, bees mainly crossbred with A. m. ligustica. Besides, for more than three decades in the last century with a purpose of selective breeding, A. m. ligustica, A. m. carnica, and A. m. caucasica were introduced to Bulgaria. A comparative analysis established that the Bulgarian honey bee is the most productive in the particular conditions of the environment which it inhabits. With the appearance of new methods, mainly molecular markers, a population structure of local honey bees was studied [16, 17, 18]. The main techniques for elucidation of the population structure of Bulgarian honey bee are divided according to the genetic markers used: 1) allozymes or protein markers; 2) microsatellite analysis; and 3) genetic markers—derived from mitochondrial genome such as 16s rDNA, COI, and ND5 gene regions.

In this paper, we will try to summarize the basic results from genotyping of the Bulgarian Apis mellifera and compare them with other bee populations to show that Bulgarian honey bees possess unique genetic features.

2. Allozyme analysis

The allozyme assay is based on gel electrophoresis results for the allelic variants of the respective enzyme loci. This is one of the first methods for characterization and application of molecular markers to study population genetics of honey bees [22, 23]. Comparative allozyme analysis showed [18] that a small number of loci were characterized for polymorphism in different honey bee populations and there were no fixed significant allelic differences between the subspecies. The importance of allozyme markers for population genetic characterization is the difference in allelic frequencies. An advantage of the allozyme genotyping method is that it is economical and easy to apply, and the results achieved are easy for interpretation.

Allozyme analysis is the most widely used approach to study the population structure of the Bulgarian honey bee. Many authors studied mainly six enzymic systems: malate dehydrogenase (MDH); malic enzyme (ME); esterase (EST); alkaline phosphatase (ALP); phosphoglucomutase (PGM) and hexokinase (HK) for genotyping of Bulgarian honey bees [24, 25, 26, 27]. The honey bee samples from Bulgaria were recognized according to morphometric analysis [2] as A. m. carnica, A. m. caucasica, A. m. ligustica, and A. m. macedonica. The summarized results of reported investigations showed the presence of the following taxonomic markers that could be used in population genetic diagnosis of honey bees subspecies: MDH-1¹²⁵ for subspecies A. m. carnica; EST-3⁸⁸ for A. m. macedonica (of non-Bulgarian origin); and PGM⁸⁰ for A. m. caucasica; HK⁸⁷ for A. m. macedonica (type rodopica). The absence of ALP⁹⁰ was specific to the bee populations of A. m. macedonica. On the other hand, the absence of the MDH-1⁸⁰ allele, a high frequency of MDH-1⁶⁵ and PGM1¹⁴, and a low frequency of EST-3⁹⁴ in the gene pool allowed researchers to make genetic comparisons within the subspecies A. m. macedonica and demonstrated specific characteristics of the Bulgarian honey bee. Moreover, the authors of reports [18, 19] observed a clear subdivision between A. m. macedonica populations from Bulgaria and the surrounding Balkan countries [28]. All honey bee populations from Bulgaria were clustered separately from the A. m. macedonica from other countries.

The local Bulgarian honey bee named by Petrov [22] as “rodopica” could be a different ecotype of A. m. macedonica (Figure 1); however, as mentioned by Ruttner [2], there is no indication of geographical variations within the A. m. macedonica subspecies.

Figure 1.

Principal component analysis (PCA) of the A. m. macedonica and A. m. carnica. PCA shows that A. m. carnica is situated in separate cluster [31].

Isozyme analyses were also used in other studies to compare honey bee population from Bulgaria and other countries of the Balkan Peninsula [17, 29, 30, 31]. The most important results of these investigations could be summarized as follows:

1. The absence of the MDH-1⁸⁰ allele from the gene pool of Bulgarian populations is an important distinctive feature, which can be used to distinguish Bulgarian honey bees from other A. m. macedonica populations of the Balkan Peninsula.

2. It is noteworthy that the EST-3⁸⁸ allele was found in the gene pool only in two of the A. m. macedonica populations and was absent in the gene pool of A. m. carnica.

In this sense, this allele could be used as a taxonomic marker for both sublevel characterization and interpopulation comparisons of the subspecies A. m. macedonica;

1. The presence of the HK87 allele was detected only in the gene pool of the Bulgarian honey bee populations and the allele HK121 was found only in A. m. carnica populations of Serbia with a low occurrence rate.

For this reason, both the HK87 and the HK121 alleles could be used for comparison between the subspecies A. m. macedonica and A. m. carnica and for interpopulation comparisons at the subspecies level.

Thus, application of allozyme markers in Bulgarian honey bees supported the differences between the populations of the two subspecies A. m. carnica and A. m. macedonica at the biochemical level; therefore, they provided a clear opportunity to distinguish between Bulgarian honey bee from the Greek and Macedonian populations of A. m. macedonica at both taxonomic and interpopulation level within the given subspecies.

3. Microsatellite analysis

Microsatellite DNA analysis is among the most preferred and currently applied genotyping approaches for studying genetic variability within populations, population structure, as well as characterization of the phylogenetic relationships. As preferred genetic markers, microsatellites (or simple sequence repeats (SSRs)) are also successfully used for characterization of specific genetic features, comparative population-genetic analyses, such as genetic mapping, as well as DNA barcoding. Bulgarian honey bees were studied using nine microsatellite loci (Table 1). The results showed that all microsatellite loci were polymorphic.

Name of the locus	Am, unified nomenclature for A. mellifera microsatellite markers Am001 to Am552	Accession no. in EBI	Size of the sequenced allele (bp)	Motifs, repeats between primers	Annealing temperature (°C)	MgCl2 (mM)
Ac011	406	AJ509637	127	(CT)19	50	1.5
A024	10	AJ509241	100	(CT)11	55	1.2
A043	25	AJ509256	140	(CT)12	55	1.5
A088	52	AJ509283	150	(CT)10(GGA)7	55	1.2
Ap226	224	AJ509455	231	(CT)8	50	1.5
Ap238	232	AJ509463	251	(AT)6(GT)3(AT)7(GA)8	50	1.5
Ap243	235	AJ509466	260	(TCC)9	50	1.5
Ap249	239	AJ509470	221	(GA)6(GA)8	50	1.5
Ap256	242	AJ509473	162	(GA)12AT(GA)3	50	1.5

Table 1.

Loci, size of allele (bp), indication of the optimal annealing temperature (°C), and concentration of MgCl₂ (mM) in local Bulgarian honey bee populations (A. m. macedonica) [16].

On the basis of different frequencies of the alleles of the given microsatellite loci, the author concluded that these markers can be used as a genetic marker to explore the genetic diversity of local Bulgarian honey bee.

Other investigation has explored the genetic characterization of Bulgarian honey bees as well as their phylogenetic relationships with all European subspecies Apis mellifera using 24 microsatellite loci. The results showed that all microsatellite loci used in the study were polymorphic with a total of 260 allelic variants (Table 2).

LocusA79	Allele	A. m. mac. Bg	A. m. mac. Gr	A. m. mac. Mac	A. m. car.	A. m. lig.	A. m. mel.	A. m. anatol.
1	91	0.000	0.000	0.011	0.190	0.081	0.673	0.000
2	97	0.188	0.167	0.091	0.024	0.000	0.010	0.379
3	99	0.000	0.000	0.034	0.048	0.000	0.010	0.000
4	101	0.104	0.208	0.159	0.048	0.068	0.010	0.229
5	103	0.250	0.500	0.295	0.143	0.324	0.077	0.271
6	105	0.188	0.042	0.091	0.357	0.203	0.019	0.057
7	107	0.104	0.042	0.159	0.143	0.149	0.048	0.050
8	109	0.083	0.042	0.114	0.024	0.054	0.038	0.007
9	111	0.021	0.000	0.023	0.024	0.054	0.058	0.000
10	113	0.000	0.000	0.011	0.000	0.027	0.000	0.007
11	115	0.042	0.000	0.011	0.000	0.014	0.038	0.000
12	117	0.021	0.000	0.000	0.000	0.000	0.010	0.000
13	119	0.000	0.000	0.000	0.000	0.000	0.010	0.000
14	121	0.000	0.000	0.000	0.000	0.027	0.000	0.000

Table 2.

Polymorphism of А79 locus: Alleles and allelic frequencies [30].

Results indicated the following: (1) low levels of genetic differentiation between the Bulgarian A. m. macedonica and the populations of A. m. macedonica originating from Greece and the Republic of Macedonia (0.009–0.017) as well as between A. m. macedonica and A. m. carnica (0.032–0.051); (2) moderate levels of genetic differentiation between population of the subspecies A. m. macedonica (0.071 for the Bulgarian population) and A. m. ligustica (0.059–0.081) and between A. m. macedonica and A. m. anatoliaca (0.087–0.110); and (3) significantly high levels of genetic differentiation between A. m. macedonica and A. m. mellifera (0.290 and 0.307, respectively).

In addition, the results are of exclusive importance to perform genetic analysis of the population of A. m. macedonica. Results help to finalize the comparison of the Bulgarian honey bee with the other populations of A. m. macedonica (Greece, Macedonia, and other Balkan countries), as well as those of the subspecies A. m. carnica, A. m. ligustica, A. m. mellifera, and A. m. anatoliaca. A comparison of honey bee species is important due to given the fact that in the modern course of development of beekeeping in Bulgaria, at certain periods (the 1960s and 1970s of the twentieth century), the introduction of bee queens of the A. m. carnica and A. m. ligustica was carried out. In that, different crossbreeding schemes between these subspecies were applied.

In this regard, Uzunov [31] investigated honey bee populations from Albania, Bulgaria, Greece, and the Republic of Macedonia in South Eastern Europe using 25 microsatellite loci. The PCA plot showed that the Slovenian bees (or known as carnica bees) were clearly distinct from all the others A. m. macedonica bees (Figure 1).

Moreover, the results based on microsatellite analyses indeed showed a certain degree of differentiation between the bees of Bulgaria and A. m. macedonica from other regions. The results strongly pointed out that the Bulgarian honey bee, belonging to the A. m. macedonica based on classical morphometry description [31], possessed some specific genetic features, which differs from the other populations of this subspecies. This significantly justifies Bulgarian honey bees as a different ecotype adapted to the specific conditions in the territory of our country. To verify this, detailed analyses of the honey bee populations in Bulgaria with all neighboring geographical regions will be needed.

The introgression of A. m. carnica alleles into Bulgarian bees has been detected. The influence of the type bees could be explained by past imports of carnica honey bee queens, given its wide commercial spread across Europe. Thus, imports of honey bee queen from other subspecies inevitably influenced genetic diversity of native Bulgarian honey bees.

In summary, results from microsatellite DNA analysis in complex with classic morphometric analysis could be used as evidence in support of an earlier hypothesis. The local honey bee of Bulgaria belongs to subspecies A. m. macedonica, but it represents a different ecotype (ecotype rodopica), adapted to the specific conditions of our country. The local honey bee of Bulgaria differs from the other populations of the same subspecies by a set of population genetic indices.

4. Mitochondrial genetic markers

4.1. PCR-RFLP analysis of 16s rDNA, COI, and ND5 gene regions

Mitochondrial DNA (mtDNA) analysis is also used for characterization of genetic diversity among A. mellifera subspecies and its ecotypes [6, 32, 33, 34, 35, 36]. The mitochondrial genome is determined as extremely suitable for population genetic studies of A. mellifera, as well as to conduct phylogenetic studies, comparisons and analyses of the entire Apis family. Because of maternal heredity manifested by mtDNA in honey bees, all bee workers and drones in the bee colony share similar mitochondrial haplotype [37].

The mtDNA variants, as well as the morphometric characteristics of honey bees, are currently associated with the characterization of their biogeographic zoning. It contributes for the discrimination of previously described evolutionary groups of Apis mellifera [38, 39].

Methods of hybridization in selective breeding may affect the distribution of mtDNA variants of the gene pool of local populations. Such changes affecting population structure (especially introgression) cannot be detected by morphometric analysis [34]. Because all the individuals of the bee family have the same haplotype, a number of authors assumed mtDNA as the suitable genetic marker for population genetic research. In the application of mitochondrial DNA analysis, however, it should be borne in mind that the results can be significantly influenced by imports of bee queens of different origins [40].

In applying mtDNA analysis, it is important to note the possibility that different haplotypes are distinguished and grouped according to classical morphological analyses [39, 41] and the geographical distribution of the subspecies A. mellifera [42, 43].

In this aspect, a PCR-RLFP analysis of three gene regions (16s rDNA, COI, and ND5) was applied in native Bulgarian honey bee A. m. macedonica on the territory of the whole country [22]. The results from this study showed that the size of PCR-amplified mtDNA products for all population of the Bulgarian honey bee is of the order of 964 bp, 1028 bp, and 822 bp for the 16s rDNA, COI, and ND5 gene, respectively.

With regard to the restriction profile of 16s rDNA gene segment, it was found that Sau3A I, Ssp I, Hinc II, and EcoR I, are useful for recognize fragment sites (Table 3).

Sau3A I	Ssp I	Dra I	Hinc II	EcoR I
B	A	A	A	A
548	628	964	598	492
416	336		366	472

Table 3.

Size (in base pairs—bp) determined by the 16s rDNA gene in the analyzed Bulgarian populations of A. m. macedonica [22].

With respect to the COI gene segment, Sau3A I has three restriction sites, but Fok I and Bcl I have two restriction sites (Table 4).

Nco I	Sau3A I	Fok I	Bcl I	Ssp I	Sty I	BstU I	Xho I
B	A	A	A	A	A	B	A
1028	371	476	465	1028	1028	1028	1028
	349	425	326
	280	127	237
	28

Table 4.

Fragment size estimates (in base pairs) of all fragment patterns observed on COI gene region among the Bulgarian populations of A. m. macedonica [22].

With respect to the ND5 gene segment, Dra I and Taq I identified two sites of recognition, Nla III—one and Ssp I—three restriction sites. The Bulgarian populations of A. m. macedonica included in this analysis did not have intra- and interpopulation variability using ND5 polymorphisms (Table 5).

Dra I	Taq I	Nla III	Hinc II	Fok I	Ssp I
A	A	A	A	A	A
429	375	585	822	822	385
285	258	237			206
108	189				124
					107

Table 5.

Fragment size estimates (in base pairs) of all fragment patterns observed on ND5 gene region among the Bulgarian populations of A. m. macedonica [22].

The results from this study [22] were compared with other similar investigations for native honey bee populations from Greece, Cyprus [24] and Crete Island [42]. It can be concluded that the studied Bulgarian bee populations of honey bees of the A. m. macedonica differ from Greek ones. The mtDNA profile of local populations of A. mellifera macedonica from Greece with regard to the COI segment after digestion with endonucleases Sty I and Nco I [36] was not found in mitochondrial profile of Bulgarian populations.

Comparing the results of Ivanova [20] with results from the studies of Bouga et al. [42], the difference between the Bulgarian and Greek populations of A. m. macedonica was reported. In that, a restriction profile of COI after Ssp I digestion, as well as ND5 the gene segment after digestion with Hinc II and Fok I restriction enzymes, was studied.

Thus, the results from PCR-RFLP analysis of three mitochondrial gene regions such as16s rDNA, COI, and ND5 did not reveal any different restriction profile in Bulgarian native honey bee A. m. macedonica populations. The difference existed in comparison with other honey bee populations of A. m. macedonica from other Balkan countries [44].

4.2. Single nucleotide polymorphism (SNP) assay followed by direct sequencing of mitochondrial COI and ND5 gene regions

Single nucleotide polymorphism (SNP) is a DNA sequence variation occurring commonly within a population (e.g. 1%) in which a single nucleotide (A, T, C or G) in the genome differs between members of a biological species or paired chromosomes. SNPs within a coding sequence do not necessarily change the amino acid sequence of the protein that is produced, due to degeneracy of the genetic code.

One of the widely used methods for searching of unknown SNPs is PCR analysis followed by direct sequencing [37]. This technique is based on the amplification of specific parts of the genome in the PCR reaction and sequencing of the obtained sequence amplification product. The great advantage of SNP polymorphism is its universality of the genome between particular species and highly efficient identification of polymorphism within the tested sequence [37]. The high density of SNPs markers in the genome leads to their extensive utilization in genetic and population analyses.

In this aspect, the first investigation examined the phylogenetic relationships and gene flow as a result of migratory beekeeping among honey bee populations from various areas of Albania, Bulgaria, Cyprus, Greece, Italy, Slovenia, and Turkey using sequencing analysis of two mitochondrial regions. In that, the ND5 and the COI gene segments were studied [43]. The results showed that among the 10 races and commercial strains studied, 7 different haplotypes were revealed for COI, 8 for ND5, and 12 for the combined dataset (Table 6). It should also be noted that the most common combined haplotype (haplotype 1) included bees from nearly all the different races examined (including Bulgarian honey bee) except A. m. ligustica. These data showed that the different races could not be discriminated, as it was known that all belonged to the East Mediterranean C lineage [8, 9, 46, 47].

Haplotype 1	A. m. cypria (PYR, KOR, DAL1 and 2), Aegean race near to A. m. adami (LMN1, 2, 3, 4, NIS 1, 2),A. m. adami (CRE4), A .m. macedonica (MAC3, 4, 5, 6 and 7, THR1, 2, 3 and 4, JIJ), commercial strain (genetically improved A. m. cypria (SB), A. m. cecropia (SAR3, MES1, 3, 4), A. m. anatoliaca (BAR), A. m. meda (OSM), A. m. carnica (LUB1 and 2)
Haplotype 2	A. m. ligustica (PER1 and 2, RAV1 and 2)
Haplotype 3	A. m. cecropia (SAR1), A. m. carnica (KEF2)
Haplotype 4	A. m. cecropia (SAR2)
Haplotype 5	A. m. macedonica (MAC1 and 2)
Haplotype 6	Aegean race near to A. m. adami (CHI1 and 2 and RHD1)
Haplotype 7	A. m. adami (CRE2 and 3)
Haplotype 8	A. m. cecropia (MES2), Aegean race near to A. m. adami (RHD2, 3, 4), A. m. ligustica (FOR)
Haplotype 9	A. m. carnica (LIT)
Haplotype 10	A. m. adami (CRE1)
Haplotype 11	A. m. carnica (GOR)
Haplotype 12	A. m. cecropia (LAR1 and 2)

Table 6.

Honey bee populations studied grouped in 12 haplotypes [43].

Abbreviations: Turkey—Osmaniye (OSM); Bartin (BAR); Albania Sarandë (SAR); Slovenia—Litija (LIT), Gorenjska (GOR), Ljubljana (LUB); Bulgaria—Jijevo (JIJ), Plovdiv (PLV); Cyprus—Pyrgas (PYR), Kornos (KOR), Dali (DAL); Italy—Ravenna (RAV), Perugia (PER), Forli (FOR); Greece–Macedonia (MAC) (including Chalkidiki), Thrace (THR), Crete island (CRE), Larissa (Central Greece) (LAR), Messinia (Peloponnese) (MES), Kefalonia (Ionian islands) (KEF), Chios (Aegean Sea Island) (CHI), Rhodes (Aegean Sea Island) (RHD), Nisyros (Aegean Sea Island) (NIS), Limnos (Aegean Sea Island) (LMN).

Another investigation provided additional information about genetic profile of local Bulgarian honey bee [46]. This research presented, for the first time, SNP analysis of a COI mitochondrial DNA (mtDNA) gene segment of the local Rhodope Mountains honey bee (Apis mellifera rodopica). The results showed, unexpectedly, that 14 haplotypes possess a second DNA fragment (named as D1) which is not part of COI gene. Its length is 768 bps from the beginning of 1343 bps of reference sequence (A. m. ligustica; Acc. No. L06178) [48]. This fragment possibly affects translation of cytochrome c oxidase I protein due to disruption of open reading frame at the C-end of protein by multiple stop codons (Figure 2).

Figure 2.

COI pseudogene organization and effects of potential translation of the COI protein on Apis mellifera rodopica [45].

Investigation of homology of mutated D1 fragment in GenBank [48] database did not show any homology, except this with a fragment of A. mellifera COI gene starting from position 264 bps to 346 bps (82 bps length). Such type of gene reorganization is a possible explanation of the duplication of this D1 fragment of A. mellifera originated from Rhodope Mountains, but the mechanism of this duplication event needs further investigation. The authors of this study [46] concluded that this specific characteristic duplication of the COI gene for the Rhodopes honey bee (Apis mellifera rodopica) can be used as a genetic marker to discriminate the local Bulgarian honey bees and to support the related conservation activities.

5. Conclusion

The conducted experiments and obtained results on the basis of various genetic methods for genotyping enabled the detailed characterization of genetic diversity among the studied bee populations of Bulgaria. The application of complex comparative genetic analysis provided genetic evidence on belonging of the Bulgarian honey bee to Apis mellifera macedonica as well as the presence of specific genetic characteristics that define them as a different local ecotype for Bulgaria, named previously as “rodopica” [26].

The genetic diversity studies provided summarized, systematized, and enriched information on genetic polymorphism in Bulgarian honey bee populations. Scientific efforts characterized the level of genetic differentiation relative to other populations of the Apis mellifera macedonica on the Balkan Peninsula and to other subspecies of Apis mellifera, in the territory of Europe.

The summarized population genetic studies from genotyping of Bulgarian honey bee also provided significant information on the availability of reliable genetic markers, applicable to the formation of a system of conservation activities for the national genetic resources of Bulgarian honey bee population.

Acknowledgments

The authors express gratitude to the National Science Fund of the Bulgarian Ministry of Education and Science, Sofia, for the financial support (grant number 06/10 17.12.2016).

Conflict of interest

No potential conflict of interest was reported by the authors.