*Article* **Further Evidence of Population Admixture in the Serbian Honey Bee Population**

**Marija Tanaskovi´c 1,\* , Pavle Eri´c <sup>1</sup> , Aleksandra Patenkovi´c <sup>1</sup> , Katarina Eri´c <sup>1</sup> , Milica Mihajlovi´c <sup>2</sup> , Vanja Tanasi´c <sup>2</sup> , Szilvia Kusza <sup>3</sup> , Andrzej Oleksa <sup>4</sup> , Ljubiša Stanisavljevi´c <sup>5</sup> and Slobodan Davidovi´c <sup>1</sup>**


**Simple Summary:** The western honey bee is one of the most ecologically and economically important pollinator species. Due to human interference, it faces serious challenges, not only in number decline and habitat loss, but also in natural subspecies diversity and distribution. The conservation of genetic diversity and perseverance of locally adapted populations and subspecies becomes a crucial task in the face of rapid environmental changes. In order to further assess present genetic variability in Serbian honey bee populations, we analyzed 14 microsatellite loci and then compared nine of them with previously published data. Our results suggest that Serbia now harbors a distinct, relatively homogenous honey bee population, although some local differences are still preserved.

**Abstract:** Socioeconomic interests and beekeeper preferences have often taken precedence over the conservation of locally native honey bee subspecies, leading to the predominance of admixture populations in human-dominated areas. To assess the genetic diversity of contemporary managed Serbian honey bee colonies, we used 14 microsatellite loci and analyzed 237 worker bees from 46 apiaries in eight localities of northern and southern Serbia. Furthermore, we compared data for nine microsatellite loci with 338 individuals from Italy, Hungary, Poland, and Spain. The standard parameters of genetic diversity in Serbian honey bee populations were in line with other analyses, although somewhat smaller. STRUCTURE analysis showed the existence of two equally distributed genetic clusters and Analysis of molecular variances could not confirm the presence of a geographically discrete population but showed local differences. Discriminant analysis of principal components showed overlapping of worker bees from different parts of Serbia. Clear genetic differentiation can be observed when comparing all populations between geographical regions and their corresponding subspecies. The absence of the *A. m. macedonica* subspecies from its historical distribution range in southern Serbia as well as the lack of distinctive geographical groups suggest that selective breeding, queen import, and migratory beekeeping practices strongly influence the genetic structure and diversity of honey bees, leading to the genetic uniformization and creation of the admixture population.

**Keywords:** honey bee; microsatellite; population genetics; genetic diversity

**Citation:** Tanaskovi´c, M.; Eri´c, P.; Patenkovi´c, A.; Eri´c, K.; Mihajlovi´c, M.; Tanasi´c, V.; Kusza, S.; Oleksa, A.; Stanisavljevi´c, L.; Davidovi´c, S. Further Evidence of Population Admixture in the Serbian Honey Bee Population. *Insects* **2022**, *13*, 180. https://doi.org/10.3390/insects 13020180

Academic Editor: Silvio Erler

Received: 30 December 2021 Accepted: 6 February 2022 Published: 9 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

The western honey bee (*Apis mellifera* Linnaeus, 1758) is one of the species that is a subject of constant human interference. Although its domestication began more than 10,000 years ago, this species has never become truly domesticated despite all efforts, mainly due to its complex mating behavior [1–3]. The species' native range of distribution in western Asia, Africa, and Europe was expanded to all other continents, except Antarctica, to regions marked by highly distinct ecological and climate conditions. In addition to largescale, transcontinental movements and expansion of the natural range, beekeeping practices and preferences for perceptively more suitable subspecies significantly changed distribution and variability in the historical range of distribution. In recent years, a substantial body of evidence has confirmed that, deliberately or not, humans shape the current diversity of honey bees worldwide [4].

Since the classical Ruttner categorization of the *A. mellifera* subspecies [5], there has been an ongoing debate about its taxonomy, number of subspecies, distribution range, and origin due to the specificity of population structure, features of biology, and resolutions of honey bee subspecies discrimination methods. Based on the results of genomic analysis [6–9], it was proposed that this species originated in northern Africa or the Middle East [10,11], but the most recent work [12] showed adaptive radiation of subspecies from Asia. From there, it colonized its native geographic range, and followed by multiple colonization waves and glacial events, it diverged into 33 existing subspecies [13] and divided into five evolutionary lineages (A, C, M, O, and Y). However, it still remains unclear how accurate this number of subspecies is, since many subspecies, due to high phenotypic plasticity, have ecotypes previously defined as subspecies [4,14–16]. The main differences between subspecies, often referred to as geographic races, are most likely the result of both local adaptations to distinct environments and geographical isolation. However, reproductive isolation does not often exist, and subspecies readily interbreed when they come into contact, although partial reproductive isolation is observed [17]. When subspecies come into contact, naturally or by human interference, it is inevitable that an admixture population will be established and that introgression of foreign genetic variants can be detected in native populations. This situation is especially prominent in Europe, which is the natural area of distribution of mitochondrial lineages A, M, and C. Beekeepers prefer subspecies *A. m. carnica,* and *A. m. ligustica*, both classified as C lineage, and their Buckfast hybrid has been imported to almost all parts of the continent [18–21]. It resulted in a considerable degree of genetic admixture among subspecies, even though in some areas the specific genetic footprint of autochthonous subspecies is still preserved [22–26]. For example, *A. m. mellifera* populations (belonging to the M lineage) have been hybridized in varying degrees in most of their native areas and in some parts, such as Germany, have been almost replaced by *A. m. carnica* because local beekeepers preferred this subspecies [27]. The loss of native subspecies and specific genetic diversity they harbor made conservationists and lawmakers in several countries establish protected areas for their native subspecies (Denmark and Great Britain) [19] as well as make laws that prohibit the breeding of nonnative bees (for example, Serbia and Croatia) [28].

Due to the spread of Varroa mites in the last decades of the twentieth century, it was believed that wild honey bee colonies became extinct in Europe, but now new evidence is emerging that some areas have thriving feral or even possibly wild honey bee colonies [29–32]. Therefore, it is not surprising that most research on the genetics of honey bees is conducted in managed colonies. Large-scale genetic comparisons show that genetic structure in any given area is heavily dependent on several factors, the most prominent of all being the level of importation of foreign queens and the presence and the type of breeding strategy implemented by beekeepers and their organizations. Extensive management by beekeepers promotes population admixture [33], which is expected when humans facilitate the movement and interbreeding of previously structured populations [34]. The general conclusion is that many lines used for contemporary beekeeping in Europe consist of a mixture of different source populations [35]. The bees from areas with frequent queen

importation show a high level of admixture and are hardly assignable to distinct subspecies anymore, but for those in areas where breeding lines were selected and maintained at their geographical origin, genetic identity was preserved and they resemble their native source populations [28,31,36–39].

It is noted that in areas that are natural contact zones of different *A. mellifera* subspecies, natural hybridization occurs [40,41] and the hybridization is inevitable in regions where human interference due to beekeeper preferences is high, which occurred in the C lineage native area of distribution [24–26,42,43]. Serbia, located in the center of the Balkan Peninsula, is geographically in the middle of the distribution range of the C lineage. Previous analysis showed that of four C lineage subspecies, Serbia harbors two (*A. m. macedonica* and *A. m. carnica*) which are clinaly distributed from the northwest (*carnica*) to the southeast (*macedonica*) with a hybrid zone between them [23,26,44,45]. Furthermore, nine described mtDNA haplotypes for *tRNAleu-cox2* of which two are novel [46] and three distinct ecotypes belonging to specific geographical regions [47] reflect significant genetic variability of *A. mellifera* in this region. In the past 30 years, the variability of Serbian honey bees has been extensively described on morphological [26,48], etiological [49,50], and genetic levels [26,44,45,47,51–53]. Microsatellite analysis of Serbian honey bees from the first decade of the 21st century showed that although substantial admixture between *A. m. carnica* and *A. m. macedonica* in the central part of Serbia can be detected, populations from the northwest and southeast retain a distinctive subspecies genetic footprint [26]. The results of microsatellite and tRNAleu-cox2 mtDNA variability [23] confirmed this clinal distribution of subspecies and their hybrids, with worker bees from northern Serbia forming a distinct genetic cluster characterized as *carnica-2* ecotype, and those from southern Serbia forming a different distinct genetic cluster characterized as the *macedonica*-1 ecotype.

However, beekeeping practices in Serbia have changed dramatically in the last decade. First, the number of managed beehives has doubled since 2009 [54], partly due to the government's financial support, and now Serbia has the largest number of beehives per capita worldwide (one beehive per six inhabitants). According to our field data, the traditional way of beekeeping is lost and the number of stationary apiaries is dwindling. The production of beekeepers who prefer *A. m. carnica* queens intensified, and the number of queen breeding institutions focused on its desired traits is growing. Serbia is one of the countries that has embedded in its legislation the preservation of autochthonous species, subspecies, and races, and as per the Law on Animal Breeding from 2009, breeding and keeping of only *A. m. carnica* are allowed in its territory [55]. Recent work based on the variability of the tRNAleu-cox2 mtDNA region [46] suggests that the composition and distribution of honey bee populations in Serbia has changed over the past decade, invoking the need for further examination of genetic variability on various levels.

Biparental inherited microsatellite loci proved to be an excellent genetic marker for inferring overall population genetic variability, deciphering the distribution of different *A. mellifera* subspecies [19,56], detecting population admixture [7,26], and determining the presence of distinct locally adapted populations [20,57]. Large-scale analysis of microsatellite loci enables a better understanding of large and fine-scale population differentiation. To better understand the genetic variability of contemporary Serbian honey bee populations, we analyzed 14 microsatellite loci in 237 worker bees from the northern and southern parts of Serbia. In addition, we compared our data for nine microsatellite loci with previously published data [58] for 338 individuals from Hungary, Poland, Spain, and Italy belonging to *A. m. mellifera, A. m. carnica*, *A. m. iberiensis*, *A. m. ligustica*, and the Buckfast hybrid.

#### **2. Materials and Methods**

#### *2.1. Sampling*

A total of 237 worker bees were collected from 46 stationary apiaries during late August and early September in 2020 and stored in 95% ethanol at −20 ◦C for further analysis. Eight localities from southern and northern parts of Serbia were chosen, four in the south (Leskovac (L), Vlasina (V), Stara Planina (SP), and Tromeda (T)) and four in the ¯ north (Subotica (S), Vršac (Vr), Deliblatska pesšˇccara (DS), and Fruška Gora (FG) (Figure 1). Approximately five worker bees from the apiary were chosen for genetic analysis, each representing one beehive. The detailed specification of sampling sites can be found in [46]. Furthermore, the DNA of 338 individuals from Hungary, Poland, Spain, and Italy belonging to *A. m. mellifera, A. m. carnica*, *A. m. iberiensis*, *A. m. ligustica,* and Buckfast hybrid from Péntek-Zakar 2015 were obtained for comparison.

**Figure 1.** Sampling localities (from [46]).
