ABSTRACT
Intense admixture of honey bee (Apis mellifera L.) populations is mostly attributed to migratory beekeeping practices and replacement of queens or colonies with non-native races or hybrids of different subspecies. These two practices are also heavily carried out in Turkey where 5 subspecies are present naturally.
Here, we carried out an analysis of population structure of Turkish honey bees sampled from six different regions (n = 250) in order to test the genetic impacts of migratory beekeeping, queen and colony trade and conservation efficacy of isolated regions. A total of 29 microsatellite markers were used in four multiplex reactions.
Direct genetic impact of migratory beekeeping was demonstrated first time based on a comparison of assignment probabilities of individuals to their geographic populations where migratory colonies showed less fidelity. The effects of queen and colony trade were revealed by the presence of very high introgression levels from the highly commercial Caucasian bees naturally limited to a narrow range. Comparison between regions that are either open to migratory beekeeping or not let us evaluate the status of isolated regions as centers of limited gene flow and showed the importance of establishing such regions.
Despite the signs of gene flow, our findings confirm high levels of geographically structured genetic diversity of four subspecies of honey bees in Turkey and emphasize the need to develop policies to maintain this diversity.
INTRODUCTION
The Western honey bee, Apis mellifera L., is a species which plays role together with other pollinators in pollination of wild and cultivated plants while the species also have significant economic importance in terms of honey and other bee products output (Morse 1991; Breeze et al. 2011). In addition to its ecological and economic importance, it is a model study organism both for evolution of eusociality and sophisticated cognitive abilities (Weinstock et al. 2006).
Natural distribution of A. mellifera includes Central Asia, Europe, Near East and sub-Saharan Africa but the species was also introduced to East and Southeast Asia, Australia and the Americas mainly on purpose 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. Studies with Single Nucleotide Polymorphisms (SNPs) in the past decade supported the hypothesis that A. mellifera have originated in the tropics or subtropics in Africa and colonize its natural range by two main routes: one through Gibraltar and one through Suez and then Bosphorus, ending up with a secondary contact between highly divergent A and C lineages around Alps (Whitfield et al. 2006).
Beekeeping is a very old tradition which dates back to 6600BC Hittites civilization and is still intensively practiced in Turkey where there are more than 6 million hives distributed all over the country, second highest in the world (Akkaya and Serhat 2007). Five different subspecies of A. mellifera -which are A. m. meda, A. m. syriaca, A. m. caucasica, A. m. anatoliaca from the O lineage and an ecotype from C subspecies group- exist in Turkey (Kandemir et al. 2005).
Major subspecies found in and around Anatolia are shown in Fig. 1a. Anatolia and Thrace, when considered together, harbors a vast diversity: at least five honey bee subspecies belonging to two different lineages meet, exchange genes and adapt to local conditions determined by diverse climatic, topographical and floristic variations available (Bouga et al. 2011). Studies concerning honey bee populations of Turkey (Bodur et al. 2007; Kence et al. 2009) demonstrated the high genetic structuring among them and confirmed presence of divergent populations pointing to different subspecies. They, all together, drew attention to the genetic diversity present in Anatolia and Thrace and to importance of its conservation.
Both the honey bees and wild pollinators are thought to be on decline (locally and/or globally depending on the species and region of concern) due to factors some of them closely related to human activities. Among them, destruction and fragmentation of natural habitats, toxicity caused by pollution and pesticides –like widely used neonicotinoids-, diseases and their spread getting easier, invasive species are leading the way (Meffe 1998; Brown & Paxton 2009; Van Engelsdorp & Meixner 2010; Blacquiere et al. 2012). Honey bees also, especially wild populations that are not managed by beekeepers (including the feral populations), take their share from the situation (Oldroyd 2007; Dietemann et al. 2009; Van Engelsdorp et al. 2009; Genersch 2010; Evans & Schwarz 2011).
Besides such negative consequences created by human activities; the genetic admixture of honey bee populations due to queen and colony trade, including complete replacement of local bees with non-natives and beekeeping practices involving movement of colonies from one region to the other impose another kind of pressure on the species: the loss and/or swamping of locally adapted gene combinations and local or global extinctions of native honey bees, not only of A. mellifera subspecies, races or ecotypes but of other species in the genus Apis –namely Apis cerena, Apis florea, Apis dorsata and other native bees of Asia (De la Rua et al. 2009).
All those factors and their interactions, including genetic and environmental ones, when combined, may have an increased effect on honey bees and can be the reasons behind continuous or discrete events of sudden colony losses with rapid depletion of worker bees while the queen continues to laying eggs accompanied by lack of dead bees in and around the hive; the syndrome called as Colony Collapse Disorder (CCD) or Colony Depopulation Syndrome (CDS) (Van Engelsdorp et al. 2009; Neumann & Carreck 2010).
Resilience of the honey bees may be lying in the adaptations they accumulated over thousands of years, and in the new potentials resided in their genetic diversity. It is highly suspected that a combination of many above mentioned factors/threats are taking their places in the recent declines by weakening the colonies step by step. Also, it is known that honey bees’ resistance or tolerance to these factors differ greatly and locally adapted variants may be encountering less stress thus standing more upright. Hence, research on honey bee diversity in the global context and at various levels (genetic, individuals, colonies, populations, ecotypes and subspecies) is of great importance for maintaining the species’ and ecosystem services they provide as well as their economic usefulness.
In the recent years’ studies conducted on honey bee population structure in European countries, it was shown that the past structure was lost or strongly disturbed (Dall’Olio et al. 2007; Canovas et al. 2011; Bouga et al. 2011). Introgression of non-native DNA was monitored in wild populations of Sudan (El-Niweiri & Moritz 2010). Among the anthropogenic effects, mainly queen and colony trade and replacement of native honey bees with non-natives as well as migratory beekeeping were the usual suspects. Therefore, the aims of this study were testing different hypotheses about recent heavy/any admixture of populations in Turkey by making use of microsatellite markers as well as i) evaluating the status of isolated regions declared as a conservation implication where migratory beekeeping is prohibited, restricted or very scarce due to lack of preference of migratory beekeepers or attitude of local beekeepers ii) acquiring and demonstrating the direct genetic outcomes of migratory beekeeping by a series of comparisons between migratory and stationary colonies iii) seeking for the effects of unregulated queen and colony trade by figuring out the origin of introgression between populations. With five subspecies dwelling within its borders and with a variety of beekeeping strategies, Turkey makes a good stage for chasing genetic evidences for the impact of anthropogenic factors on one of the most important crop and wild plant pollinators.
METHODS
Sampling
We sampled a total of 250 honey bees each from different colonies from 18 provinces during the period of March 2010 and August 2012. Of those 250 honey bees, 174 were from apiaries that were stationary and 76 were from migratory ones. We grouped samples from provinces with a small sample size with nearby provinces to form 10 major localities: Kırklareli, Edirne+ (Edirne and Tekirdağ), Muğla, Eskişehir+ (Eskişehir, Kütahya and Bilecik), Düzce+ (Düzce, Zonguldak and Bolu), Ankara, Hatay, Bitlis+ (Bitlis, Elazığ, Erzurum and Ordu), Ardahan, and Artvin. We carried out combinations according to geographical proximity; similarity in terms of climatic, topographic and floral variables; results of previous studies as well as results of preliminary analysis of this study. Sampling sites and sample sizes can be seen in Fig. 1b. The samples were kept in −80°C until genetic analysis.
Genotyping
We isolated DNA from bee heads by QIAGEN DNeasy Blood and Tissue Kit following the procedure of the producer for insect samples with slight modifications. We grouped a set of 29 microsatellite loci into four clusters for three 7-plex (AP218, A113, AB024, AP249, A088, AP001, AP043 and AP019, AB124, A043, AP273, AP289, HBC1605, A028 and AP049, AP238, AC006, AP243, AP288, HBC1602, A107) and one 8-plex (A079, AC306, AP226, A007, HBC1601, AP068, A014, AP223) polymerase chain reactions. Detection of microsatellite allele sizes was achieved by capillary electrophoresis with ABI 3730XL sequencing machines.
Population Structure
We calculated pairwise FST values by Arlequin 3.5 (Excoffier et al. 2005) and the pairwise population distances by the Populations 1.2.32 software which were visualized by the same program (Langella 2011). We used PCA-gen software to plot populations on a two-dimensional space (Goudet 1999). Population structure was estimated by Structure 2.3.3 (Pritchard et al. 2000), K values of distinct populations were analyzed by Structure Havester software (Earl & von Holdt 2012), and we used Clumpp software to permute the membership coefficients of individuals determined by Structure 2.3.3 (Jakobsson & Rosenberg 2007) and Distruct software to visualize the results obtained by Clumpp (Rosenberg 2004).
Statistical analyses
We then used membership coefficients obtained to test hypotheses about beekeeping practices, isolated regions and queen/colony trade. For the analysis, we first arcsine transformed the coefficients. Then we carried out Mann-Whitney U and t tests to compare mean membership coefficients.
Beekeeping practice: migratory vs stationary
For the first hypothesis to be tested, we compared assignment probabilities of migratory and stationary colonies in Ankara, Muğla and Hatay separately, for the three provinces combined and for the total data set. If the migratory colonies acted as a potential vector of foreign alleles then they would have much lower probabilities of being assigned to their own clusters.
Isolated regions as a conservation practice
The second hypothesis was about isolated regions. If the isolated regions were efficient in preserving genetic diversity by preventing gene flow between different clusters then one would expect to see a higher assignment probability for stationary individuals belonging to these regions and lower for stationary individuals that belong to regions open to migratory beekeeping.
Kırklareli is a province that is declared officially as an isolated region where migratory beekeepers could not visit for years at first thanks to local beekeepers’ negative attitude towards them. The region is home to a Carniolan ecotype carefully maintained by local beekeepers. Ardahan is legally declared a conservation and breeding area for A. m. caucasica so migratory beekeepers cannot enter the province and queen import from other subspecies is forbidden. A part of Artvin is also officially declared as an isolated region for conservation of A. m. caucasica as a pure race. The province in general is rarely visited by migratory beekeepers for geographical reasons and beekeepers there, dealing with mass queen breeding, do not use non-native queens. We compared these three provinces 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
Third set of tests were about the impacts of queen and colony trade. We compared the probabilities of being assigned to a different cluster than the native cluster among individuals of the total data set to find out which cluster contributed most to other clusters’ gene pools. Ardahan and Artvin provinces host the A. m. caucasica subspecies which is also widely used for commercial purposes and the caucasica queens and their hybrids are sold all over the country. But these provinces are also isolated regions so a possible high introgression of their alleles would mostly, if not completely, be due to queen and colony trade.
RESULTS
We calculated FST values by using both the obtained frequencies in this study and by using the null allele corrected frequencies. We calculated for the stationary (n = 174) colonies an overall FST of 0.065 and an FST of 0.067 after null allele corrections. For migratory colonies the values were 0.011 and 0.015 respectively and for all the 250 samples the values were 0.046 and 0.047. We used pairwise population distances to construct phylogenetic trees for the populations in the regions of sampling (Fig. 2) where 4 distinct branches were easily resolved. We plotted stationary colonies, migratory colonies and the overall data for the regions of sampling on 2D spaces by carrying out Principle Component Analysis (Fig. 3) which showed a similar pattern with the UPGMA tree.
The best K values were selected by the Structure Harvester program as 2 and 4 with close results, K=2 slightly being higher than K=4. We calculated membership coefficients of individuals to the observed clusters in K=4 and used them for further hypothesis testing. Clustering analyses showed no population structuring for migratory colonies (Fig. 4b) in contrast to stationary colonies and the overall data (Fig. 4c and 4d).
We compared individuals from stationary and migratory colonies according to their membership coefficients to their native clusters (or natural populations as similar). Stationary colonies from Muğla and Hatay were quite more likely to be assigned to their own clusters than the migratory colonies from these provinces (p<0.01), the same held when we compared the data from the three provinces (p<0.01) or all the migratory and stationary colonies (p<0.001), however the situation was the reverse in Ankara (Table 1). When we compared isolated regions (Kırklareli, Ardahan, Artvin) and regions open to migratory beekeeping (Edirne+, Muğla, Ankara, Düzce+, Eskişehir+, Hatay) in terms of their membership coefficients to native clusters we have seen that (p<0.001) stationary colonies within isolated regions showed higher fidelity to native clusters (Table 1).
Even if the individuals are assigned with high probability to their own clusters, let’s say with a 90% of probability, this means that 10% of their genome still belongs to other clusters. Given there are four clusters, we investigated if any of these “mis”assigned genome parts were over represented by any of the clusters. This comparison (Table 2) showed that “mis”assignments to A. m. caucasica and A. m. anatoliaca clusters were significantly more frequent than the others (p<0.05). Despite observation of the highest values in A. m. caucasica “mis”assignments, the results were not significant between A. m. caucasica and A. m. anatoliaca clusters.
DISCUSSION
FST values obtained were highly significant but they were lower than what Bodur et al. (2007) estimated -a total FST of 0.077 and also higher values for pairwise comparisons among populations- with the samples collected ten years before ours. This may indicate a recent increased gene flow and can be an alarm signal for a trend. Similar studies should be actualized in the future to see if there is such a trend really. The high degree of structuring in stationary colonies according to FST results was lost in migratory ones, meaning they are less differentiated from each other due to high degree of gene flow.
Phylogenetic tree clearly showed that Thracian samples were completely distinct from others pointing to an early division of populations and limited gene flow. This supports the hypothesis for a Carniolan (C-lineage) descent of Thracian bees in Turkey. Making use of samples from the major C-lineage subspecies would confirm the subspecies of these bees highly differentiated from Anatolian samples. West Anatolian, Hatay and Caucasian populations did also form separate clusters in the tree. PCA results confirmed those 4 different clusters inferred from tree topology. Bitlis+ samples resided with Central and West Anatolian populations in both phylogenetic tree and PCA results.
The two most possible K values in structure analysis for the whole sample and the stationary colonies were K = 2 and K = 4, both results supporting the hypotheses of populations belonging to 2 separate lineages (C and O) and 4 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 regions). Results for migratory beekeepers’ samples lacked any population structuring in the cluster analysis further clarifying the highly hybridized status of migratory apiaries.
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 Marmara Sea and Aegean. This may be a hybrid zone between the C and O lineages like the ones identified before between M and C lineages in Alps and Apennine Peninsula and between A and M lineages at the Iberian Peninsula and Mediterranean islands. When K was considered as 4, all four subspecies were easily differentiated from each other, in accordance with the expectances. The significance of two distinct clusters (K = 2) was higher than four (K = 4) which means that the differences between the populations belonging to C (in Thrace) and O (in Anatolia) lineages are more clear-cut than differences between the populations of four different subspecies. A. m. anatoliaca samples fell in the middle of the other subspecies in ordinations, being similar to all other populations according to FST values despite being a distinct cluster in structure analysis which may point to a significant historical contribution to A. m. anatoliaca populations from the neighboring regions. However this was quite different than what was observed in let’s say, all-migratory Bitlis+ samples where we observed a mixture of different gene pools instead of a distinct identity.
A better understanding in terms of phylogenetic relationships between the populations in Turkey would be developed if populations neighboring Anatolia and Thrace were also sampled. This can be a direction for future research too for shedding light on the complicated taxonomic status within and between the C and O lineages and for drawing edges and transition zones of the subspecies present in the region.
Results from different analyses conducted here confirmed the presence of clusters but also they all together pointed to the status of migratory colonies: they might be acting as a hybrid zone mobile in space and time, being at one region in spring and at others in summer and fall, becoming vectors of otherwise local gene combinations. Statistical results concerning a comparison between migratory and stationary colonies confirmed the significant gene flow towards the migrants from local bees. A significant gene flow towards local bees was also observed in the comparison between isolated and not-isolated regions. This result is pointing to the vitality of establishing areas away and free from migratory beekeeping for preservation of honey bee genetic diversity.
One interesting point in those was that the trend of the stationary colonies in Ankara. They had a significantly lower probability of being assigned to their own clusters. This may be related with the regions migratory apiaries visit during their migratory cycle or due to preference of using native queen bees by migratory beekeepers. The low assignment degree of stationary colonies in Ankara may also be related with Kazan apiary of TKV (Development Foundation of Turkey) placed there where hundreds of colonies of Caucasian bees are raised and sold around for more than 30 years. The same practice is also carried out by many queen bee breeders in Kazan region. Gene flow through these apiaries and queen bees distributed locally by trade may contribute quietly to such an admixture observed in stationary colonies in Ankara. The high “mis”assignment probability of colonies in Ankara to the Caucasian cluster also revealed such a process as probable.
Honeybees from stationary colonies were assigned more often to their native clusters but they were also assigned to other clusters with lower probabilities. Samples in the whole range “mis”assigned to Caucasian cluster more often than they were “mis”assigned to others. This is probably due to wide distribution of Caucasian queen bees by trade. Migratory beekeeping is not practiced in Ardahan and Artvin where highly commercial Caucasian bees are native. So, the observed introgression of Caucasian alleles to the stationary colonies elsewhere whose beekeepers let them change their queens on their own rather than purchasing queens of different origins, could mainly be attributed to frequent queen bee and colony replacements within those regions. It is shown here that practicing of queen and colony trade is increasing the level of introgression within the gene pool of honey bees of Turkey. As previously discussed, a very high level of Caucasian introgression was observed in Ankara. A. m. anatoliaca alleles also showed high introgression especially in Edirne+ of Thrace region but also at average levels in other regions. These high levels may be related to this subspecies’ geographical proximity to other populations which might have led to historical and recent gene exchange, or the widespread practicing of migratory beekeeping by Western and Central Anatolian beekeepers throughout Turkey, rather than queen or colony replacements since there are very few commercial queen breeders within the distribution range of A. m. anatoliaca.
Results of the various statistical tests carried out and analysis applied in this study clearly showed that the genetic structure of honeybee populations in Turkey were highly conserved and still maintained. But it doesn’t mean that the structure and diversity observed is secure. Rather it should be considered under threat since the anthropogenic factors leading to gene flows are still functioning and admixing the populations. A quiet interesting point was that, the preservation of population structure was achieved despite a very high number of colonies moved from one location to the other by migratory beekeeping practice and despite unregulated and frequent queen and colony sales. Future research may also need to focus on how this biodiversity and structuring were preserved and its relation to natural selection.
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. Such regions were proven to be effective in conservation of unique diversity present within. In the light of this study we propose massive establishing of such regions for conserving locally adapted native bees throughout the natural distribution of the species. In such isolated regions migratory beekeeping must be strictly prohibited as well as replacement of queen bees with non-native ones. However these isolated regions should also be wide enough involving additional buffer zones where restrictions on migratory beekeeping and bee trade are applied for efficient isolation and for fulfilling sufficient effective population size.
Queen bee trade is not currently subjected to any restrictions or regulations in Turkey and to our knowledge there are very few pioneering measures within the natural distribution range. Such measures should be applied from a conservation perspective to avoid extinction of native races, ecotypes and diversity present in these populations. Genetic similarity of donor and recipient populations should be considered while determining migration routes for migratory beekeepers and determining permissions for queen and colony trade.
Central and western Anatolian populations suffer from significant gene flow from Caucasian populations. Despite its wide range of distribution spanning Anatolia from one side to other, special consideration should be taken for preserving A. mellifera anatoliaca subspecies. The case with Hatay populations too, can get worse and worse since the migratory beekeeping practice is heavily carried out in the region and queen bee replacement with non-native races was frequent in the last years. In the future this may end up in A. m. syriaca colonies getting limited to a few localities and apiaries since the range of the subspecies in Turkey is very narrow. A conservation program in the light of these findings should be actualized immediately in this region too. Thracian populations show a significant differentiation from the rest of the bees in Anatolia but the subspecies which they belong to is not characterized on a strong basis yet and this unique population is not registered officially like the case with A. m. syriaca of Hatay. Only subspecies officially recognized in Turkey is A. m. caucasica so identification and registration procedures should be put into practice as soon as possible. An improvement based on molecular genetic techniques can be applied to the ongoing conservation programs for the A. m. caucasica subspecies. This proposal holds for other subspecies too. More attention should be paid to genetically characterize A. m. meda subspecies that was out of the reach of this study and which can be threatened by anthropogenic factors listed and studied here.
Rather than queen bee replacement, use of locally adapted native bees improved for desired characters should be preferred. Such improved breeds would be used locally and not distributed in a country-wide manner so that local adaptations can be preserved while bees are selected for resistance to pests and pathogens, hygienic behavior, reduced aggressiveness, reduced tendency for swarming, higher winter survival, higher productivity or for increased pollination. For obtaining better results, research concerning the smoothing, development and extension of breeding locally adapted native bees and artificial insemination techniques should be given higher priority, globally throughout the natural distribution range of local subspecies.
Acknowledgments
We dedicate this paper to the memory of late Professor Aykut Kence who spent His last years studying honey bee diversity in Turkey together with struggling to train many valuable young evolutionary biologists under very harsh conditions and tried to draw attention of researchers, decision makers and beekeepers to the conservation of locally adapted native bees as a prescious legacy of our kind. He passed away after the completion of the study and at the very beginning of manuscript formation. This study was funded by The Scientific and Technological Research Council of Turkey (project no: 109T547).
Footnotes
↵2 Deceased