Honey Bee Diversity Is Swayed by Migratory Beekeeping and Trade Despite Conservation Practices: Genetic Evidence for the Impact of Anthropogenic Factors on Population Structure

Kükrer, Mert; Kence, Meral; Kence, Aykut

doi:10.3389/fevo.2021.556816

ORIGINAL RESEARCH article

Front. Ecol. Evol., 15 April 2021

Sec. Biogeography and Macroecology

Volume 9 - 2021 | https://doi.org/10.3389/fevo.2021.556816

This article is part of the Research Topic Adaptation of Invasive Species to Islands and the Puerto Rican Honey Bee View all 14 articles

Honey Bee Diversity Is Swayed by Migratory Beekeeping and Trade Despite Conservation Practices: Genetic Evidence for the Impact of Anthropogenic Factors on Population Structure

$\r\nMert Kükrer,*$ Mert Kükrer^1,2*

Meral Kence¹

Aykut Kence^1†

¹Biology Department, Middle East Technical University, Ankara, Turkey
²Molecular Biology and Genetics Department, Kilis 7 Aralik University, Kilis, Turkey

The intense admixture of honey bee (Apis mellifera L.) populations at a global scale is mostly attributed to the widespread migratory beekeeping practices and replacement of queens and colonies with non-native races or hybrids of different subspecies. These practices are also common in Anatolia and Thrace, but their influence on the genetic make-up of the five native subspecies of honey bees has not been explored. Here, we present an analysis of 30 microsatellite markers from honey bees from six different regions in Anatolia and Thrace (N = 250 samples), with the aim of comparing the impact of: (1) migratory beekeeping, (2) queen and colony trade, and (3) conservation efforts on the genetic structure of native populations. Populations exposed to migratory beekeeping showed less allegiance than stationary ones. We found genetic evidence for migratory colonies, acting as a hybrid zone mobile in space and time, becoming vectors of otherwise local gene combinations. The effect of honey bee trade leaves very high introgression levels in native honey bees. Despite their narrow geographic range, introgression occurs mainly with the highly commercial Caucasian bees. We also measured the direction and magnitude of gene flow associated with bee trade. A comparison between regions that are open and those closed to migratory beekeeping allowed the evaluation of conservation sites as centers with limited gene flow and demonstrated the importance of establishing such isolated regions. Despite evidence of gene flow, our findings confirm high levels of geographically structured genetic diversity in four subspecies of honey bees in Turkey and emphasize the need to develop policies to maintain this diversity. Our overall results are of interest to the wider scientific community studying anthropogenic effects on the population diversity of honey bees and other insects. Our findings on the effects of migratory beekeeping, replacement of queens and colonies have implications for the conservation of honey bees, other pollinators, and invertebrates, in general, and are informative for policy-makers and other stakeholders in Europe and beyond.

Introduction

The Western honey bee, Apis mellifera L., plays an important role, together with other pollinators, in the pollination of wild and cultivated plants. Likewise, honey bees have significant economic importance due to their production of honey and other products (Morse, 1991; Breeze et al., 2011). In addition to their ecological and economic importance, honey bees serve as model organisms for the study of fundamental questions on sociality and on cognition (Weinstock et al., 2006).

The natural distribution of A. mellifera includes Central and Southwest Asia, Europe, and Africa but the species was also introduced to East and Southeast Asia, Australia, and the Americas, mainly for its economic benefits (Ruttner, 1988). Morphological and molecular studies point to four major lineages of numerous—more than 20—subspecies (Ruttner, 1988; Whitfield et al., 2006). The four widely recognized lineages are A (Africa), M (western and northern Europe), O (Near East and Central Asia), and C (Eastern Europe) lineages.

In the past decade, various molecular-based studies have garnered support for the hypothesis that A. mellifera originated in the African tropics or subtropics and colonized its present European range by two main routes: through the strait of Gibraltar and through the Suez and then Bosporus regions, with a subsequent contact between the highly divergent M and C lineages in the region surrounding the Alps (Whitfield et al., 2006; Han et al., 2012; Harpur et al., 2014; Wallberg et al., 2014; Cridland et al., 2017).

Honey bees and wild pollinators are experiencing a worldwide decline due to factors closely related to human activities. Levels of decline vary and are related to species and geographic region. Some of the anthropogenic factors implicated in the decline are the destruction and fragmentation of natural habitats, toxicity caused by pollution and pesticides—such as the widely used neonicotinoids—and diseases. The latter is being facilitated by the spread of invasive species (Meffe, 1998; Brown and Paxton, 2009; Van Engelsdorp and Meixner, 2010; Blacquiere et al., 2012). Wild populations of honey bee species including feral populations in the genus Apis also have been negatively affected, namely, Apis cerana, Apis florea, Apis dorsata, and other native bees of Asia (Oldroyd, 2007; Dietemann et al., 2009; Van Engelsdorp et al., 2009; Genersch, 2010; Evans and Schwarz, 2011).

In addition, native honey bees are experiencing local losses, extinction, and/or genetic swamping as a result of genetic admixture due to bee trade, including the replacement of local bees with non-native strains and the beekeeping practice of moving colonies between geographic areas (De la Rúa et al., 2009).

The above genetic and environmental factors, and their interaction, have a cumulative adverse effect on honey bees and likely contribute to continuous or discrete events of sudden colony losses. This Colony Collapse Disorder (CCD) or Colony Depopulation Syndrome (CDS) (Van Engelsdorp et al., 2009; Neumann and Carreck, 2010) as it is referred to, is characterized by rapid depletion of worker bees while the queen continues egg laying and is accompanied by a lack of dead bees in and around the hive.

Honey bees may be able to adapt to this new challenge by relying on the adaptations and genetic diversity they accumulated over their evolutionary history. Honey bee subspecies perform differently in different environments and some locally adapted populations may display greater resistance to anthropogenic threats (Büchler et al., 2014). Hence, research on honey bee genetic diversity at the individual, colony, population, ecotype, and subspecies level is of great importance for safeguarding the species, the ecosystem, and the economic services they provide.

Recent research conducted on European honey bee population structure has shown that past distribution patterns have been disturbed (Dall’Olio et al., 2007; Bouga et al., 2011; Cánovas et al., 2011). In Africa, introgression of non-native DNA was detected in wild populations of Sudan (El-Niweiri and Moritz, 2010). The causes for these disturbances have been attributed mainly to queen and colony trade, replacement of native honey bees with non-natives as well as migratory beekeeping. However, there are very few studies on the direct genetic consequences of human practices on honey bee diversity.

Turkey has five subspecies of honey bees within its borders and beekeepers practice a variety of strategies, thus it provides an ideal environment to test the impact of anthropogenic factors on one of the most important pollinators of crops and wild plants.

Beekeeping in the region of Anatolia is a practice dating back to 6600 BC when the Hittite civilization presided over this region (Akkaya and Alkan, 2007). Beekeeping has been persevered and continues to be intensively practiced in Turkey where there are more than eight million hives distributed throughout the country. This is the third-highest number of hives in a single country. It is three times higher than the number of hives in the United States and reaches half of the EU countries total (European Parliament, 2017; USDA NASS, 2019).

As a reflection of the long association of the genus Apis with the region of Anatolia, one-fourth to one-fifth of the recognized subspecies of A. mellifera, namely, A. m. meda, A. m. syriaca, A. m. caucasica, and A. m. anatoliaca from the O-lineage and an ecotype from the C subspecies group occur in Turkey (Kandemir et al., 2005). In addition, A-lineage genetic material was also detected in native bees from the Levantine coast of Turkey (Kandemir et al., 2006) bringing together genetic elements from three continents, Africa, Europe, and Asia. The major subspecies found in and around Anatolia are shown in Figure 1A.

FIGURE 1

Figure 1. Geographic distribution of (A) major honey bee (A. mellifera) subspecies in and around Anatolia (B) sampling sites and sample sizes.

Together, Anatolia and Thrace harbor a vast diversity of honey bees belonging to three different lineages. In this region, they meet, exchange genes, and adapt to local conditions determined by local diverse climatic, topographical, and floristic variations (Bouga et al., 2011). The refugial status of Anatolia during the ice ages contributed to the enhanced levels of the present floral and faunal biodiversity (Hewitt, 1999). Studies of Turkish honey bee populations (Bodur et al., 2007; Kence et al., 2009) demonstrated high genetic structuring and confirmed the presence of divergent populations pointing to different subspecies. These researchers pointed to the rich diversity of honey bee populations in Anatolia and Thrace, and highlighted that they are under threat and that there is an urgent need to take steps for their conservation.

However, despite the above research, arguments prevail in the beekeeping environment that locally adapted honey bee ecotypes have been irrevocably lost, due to gene flow, and thus steps to safeguard locally adapted honey bee variants do not have merit. This argument is further strengthened since queen bee trade is not currently subject to any restrictions or regulations in Turkey. There are few pioneering measures of conservation within the natural distribution range of the subspecies, likely not enough to guarantee the preservation of genetic structure in the next decades.

The aim of this study is to investigate the impact of anthropogenic factors and conservation efforts on the current pattern of genetic diversity of honey bee populations in Turkey. Should genetic structure be identified, this could inform policies. Conservation measures could avoid extinction of native races, ecotypes, and diversity to be found in these populations. Genetic similarity of donor and recipient populations may be considered in recommending migration routes for migratory beekeepers and bee sales.

The research herein tested hypotheses regarding the occurrence of recent admixture in Turkish honey bee populations across the subspecies of A. m. syriaca, A. m. caucasica, A. m. anatoliaca, and the C-lineage ecotype in Thrace using 30 microsatellite markers. In addition, we: (i) investigated the robustness of genetic diversity of honey bees in geographic areas where migratory beekeeping is restricted for varied reasons; (ii) compared patterns of genetic diversity of honey bees between migratory and stationary colonies; and (iii) determined the degree, origin, and direction of introgression in the Turkish honey bee populations to assess the consequences of unregulated queen and colony trade.

Materials and Methods

Sampling

We sampled a single honey bee each from 250 colonies located in 18 Turkish provinces during the period of March 2010 through August 2012. Of the 250 honey bees sampled, 174 were from stationary and 76 from migratory colonies. Beekeepers who participated in this study declared that they used honey bees from stocks native to their area and that they had not purchased non-native queens or colonies in the last 10 years. Honey bee samples were stored at −80°C prior to genetic analysis.

We grouped samples from provinces with small sample sizes with nearby provinces to form 10 major localities: (1) Kırklareli; (2) Edirne + (Edirne and Tekirdağ); (3) Muğla; (4) Eskişehir + (Eskişehir, Kütahya, and Bilecik); (5) Düzce + (Düzce, Zonguldak, and Bolu); (6) Ankara; (7) Hatay; (8) Bitlis + (Bitlis, Elazığ, Erzurum, and Ordu); (9) Ardahan; and (10) Artvin. The localities sampled correspond to the natural distribution range of the five subspecies that occur in Turkey: A. m. syriaca in Hatay, A. m. caucasica in Ardahan and Artvin, A. m. anatoliaca in Düzce, Eskişehir +, Muğla, and Ankara from the O lineage as well as the ecotype from the C subspecies group that occurs in Kırklareli and Edirne + and A. m. meda. Geographic locations were considered based on geographical proximity, and similarities in climate, topography, and floral profiles as well as preliminary data from previous studies. Sampling sites and sample sizes are listed in Figure 1B.

Genotyping

We isolated DNA from bee heads using the DNeasy Blood and Tissue Kit (Qiagen, Ankara) following the manufacturer’s instructions, with slight modifications for insect samples. For polymerase chain reactions (PCR), we grouped a set of 30 microsatellite loci into four clusters for two 7-plex (set 1: AP218, A113, AB024, AP249, A088, AP001, AP043; set 2: AP049, AP238, AC006, AP243, AP288, HBC1602, A107) and two 8-plex (set 3: A079, AC306, AP226, A007, HBC1601, AP068, A014, AP223; set 4: AP019, AB124, A043, A076, AP273, AP289, HBC1605, A028) (Estoup et al., 1995; Solignac et al., 2003; Bodur et al., 2007; Shaibi et al., 2008; Tunca, 2009). The program, Multiplex Manager 1.2 (Holleley and Geerts, 2009) was used to determine the multiplex groups. Information on primer pairs, fluorescent dyes, and PCR conditions are provided in Supplementary Tables 1, 2.

The microsatellite allele sizes were determined by capillary electrophoresis with the ABI 3730XL sequencing machine (Applied Biosystems, Foster City CA). Locus A076 did not consistently amplify across samples; thus, it was excluded from the data and the downstream analysis.

Population Structure

We calculated pairwise F_ST values using Arlequin 3.5 (Excoffier and Lischer, 2010); the Mantel test with 10,000 permutations was used to test for isolation by distance. Pairwise population distances were calculated (Reynolds et al., 1983) using Populations 1.2.32 software (Langella, 2011) and visualized with the online tool Interactive Tree of Life v4 (Letunic and Bork, 2019). We used PAST4 and PCAgen software to plot relationships of populations on a two-dimensional space using a correlation matrix between groups (Goudet, 1999; Hammer et al., 2001).

Population structure was estimated by Structure 2.3.3 (Pritchard et al., 2000), K-values of distinct populations were analyzed by Structure Harvester software (Earl and von Holdt, 2012), and we used the Clumpp software (Jakobsson and Rosenberg, 2007) to permute the membership coefficients of individuals determined by Structure 2.3.3 and Distruct software (Rosenberg, 2004) to visualize the results obtained by Clumpp.

Other population genetic parameters and diversity indicators were calculated and include the frequency of null alleles, allelic richness and diversities, inbreeding and prevalence of close relatives, number of effective alleles, levels of heterozygosity, deviations from Hardy–Weinberg, linkage disequilibrium, bottlenecks, effective population sizes, and microsatellite information index (Supplementary Tables 3–10).

Statistical Analyses

To test the hypotheses regarding beekeeping practices, conservation sites, and queen/colony trade, we used membership coefficients. We first applied the arcsine square-root (angular) transformation to the coefficients since the data were composed of proportions and not normally distributed (Rohlf and Sokal, 1995). Then we performed Shapiro, Mann–Whitney U, Kruskal–Wallis, Dunn’s, F, ANOVA, Tukey’s, and t-tests wherever necessary and applicable to compare mean membership coefficients and estimated Cohen’s d to determine effect sizes. The above tests were carried out with R statistical software using packages pwr, effsize, dunn.test, and dabestr (R Core Team, 2013; Torchiano, 2016; Dinno, 2017; Champely et al., 2018; Ho et al., 2019). The associated code is provided as Supplementary Material (R code).

Estimation plots were used to visualize untransformed data for membership coefficients and the impact of experimental factors. This is a less conventional method than bar or boxplots and the reporting of significance tests but more convenient and powerful to summarize all the data in an unbiased manner by displaying all measurements and effect sizes as well as the precision of estimates and distribution of mean differences (Ho et al., 2019).

Beekeeping Practice: Migratory vs. Stationary

To test the hypothesis of whether beekeeping practice affects population structure and subspecies identity, we compared membership coefficients of migratory and stationary colonies in Ankara, Muğla, and Hatay separately, combined, and for the total data set. We propose that if migratory colonies acted as a potential vector of foreign honey bee alleles, then samples would have much lower probabilities of being assigned to the clusters of origin.

We used all samples (N = 250) to quantify the differences in membership coefficients for migratory and stationary colonies. For the remaining analysis, we used a subset of the samples from stationary colonies (n = 174) since this can give a better perspective of the population structure.

Isolated Regions as Conservation Sites

If isolated regions preserve genetic diversity by preventing gene flow, we predict higher membership coefficients for samples that originate from isolated regions compared with those from regions exposed to migratory beekeeping.

Kırklareli is officially declared as an isolated region. This is due to local beekeepers’ long-standing negative attitude and resistance to migratory beekeepers. As a result, they have not accessed this region for many years. This region is home to a C-lineage honey bee ecotype, carefully maintained by local beekeepers. Ardahan is legally declared as a conservation and breeding area for A. m. caucasica, therefore migratory beekeepers cannot enter the province, and queen import from other subspecies is forbidden. Parts of Artvin province are also officially declared as isolated regions for the conservation of A. m. caucasica as a pure race. The province, in general, is rarely visited by migratory beekeepers because of the difficulties in transportation in the rough terrain. Moreover, local beekeepers often engage in commercial queen sales so, they only use native-bred queens. We compared the above three provinces with restricted inflow of migratory beekeepers with the other six regions (Edirne +, Muğla, Düzce +, Eskişehir +, Ankara, and Hatay) where migratory beekeeping and bee trade are freely exercised.

Effect of Queen and Colony Trade

Using all samples, we compared membership coefficients in non-native clusters between each other to determine which groups contributed most to other populations’ gene pools.

Ardahan and Artvin provinces host the A. m. caucasica subspecies, which is widely used for commercial purposes. A. m. caucasica queens and their hybrids are sold throughout Turkey. However, these provinces are limited to a very narrow range in the Northeast of the country and are declared isolated regions. Therefore, a high introgression of caucasica alleles from these regions would mostly, if not completely, be due to the replacement of queens and colonies.

We also tested for the presence of other genetic patterns within the Turkish honey bee population to understand the magnitude and direction of gene flow within and across the sampled localities.

Results

Population Genetic Structure

We calculated F_ST values by using the frequencies obtained in the study and the null allele corrected frequencies. We calculated an F_ST of 0.065 for all samples and an F_ST of 0.067 after correction for the stationary colonies (n = 174). The F_ST values for migratory colonies were 0.011 and 0.015, respectively, and for all the 250 samples, they were 0.046 and 0.047. We plotted stationary colonies on 2D space by carrying out the principal component analysis (PCA) (Figure 2A). The x and y axes explained, respectively, 41.8 and 32.1% of the variance within the samples. For stationary colonies, the phylogenetic tree constructed using pairwise population distances resolved four distinct branches (Figure 2B). Using the Structure Harvester clustering program, we determined that K = 2 and K = 4 gave similar outcomes with the latter being more likely as this mirrors the number of subspecies present in the regions sampled.

FIGURE 2

Figure 2. (A) Principal component analysis (PCA) of stationary colonies with 66% concentration ellipses shown. Component 1 and Component 2 explain 41.8 and 32.1% of the variance within the samples, respectively. The first axis differentiated samples in Thrace indicating strong divergence between those and others whereas the second axis differentiated subspecies throughout Anatolia. Each dot represents an individual (orange: C-lineage subspecies in Thrace; O-lineage subspecies: yellow: anatoliaca, blue: caucasica, violet: syriaca). (B) UPGMA tree of honey bee populations based on Reynolds, Weir, and Cockerham’s genetic distances. Tree resolves four distinct branches corresponding to four subspecies. Thracian populations constitute the extreme end of the unrooted tree. The other end is divided into three almost equidistant branches of Caucasian, Levantine, and Anatolian populations.

We calculated membership coefficients of individuals to the observed clusters in K = 4 and found no population structure for migratory colonies (Figure 3A) in contrast to samples from stationary colonies and the entire data set (Figures 3B,C).

FIGURE 3

Figure 3. Estimated population structure and clustering of honeybees in Anatolia and Thrace for (A) migratory colonies, (B) stationary colonies, and (C) the whole sample. Structure analysis is based on microsatellite data and suppose either K = 2 (orange: C-lineage, blue: O-lineage) or K = 4 (orange: Thracian, yellow: Anatolian, blue: Caucasian, violet: Levantine) hypothetical populations. Each individual is represented by a vertical bar and colored according to membership coefficients belonging to each cluster. K = 2 is found to be slightly likelier than K = 4. For stationary colonies at K = 2, the transition from C-lineage to O-lineage is gradual. No population structure is observed in migratory colonies in contrast to stationary ones and the overall data where four different subspecies are evident. Note the higher admixture levels in the overall data in comparison to stationary colonies.

Effects of Beekeeping Practices and Conservation

Results from the Mantel test showed a significant correlation with geographic distance between populations (r = 0.60, p < 0.01) for stationary but not migratory colonies. Distance matrices and test results are provided in Supplementary Tables 11, 12. A significant difference was detected in a comparison of membership coefficients of individuals from stationary and migratory colonies (Figure 4A). Stationary colonies from Muğla and Hatay had a higher likelihood to be assigned to their own clusters than migratory colonies sampled from these provinces (p < 0.05 and p < 0.001, respectively, Mann–Whitney U and t-tests). The same pattern was observed when the combined data from the three provinces (p < 0.01), or all the migratory and stationary colonies (p < 0.001) were considered. However, the situation was reversed in Ankara (p < 0.05). In all the comparisons but one, 95% CI of the mean differences between the membership coefficients of migratory and stationary colonies lie below the zero-line (Figure 4B). The mean values, effect size, and significance level of the differences are summarized in Table 1.

FIGURE 4

Figure 4. Comparison between stationary (Sta_) and migratory (Mig_) colonies in Ankara (p < 0.05), Muğla (p < 0.05), and Hatay (p < 0.001), as well as these three provinces combined (p < 0.01) and the whole data set (p < 0.001). (A) Boxplot display of arcsine square-root transformed membership coefficients used in significance testing of comparisons. (B) Scatter plot with estimations of mean differences based on raw individual membership coefficients (yellow: Ankara and Muğla belonging to the Anatolian cluster, violet: Levantine cluster, coral: for a combination of three provinces, firebrick: whole data). Stationary colonies are annotated as < Group name > 0 and migratory colonies as < Group name > 1. Bars right to the data points refer to the 25 and 75% quartiles and the gap between them is the median value for the sample. The zero-line below corresponds to the mean membership coefficients of stationary colonies in each pairwise comparison. The Euclidean distances from those means for the migratory colonies are shown as dots with a 95% confidence interval bar around. Also, distributions of the estimation statistics are included to comprehensively compare the strength of the drift for different populations and subsets of the data. Stationary colonies exhibit higher mean membership coefficients than migratory ones except for Ankara where the vice versa is true.

TABLE 1

Table 1. Genetic impact of beekeeping and conservation practices on (arcsine square-root transformed) membership coefficients to native clusters (*p < 0.05, **p < 0.01, and ***p < 0.001).

The comparison of isolated regions with those open to migratory beekeeping (Table 1 and Figure 5A) showed that stationary colonies within isolated regions have significantly higher fidelity to their original clusters (p < 0.001, Mann–Whitney U and t-tests). This can also be seen in the estimation plot (Figure 5C) where the mean membership coefficients of samples that are from regions open to migratory beekeeping are lower and fall beyond the 95% confidence interval of the estimated mean of the difference between the two groups. In addition, despite the lack of conservation efforts, samples from Hatay and Düzce + showed membership coefficients comparable with those of Kırklareli, Ardahan, and Artvin (Figure 5B).

FIGURE 5

Figure 5. Samples within isolated regions assigned to their clusters with higher probabilities in contrast to samples from regions open to migratory beekeeping. (A) The first boxplot displays the arcsine square-root transformed membership coefficients for nine populations, whereas the second one presents a comparison of samples within isolated regions and those are not (p < 0.001). (B) Scatter plot based on a comparison of raw individual membership coefficients to their native clusters for nine populations against Thracian samples. Bars right to the data points refer to the 25 and 75% quartiles and the gap between them is the median value for the sample. Note that despite lacking a conservation status, samples from Hatay and Düzce + have membership coefficients comparable to those of Kırklareli, Ardahan, and Artvin which are isolated regions. (C) Scatter plot contrasting individual raw membership coefficients with an estimation of the mean difference between isolated regions and those are not (orange: Thracian, yellow: Anatolian, blue: Caucasian, violet: Levantine clusters; orchid and “1”: isolated regions, green and “0”: regions open to migratory beekeeping). The zero-line corresponds to the mean membership coefficient of colonies in regions open to migratory beekeeping. The Euclidean distance of the colonies in isolated regions from that mean is shown as a dot with a 95% confidence interval bar around. The distribution of the estimation statistic is included to account for the precision.

Impact of Queen and Colony Trade

If an individual is assigned with high probability to its own cluster, i.e., 90% probability, there remains a 10% chance that it can be assigned to other clusters. Given four clusters, we investigated whether these mis-assignment probabilities were enriched for any particular cluster. The mean transformed values of cluster mis-assignments among individuals of other populations were as follows: Thracian 0.16, Anatolian 0.25, Caucasian 0.26, and Levantine 0.20 (Figure 6A).

FIGURE 6

Figure 6. Mis-assignment of individuals to caucasica and anatoliaca clusters were significantly more frequent than the others (p < 0.001 for both subspecies against C-lineage Thracian bees and p < 0.05 against syriaca group). (A) Boxplot display of arcsine square-root transformed membership coefficients mis-assigned to each cluster. (B) Scatter plot with estimations of mean differences against Thracian mis-assignments based on raw individual membership coefficients (orange: C-lineage Thracian cluster, yellow: anatoliaca, blue: caucasica, violet: syriaca). Bars right to the data points refer to the 25 and 75% quartiles and the gap between them is the median value.

A significant Kruskal–Wallis test (p < 0.001) and a post hoc Dunn’s test, accompanied by a significant ANOVA result (p < 0.001) followed by a Tukey’s test, showed that mis-assignments to A. m. caucasica and A. m. anatoliaca clusters were significantly more frequent than to the other subspecies (p < 0.001 for both subspecies against C-lineage Thracian bees and p < 0.05 against syriaca group). The effect sizes according to Cohen’s d varied from 0.34 to 0.54 with estimation plots verifying the precision of the difference observed (Figure 6B). Despite the observation of the highest values in A. m. caucasica mis-assignments, the results between A. m. caucasica and A. m. anatoliaca clusters were not significant. We tested whether these differences were due to many individuals with high admixture levels, but such data only constituted 7.5% of all the observations. This figure is obtained by a threshold of 0.5 for the transformed values, which corresponds to a second hybrid, implying a 25% contribution of non-native origin.

We also investigated if these small drifts in admixture proportions were more prominent in some localities and if populations differed as to the identity of the subspecies from which they receive gene flow. This led us to learn the extent, magnitude, and direction of gene flow among the subspecies with a particular sensitivity to specific populations (Figure 7). Results of Dunn’s test for each pairwise comparison between populations are in Supplementary Tables 13–16 (12 significant differences out of 80 comparisons in total).

FIGURE 7

Figure 7. Patterns of gene flow between populations. (A) Boxplot displays of arcsine square-root transformed membership coefficients mis-assigned to each cluster (12 significant differences out of a total of 80 comparisons are provided in the Supplementary Material). (B) Scatter plots with estimations of mean differences based on raw individual membership coefficients to each cluster (orange: Thracian, yellow: Anatolian, blue: Caucasian, violet: Levantine clusters) contrasted against Kırklareli, Düzce +, Ardahan, and Hatay populations representative of four subspecies. Bars right to the data points refer to the 25 and 75% quartiles and the gap between them is the median value.

Discussion

Given the promiscuous nature of the honey bee mating system, it has been suggested that large-scale migratory beekeeping and bee trade have exposed local populations to introgression (De la Rúa et al., 2009). Although there is evidence that management actually increased genetic diversity (Harpur et al., 2012), admixture can also drive the loss of valuable local adaptations (De la Rúa et al., 2013). Since the global environment alters with an increasing pace, honey bees face new challenges in which they need to rely on adaptations and genetic diversity they accumulated over the course of their evolutionary history (Kükrer and Bilgin, 2020).

The main finding of this study is that there are distinct populations of subspecies of bees, isolated by distance, yet migratory colonies and bee trade likely cause gene flow across these populations in Turkey. The differences in F_ST values between stationary and migratory colonies indicate that the latter, with lower F_ST values, experience a high degree of gene flow. This conclusion is also reflected by the absence of positive correlation between genetic and geographic distances in migratory colonies in contrast to stationary colonies where an isolation by distance pattern was observed (Supplementary Tables 11, 12). Overall, F_ST values obtained were highly significant but lower than those from Bodur et al. (2007), estimated for samples collected 10 years prior to the study herein, that showed total levels of F_ST of 0.077 together with higher values for pairwise comparisons among populations. This may indicate recent increased gene flow and may signal an alarming trend toward greater movement of honey bees in the regions sampled. Long-term studies are needed to determine if this is a persistent trend.

Structure of Bee Populations in Turkey

PCA results confirmed the four different clusters inferred from the UPGMA tree topology (Figure 2). The first axis designating the first principle component differentiated Thracian samples, whereas the second axis, corresponding to the second component differentiated subspecies in Anatolia (syriaca, anatoliaca, and caucasica). Bitlis + samples clustered with Central and West Anatolian populations in both phylogenetic tree and PCA results (Supplementary Figures 1, 3, 4). However, all samples from this locality were from migratory colonies thus resampling this area with the inclusion of stationary colonies from East Anatolia would render a clearer picture of the phylogenetic relationship of these populations.

The two most likely K-values in structure analysis for all samples and the stationary colonies were K = 2 and K = 4, both results support the hypotheses of the sampled populations belonging to two separate lineages (C and O) in line with (Kandemir et al., 2005) and four distinct subspecies (a Carniolan ecotype in Thrace, A. m. caucasica in Artvin and Ardahan, A. m. syriaca in Hatay and A. m. anatoliaca, widely distributed, covering the rest of the country) (Figure 3). In contrast to the belief that migratory beekeepers make use of native stocks for their operations, our results showed the absence of structuring in these samples, and support the conclusion that migratory apiaries are highly hybridized.

Distinct Phylogeographic Patterns in Stationary Bees

Stationary apiaries, as expected, yielded highly structured groups where all the subspecies could be detected. When K was 2, the structure analysis of two distinct clusters showed that there was a transition zone between Thracian and Anatolian samples around the Marmara Sea and Aegean. Ruttner’s analysis based on morphometry (1988) distinguished bees in Western Anatolia from the rest of the anatoliaca group. Contributions from the Thracian cluster are significantly high in Düzce + and Eskişehir + located southeast of the Sea of Marmara across the Bosporus. Also, there are some non-significant overabundant Thracian contributions in Muğla province on the Aegean coast (Figure 3). This may constitute a hybrid zone between the C and O lineages and resemble the hybrid zones identified between M and C lineages in the Alps and the Apennine Peninsula and between A and M lineages in the Iberian Peninsula and Mediterranean islands (De la Rúa et al., 2009). An expected symmetrical introgression might be the reason behind the East–West cline observed by Muñoz and De la Rúa (2020) in four distinct ecotypes of A. m. carnica and A. m. macedonica in the Balkan Peninsula.

When K was considered as 4, all four subspecies were easily differentiated from each other. The significance of two distinct clusters (K = 2) was higher than four (K = 4) indicating the evolutionarily greater differences between the lineages belonging to C (in Thrace) and O (in Anatolia).

Thracian samples form a clade in the unrooted phylogenetic tree while the other three populations, Caucasian, Levantine, and Anatolian are equidistant from each other and form a separate clade. These results indicate that the Thracian population is distinct from the others and likely has experienced limited gene flow in allopatry, supporting the hypothesis for a Carniolan (C-lineage) descent of Thracian bees in Turkey. A direct comparison with honey bee samples of the major C-lineage subspecies would confirm the subspecies of these bees which are highly differentiated from Anatolian samples. This finding is in contrast to the conclusions of Ruttner (1988) that Thracian bees are part of the anatoliaca subspecies groups and merits further investigation.

A. m. anatoliaca samples formed a distinct cluster in structure analysis, yet fell in the middle of the other subspecies in ordinations according to F_ST values. This similarity may point to a significant historical contribution to A. m. anatoliaca populations from the neighboring regions. Another explanation is that the putative basal position of anatoliaca for O-lineage honey bees places this group at the center of genetic diversity. In contrast to anatoliaca bees forming a distinct group, all-migratory Bitlis + samples were a mixture of different clusters and did not form a separate group.

A greater understanding of phylogenetic relationships of the populations of bees in Turkey can be achieved only if neighboring populations in the Balkans, Iran, Caucasus, and Southwest Asia are also sampled. This future research direction may clarify the complex taxonomic relations within and between the C and O lineages, and delimit distributions and transition zones of the subspecies in this region.

Homogenizing Effect of Migratory Beekeeping

Migratory colonies are acting as a hybrid zone mobile in space and time. The colonies are in one region in spring and in others in summer and fall. As such, these bees serve as vectors of otherwise local gene combinations. Statistical comparison of migratory and stationary colonies confirms the significant gene flow toward the migrants from local bees (Figure 4). Likewise, a significant gene flow toward local stationary bees was also observed outside the conservation sites. These results, derived from direct comparison of two distinct contrasts, demonstrate the vitality of establishing areas away from migratory beekeeping for the preservation of honey bee genetic diversity. This conclusion is in agreement with other studies on conservation practices (Oleksa et al., 2011; Pinto et al., 2014).

An exception that proves the point is the lower assignment probability of bees sampled from Ankara to their province, even in comparison to migratory bees in the same location. There are two factors: First, the region’s beekeepers prefer to use queen bees native to the region. The second factor is that this region is a principal queen breeding area. The Kazan apiary of TKV (Development Foundation of Turkey) uses hundreds of colonies of Caucasian bees and raised queens are sold around the country for over 30 years. Many independent queen bee breeders in the Kazan region continue the same practice. Gene flow from these queen breeders’ apiaries may contribute to the admixture observed in stationary colonies in Ankara. The high mis-assignment probability of colonies in Ankara to the Caucasian cluster supports this hypothesis.

Direction and Magnitude of Introgression Determined by Bee Trade

It is hard to directly quantify the effect of queen and colony trade on genetic mixing. The availability of several naturally occurring subspecies in Anatolia and Thrace helps in understanding the relative role of queen and colony trade in gene flow. Honey bees from stationary colonies were assigned more often to their native clusters, yet they were also assigned to other clusters with lower probabilities. Samples in the whole range of the study mis-assigned to the Caucasian cluster more often than they were mis-assigned to others (Figure 6). This is most likely due to the wide distribution of Caucasian queen bees by queen trade.

Migratory beekeeping is not practiced in Ardahan and Artvin where highly commercial Caucasian bees are native. Hence, no bees go in or leave out the region as migratory colonies. We infer that the observed introgression of Caucasian alleles to the stationary colonies elsewhere could mainly be attributed to the frequent purchase of caucasica queen bees and colony replacements in neighboring apiaries within those regions. Practices of neighboring beekeepers become important because even if beekeepers included in this study within a region do not purchase caucasica queens, their colonies may be subject to queen supersedure and natural mating in that region.

Central and western Anatolian populations suffer heavily from gene flow from Caucasian populations as demonstrated by our results (Figure 7). Muğla, which receives millions of migratory colonies during the honeydew season, and Ankara showed high levels of significant gene flow from other subspecies, especially the caucasica. This is especially alarming because Muğla (in the southwest) and the Caucasus region (in the northeast) lie at the diagonal extremes of the country, some 1,500 km apart.

A. m. anatoliaca alleles also showed high introgression especially in the Thrace region but also at average levels in other regions. These high levels may be related to the geographical proximity of this subspecies to other populations. The proximity may explain historical and recent gene exchanges. Alternatively, widespread practice of migratory beekeeping by Western and Central Anatolian beekeepers throughout Turkey may have contributed to observed introgression. In this case, queen replacement could be a minor contributor since there are very few commercial queen breeders within the distribution range of A. m. anatoliaca.

Conservation Sites

The importance of establishing isolated regions was highlighted with genetic data. The results of the statistical tests showed a significant difference between the conservation of identity in and out of isolated regions with isolated regions staying purer in terms of subspecies composition (Figure 5). Such regions were proven to be effective in the conservation of unique diversity present within (Requier et al., 2019).

In the light of this study, we propose a renewed effort to address the need for massive establishing of such regions for conserving locally adapted native bees throughout the whole natural distribution of the species. This especially holds for underrepresented regions in terms of local diversity hotspots. A gap analysis aiming for complementarity in the planning of systematic conservation efforts is urgently needed globally.

In such isolated regions, naturally, migratory beekeeping, as well as replacement of queen bees with non-native ones, must be strictly prohibited and checked by relevant molecular monitoring techniques. However, these isolated regions should also be wide enough involving additional buffer zones where further restrictions on migratory beekeeping and bee trade are applied for efficient isolation and for fulfilling sufficient effective population sizes.

Thanks to increasing awareness in the last decade within the industry, now there are at least 11 isolated regions in service or being established in Turkey. These conservation sites make ideal places for breeding purposes. The establishment of such sites is achieved through the significant efforts of scientists and their collaboration with the Turkish Beekeepers Association (Kükrer and Bilgin, 2020). There is an ever-growing need for establishing closer links with decision makers and stakeholders and the necessity of investing more efforts in communicating the results of scientific studies to all involved.

Conclusion

Overall results of this study clearly show that the genetic structure of honeybee populations in Turkey is highly conserved. This, however, does not mean that the structure and the diversity observed are secure. Rather the honey bee genetic diversity in Turkey should be considered under threat. We demonstrated continued gene flow and admixing of populations, likely due to anthropogenic factors.

The preservation of population structure despite movement of the high number of colonies and unregulated and frequent queen and colony sales is biologically interesting. Future research may also need to focus on how this biodiversity and its structuring were preserved and its relation to natural selection. The relative effects of natural selection and gene flow should be compared; the former could significantly counterbalance the latter.

Genetic variation eventually leading to local adaptations with such a significant outweighing effect can be considered as a valuable resource for honey bee populations in the global context at this time of unusual bee losses as well as global climate change. A better understanding of present adaptation to both local climate and geographic conditions as well as adaptive capacity to future changes is important for bees and stakeholders. A fair amount of effort should be invested in more studies focusing on candidate functional variants at the genome level that play role in due process in different parts of the world. Novel and innovative ways of coping with environmental and climatic stressors developed by honey bee populations or exploration of interesting patterns of convergent evolution are waiting ahead to be yet discovered.

Our overall results are restricted to the present situation of honey bee subspecies in Turkey, yet they highlight the significance of local populations and provide a preliminary quantification of human impact. We expect our findings on migratory beekeeping, trading of queens and colonies as well as conservation implications to be of use for the decision makers and other stakeholders.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

MKü carried out the experiments and statistical analyses and took the lead in writing the manuscript. AK was in charge of overall direction and planning. All authors conceived and planned the experiments and contributed to the fieldwork, contributed to the interpretation of the results, provided critical feedback and helped shape the research, analysis, and manuscript.

Funding

This study was funded by the Scientific and Technological Research Council of Turkey (project no: 109T547), Republic of Turkey Ministry of Agriculture and Forestry (project no: TAGEM/11/AR-GE/13), and Middle East Technical University Revolving Funds (project no: 2007-16-12-00-3008).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Dedication

This article is dedicated to the memory of the late AK who spent his last 20 years studying honey bee diversity in Turkey trying strongly to draw the attention of researchers, decision makers, and beekeepers to the conservation of locally adapted native bees as a precious legacy of our kind. He will be remembered also for defending the theory of evolution to be taught in the science curriculum and for training many valuable young evolutionary biologists under very harsh conditions. He passed away on February 1, 2014 after the completion of the study and at the very beginning of the manuscript preparation.

Acknowledgments

We would like to thank the numerous beekeepers who provided samples to the study and the Turkish Beekeepers Association. We thank to Devrim Oskay for his contribution in the form of consumables. We would also like to thank Mustafa Nail Cırık, Mehmet Ali Döke, Okan Can Arslan, Cansu Özge Tozkar, Mehmet Kayım, and Eda Gazel Karakaş for the fieldwork and Esin Öztürk and Ezgi Ersin as well as Batuhan Çağrı Yapan, Ayshin Ghalici, Gizem Kars, and Batuhan Elçin for contributing lab sessions. We thank to Çiğdem Akın Pekşen, Cemal Can Bilgin, Tuğrul Giray, and Felipe Soto for their advice on the manuscript. This manuscript has been released as a pre-print at https://www.biorxiv.org/content/10.1101/154195v2 (Kükrer et al., 2020) and includes elements from the corresponding author’s thesis (Kükrer, 2013).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021.556816/full#supplementary-material

References

Akkaya, H., and Alkan, S. (2007). Beekeeping in Anatolia from the Hittites to the present days. J. Apicultural Res. 46, 120–124. doi: 10.3896/ibra.1.46.2.10