Microsatellite Analysis of Apis mellifera from Northern and Southern Parts of Serbia

: Practice of commercial honey bee breeding and selection for desired traits, intensification of queen importation and migration of once stationary apiaries significantly influences distribution and genetic diversity of local subspecies, populations and ecotypes. Bee colonies worldwide are facing serious declines resulting in colonies loss and reduction of genetic diversity. Thus, reassessing the genetic status of native honey bee populations becomes imperative. The latest reports, which include samples from nine years ago, suggest the presence of both Apis mellifera carnica in north and A. m. macedonica in the south of Serbia and significant hybridization between two subspecies.


Introduction
The western honey bee (Apis mellifera Linnaeus, 1758) is considered to be the most important managed pollinator species currently facing serious challenges in its whole area of distribution [1].It is a highly polymorphic species with approximately 30 described subspecies [2] sometimes referred to as geographic subspecies due to their distribution in specific geographic areas.Additionally, several ecotypes and breeding lines exist within subspecies reflecting significant variation and adaptation to regional climate factors, vegetation, pests and pathogens.Based on morphological, genetic and ecological traits, all subspecies are divided into five evolutionary lineages of which Europe harbours four: the African lineage (A), the Southeast Europe lineage (C), the West and North European  lineage (M) and the Near and Middle Eastern lineage (O) [3].Commercial honey bee breeding and selection for desired traits, intensification of queen importation and migration of once stationary apiaries considerably influence this natural diversity and distribution range which may result in loss of both genetic variability and local adaptations as well as genetic homogenization of admixed populations [4][5].Serbia, geographically in the middle of the natural area of distribution of the Southeast Europe lineage C, has a long tradition of beekeeping and strong breeding and beekeeper organizations.Serbian honey bees are thoroughly characterized on morphological, genetic, etiological and ecological levels revealing the presence of A. m. carnica in the north-west, A. m. macedonica in the south-east [6][7][8], their hybrids [6,9], three ecotypes of A. m. carnica (Banatski, Sjeničko-Pešterski, and Timočki) [10] and locally adapted populations.Both mtDNA and microsatellite variability of Serbian honey bees is well described.Widely accepted, highly variable mtDNA COI-COII was used in several studies [7,9,[11][12] and more than seven haplotypes were reported showing significant genetic variation and possible ways of introduction of different mitotypes.Investigation of microsatellite variability [6,8] showed clinal south-north and west-east distribution of two distinct honey bee populations that corresponds to A. m. carnica and A. m. macedonica and their hybrids.
However, all molecular investigations both on microsatellite and mtDNA loci, even the most recently published, are based on samples from the first decade of the 21 st century.For example, the samples used in previously mentioned studies were gathered from 2007-2012.[6,8] and don't reflect the current situation in Serbia.Since then, beekeeping practice in Serbia has significantly changed and the number of beekeepers has vastly increased while the traditional way of beekeeping has been abandoned and modern beekeeping practices have been utilized.Foremost, a number of queen breeding institutions focused on desired traits of A. m. carnica were established in concordance with Serbian legislation which allows breeding of this subspecies only [13].Long-distance migratory beekeeping became more prevalent and importation of Carnolian queens intensified.
To shed a light on the current status of genetic diversity in Serbian managed honey bees we sampled apiaries from northern and southern parts of Serbia.

Methods
For this study, 227 beehives originating from 46 apiaries and the same number of different locations were sampled (Figure 1.).These locations were grouped in eight localities from southern and northern parts of Serbia: four in South (Leskovac (L), Vlasina (V), Stara Planina (SP) and locality between L-V-SP hereby named as Tromedja (T)) and four in North (Subotica (S), Vršac (VR), Deliblatska peščara (DP) and Fruška gora (FG)).Grouping into eight localities was performed based on the geographical proximity of the sampled apiaries and regional division into districts.Apiaries were chosen according to the following criteria: they needed to be stationary, desirably old and with no or minimal queen importation.
Whole-genomic DNA was extracted according to the protocol described in [14].To distinguish A. m. carnica from A. m. macedonica PCR-RFLP method described in [5] was used.Primers 5'GATTACTTCCTCCCTCATTA3' and 5'AATCTGGATAGTCTGAATAA3' were used for amplification of mtDNA COI segment, PCR products purified according to [15] digested with restriction enzymes NcoI and StyI and then visualised on 2% agarose gel stained with ethidium bromide under UV light.For fragment analysis, a total of 14 microsatellite loci was used (A7, A8, A14, A24, A28, A35, A43, A79, A88, A107, A113, Ap43, Ap249, B124).Microsatellite loci were amplified in four multiplex PCR reactions with different annealing temperatures in which forward primers were labelled with a fluorescent dye (FAM, NED, VIC and PET).Parameters for assessing genetic diversity (number of alleles, mean allelic range, average gene diversity over loci, observed (HO) and expected (HE) heterozygosity, random match probability and mean number of pairwise differences), pairwise population, overall FST values and Analysis of molecular variances (AMOVA) were calculated in Arlequin ver.3.5.2.2 [15].The matrix of pairwise population FST values was visualized both by R functions and presented in two dimensions by the means of multidimensional scaling (MDS) using PAST 3.25 [17].The number of genetic clusters represented in our sample was estimated in the STRUCTURE v 2.3.4 software [18] using the admixture model with a burn length of 10,000 and a Markov Chain Monte Carlo (MCMC) of 100,000 randomizations.STRUCTURE harvester [19] was used to analyze the results obtained by STRUCTURE.Clustering of the analyzed populations and the observed distances among samples were evaluated by discriminant analysis of principal components (DAPC) [20] using the R package adegenet [21].

PCR-RFLP
The size of the PCR amplified COI segment for RFLP analysis was 1029bp.Digestion with both NcoI and StyI showed no restriction pattern characteristic for mtDNA lineage found in A. m. macedonica.Since no restriction sites were observed after RFLP analysis, we presume that all individuals in our sample belong to A. m. carnica.

Genetic diversity analysis
Standard diversity parameters obtained from Arlequin are presented in Table 1.Average numbers of alleles and observed heterozygosity are the highest in L, while T locality showed the lowest average number of alleles.The average number of pairwise differences between and within populations together with Nei's distances is visualized in Figure 2a and pairwise population FST is visualized in Figure 2b.Differences between some pairs of localities are consistent in all analyses with T locality showing the highest number of significant pairwise FST differences, both compared with geographically proximate south and distant north localities.
The AMOVA performed across all 14 loci showed a very low value of genetic variance among localities (0.01) and only 0.64% of the genetic variance can be attributed to the variation among the localities with no statistical significance (p=0.89).When localities were grouped according to their geographical region, the percentage of variation is higher among regions than within regions.In addition, when localities are grouped according to region, the percentage of variation among localities decreases reflecting similarity within the region.However, these results are not statistically significant (Table 2).Results of analysis performed by the DAPC method are shown in Figure 3a and visualization by MDS plot shows the positioning of populations in two dimensions (Figure 3 b).Although individuals from T, SP and DP appear to cluster separately from others, DAPC analysis showed that individuals from geographically remote localities are grouped in cluster overlap indicating similarity between them.However, the MDS plot positioned localities from north and south on the opposite sides of the first axis, which is in correlation with AMOVA that suggests the existence of two groups.STRUCTURE Harvester showed that K = 3 is the most likely scenario.All three clusters inferred by STRUCTURE are equally distributed in all sampled localities.Additionally, the number of clusters inferred with the DAPC method was 3 and all clusters were presented in all localities.

Discussion
Genetic variability of managed honey bees is under the strong anthropogenic influence due to selection for desired traits, queen importation and migratory beekeeping.Previous work [6][7][8][9][10][11][12] showed that Serbia harbours substantial genetic diversity of this species reflected in presence of two subspecies A. m. carnica and A. m. macedonica, their hybrids and three ecotypes.In this study, we tried to assess existing genetic variability based on 14 microsatellite loci and PCR-RFLP analysis of mtDNA COI gene region in managed honeybees in Serbia.
Based on the results of PCR-RFLP analysis, the presence of A. m. macedonica and its specific mtDNA lineage couldn't be confirmed due to the absence of a distinct restriction pattern previously described for this subspecies [5,7] suggesting that all sampled individuals belong to A. m. carnica subspecies.
Parameters of genetic diversity are relatively high and in range with previously published data and although significant FST values were obtained between some pairs of localities there is no clear pattern that indicates a geographical distribution of microsatellite loci since statistical significance was observed for pairs of south populations, north populations and north/south populations.Together with equal distribution of all three clusters inferred by STRUCTURE analysis in all localities and cluster overlapping inferred by DAPC analysis, presented results may indicate population admixture.However, results of AMOVA analysis suggest that grouping according to region, although not statistically significant, still may be observed.
Our results are partially in contrast with previously published data due to the absence of A. m. macedonica and still observable but very weak north/south differences between individuals from different parts of the country.The formation of hybrid populations may be one of the reasons behind the presented results.Extensive hybridization between A. m. carnica and A. m. macedonica subspecies was previously described in [6,[8][9] and it is possible that the hybridization zone expanded reflecting the close connections among beekeepers in the form of intensive queen trade and introductions of queens of different origins.The implementation of Serbian legislation that only A. m. carnica subspecies can be present in apiaries [13], which is strongly encouraged by a leading beekeeper organization as well, could be the reason behind the observed loss of A. m. macedonica and weak genetic differentiation between localities.Although A. m. carnica characterization by beekeepers is mainly morphological (our observation from the terrain concluded after the conversation with the beekeepers); it is reasonable to assume that beekeepers hold to this practice and this subspecies become predominant in apiaries.Overlapping of individual bees from South and North localities can also be partially contributed to the intensification of migratory beekeeping since apiaries from the South are transported to the North during the flowering season of agricultural plants.This can lead to further admixture once the migratory apiaries are returned in the region where stationary apiaries (sampled in this study) are located since there can be no control of mating between the queens and drones from different apiaries.AMOVA results are an indication that some weak regional differences between north and south localities still exist, but that difference is not statistically significant.
After a thorough examination of our results, we can conclude that the exclusion of autochthonous A. m. macedonica subspecies from the legislation is questionable since its implementation may have caused a disappearance of this autochthonous subspecies from many localities in which it was previously reported.Maintaining the high biological and genetic diversity of pollinators is essential for the wellbeing of ecosystems and is recognized as an important investment for the future of food production [1].Thus the current Serbian legislation needs to be updated and it needs to include other autochthonous honey bee subspecies and other bee species as well.In addition, the legislation has to be updated in the terms of honey bees management which will clearly specify the periods when the migrations of apiaries can be performed and in what regions in order to prevent the admixture of different honey bee populations that are adapted to the local environmental conditions.We further recognize the need to limit the distribution of commercially bred queens from different regions and their importation into the regions with different environmental conditions from the ones in which they are bred.In this way, the local genetic variability can be preserved and the admixture and uniformisation can be further prevented.The conservation of the biological diversity of honey bees is already recognized by the Serbian legislation [13] although some key aspects for the conservation of genetic diversity are lacking because of which the implementation of this legislation may have caused the observed genetic uniformization and loss of distinctive local populations.
However, our samples were not yet compared with reference and other honey be populations in the neighbouring countries.This comparison is under work and may yet reveal different positioning of investigated localities compared to others, especially when they are grouped according to the geographical region rather than locality.Nevertheless, present results indicate that change in beekeeping practice significantly influences genetic variability of managed honey bees which may result in population admixture and loss of distinctive local populations.Beekeeping is a dynamic process under constant change that requires a continuous investigation of the biological and genetic diversity of honey bees, an ecologically and economically important species.Only by constant monitoring of the honey bee population the proper measures for their protection and conservation can be implemented on time, thus preventing further losses of genetic and biological diversity which can also help to further improve the beekeeping practice in general.
Copyright: © 2021 by the authors.Submitted for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses /by/4.0/).

Figure 2 .
Figure 2. Matrixes of the average number of pairwise, Nei's (a) and FST (b) distances based on the analysis of microsatellite loci.a) The average number of pairwise differences between populations is presented above diagonal, the average number of pairwise differences within the population is presented diagonal and Nei's distances are presented below diagonal.b) Statistically significant FST values are marked with an asterisk (*).

Figure 3a
Figure 3a Discriminant analysis of principal components in which LDA was performed on the Figure 35.The first and second linear discriminants are presented in the plot.Figure 3b Non-metric multidimensional scaling plot of FST distances between localities in Se.rbia E. and S.D.; data curation, M.T. and S.D.; writing-original draft preparation, M.T. and S.D.; writing-review and editing, M.T., A.P., K.E., P.E., Lj.S. and S.D.; visualization, M.T. and S.D.; supervision, S.D. and Lj.S.; project administration, S.D.; funding acquisition, S.D.All authors have read and agreed to the published version of the manuscript."Funding: This research was funded by the Science Fund of the Republic of Serbia, PROMIS, #GRANT No 6066205, SERBHIWE for M.T., A.P., P.E., K.E., LJ.S., S.D. and Ministry of Education, Science and Technological Development of the Republic of Serbia, grant number 451-03-9/2021-14/ 200007 for M.T., P.E.A.P., K.E. and S.D.Informed Consent Statement: Not applicableData Availability Statement: Data are available on request from the corresponding authorConflicts of Interest: "The authors declare no conflict of interest.""The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results".
N -number of genotyped individuals, A -mean number of alleles, Ar -mean allelic range, Agd -average gene diversity over loci, Ho -observed heterozygosity, HE -expected heterozygosity G-W Garza-Williamson index, MPD -mean number of pairwise differences, RMP -random match probabili.ty

Table 2 .
Outcomes of AMOVA when localities are grouped according to their geographical region.
d.f.-degrees of freedom, SS-sum of squares, p-statistical significance