Genomic divergence and a lack of recent introgression between commercial and wild bumblebees (Bombus terrestris)

Abstract The global movement of bees for agricultural pollination services can affect local pollinator populations via hybridization. When commercial bumblebees are of the same species but of different geographic origin, intraspecific hybridization may result in beneficial integration of new genetic variation, or alternatively may disrupt locally adapted gene complexes. However, neither the existence nor the extent of genomic introgression and evolutionary divergence between wild and commercial bumblebees is fully understood. We obtained whole‐genome sequencing data from wild and commercial Bombus terrestris collected from sites in Southern Sweden with and without long‐term use of commercially imported B. terrestris. We search for evidence of introgression, dispersal and genome‐wide differentiation in a comparative genomic analysis of wild and commercial bumblebees. Commercial B. terrestris were found in natural environments near sites where commercial bumblebees were used, as well as drifting wild B. terrestris in commercial bumblebee colonies. However, we found no evidence for widespread, recent genomic introgression of commercial B. terrestris into local wild conspecific populations. We found that wild B. terrestris had significantly higher nucleotide diversity (Nei's pi, π), while the number of segregating sites (Watterson's theta, θw) was higher in commercial B. terrestris. A highly divergent region on chromosome 11 was identified in commercial B. terrestris and found to be enriched with structural variants. The genes present in this region are involved in flight muscle contraction and structure and pathogen immune response, providing evidence for differing evolutionary processes operating in wild and commercial B. terrestris. We did not find evidence for recent introgression, suggesting that co‐occurring commercial B. terrestris have not disrupted evolutionary processes in wild B. terrestris populations.


| INTRODUC TI ON
Deliberate or inadvertent human-assisted movement of non-native species into new locations creates opportunities for novel ecological interactions (Bartz & Kowarik, 2019;Crispo et al., 2011;Keller et al., 2011). However, the movement of a species to a new region where a closely related species co-occurs may also result in genetic exchanges and altered evolutionary outcomes (Crispo et al., 2011;Kanbe et al., 2008;Yoon et al., 2009Yoon et al., , 2011. Genomic introgression, the exchange of genetic material between species or subspecies through hybridization and repeated backcrossing (McFarlane & Pemberton, 2019), can be either beneficial or detrimental to local populations. Beneficial effects include 'genetic rescue' whereby novel genetic variation is introduced into genetically depauperate populations, or through 'adaptive introgression' by creating novel selection pathways that may reduce extinction risk (McFarlane & Pemberton, 2019) or simply a general increase in genetic diversity to better adapt to changing environments (Nelson et al., 2017;Suni et al., 2017). Alternatively, detrimental effects include reduced genetic diversity and disruption of local co-adapted gene complexes, with novel alleles replacing locally adapted ones (Roesti, 2018;Schumer et al., 2018). Such processes may lead to population declines and even local extinctions (Jensen et al., 2005;Keller et al., 2014;Todesco et al., 2016). Introgression is therefore recognized as a process of significance for wildlife conservation (Laikre et al., 2010; McFarlane & Pemberton, 2019; Todesco et al., 2016).
There are many historical examples where exotic species have been deliberately or accidentally introduced beyond their natural range, leading to hybridization or introgression between closely related species (Biedrzycka et al., 2012;Chazara et al., 2010;Escalante et al., 2014). The use of exotic bee species for crop pollination and honey production has played a major role in the global movement of different bee species, repeatedly creating conditions for hybridization and introgression among native and non-native bees (Byatt et al., 2015;Goulson, 2010;Kraus et al., 2007), which in some cases, has resulted in inviable hybrids between species (Kanbe et al., 2008;Tsuchida et al., 2010). More recently, attention has turned toward the ecological and evolutionary consequences of commercially produced bumblebees for native pollinator health and genetic diversity, due to the industrial-scale global transport and use of commercial bumblebees for agricultural crop pollination (Dafni et al., 2010;Goulson, 2010;Velthuis & Van Doorn, 2006).
Bumblebees have several characteristics that make them highly effective agricultural pollinators, such as varied tongue lengths, tolerance to a wide range of weather conditions, and large hairy bodies providing ample surface area for pollen attachment. In particular, their ability to buzz-pollinate by shaking floral anthers with high frequency to release pollen (Nayak et al., 2020;Nielsen et al., 2017;Wahengbam et al., 2019) makes bumblebees very effective pollinators for some crops (e.g., tomatoes : Cooley & Vallejo-Marín, 2021) and therefore difficult to replace in agricultural settings. Commercial and introduced bumblebees, however, may be detrimental to local wild pollinators. For example, they are known to compete for resources with local native bees (Hingston & McQuillan, 1998;Ings et al., 2005Ings et al., , 2006; Morales et al., 2013) and spread pathogens (Dafni et al., 2010;Evans, 2017;Meeus et al., 2011) (but see Trillo et al., 2021), thereby disrupting native plant-pollinator relationships (Aizen et al., 2019;Hingston & McQuillan, 1998). However, the implications of commercial and introduced bumblebees for the genetic integrity and evolutionary trajectory of native pollinators is little understood (Seabra et al., 2019).
Since the commercialization of bumblebees in the late 1980s for the purpose of tomato pollination, the industry has grown rapidly (Owen, 2018). By 2006, it was estimated that more than one million reared colonies had been transported for introduction across the world (Velthuis & Van Doorn, 2006). While the subspecies B. terrestris terrestris was initially used for domestication (Velthuis & Van Doorn, 2006), several of the other nine B. terrestris subspecies have since been domesticated (Rasmont et al., 2008) and introduced outside of their native ranges (Inari et al., 2005;Kraus et al., 2011;Schmid-Hempel et al., 2014). These nine subspecies are assumed to be adapted to the environmental conditions of their native geographic ranges and differ in phenology, colour patterning, behavioural traits and parasite resistance (Rasmont et al., 2008). Since only a few countries (e.g., Canada, United States and Japan) have implemented trade regulations for the importation of commercial B. terrestris (or for certain subspecies; Lecocq et al., 2016;Moreira et al., 2015;Velthuis & Van Doorn, 2006), introgressive hybridization among nonlocal or non-native species or subspecies poses a threat to the genetic integrity of native or local bumblebee populations (Ings et al., 2006;Kanbe et al., 2008;Tsuchida et al., 2010;Yoon et al., 2009Yoon et al., , 2011. Several studies have found evidence of introgression between native and commercial subspecies of bumblebees in the Iberian Peninsula, and Western and Eastern Europe (Bartomeus et al., 2020;Cejas et al., 2018Cejas et al., , 2020Moreira et al., 2015;Seabra et al., 2019). In other cases, introgression between commercial and native subspecies of bumblebees has not been detected, such as in the UK (Hart et al., 2021) and New England, USA (Suni et al., 2017).
Genetic differentiation between wild and commercial bumblebees is of further importance to investigate since different selection pressures may arise from the process of domesticating commercial bumblebees. Although there are no artificial selective breeding programs specifically for bumblebees (Lecocq, 2019), and genetic divergence between commercial and wild bumblebees is not necessarily expected, only the fact that they are raised in an artificial environment may result in divergence. Accordingly, studies have found significant pairwise genetic differentiation (F ST ) between commercial and wild-caught Bombus impatiens and B. terrestris within New England, USA, and within the Iberian Peninsula (Seabra et al., 2019;Suni et al., 2017). Although the drivers of this differentiation have not been thoroughly investigated, this suggests that demographic processes and natural selection may be operating differently on commercial and wild bumblebees.
Thus, the identification of introgression among differentially selected bumblebees could result in novel alleles replacing locally adapted ones, interfering with processes of local adaptation. To our knowledge, no studies have investigated whether there is introgression between commercial and wild B. terrestris in Northern Europe. Additionally, to date, studies investigating introgression in bumblebees have used either microsatellites or reduced representation approaches (i.e., Restriction-site Associated DNA sequencing-'RADseq' or Genotype By Sequencing-'GBS') and have therefore been limited in their ability to detect introgressive hybridization and its role in selection processes in wild and commercial bumblebees (Bartomeus et al., 2020;Cejas et al., 2018Cejas et al., , 2020Moreira et al., 2015;Seabra et al., 2019;Suni et al., 2017).
Here, we use whole-genome sequencing (WGS) data, which represents an opportunity to examine introgression, selection and evolutionary divergence at high resolution, using high-density genetic markers across the entire genome, from commercial (most likely subspecies B. t. terrestris and/or B. t. dalmatinus;Goulson, 2010) and wild-caught B. terrestris (dominant subspecies B.t. terrestris;Rasmont et al., 2008) sampled across the southernmost region of Sweden. This area is dominated by agriculture and includes a widespread but localized long-standing use of commercial bumblebees for pollination, for example fruit, berries and tomatoes. We ask: (1) Is there evidence for genomic introgression between wild and commercial populations of B. terrestris? And (2) do genome-wide selection signatures differ between commercial and native bumblebees? Our study has implications for commercial practices involving the distribution and containment of bumblebees, understanding adaptive genetic variation in commercial and native bumblebees and, therefore, monitoring the evolutionary resilience of pollinators in a rapidly changing world.

| Experimental design
Ideally, introgression should be investigated before and after the introduction of commercial bumblebees in a landscape, which however is logistically infeasible as detectable signals are only expected to develop over time. Instead, we used an experimental landscape design, which assumes that contemporary spatial ecological patterns are roughly equivalent to changes over time (cf. Pickett, 1989), to determine the occurrence and extent of genomic introgression of commercial B. terrestris genetic material into wild B. terrestris genomes.
This was done by comparing the genomes of commercial B. terrestris with those of (wild) B. terrestris caught at sites with long-term use of commercial B. terrestris used for pollination services (hereafter 'experimental sites') with those of wild B. terrestris caught at sites with no history of commercial bumblebee use (hereafter 'control sites'). Control sites were at least 1500 m from the nearest commercial bumblebee colony. Commercial B. terrestris were collected at the experimental sites and sites located in our general study area. To account for potential environmental variability between the experimental and control sites, we chose control sites that were located in areas with similar high agricultural land cover as the experimental sites. All sampling sites were located on the same latitude, and hence, temperature and precipitation did not vary markedly between sites (Table S1; data from WorldClim v2.1, Fick & Hijmans, 2017).

| Wild-caught bumblebees
Seventy-eight free-flying 'wild' female diploid worker B. terrestris (WB) were collected across ten sites in Southern Sweden during July 2018 (nine sites) and July 2019 (one site) ( Figure 1; Table 1). At five 'experimental' sites, commercially reared bumblebees had been used for 12-27 years (experimental sites) for pollination in either greenhouse cultivation of tomatoes (three sites), tunnel cultivation of berries (one site) or a combination of greenhouse (tomato), tunnel (berries) and free land (berries and apples) cultivation (one site), at a stocking rate of 50 to 350 colonies/year. At these experimental sites, wild free-flying bumblebees (hereafter, called wild experimental 'WE') were collected between 700-1000 m from the closest greenhouse, tunnel and/or free land with commercial bumblebee colonies. At the sixth experimental site, commercially reared bumblebees have been used since 2013 (except during the years of 2015 and 2017) at a stocking rate of 12 colonies/year. At this site, wild free-flying bumblebees were collected at 3000 m from the closest free land with commercial bumblebee colonies (hereafter, called wild experimental 'WE'). The site was added later to increase the total sample size, but only three individual samples from this site were included in the dataset used for analyses due to quality control of the dataset and identification of full siblings. The selection of sampling locations for all experimental sites allows for the detection of introgression while minimizing the incidence of commercial 'escapees' since the mean foraging distance for B. terrestris workers are <300 m from the colony (Wolf & Moritz, 2008), although foraging trips may be longer than 700 m (Goulson, 2010).
Free-flying B. terrestris workers (hereafter, called wild control 'WC') were also collected at four 'control sites' at least 15-20 km from the nearest location of commercial bumblebee use (data collection summarized in Table 1). These sites were selected to maximize the distance between where wild 'control' B. terrestris were sampled, and the location of the closest commercial bumblebee colony(ies).
Out of necessity, the sites were relatively close to each other since it proved difficult to find locations in Southern Sweden where we could safely assume that commercial bumblebees have not been used or are currently being used. We ensured the absence of use of commercial bees by avoiding known locations using bumblebees and by collecting customer postcodes of sales records from the major bumblebee importers to avoid these areas.
The locations of commercial bumblebees (experimental sites) were identified by contacting growers producing bumblebeedependent crops in greenhouses or tunnels; the names of enterprises and their exact location remain confidential by agreement.
We collected all encountered foraging or flying wild bumblebees using sweep nets and brought back specimens to the laboratory in a cooling box, preserved in 70% ethanol and stored at −20°C.

| Commercial bumblebees
We collected 18 commercial bumblebees (hereafter 'CB') directly from 18 different colonies, one bee per colony (Biobest Group NV, Belgium) ( Table 1). Four of these colonies came from the greenhouse, tunnel and free-land cultivations where wild experimental (WE) samples were also collected. The remaining 14 colonies were from the same general study area as our wild experimental sites (WE) and at least 15 km from our wild control (WC) sampling sites ( Figure 1). All commercial bumblebees originated from the same supplier. We preserved all bumblebees in 70% ethanol and stored them at −20°C.

| DNA extraction, barcoding and sequencing
Genomic DNA was extracted from the head and two legs of each bumblebee using a Qiagen Blood & Tissue Extraction Kit (QIAGEN GmbH) following a modified version of the manufacturers Supplementary Protocol (Text S1). To confirm species identification of the free-flying bumblebees, the COI mitochondrial gene was amplified in all samples according to Wahlberg and Wheat (2008) (Text S2). COI sequences were compared with both complete and partial mitochondrial genome sequences from B. terrestris in GenBank F I G U R E 1 Sampling locations for all free-flying wild collected bumblebees and locations of the commercial hives in Southern Sweden. Free-flying female workers collected at experimental sites (WE) in 2018 and 2019, where commercial bumblebees have been used for 12-27 years are shown as dark purple triangles (site from 2019 has an '*' next to the symbol) and were collected 700-3000 m from the closest location with commercial bumblebees. Free-flying female workers collected at control sites (WC) at least 15 km from the nearest location with commercial bumblebees are shown as dark purple diamonds. The locations of the commercial hives (CB) are shown as green circles (commercial, n = 18). Six of the zoomed-in circles illustrate the 'experimental sites' and the three different agricultural practices: greenhouse; open tunnel-cultivation and free land, which use commercial bumblebees for pollination services. The final zoomed-in circle illustrates the 'control sites'

27-07-2018
Note: The name of the sampling sites and the name of the three sampling groups are illustrated in the table. Indication if commercial colonies are present in the agricultural practices located in our sampling sites, the cultivation type, the number of commercial bumblebee (CB) colonies used on a yearly basis, the number of years with CB colony usage for the agricultural practices used in the study, the number of free-flying wild bumblebees (WBs) caught at each site and the distance in meters from each sampling site to the nearest CB colony is also illustrated in the table.
We collected the free-flying and commercial bumblebees between 17/07/2018 and 03/08/2018. Abbreviations: WE, wild experimental bumblebees (n = 57); CB, commercial bumblebees (n = 18) and WC, wild control bumblebees (n = 21). NA*, the exact number of CB colonies used per year is not known but the greenhouse tomato cultivation is of large scale and has been operating since 2000.  (Purcell et al., 2007). Analysis was performed on the final set of 652,002 SNPs from the 18 assembled B. terrestris chromosomes (excluding scaffold sequences).

| Genetic structure and introgression
We first performed a principal component analysis (PCA) on the whole dataset using the linkage-pruned SNP dataset. Resulting eigenvectors were plotted in R using ggplot2 (Wickham, 2016). We also performed two separate PCAs, one on the CB group and one on the WE and WC together (hereafter wild bumblebees, 'WB') to identify any further genetic sub-structuring within each of the two groups. Before grouping the WE and WC groups together, we calculated pairwise F ST (Weir & Cockerham, 1984)  the optimal number of clusters within the dataset), we ran a crossvalidation procedure with K values of 1-5 and with a 10-fold crossvalidation (CV). We performed 10 separate runs of the algorithm, using a different random seed for each run, taking the average of the 10 runs for our final result ( Figure S3). The same cross-validation procedure was run separately for the CB group and the WB group to identify any further sub-structure within groups. The CV error rates are reported as an output logfile for each value of K ( Figure   S3). For fineSTRUCTURE, we first phased the SNP-dataset using BEAGLE v.3.1.2 (Browning & Browning, 2009) and then SHAPEIT v2.r904 (Delaneau et al., 2012) to create a recombination map in

| Selection detection in commercial and wild bumblebees
To explore differentiation along the genome between wild and commercial bumblebees, we calculated pairwise F ST (Weir & Cockerham, 1984) and mean F ST per chromosome in 10 kb nonoverlapping windows using PopGenome v.2.7.5 (Pfeifer et al., 2014). As F ST is only a relative measure of divergence and is strongly influenced by the level of within-population diversity, we also used the cross-population extended haplotype homozygosity statistical test (XP-EHH) in selscan v1.2.0 (Szpiech & Hernandez, 2014). The method is used to detect recent directional selection via selective sweeps and compares haplotype lengths between two populations (Sabeti et al., 2007).
To identify stretches of extended homozygotic haplotypes in each group (CB and WB), we computed XP-EHH scores using the other group as the diversity reference. XP-EHH scores were standardized across 10 kb windows. To identify regions of the genome potentially under selection, we identified regions in which 10 kb windows over-

| Structural variant detection
To detect structural variants (SV) from our short-read whole-genome sequencing data, we used two SV detection programs: smoove v.0.2.5, which uses lumpy (Layer et al., 2014) and Breakdancer v.1.4.5 (Fan et al., 2014). Both SV callers were run for the CB and WB groups separately. For smoove outputs, we did not consider SVs that were marked as 'imprecise', <=1kb in length, and which had read pair support lower than the median for each group (CB: 39; WB:35). For Breakdancer outputs, we removed SVs that had a quality score <99, <=1kb in length, and which had a fewer number of reads than the median for each group (CB: 50, WB: 334).

| Gene annotation
To explore candidate genes within identified regions of high divergence between wild and commercial bumblebees, we extracted the gene ID, start and end position, accession IDs and gene name identified from the B. terrestris reference genome. Identified genes were used in gene ontology enrichment analysis (GO) using the R package biomaRt (Durinck et al., 2009) via Bioconductor v. 3.12 (Yu, 2014).
The reference list of genes was also matched with the KEGG (Kyoto Encyclopaedia of Genes and Genomes) ENZYME database to identify potential enzymes and their function using biomaRt.

| Data processing and quality control
We obtained a total of 5,672,780,248 reads across all genotyped individuals and the mean number of reads per individual was 5,909,1461 (Table S3). Of the total 96 sequenced individuals, three had >70% missing data and were therefore removed from the dataset. We identified five wild-caught full siblings (one sibling pair and one sibling trio). Based on the estimated kinship coefficient (1 st degree sibling: 0.177-0.354; Manichaikul et al., 2010), we retained two unrelated individuals of the five in the dataset. We also identified one additional individual with a highly negative estimated kinship coefficient of −0.9 in all pairwise tests, indicative of much lower relatedness than expected by chance, and increased genetic divergence (Manichaikul et al., 2010). We suspect that this individual was erroneously identified as our target species during COI barcoding and it was therefore removed. The final dataset consisted of 89 individuals for analysis.

| Genetic structure of commercial and wild bumblebees
The PCA showed two defined clusters, where the wild bumblebees (WB) cluster together and the commercial bees (CB) formed a separate cluster along PC1 (12.7% of the variance) (Figure 2a). There was no further genetic substructure observed within either the CB or the WB group when these were analysed separately using PCA ( Figures   S1 and S2). Pairwise F ST between the WE and WC groups was low (F ST = 0.0004), and we therefore grouped the WE and WC groups together in all following analyses. One CB individual (Sample 25) showed genetic distinctiveness along PC2 (6.9% variance). We did not find evidence that the underlying data for this individual were biased in any way, that is in the number of reads, coverage, alignment rate, or missing data. We did however find that this individual had considerably higher observed homozygosity (Ho: 0.92; He: 0.68) and F IS (0.980), suggesting that the low genetic diversity of this sample skewed its placement in PCA space. The PCA showed no clear evidence for shared genetic assignment between the wild and commercial bees. Notably, two of the WE samples (samples 3 and 71) from two different experimental sites (V2 and T6, see Table 1) clustered with the CB group. We believe that these individuals represent commercial 'escapees' that were caught foraging outside of the cultivation (greenhouse, open tunnel cultivations and free land). They also grouped within the CB cluster (Figure 2b,c). Additionally, three of the CB samples (samples 22, 30 and 37) were genetically assigned to the WB group (Figure 2b,c), which is most likely due to wild bees entering the commercial bumblebee nest, either accidentally due to disorientation ('drifting') or deliberately, for theft of resources or shelter, which is a documented phenomenon among both wild and commercial bumblebees where bees visit non-natal colonies (Zanette et al., 2014). We did not find any evidence for genetic admixture in the identified drifters or escapees.
The cross-validation procedure of ADMIXTURE to identify potentially introgressed individuals performed on the two groups (CB and WB) gave the highest support for one genetic cluster, K = 1 (0.488) ( Figure S3, see Figure S4 for variation across runs). However, two genetic clusters, K = 2, also exhibited a low cross-validation error rate (0.501). Arguably, K = 1 and K = 2 often cannot be distinguished and it is recommended to explore multiple K values of the dataset (Janes et al., 2017) (see Figure S5 for K = 3). Apparent from the K = 2 result (and K = 3, see Figure S5) was a lack of admixed samples (Figure 2b), since none of the WB samples had high ancestry proportions derived from the CB group (see Figure S6 for bootstrapped admixture proportions). Again, this was also evident for the drifters and escapees. The cross-validation procedure run with the CB group and the WB group separately exhibited the highest support for K = 1 for both groups (CB = 0.64, WB = 0.48) ( Figure S3).
The fineSTRUCTURE coancestry matrix showed a clear separation between the CB and the WB groups, with members of the former tending to share the highest ancestry with each other (Figure 2c).
Some genetic substructuring was evident within the two groups.
Sample CB25 (also observed as an outlier in the PCA, Figure 2a

| Genetic diversity
All our quantified diversity statistics (Watterson's theta (θw), Nei's pi (π), Tajima's D and F IS ) were statistically different between the WB and CB groups. Although genome-wide π distributions were similar (Figure 2d), the WB group had a higher median π (0.00082, 95% CIs: 0.000-0.0016) compared to the CB group (0.00077, 95% CIs: 0.000-0.0016) (p < 0.0001). Higher median π in the WB group was also reflected in the per-chromosome calculations (Table S4). In contrast, the median θw was significantly lower in the WB group (5.85, 95% CIs: 0.549-10.999) compared to the CB group (6.86, 95% CIs: 0.527-13.272) (p < 0.0001), also apparent from the genome-wide distribution of the statistic (Figure 2e). The number of segregating sites (θw) per chromosome was also marginally higher in the CB group (Table S4). The discrepancy between π and θw was reflected in

| Genome-wide differentiation
The median global F ST across the whole genome was moderate to low between the CB and WB groups, 0.034 (99% CIs: 0-0.182).
When analysing individual chromosomes, the median F ST ranged from 0.028-0.041 (Table S5). When we visualized F ST across the whole genome, chromosomes 10 and 11 both had notable peaks of high F ST (Figure 3a). Compared to the rest of the genome, chromosome 11 also showed positive XP-EHH in the CB group, and respectively negative XP-EHH in the WB group (Figure 3b and 3c). Across

| Gene functions under potential selection on chromosome 11
The 177 genes present in the candidate region on chromosome 11 ( Figure 4f) were associated with 102 unique gene ontology (GO) terms (Table S10). Five GO terms were significantly enriched at p < 0.05, although these were nonsignificant after correction for the false discovery rate using multiple comparisons: calcium ion binding

| DISCUSS ION
Using a whole-genome sequencing approach, we investigated evidence for introgression and genomic divergence between commercial and wild-caught B. terrestris in Southern Sweden. Despite the discovery of escaped commercial bumblebees among the wild-caught group, as well as wild B. terrestris 'drifters' inside the commercial bumblebee colonies, we did not find evidence for recent genomic introgression among our samples. Although wild bumblebees had significantly higher nucleotide diversity (π) than commercial ones, the number of segregating sites (θw) was higher in commercial bumblebees. This discrepancy was supported by positive Tajima's D values in wild bumblebees. Selection scans across the genome revealed a highly divergent region on chromosome 11 in commercial bumblebees, as well as a high degree of linkage and chromosomal deviations within this region. This is indicative of differential selection processes operating in this region. Although we found no clear evidence for a single large structural variant (SV) associated with this divergent region, we found an overall enrichment of shorter SVs associated with this region. The divergent region contained genes involved in cellular processes that have previously been associated with flight muscle contraction and structure in insects (Bullard & Pastore, 2011;Cao & Jin, 2020;Kržič et al., 2010;Rusu et al., 2017), associated with pathogen immune response (Ramsey et al., 2017). The lack of introgression indicates that the extensive and historical use of commercial bumblebees in Southern Sweden has so far not affected the evolutionary processes within wild populations of B. terrestris via introgression.

F I G U R E 3
Genome scans run across sliding windows of 10 kb, where the blue highlights the region on chromosome 11 containing SNPs with increased differentiation that are under putative selection; red SNPs indicate the upper 99% and lower 1% confidence intervals. (a) The genetic differentiation (pairwise F st ) between the two groups; (b) standard mean XP-EHH scores for the WB group; (c) standard mean XP-EHH scores for the CB group

| Genetic structure and a lack of introgression
We identified two distinct genetic groups that separated wild bumblebees and commercial bumblebees. The groups showed no substantial evidence for shared recent genetic ancestry, indicating a distinct lack of introgression within Southern Sweden. We found subtle variation in genetic co-ancestry among individuals within the two clusters, which was likely due to the presence of 2 nd and 3 rd degree siblings (Figure 2c). The dendrograms for each group (Figure 2c) showed that the commercial bumblebees consisted of nine lineages and the wild bumblebees had four lineages. This may reflect the use of B. terrestris queens originating from multiple locations in commercial breeding programs (Velthuis & Van Doorn, 2006). Notably, there was no genetic structure detected among the wild-caught bumblebees, including bees sampled from both experimental (WE) and control sites (WC), which may be explained by the capacity of bumblebees to disperse distances of several kilometres (Kraus et al., 2009). The fact that we did not find any genetic structuring is consistent with previous studies which have evidenced (even at large spatial scales) low to moderate genetic differentiation among wild B. terrestris populations in mainland Europe (Estoup et al., 1996;Seabra et al., 2019;Silva et al., 2020) (Table 1).  (Pedersen et al., 2020). However, we did not find any significant genetic divergence within our CB samples, indicating that distinct subspecies were either absent or highly admixed, likely because of interbreeding in the breeding facility. Historical or ongoing admixture among different subspecies might also partially explain why we detected more than twice the number of lineages in the CB group compared to the WB group (Figure 2c) (Rasmont et al., 2008). Unfortunately, we cannot be sure which subspecies are present in the CB group used in this study.
We observed that two of the wild-caught B. terrestris (from the experimental sites V2 and T6 sites) were genetically assigned to the commercial B. terrestris group. These two individuals were most likely foragers from commercial colonies, since foraging distances for B. terrestris at the distance they were caught (~700 m) are not exceptional (Wolf & Moritz, 2008). Importantly, these individuals do

| Genetic diversity and divergence between commercial and wild bumblebees
Both groups showed moderate genetic diversity. However, nucleotide diversity (π) was higher in wild bumblebees, while the number of segregating sites (θw) was higher in commercial bumblebees. The discrepancy between the pattern for these two diversity statistics was supported by significantly positive Tajima's D values in wild bumblebees, which could be indicative of a sudden population contraction. Numerous studies document population contractions and wild bee declines in both abundance and diversity due to anthropogenic change (Goulson et al., 2008;Goulson & Hughes, 2015;Kerr et al., 2012;Kosior et al., 2007;Marshman et al., 2019). However, in our study area of Southern Sweden, the relative abundance of in Southern Sweden. Given the effect of pathogens on wild bee populations globally, and also in commercial colonies (Cameron & Sadd, 2020), studying the role of balancing selection in conferring potential immunity differences between these groups represents interesting future research.
Other studies have found moderately higher genetic diversity (observed heterozygosity and mean allelic richness) in wild compared to commercial bumblebees in the UK (Hart et al., 2021) and New England, USA (Suni et al., 2017), and in mainland Europe (Moreira et al., 2015). Seabra et al. (2019), however, found no difference in genetic diversity (expected heterozygosity) between wild and commercial B. terrestris. Interestingly, studies that report higher genetic diversity used microsatellite markers, while Silva et al. (2020) and Seabra et al. (2019), who reported similar values of genetic diversity to those detailed here, used RADseq. RADseq and WGS provide a much higher genetic data resolution compared to microsatellitebased analyses and might provide a more accurate estimate of bumblebee genetic diversity, as discussed by Lozier (2014).
We did not detect significant inbreeding within either the commercial or wild bumblebee groups. However, although inbreeding was low in both groups, there was overall a higher level of global

| Identification of a highly divergent region on chromosome 11
We found moderate to low global genetic differentiation between commercial and wild B. terrestris (median F ST = 0.034), but discovered a notable region of high pairwise F ST at the tail-end of chromosome 11 (median F ST = 0.04) (Figure 3a). This high peak of F ST , together with evidence of an excess of extended homozygous haplotypes ( Figure 4b) and reduced genetic diversity within this region in the commercial B. terrestris group, are indicative of differential selection processes operating in wild and commercial B. terrestris.
We observed evidence for reduced recombination in this region in the form of elevated LD (Figure 4e) and shifts in localized heterogeneity ( Figure 4d). Based on a linkage map of B. terrestris (Stolle et al., 2011), chromosome 11 did not show any evidence of having unusual recombination rates and also showed high homology to A. mellifera.
B. terrestris is reported to have higher than average genome-wide recombination rates compared to other insects (Stolle et al., 2011;Wilfert et al., 2007), which lends further support to the hypothesis that chromosome 11 may be under selection in commercial bees.
While we did not observe evidence for large structural variants encompassing the entire candidate region on chromosome 11, we did see an increase in the number of structural variants (SVs). However, to determine the role of these SVs and to assess if they are being maintained via linkage, long-read sequencing would be required as opposed to the short-read data we analyse in our study. Additionally, we observed an increase in gene density in the diverged region of chromosome 11 ( Figure S11), which can also cause signals of chromosomal deviations. Future studies should aim to fully characterize this region and the underlying processes maintaining the differentiation between wild and commercial bees.
We identified 177 genes in the candidate region on chromosome 11. However, since none of the GO terms were significant after multiple-test correction, we cannot be certain about which genes, or potential pathways may be under selection. We will therefore only briefly discuss the potential relevance of some genes which could be of functional importance in driving differences between wild and commercial bumblebees. Two of the seven genes with the 'calcium ion binding' GO term were identified as Troponin C (TnC) and are associated with regulating flight muscle contraction in insects (Bullard & Pastore, 2011;Cao & Jin, 2020;Kržič et al., 2010;Rusu et al., 2017). Two additional genes (sarcomeric α-actinin and sarcoplasmic calcium-binding protein 1) were also identified. Sarcomeric α-actinin acts as structural support in the Z-line in the flight muscles of insects (Cao & Jin, 2020;Kržič et al., 2010), while the sarcoplasmic calciumbinding protein 1 has several functions commonly found in muscle and neuronal tissues exclusively in invertebrates (Hermann & Cox, 1995;Pauls et al., 1993). Arguably, recent selection on commercial B. terrestris workers has taken place under the laboratory conditions within which they are bred. Under such conditions, flight is highly constrained, suggesting that the observed selection on flight muscle genes may actually be for decreased use of the flight muscles.
Moreover, two out of the four genes associated with 'chitin binding protein' GO terms were identified as mucin-5AC and chitinase domain-containing protein 1 (CHID1). Mucin-5AC is a glycoprotein (gel-forming) that is a major component in the mucus lining in both the respiratory tract and stomach in humans (Cornick et al., 2015;Lang et al., 2016), but is also present in the insect gut lining (Ramsey et al., 2017). The insect midgut epithelium is lined by the peritrophic membrane (PM), which primarily consists of chitin and glycoproteins (Wang & Granado, 1997). This acts as a protective barrier against ingested pathogens and aids food digestion (Kuraishi et al., 2011).
CHID1 was the second gene found in the candidate region that is involved in chitin binding and has been suggested to play a role in the innate immune response against chitin-binding pathogens (Ab et al., 2020). Nosema bombi and Crithidia bombi, two of the most common gut pathogens that infect bumblebees (Geslin et al., 2017), infect via digestion of spore-infected food or via contaminated bee faeces and develop in the midgut epithelial lining (Paris et al., 2018;Pham & Schneider, 2008). It is possible that the eusocial behaviour of bumblebees selects for increased gut pathogen immunity since it creates a favourable environment for pathogen spread to genetically similar individuals (Simone-Finstrom, 2017). However, further research is needed to better understand the function of genes that may be selected for in commercial versus wild bumblebees.

| CON CLUS IONS
Human-mediated introduction of non-native commercial bumblebees for pollination services is known to have negative ecological effects on some local pollinator fauna and plant-pollinator relationships (Aizen et al., 2019;Evans, 2017;Hingston & McQuillan, 1998;Morales et al., 2013). However, the evolutionary outcome of this remains uncertain. To our knowledge, this is the first study to use whole-genome sequencing data for the purpose of detecting introgression between wild and commercial B. terrestris, which allows for higher genetic data resolution compared to reduced marker techniques. We found individuals from the commercial hives foraging up to 700 m from the colony where we collected wild bumblebees, suggesting that hybridization is a possibility if queens and drones also escape into wild environments. The lack of recent introgression in our study suggests that the historical and current extensive use of commercial B. terrestris in Southern Sweden is yet to affect evolutionary processes within local wild populations of B. terrestris.
However, given our relatively low sample size, we cannot exclude the possibility of introgression between wild and commercial B. terrestris since commercial individuals were identified among wild bumblebees in the natural environment. Since different subspecies of B. terrestris differ in traits related to phenology, foraging efficiency, colony size and parasite resistance (Rasmont et al., 2008), divergent selection processes can act on wild and commercial B. terrestris.
Additionally, the discovery of a highly divergent region on chromo-

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
Whole-genome sequencing data are available at the European Nucleotide Archive (ENA) under the Study Accession PRJEB49221.