Winter activity unrelated to introgression in British bumblebee Bombus terrestris audax

Bombus terrestris is a bumblebee with a wide geographic range, with subspecies showing a variety of local adaptations. Global export of commercially-reared B. terrestris started in the 1980s; the bees are a mixture of subspecies bred for ease of rearing, bivoltinism and large nests. This paper investigated whether the increase in bivoltinism in UK resident B. terrestris audax populations was related to introgression with imported foreign subspecies. Workers were collected from wild populations in London and Bristol, as well as two commercial suppliers. Fourteen microsatellite loci were used to study population structure, hybridisation and introgression. No introgression with commercial B. t. dalmatinus was detected in wild populations. Hence, the increase in winter activity appears unrelated to introgression.


INTRODUCTION
Bombus terrestris (Linnaeus 1758) is a large, highly abundant bumblebee with a wide distribution across Europe (Figure 1), where it can be found in a range of habitats, excluding arctic, alpine and desert regions (Rasmont et al. 2008). Nine distinct subspecies have been described, with physiological and behavioural adaptations appropriate to their local environment, such as differences in colour pattern, pheromones and aggression (Rasmont et al. 2008;Coppée 2010). Despite such variation, the vast majority of B. terrestris populations are univoltine, with queens emerging from diapause in the spring, producing one colony over the summer, and with new queens initiating diapause at the end of summer (Løken 1973;Rasmont et al. 2008). A few exceptions do exist; B. terrestris exhibits facultative bivoltinism in areas where climatic conditions can support colonies over winter, such as the Mediterranean Basin (Beekman et al. 1999;Rasmont et al. 2008). In warmer Mediterranean regions, B. terrestris queens aestivate in response to the summer heat (Gurel et al. 2008).
Chosen for their ability to pollinate a wide range of plants and to thrive in a variety of environmental conditions, B. terrestris has been bred commercially since the 1980s for agricultural pollination services (Velthuis and Van Doorn 2006;Rasmont et al. 2008). Through a behaviour known as buzz pollination, bumblebees are capable of pollinating several crops (e.g. aubergine, tomato and blueberry) which cannot be effectively pollinated by honeybees, and are thus a cheaper alternative to mechanical pollination in which humans use mechanical, hand-held devices to manually pollinate such crops (Velthuis and Van Doorn 2006). Commercial breeding practices initially relied on the collection of tens of thousands of field-caught queens (Velthuis and Van Doorn 2006). Although fewer field-caught queens are caught today, the numbers are still estimated to be in the hundreds (Velthuis and Van Doorn 2006). Initial stocks are thought to have been built from a mixture of central European B. t. terrestris , eastern European B. t. dalmatinus , and Iberian B. t. lusitanicus . Over time, B. t. dalmatinus has become the market favourite subspecies due to its ease of rearing and development of large colonies (Velthuis and Van Doorn 2006). As well as colony size and domesticity, commercial populations were also bred to discourage diapause, leading to colony initiation immediately after mating, which has a known genetic component (Beekman et al. 1999;Velthuis and Van Doorn 2006). These mixed, domesticated colonies have been distributed worldwide, often with no regard for the impact on local, native bumblebee species. For example, escaped bumblebees have established feral, invasive B. terrestris populations in New Zealand, Chile and Japan (Donovan and Wier 1978;Ruz and Herrera 2001;Inoue et al. 2008;Rasmont et al. 2008;Aizen et al. 2018). Colony import has also been linked to the introduction and spread of parasites to wild bee populations (Meeus et al. 2011;Graystock et al. 2013;Graystock et al. 2016;Chandler et al. 2019). In the UK, regulations controlling B. terrestris subspecies import were only implemented in 2015, but non-native subspecies may still be imported with a licence (England. Natural England 2015). These domesticated bumblebee species could cause negative effects on local species and subspecies within their native range such as spreading disease, disrupting local genetic adaptation or outcompeting local populations (Ings et Rasmont et al. (2008) and Lecocq et al. (2016). Hybrid zones are indicated by striped regions. Note that one subspecies, B. t. canariensis , is not indicated on the map, but can be found solely in the Canary Islands. The blank map of Europe was created by Lvcvlvs and edited under the Creative Commons BY-SA 3.0 license.
A. F. Hart et al. Bombus terrestris audax is the resident subspecies of the British Isles, and like most B. terrestris populations, has historically been considered univoltine (Alford 1969). However, records of winter active colonies can be found as early as 1991, when queens would typically undergo diapause, especially in the South and Midlands of England (Robertson 1991;Stelzer et al. 2010). Reports of winter active bumblebees have been increasing in frequency, with a number of proposed mechanisms (Stelzer et al. 2010). Given the dependence of diapause induction and termination on temperature cues, the warmer winters associated with climate change (and an abundance of introduced winterflowering plants in urban areas) are thought to be facilitating bivoltinism (Alford 1969;Stelzer et al. 2010). Owen et al. (2013) suggested that imported B. t. dalmatinus , which are more fecund and more inclined to bivoltinism, may be displacing native B. t. audax (Ings et al. 2006). Additionally, Kraus et al. (2011) found that in Poland, commercial colonies frequently escape containment and interbreed with local populations. In Spain and Portugal, Seabra et al. (2018) and Cejas et al. (2020) both found evidence of commercial B. terrestris introgression with wild populations of local subspecies B. t. lusitanicus . Assuming the same is true for the UK, introduction of more bivoltine-inclined bumblebees may encourage winter activity in hybrid populations.
In this paper, we aim to investigate whether hybridisation and consequent introgression with imported, non-native subspecies is underpinning the recent trend towards bivoltinism in British B. terrestris , by using microsatellite analysis to detect differences and similarities between wild winter and summer populations, and commercial, imported bumblebees. To do so, wild B. t. audax workers were collected from London and Bristol across consecutive summers and winters (2009)(2010)(2011) and were genotyped with 14 microsatellite loci. Genetic profiles of wild bumblebees were compared with mass-reared workers of B. t. audax and B. t. dalmatinus .
Several different scenarios may be detectable through this method. If introgression is absent, bivoltinism will presumably be the consequence of behavioural plasticity. If this is true, then the summer and winter populations would look genetically similar, wild populations would be strongly differentiated to commercial populations, and no signs of introgression would be present. If introgression is present but not linked to winter activity, then summer and winter populations would not be distinct, and both should show medium to high similarity to commercial stocks. Finally, if winter activity is induced by introgression, then we would expect to see one bivoltine population present in both summer and winter, and a distinct, un-introgressed summer population, with the winter population showing higher similarity to commercial stocks.

Specimen collection
Between 2009 and 2011, 317 presumably wild B. terrestris audax workers were collected during summer and winter from London and Bristol areas ( Figure 2). Other B. terrestris subspecies are not expected to be present at these urban locations, as they are associated with agricultural use of colonies (Kraus et al. 2011;Moreira et al. 2015;Bartomeus et al. 2020). Individuals were later genotyped.
Samples were frozen whole and stored in individual vials at − 20°C. Two specimens were identified as B. lucorum by RFLP analysis, as workers of these two species cannot be easily distinguished (Murray et al. 2008 Table 1). In addition to the wild bees, individual workers were also sampled from different commercially-reared colonies: 29 B. t. dalmatinus from Syngenta (Weert, the Netherlands) and 14 B. t. audax from Agralan Ltd. (Swindon, UK). OpenStreetMap data were used to indicate sampling site locations (OpenStreetMap 2017).

DNA extraction and microsatellite protocol
DNA extractions were performed using one hind leg per sample, using the Promega Wizard Has introgression caused winter activity? SV 96 Genomic DNA Purification System, and automated using a Biomek Liquid Handling Platform. Quality control was performed using a Nanodrop spectrophotometer.

Data validation
Locus BTMS0106 was found to be monomorphic and excluded from further analysis. Individuals with data missing at three or more loci were also excluded. Colony version 2.0.6.5 (Wang 2004) was used to detect siblings via the combined full likelihood and pairwise likelihood score (FPLS) algorithm. Upon full sibling detection, the individual with the most missing data was removed; if equal, one was randomly chosen. After validation, 248 individuals with data for 13 loci remained (Supplementary Table 3).
The presence of null alleles was determined using Microchecker version 2.2.3 (Van Oosterhout et al. 2004) using the Bonferroni 1000× analysis option. GenAlEx 6.503 (Peakall and Smouse 2006) was used to check for deviations from the Hardy-Weinberg equilibrium, after which Bonferroni correction was manually applied. FSTAT 2.9.4 (Goudet 1995) was used to check for linkage disequilibrium between loci, which also underwent Bonferroni correction. Reported p values are post-Bonferroni correction.

Genetic diversity and population structure
GenAlEx 6.503 was also used to calculate expected and observed heterozygosity (H E and H O ), as well as allelic differentiation measure D EST , and fixation measures F ST and G ST (Jost 2008). Rarification and calculation of allelic richness (normalised to 10 genes) was performed with HP-Rare version June-6-2006 (Kalinowski 2004). Averages of genetic diversity (but not allelic differentiation or fixation) were weighted according to the number of individuals in each population.
BAPS 5.4 (Corander et al. 2003) was used to indirectly infer the degree of introgression using the 'population mixture analysis trained clustering' method. Using the two commercial populations as clusters of known origin, the algorithm assigns individuals of unknown origin (here, the wild samples) into clusters. Wild bees with a high degree of similarity to the commercial populations will be assigned in the same cluster. Fifty replicates of this algorithm were used to generate a mean and standard deviation to infer introgression, a method previously used for detecting introgression in B. terrestris by Kraus et al. 2011. STRUCTURE 2.3.4 (Pritchard et al. 2000) was used to find evidence of population structure by geography, season and origin (i.e. wild or commercial). K was run from 1 to 10, with 3 iterations for each with 250,000 MCMC burn-in reps and 750,000 after burn-in. STRUCTURE harvester 0.6.94 (Earl and vonHoldt 2012) was used to find the most likely value of K. CLUMPP mac version 1.1.2 (Jakobsson and Rosenberg 2007), and DISTRUCT version 1.1 (Rosenberg 2004) were used to visualise STRUCTURE results.

Data validation
Twenty-four individuals were removed from the dataset due to excessive missing data (Supplementary Table 3). Colony detected the presence of multiple full sibling pairs, leading to the exclusion of 12 bees from further analysis (Supplementary Table 3). In some cases, three siblings were collected from the same colony, thus two were removed. Microchecker detected three null alleles at allele frequencies above 10%, leading to the exclusion of three further loci (BT119, BTMS0067, BT132). FSTAT detected no significant (p = 0.0014) disequilibrium between loci over 78 pairwise comparisons when using Bonferroni correction for multiple comparisons. Five significant (p = 0.0007) deviations from the Hardy-Weinberg equilibrium were detected after adjusting for multiple comparisons, four in the summer 2010 population (BT06, BT09, BT20 and BTMS0014) and one in the summer 2011 population (BT30).

Genetic diversity
Overall, wild populations showed a marginally higher average allelic richness (A R ) and observed heterozygosity (H O ), but lower private allelic richness, than the commercial populations (Table I). Average A R for the commercial populations was 4.34, but 4.64 for wild populations. For private allelic richness, wild populations averaged 0.40, while commercial populations averaged 0.45. H O was very similar for the wild (0.680 ± 0.070) and commercial (0.660 ± 0.080) populations. Expected heterozygosity was also very similar for both wild (0.700 ± 0.130) and commercial (0.650 ± 0.080) populations.

Population differentiation by origin
Three different metrics of population differentiation were tested, namely fixation indexes F ST and G ST , and the allelic differentiation measure D EST , comparing the commercial B. t. audax , commercial B. t. dalmatinus and wild B. t. audax populations (Table II)  dalmatinus, p = 0.044 and 0.001). Not surprisingly, the strongest differentiation was observed between the two commercial populations, with an F ST of 0.037, a G ST of 0.018 and a D EST of 0.089 (p = 0.002 for each). STRUCTURE and STRUCTURE harvester were also used to study and visualise population differences (Figure 3). STRUCTURE harvester suggested that the most likely number of populations (K) when separated by origin was 4. However, this result is probably a consequence of a limitation of the program, which cannot detect the value of K if it is equal to 1 (Earl and vonHoldt 2012). Visualisation of the structure results indicates the approximately equal contribution of each population (K) in Figure 3, implying that the true value of K is indeed 1, with no structuring over origin.

Population differentiation by season
F ST , G ST , and D EST were also tested between commercial populations and wild populations sorted by season (Table III). F ST -values between seasons were very low, with an average 0.010, Figure 3. Population structure results obtained in STRUCTURE and visualised using CLUMPP and DISTRUCT. Each plot shows the proportion of the genome, for each individual, likely belonging to each cluster (K ). When data was analysed by season, the results of K = 2 are seen in A, and K = 5 in B. When analysed by origin, K = 2 in C, and K = 4 in D. Labels of 'Winter' or 'Summer' are referring to wild populations.
Has introgression caused winter activity? with varying levels of significance (Table III). Although the average genetic difference between wild B. terrestris grouped by season and both commercial populations was slightly higher (average F ST of 0.030), these pairwise measurements showed a low but significant genetic differentiation (p < 0.05; Table III). D EST and G ST were in agreement with the observed pattern of F ST ; seasons were not differentiated, and medium (D EST ) or low (G ST ) but significant differentiation was present between season and both commercial B. t. audax and B. t. dalmatinus populations.
Here, STRUCTURE harvester suggested that the most likely number of populations (K ) when separated by season was 5 (Figure 3). However, the structure results indicated, by the approximately equal contribution of each population in all bees (see Figure 3), that K is 1. Similar to the results of separation by origin, no structuring was thus observed for season by STRUCTURE.

Introgression
BAPS was used to infer the presence of introgression by attempting to assign the 216 wild individuals to commercial B. t. dalmatinus , commercial B. t. audax or their own cluster. After 50 repeats, BAPS consistently assigned all wild individuals to the same cluster as the commercial B. t. audax , with identical results for all 50 repeats.

Seasonal population structure
Despite the known genetic component of bivoltinism (Beekman et al. 1999), no population structure was apparent between different seasons of wild bees, implying there is little genetic difference between summer and winter populations (if any), in contrast to the significant, low or medium differentiation between seasons and commercial populations. These results, including the STRUCTURE results, are in agreement with Moreira et al. (2015), who found that wild B. t. audax across the island of Great Britain were one population. Furthermore, the BAPS analysis did not assign wild B. t. audax of any season to the commercial B. t. dalmatinus cluster, which excludes the scenario in which winter activity is due to introgression, as the winter population showed no higher similarity to commercial stocks than the summer populations.

Introgression
The BAPS results showed no evidence of hybridisation or introgression between the commercial B. t. dalmatinus and wild B. t. audax . However, this does not necessarily dictate that introgression is not occurring between these two populations, as introgression with imported commercial bumblebees and local subspecies appears widespread across other parts of Europe (Kraus et al. 2011;Seabra et al. 2018;Bartomeus et al. 2020;Cejas et al. 2020). This may be a consequence of the limitations of microsatellites, as discussed below.

Explaining winter activity
In general, our results support the scenario in which winter activity is not due to a genetic similarity to a bivoltinism-inclined commercial population. Indeed, given that the winter bees are not a subset of summer bees, or a consequence of introgression, it seems likely that the winter active populations are the same wild populations exhibiting behavioural plasticity, taking advantage of the warmer weather and abundant urban winter resources (Stelzer et al. 2010). This is supported by the incidence of winter activity occurring mostly in the South of England and climate data showing that the southern regions are also more intensely experiencing climate change, with a severe North-South gradient in the temperature of the hottest days of summer and winter (Stelzer et al. 2010;Stainforth et al. 2013).
However, it could also be the case that the wild and commercial populations are extremely highly introgressed, and have entirely replaced the native wild B. t. audax populations. STRUCTURE and STRUCTURE harvester indicated that the wild and both commercial populations were actually one population, due to high genetic similarity. If extreme introgression was actually the case, we would not be able to tell from our data whether the presence of introgression was related to winter activity, as all populations were bivoltine.

Comparison of commercial and wild populations
Given that domesticated populations are known to be under strong selective pressure for commercially-valuable traits, the level of differentiation between commercial and wild populations is surprisingly low (Velthuis and Van Doorn 2006;Abadía-Cardoso et al. 2016). Levels of observed and expected heterozygosity were similar to those seen in both British and continental European populations of B. terrestris (Estoup et al. 1996;Moreira et al. 2015), but our own data showed lower allelic richness.
Although the commercial B. t. dalmatinus population showed higher genetic diversity than the commercial B. t. audax , and similar levels to wild audax populations, the level of differentiation was unexpectedly similar to the wild audax population, showing only marginal difference in F ST and G ST to the values of the commercial audax (Tables I; II). This may be attributable to a number of different causes; a 25-year period of importing B. t. dalmatinus into the UK for commercial reasons is likely to have established some level of introgression, a pattern seen elsewhere in Europe (Kraus et al. 2011;Seabra et al. 2018;Bartomeus et al. 2020;Cejas et al. 2020). The 1980s practices of crossing different subspecies may mean that commercial stocks are hybridised, and even populations labelled as B. t. dalmatinus may contain genetic material from other subspecies such as B. t. audax , contributing to the lack of differentiation (Velthuis and Van Doorn 2006). B. t. audax was also domesticated much more recently than B. t. dalmatinus , and thus the commercial stock has had a shorter time period to diverge from its parent population (Velthuis and Van Doorn 2006). Finally, the intermittent addition of field-caught B. t. audax queens, which may be hybrid, to commercial B. t. audax stocks could also have potentially contributed to the lack of differentiation between commercial B. t. dalmatinus and wild population (Velthuis and Van Doorn 2006).

Alternative approaches to microsatellites
Microsatellites are a useful and dependable tool for population genetics studies, with the ability to reliably resolve population structure, identify inbreeding and genetic diversity parameters (Maebe et al. 2016(Maebe et al. , 2018(Maebe et al. , 2019a and have previously been used to find evidence of introgression in B. terrestris populations (Kraus et al. 2011;Lemopoulos et al. 2019). However, no technique is without limitation. Microsatellites may struggle to find siblings, or to identify the nuances of poorly differentiated populations (Lemopoulos et al. 2019).
In this study, a locus density of less than one microsatellite per chromosome may have limited our ability to uncover fine-scale clues in the dataset. Introgression or hybridisation can be hard to detect using microsatellites as large introgressed regions may lie undetected between sparse markers, which may help to explain why no evidence of introgression was detected through BAPS analysis (Jeffries et al. 2016). Additionally, the genes related to winter activity may not be linked with any of the microsatellite markers used. Such a gene must be in linkage (i.e. physically proximal) with one of the microsatellite loci in order to be detected; increasing the number of microsatellite loci used increases the resolution, thus increasing the likelihood that a gene is linked to a microsatellite locus.
A more thorough genomic surveying technique such as Restriction-site Associated DNA Sequencing (RADSeq), which can provide tens of thousands of informative SNP loci, may answer questions where microsatellites have failed to do so. In particular, RADSeq may help to understand why commercial and wild populations so closely resemble each other.

CONCLUSION
Microsatellite analysis of wild and commercial populations of B. terrestris has not revealed the presence of introgression; thus, introgression is unlikely to be the cause of the observed winter activity in wild B. terrestris populations in the UK. Furthermore, we found a lack of population structure and differentiation, which is consistent with previous literature. Further investigation is however necessary to determine if the import of foreign B. terrestris has had negative effects on the genetic diversity and fitness of native, wild populations.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.