Level of Genetic Diversity in European Bumblebees is Not Determined by Local Species Abundance

Bumblebee species with declining population trends tend to show lower genetic diversity levels than stable species. The observed difference might be explained by abundance differences, with declining bumblebee species having lower genetic diversity levels simply due to their lower local species abundances. However, whether this holds true is not known. Here, we investigated whether bumblebee local abundances determines population genetic diversity levels. Therefore, local species abundances were measured for bumblebee species at four locations in Belgium and two locations in Estonia during bumblebee foraging years 2013–2017. These locations and countries were chosen to ensure the greatest possible variance in both local abundances and population trends for these species. Hence, genetic diversity levels were obtained for seven species by genotyping collected specimens with 16 microsatellites. Our results showed that the observed patterns in genetic diversity did not depend on local species abundance. So, although declining bumblebee species can become locally abundant, they will still show lower genetic diversity levels than stable species. This result implies that measuring bumblebees’ local abundance cannot be used to directly determine the health status of a population. Furthermore, this result has also major impact on future conservation strategies as increasing the genetic diversity levels of declining species will be very difficult, and habitat quality should be high to maintain their populations, otherwise these species are doomed to disappear first.


INTRODUCTION
An essential factor for species survival is the level of genetic diversity present within its populations. Species which have a low level of intraspecific genetic diversity will have more limited potential to adapt to current and future changes in the environment (Frankham, 2005;Zayed, 2009;Habel et al., 2014;Koch et al., 2017). Small populations are more likely to have lower genetic diversity levels than large populations due to the impact of genetic drift, the random loss of genetic diversity over time (Reed and Frankham, 2003;Frankham, 2005;Zayed, 2009). However, high levels of gene flow and improved dispersal capacities can compensate for loss of allelic richness due to drift. In smaller less-connected populations, gene flow may be more limited, decreasing its buffering effect, which leads to lower genetic diversity and increase brother-sister mating chances, and in turn to inbreeding and inbreeding depression (Frankham, 2005;Zayed, 2009;Habel et al., 2014). The latter dynamics can further diminish genetic variability within small populations, potentially creating a vicious circle, known as the extinction vortex, which could ultimately lead to extinction (Frankham, 2005;Zayed, 2009;Habel et al., 2014).
Bumblebees are eusocial insects, living in large colonies mostly founded by one queen and many workers (Alford, 1975;Goulson, 2010). In most bumble bee species, queens are monoandrous with their offspring being fullsibs, having 75% or more genetic similarity between them (Estoup et al., 1995;Schmid-Hempel and Schmid-Hempel, 2000;Goulson, 2010). The high relatedness of sisters can be explained by bumblebees' haplodiploid sex-determination system in which unfertilized eggs, being hemizygous at the single sexdetermining locus, will develop in haploid males or drones, and at which fertilized eggs will develop in diploid female offspring (workers and daughter queens) when being heterozygous at the sex locus. However, fertilized eggs can also develop into diploid males if being homozygous (Duchateau et al., 1994;Cook and Crozier, 1995;Whitehorn et al., 2009). The probability of such matched-pair matings at the sex locus is normally low due to the rather high number of alleles at this loci (46 alleles in B. terrestris; Duchateau et al., 1994). However, the presence of diploid drones will cause negative effects on population growth and survival, particularly in populations with low variation at the sex locus, which are typically small and/or inbred populations (Cook and Crozier, 1995;Duchateau and Marien, 1995;Gerloff and Schmid-Hempel, 2005;Whitehorn et al., 2009).
Multiple studies have shown that bumblebee species with declining population trends tend to have lower genetic diversity levels compared with stable bumblebee species Ellis et al., 2006;Goulson et al., 2008;Charman et al., 2010;Cameron et al., 2011). This observation has been often explained as a reduction of genetic diversity over time (Goulson et al., 2008;Charman et al., 2010) due to the impact of one or multiple possible drivers of bee decline, e.g. agricultural intensification, new pathogens, and climate changes, lowering population sizes (Potts et al., 2010;Meeus et al., 2011;Vanbergen et al., 2013;Rasmont et al., 2015). However, studies using historical populations found a similar difference in genetic variation between declining and stable bumblebee species (Lozier and Cameron, 2009;Lozier et al., 2011;Maebe et al., 2015;Maebe et al., 2016), detecting no major drop in genetic diversity over one century in Belgium, possibly due to dispersal countering drift effects (Maebe et al., 2016). The latter result indicates that, at least for Belgian populations, species abundance, and genetic diversity are not linked. This is in conflict with the rather general theory, that when a certain species is locally abundant, its large population should have a higher amount of genetic variation than observed at a location where this species is less abundantly present (Reed and Frankham, 2003;Frankham, 2005;Zayed, 2009). However, for social insects, such as bumblebees, this might thus not be always true.
Here, we investigated if species abundance determines genetic diversity levels in eusocial bumblebee species, and this on both local and larger scale. Hence, we searched whether genetic diversity is linked with species IUCN status (or European pattern of occurrence), and compared this with local abundance measures. Therefore, specimens of seven bumblebee species (B. ruderarius, B. soroeensis, B. sylvarum, B. hortorum, B. hypnorum, B. lapidarius, and B. pascuorum) were collected from two countries: Belgium and Estonia. Three species have restricted distribution and declining population trends in Belgium (B. ruderarius, B. soroeensis, and B. sylvarum), while the remaining four species have a more widespread distribution and showed fairly stable population trends (B. hortorum, B. hypnorum, B. lapidarius, and B. pascuorum) (division based on Maebe et al., 2016 and refs herein). Although the latter bumblebee species are also dominantly present in Estonia, these species occur in different abundances, while the "declining" species are more abundantly present and even showing increasing population trends ( Table 1). Bumblebee workers were collected from six locations, four in Belgium (Francorchamps, Moorsel, Torgny, and Trivières) and two in Estonia (Harjumaa and Põlvamaa) during bumblebee foraging years from 2013 to 2017. For each species, genetic diversity was estimated from each sampling location using 16 microsatellite markers. Furthermore, species local abundance was determined as its relative abundance (= its present abundance compared to the total species abundance presented at that location). This approach allowed us to investigate the link between bumblebee local abundance and genetic diversity.

Sampling and Proportional Abundance Measurement
A sampling effort for specimens from seven bumblebee species (B. hortorum, B. hypnorum, B. pascuorum, B. ruderarius, B. soroeensis, B. sylvarum, and B. hypnorum) was performed at four sampling locations in Belgium and two in Estonia during bumblebee foraging seasons from 2013 to 2017. Selected sampling locations within both countries were around 56-210 km apart, while the distance between both countries ranged 1,650 km ( Figure 1). Specimens were collected from each location for three species, B. hortorum, B. hypnorum, and B. pascuorum. Unfortunately, the other species were only present at three out of seven sampling locations (one Estonian and two Belgian locations for B. hypnorum, one Belgian and two Estonian locations for B. ruderarius, B. soroeensis, and B. sylvarum; Table 2). 392 bumblebee specimens were collected from the four selected Belgian locations: Francorchamps, Moorsel, Torgny, and Trivières, and were already genotyped at 16 microsatellite loci within a previous study (Maebe et al., 2016; Table 2). These specimens were part of an intensive sampling performed within a national bumblebee project (Belbees) during 3 days in the bumblebee foraging season of 2013 and 5 days in 2015. At pollinator-suitable weather conditions, all encountered bumblebees were sampled during straight ahead transect walks using a net. Species abundances were calculated as the proportion of a species abundance versus the total bumblebee abundance collected at that particular location ( Table 3).
In Estonia, Harjumaa and Põlvamaa were selected as sampling locations ( Table 2). Sampling was performed during the bumblebee foraging seasons of 2015 and 2017. Individual specimens were collected with small glass jars from a flower, killed with chloroform, and stored in a −20°C freezer. The Estonian bumblebee abundance data were derived from Marja et al. (2014), from which the North and South Estonia abundances were used as estimates for Harjumaa and Põlvamaa, respectively. Here, per species local abundances were also calculated ( Table 3).

DNA Extraction and Microsatellite Protocol
DNA extraction and analysis were performed only on the Estonian specimens, with the same method as was done in a previous study for the Belgian specimens (see Maebe et al., 2016). One middle leg of an individual bumblebee worker was used for DNA extraction using a Chelex DNA extraction protocol (see methods described in Maebe et al., 2015). Specimens were genotyped with 16 microsatellite (MS) loci which gave reliable signals in previous research using different bumblebee species (Supplementary Table 1; Maebe et al., 2015;Maebe et al., 2016;Maebe et al., 2018). MS were amplified using the Type-it QIAGEN PCR kit using four multiplex mixes as described in Maebe et al. (2016). The PCR protocol, capillary electrophoreses, and fragment scoring were performed with the method as described in Maebe et al. (2013).

Impact of Local Abundance on Level of Genetic Diversity
To be able to detect if species abundance influences the observed differences in species genetic diversity, Linear Mixed Models (LMMs) were performed for each species and both genetic diversity parameters (A R and H E ) in RStudio (R Development Core Team, 2008) with R package lme4 version 1.1-10 (Bates et al., 2015) as described in Maebe et al., 2018. In short: species, local abundance, and location were set as fixed factors and microsatellite loci as a random factor. The best model fitting our data was selected based on the Akaike's Information Criterion (AIC) by using the "dredge" command within the MUMIn package (Barton 2015;Maebe et al., 2016;Maebe et al., 2018). The main effect of the factor of interest was analyzed for each selected LMM by performing Likelihood Ratio Tests (LTR) in which the model with factors was compared with a "null" model without these factors (Maebe et al., 2018). After computing the marginal and conditional coefficient of determination, Tukey HSD post hoc comparisons were performed using the R package multcomp to find differences in genetic diversity due to species local abundance, location and/or species (Hothorn et al., 2008;described in Maebe et al., 2018, and references therein).

Genetic Diversity and IUCN Status
Based on the IUCN bumblebee red list status, bumblebee species were organized in two groups, either "stable" or "declining" species ( Table 1). The group of the "stable" species consisted out of four species (B. hortorum, B. hypnorum, B. lapidarius, and B. pascuorum) of which all species had an IUCN status of "Least Concern" (LC), except for B. hortorum which have a "Near Threatened" (NT) status (Rasmont et al., 2015). The three remaining species, B. ruderarius, B. soroeensis, and B. sylvarum, were grouped as "declining" species, as they all have a threatened IUCN status ("Endangered" (EN), "Vulnerable" (VU), and "Critically Endangered" (CE), respectively; Rasmont et al., 2015; Table 1). Secondly, for each country, bumblebee species were divided in two groups based on their population trends, being "declining" or "not-declining" species. Such a population trend was made for each species based on the difference in relative abundance of that species between two time periods: for Belgium between 1910-1930and 1990-2016(Vray et al., 2019; Table 3), while for Estonia between 1955-1967-2018(Kotkas, 1968Viik, 2019; Table 3). Genetic diversity parameters (A R and H E ) were compared between the different groups by LMMs in RStudio (R Development Core Team, 2008) with R package lme4 version 1.1-10 ( Bates et al., 2015). Here, models were run with microsatellite "loci" as random factor, and "group" as fixed factor. Tukey post hoc tests were performed as described above.

RESULTS
All 16 microsatellites could be amplified and scored reliably in the seven Bombus species. In total, 628 specimens (309 Belgian and 319 Estonian specimens) remained for further genetic analysis out of 729 specimens, due to the removal of full-sibs which were detected by Colony 2.0 and Kinalyzer analyses, and discarding of specimens which have more than five loci of missing data from both Belgian and Estonian dataset. Furthermore, no significant linkage disequilibrium between loci were detected, and Hardy-Weinberg equilibrium tests displayed no or only limited heterozygote deficits.

Correlation Between Local Abundance and Genetic Diversity
Overall for populations and species, allelic richness (A R ) ranged from 1.87 to 3.66, with a mean A R of 2.70. Mean H E was 0.488, with individual population values ranging from 0.269 to 0.739 ( Table 2). Linear regression showed low and non-significant correlations between local abundance and A R or H E (R 2 = 0.0261, p = 0.393; R 2 = 0.0419, p = 0.278; respectively; Figure 2) suggesting that only a low amount of variation of genetic diversity is explained by local abundance. By comparing the AIC scores of the different LMM models, the best model explaining the observed patterns of both genetic diversity parameters was not the model with "abundance", but with "species" as the only fixed factor ( Table 4). The importance of the factor "species" for both A R and H E was clearly shown by comparing the models with and without "species" as fixed factor (LRT, c 2 = 154.07, d.f. = 6, p < 0.001; and c 2 = 140.72, d.f. = 6, p < 0.001, respectively). However, addition of abundance as an extra fixed factor had no significant impact on the model, as became clear after comparing the models with both "species" and "abundance" as fixed factors to the models without abundance as additional fixed factor (LRT, A R and H E : c 2 = 0.1307, d.f. = 1 p = 0.718; and c 2 = 0.0598, d.f. = 1, p = 0.807, respectively). Marginal and conditional R 2 were similar for both A R and H E , 16.4% and 57.1% versus 15.2% and 55.7%, respectively (Nakagawa and Schielzeth, 2013). Post hoc tests were performed and showed significant differences in genetic diversity levels between the seven bumblebee species (Table 5 and Supplementary Tables 2 and 3).

Genetic Diversity Differences Between Species Groups?
Comparing both A R and H E between the two IUCN groups showed a significantly higher genetic diversity level in "stable" bumblebee species compared to the contemporary "declining" species (post hoc tests: A R , t = 5.874, p < 0.001, and H E , t = 6.226, p < 0.001; Table 6). However, genetic diversity was similar when comparing populations with "not-declining" compared with those with "declining" population trends (post hoc tests; A R , t = −0.489, p = 0.625, and H E , t = 0.172, p = 0.864; Table 6).

DISCUSSION
In this study, we investigated whether genetic diversity is determined by species local abundance. By genotyping specimens of seven bumblebee species from Belgium and Estonia with 16 microsatellites, we were able to estimate genetic diversity and to compare this with local species abundance. This approach allowed us to conclude that observed genetic diversity patterns were not significantly dependent on species local abundance. Indeed, the best models explaining genetic diversity patterns did not include local abundance as a factor, and both A R and H E were not significantly correlated with abundance ( Figure 2). This is rather unexpected as one might think that when a certain species is locally abundant, this large population should  Of all possible models run under MUMIn using location (= Loc.), species abundance (= Abn), species (= Spss), and time periods (= Per.) as fixed effects, and locus as a random effect, the best fitting linear mixed-effect models (with a delta <10) are given. Based on their high (negative or positive) Akaike's Information Criterion (AIC), the best models for A R (A) and H E (B) with + = parameters included in the model, and NA, not included parameters. Selected best models for AR and HE are given in bold.
normally have a higher amount of genetic variation than may be observed at a location where this species is less abundantly present. However, this is not always true. Indeed, when the local environment has enough food, resources, and is well suited to harbor a large population of bumblebees, then a bumblebee population can become locally large without increasing its genetic diversity levels. This may be partly explained by bumblebees' haplodiploid sex-determination system in which unfertilized eggs develop into haploid males and fertilized eggs into diploid females, where workers do not mate and normally do not contribute to the next generation, and where each colony consists of a single mated queen and her offspring (Alford, 1975;Goulson, 2010). Such colonies can have large numbers of specimens which are in fact all "clones" with 75% or more genetic similarity between them (Goulson, 2010). A population consisting of large colonies founded by closely related queens will harbor many closely related workers and thus can have large population numbers but lower genetic variation than one might expect. Furthermore, several studies have shown that bumblebees are able to recover from genetic bottlenecks by increasing population sizes without increasing their genetic diversity levels (Schmid-Hempel et al., 2007). Although these populations are characterized with having many specimens, and thus are at first glance in a healthy condition, these populations are under severe fitness risks, due to the low genetic variation which makes these populations vulnerable to additional future environmental changes (Frankham, 2005;Zayed, 2009;Habel et al., 2014;Koch et al., 2017). Indeed, populations with low genetic variation show less population growth and survival (Whitehorn et al., 2009), and are more susceptible to diseases (Whitehorn et al., 2011;Cameron et al., 2011;Whitehorn et al., 2014). Hence, when neighboring populations have rather similar levels of genetic diversity, then increasing a populations' genetic diversity through high amounts of gene flow will not be possible within a short time span, making these populations even at higher risk. The latter, implying that neighboring populations have similar genetic diversity cannot be proven in our study, as no in depth sampling of multiple locations was performed and much variation may thus be missed from unsampled populations. However, the similar genetic diversity levels we observed within the populations of a certain bumblebee species ( Table 2), suggest that genetic diversity levels in neighboring populations may be very similar. This in turn would mean that for conservational purposes, increasing genetic diversity levels may not so easily be achieved.
Here, grouping the bumblebee species by IUCN status into "declining" and "stable" ( Table 6) showed similar results as described by Maebe et al. (2016), confirming lower genetic diversity levels within the populations of declining species than observed within widespread, stable species (as been observed and/or discussed in: Darvill et al., 2006;Ellis et al., 2006;Goulson et al., 2008;Charman et al., 2010;Cameron et al., 2011).
*** NS *** *** NS *** − Impact of the factor lrdquo;species" in the model on A R and H E . Significant levels p < 0.001 = ***, p < 0.01 = **, p < 0.05 = *, and NS = not significant, with in bold and italic negative and positive interactions, respectively. With the estimate, standard error (SE) and p-value of each factor in the model. Significant factors are indicated in bold, with significant levels p < 0.001 = ***, p < 0.01 = **, and p < 0.05 = *. Furthermore, our results showed clear differentiation of genetic diversity based on species. Indeed, the best models explaining observed A R and H E included species as fixed factor (Table 4), and post hoc tests showed multiple significant differences in genetic diversity between species (Table 5). Those results assume that the maximum level of genetic diversity which a species could achieve might be always lower than in other species. The latter might be due to species-specific population characteristics such as lower queen production, smaller colonies, low dispersal capacity, etc., which all could limit the generation of higher genetic diversity levels within the populations of a certain species, and/or due to past population dynamics which have altered diversity levels within species' populations in the past (as hypothesized in Maebe et al., 2016). The hypothesis of past population dynamics influencing genetic variability has been raised attention as a recent study showed a clear drop in genetic diversity over time in south-Brazilian bumblebee species (Maebe et al., 2018) which contradicted the temporal stability observed within the European bumblebees (Maebe et al., 2016). The authors explained this difference due to deforestation differences between both continents, hypothesizing that in Europe major reductions in genetic diversity occurred earlier in time (Maebe et al., 2018). Whether this holds true must be further investigated, but the lower levels in genetic diversity observed in certain species might thus be due to different impact of land-use changes, such as deforestation, on particular bumblebee species. To conclude, our results showed that measuring bumblebees' local abundance cannot be used to directly determine the health status of a population, especially on a longer term. Furthermore, as bumblebee species with less genetic diversity have lower fitness (e.g. Cameron et al., 2011;Whitehorn et al., 2011;Whitehorn et al., 2014), it seems that for conservation purposes, measuring genetic diversity parameters can give a good prediction value whether a certain population or bumblebee species is more or less vulnerable towards future population decline.

AUTHOR CONTRIBUTIONS
KM and IM designed the research. RK provided specimens. KM analyzed the data. KM wrote the paper, while all co-authors (RK, IM, MM, and GS) helped in improving the manuscript.