The tragedy of the common? A comparative population genomic study of two bumblebee species

Within the theoretical framework of the small population paradigm, we investigated the population genomics and parasite load of two bumblebee species across the UK and Ireland. Bombus pratorum is widespread and common throughout its range while Bombus monticola is restricted to higher altitudes and shows a more fragmented distribution. Bombus monticola showed stronger population structuring, isolation‐by‐distance, and a deficit of heterozygotes in the most isolated population in the south of its range (Dartmoor). Heterozygosity and inbreeding coefficients (FIS) were comparable between both species, but the proportion of polymorphic sites was much greater in B. pratorum. Notably, both species have suffered significant declines in Ne over the last 100 generations and estimates and declines for both species were of similar orders of magnitude. No pattern of increased parasite prevalence in populations of lower heterozygosity was observed. Instead, ecological and demographic factors (age, latitude, date, habitat suitability) were the main drivers of parasite prevalence. Distinct patterns of selection were observed in both species in regions involved in regulation of transcription and neurotransmission and in particular pathways targeted by neonicotinoid insecticides. Our results highlight the pressing need for monitoring to include common as well as rare species. This should not focus solely on census population counts, but include estimates of Ne. We also highlight the need for further work to establish adaptive shifts in globally important pollinator communities.


INTRODUCTION
Conservation efforts and priorities are generally focused on species that exist in small and/or declining populations since they are the most vulnerable to extinction.This is underpinned by the 'small population paradigm' (SPP, Gilpin & Soulé, 1986;Caughley, 1994), which predicts that small populations are vulnerable to the extinction vortex.
The vortex is created when populations decline to a small effective size (Ne), leading to high drift load (Willi et al., 2013) and decreasing adaptive potential.The consequent impacts on individual fitness become increasingly acute through each generation and without intervention, the interaction of demographic and genetic stochasticity can lead to extinction.Conservation genetics has since been dominated by efforts to measure standing genetic diversity, drift, inbreeding, gene flow and effective population size (Ne) (Ouborg et al., 2010).These parameters are quantified by screening neutral genetic markers since they are free from the complex potential influences of selection and are expected to behave in an evolutionarily predictable way in populations through time.Typically, such studies are restricted to single species and often in very limited geographic locations (i.e., those populations/species of concern that are rare and localised).
While this framework has led to a wealth of insightful research, criticisms have followed.First, the single rare species approach often lacks taxonomically relevant and directly comparable datasets for context, hampering the ability to correctly interpret and identify key conservation threats in wild populations.Taxonomic groups vary markedly in levels of standing genetic diversity (Ellegren & Galtier, 2016) so population genetic metrics are best interpreted in the context of comparative data from related species with similar life histories (Leffler et al., 2012).Such comparative studies are rare, however, and while such taxonomic context may be available through other independent studies, the use of different marker sets can often make direct comparisons difficult.A single focus on rarity can also obfuscate trends in species that have historically been (and may appear to remain) relatively common.From some conservation perspectives, such species are (far more significant than those that are rare or localised since the impact of their loss or decline may be more ecologically severe (Gaston & Fuller, 2008)).Despite the fate of the passenger pigeon, among many (see e.g., Roycroft et al., 2021), and regular calls for more widespread monitoring to avoid catastrophic and rapid extirpations of more common organisms (e.g., Gaston & Fuller, 2008), the conservation lens remains strongly focused on rare species or populations.
Secondly, there are a number of scenarios where populations may not behave in accordance with the theoretical predictions of the SPP.Habel and Schmitt (2012), for example, argue that highly localised 'specialist' species occurring in relatively low abundances may carry much lower genetic load due to purging or other strong selection.There is certainly good evidence that purging through selection can remove genetic load (e.g., Gloag et al., 2016;Robinson et al., 2018;Khan et al., 2021;Xie et al., 2022, although see also Wang et al., 2021) but how effective purging will be will depend strongly on how rapid and severe declines in population sizes are (i.e., historical demography): selection driven purging processes can be over-powered by the strength of drift if that reduction is dramatic and fast (see e.g., L opez- Cortegano et al., 2016;Pérez-Pereira et al., 2022).
Thirdly, whether genetic factors present a real extinction threat to small in situ populations has been the subject of debate since the inception of the theoretical framework.Several authors have suggested that demographic and environmental factors are more likely to drive wild populations to extinction before these genetic effects become so extreme that they directly cause extirpation (see e.g., Caughley, 1994;Lande, 1988).There is now compelling evidence from a wide variety of studies that the increased mutational (genetic) load created by drift and inbreeding in small populations can be linked directly to extinction risk in wild populations (Nonaka et al., 2019;O'Grady et al., 2008;Saccheri et al., 1998;Willi et al., 2013).Of course, these studies do not establish a general rule, but they do provide examples of where genetic factors are likely to have been important.
While there has been an emphasis on studies assessing the effects of genetic load on fitness (inbreeding depression), the parallel issue of maintaining genetic diversity to enable adaptive potential has been less controversial.Although adaptive variation is more difficult to study given our limited toolbox for identifying, functional genetic variation, clear links have been established between low genetic diversity and extinction risk (e.g., Evans & Sheldon, 2008;Spielman, Brook, Briscoe, & Frankham, 2004;Spielman, Brook, & Frankham, 2004; see also Lanfear et al., 2010).Thus, while the ability to assess adaptive genetic diversity, particularly in non-model systems, is still in its infancy, the pressing calls (e.g., Ralls et al., 2018;Razgour et al., 2019) for maintaining genetic diversity are uncontentious, particularly in the face of current climate change predictions (Bush et al., 2016).
A final aspect of the SPP that has not been without controversy is which component of genetic diversity is most appropriate to assess.
Recently, Teixeira and Huber (2021) asserted that variation at neutral genetic markers is at best only loosely linked to significant population declines and instead advocate a focus on functional genomic regions.DeWoody et al. (2021) and Kardos et al. (2021), however, present compelling evidence to the contrary and, as noted above, from a practical perspective we currently lack the sophisticated level of understanding of functional genetic diversity to enable an exclusive functional approach.There is also good evidence however, that, in some circumstances, differences in functional markers can arise between populations without notable patterns in neutral diversity being observed (see e.g., Allendorf et al., 2010;Krohn et al., 2019).Fortunately, much more comprehensive genome-wide screening has recently become accessible and affordable for population-level studies.This allows both neutral and functional markers to be screened simultaneously, and at much greater scale, significantly enhancing the resolving power of neutral markers while also giving important initial insights into variation at functional markers in non-model organisms (see e.g., Hoelzel et al., 2019).
Within this framework, we took a comparative, conservation genomic approach to assess paired samples from populations of two congeneric species of bumblebee (Bombus spp.) across their British and Irish ranges, relating aspects of both neutral and functional genomic diversity to measures of fitness (parasite load) and habitat suitability.While data remain surprisingly limited for many species, there is strong evidence for global declines in bumblebee species (e.g., Cameron & Sadd, 2020;Goulson et al., 2005;Vray et al., 2019) along with indications that specialised species may be most vulnerable (e.g., Casey et al., 2015), and that declines may be phylogenetically structured (Arbetman et al., 2017;Cameron et al., 2007Cameron et al., , 2011)).Furthermore, as eusocial, haplodiploid organisms, bumblebees are especially prone to extinction (Zayed & Packer, 2005) and census-based abundance estimates can be particularly misrepresentative of genetic effective population sizes (Ne).As key ecosystem service providers pollinating both crops and wildflowers, this is a group of particular conservation concern.Bumblebees also provide ideal models for testing predictions of the small population paradigm since species are variably threatened and isolated, ranging from the very rare and highly localised to those that are widespread and common.
The two congeners studied belong to the subgenus Pyrobombus, a group showing resilience to overall patterns of decline (Arbetman et al., 2017).Bombus monticola is highly localised in its distribution and included on English Nature's Programme for Recovery due to recent significant declines (BWARS, 2022;BCT, 2022).In contrast, Bombus pratorum is one of the most ubiquitous species found across the UK (BWARS, 2022;BCT, 2022).Several species in this subgenus have recently expanded their range (Biella et al., 2021;Huml et al., 2021;Rasmont et al., 2015;Richardson et al., 2019) including B. monticola and B. pratorum, which were first recorded in Ireland in 1974and 1947, respectively (Fitzpatrick et al., 2007;Speight & MCD, 1974).Despite these range expansions, B. monticola has historically been restricted to areas of higher altitude across Europe, while B. pratorum represents one of the most widespread and abundant species found across a range of habitats in Europe (Rasmont et al., 2015).Neither has been the subject of landscape-scale population genetic studies to date.Specifically, our aims were: 1.To assess and interpret genetic diversity (heterozygosity and % polymorphic sites), drift, inbreeding, gene flow and effective population size (Ne); 2. To comparatively assess any associations between population genetic parameters and environmental (habitat) and biological (parasite) factors; 3. To identify whether particular and/or consistent functional regions of the genome are under strong selection; 4. To interpret the results in the context of the factors potentially driving declines and key conservation priorities.

Sample collection
Bombus monticola and B. pratorum workers (females) were sampled from nine and eight locations, respectively (Figure 1) between June October 2011 (B.pratorum) (at some sites the species did not cooccur).Precautions were taken to minimise the probability of collecting sisters within sampling sites (Goulson et al., 2011), although putative sisters were also filtered from the subsequent analysis (see below).Individuals were frozen until processing using a portable freezer.

Wing wear
As immunocompetence can decline with age in Bombus (Doums & Schmid-Hempel, 2000), the relative age of individual bees was estimated by assessing the level of wing wear, following a modified version of Mueller and Wolf-Mueller's (1993) method where 0 was the lowest and 4+ the highest level of wing wear (as in Goulson et al., 2012).Wing wear is frequently used as a correlate of age in bumblebees in the absence of other reliable individual age estimators (e.g., Colla et al., 2006;Allen et al., 2007;Whitehorn et al., 2011) however, it is important to note that other factors, such as individual foraging behaviour or predation pressure also affect wing wear (Foster & Cartar, 2011).

Measuring parasite load
The tracheal tubes and tissues of the abdomens of all samples were inspected for the presence or absence of the tracheal mite Locustacharus buchneri.
All samples were also screened for the presence or absence of the microparasites Apicystis bombi, Crithidia bombi and Vairimorpha bombi using species-specific molecular markers (SSUrRNA, Klee et al., 2006;Neo, Meeus et al., 2010;ITS1, Schmid-Hempel & Tognazzo, 2010).The contents of the abdomen were added to 500 μl of 10% Chelex ® (10% volume Chelex ® resin, 50 mM TRIS) and homogenised using a sterile micropestle.These were heated at 56 C for 30 min, then, boiled for 8 min at 98 C, before being centrifuged at 13,000g for 5 min.The supernatant was removed and 3 μl was used for subsequent PCR reactions.The PCR conditions were: 3 μl template, 1.25 μM each primer (0.625 μM for C. bombi primers), 2X Biomix (Bioline), made up to 10 μl using molecular grade water.The PCR cycle was: 98 C for 2 min; 50 cycles of 98 C for 30 s, 45 s annealing at 48 C for SSUrRNA, 57 C for Neo and 55 C for ITS1, 72 C for 1 min; 72 C for 5 min.After PCR, 5 μl of the reaction was run for 1 h at 110 V on a 1% agarose gel stained with 5 μl of SYBRSafe (Life Technologies).

RAD library preparation
DNA was extracted from individuals from a sub-sample of six and five locations for B. monticola and B. pratorum, respectively (Figure 1) from thorax and leg tissue using a CTAB ammonium acetate protocol (Nicholls et al., 2000).Not all originally sampled locations could be genetically tested due to issues with sample quality.Extractions were quantified on a Qubit 3.0 Fluorometer (Thermo Fisher Scientific) after RNAse treatment.DNA extracted from individual bees was digested using either 100-200 ng of DNA and 5 Units of MluCl (New England Biolabs), 10 Units of MspI (New England Biolabs) in a 10 μl volume, or 200-400 ng of DNA, 10 Units of MluCl (New England Biolabs), 20 Units of MspI (New England Biolabs) in a 20 μl volume at 37 C for 3 h.Library preparation followed Peterson et al. (2012)

Data processing and variant calling
Raw reads were de-multiplexed, trimmed and quality filtered using the process_radtags.plscript in Stacks 2.48 (Rochette et al., 2019).
Reads were filtered and removed using default values in Stacks 2.48 (Rochette et al., 2019;i.e., those without intact restriction sites or an average phred score below 10 in any 22 bp window).All reads were trimmed to a uniform length of 140 bp to remove low-quality base calls at the end of the reads as standard (O'Leary et al., 2018).Individuals with low sequencing coverage and a representation below 1 M paired reads were excluded from the analysis (N = 54).
Filtered reads were assembled de novo using the denovo.plscript in Stacks 2.48 (Rochette et al., 2019).Assembly parameters were optimised following the recommendations of Paris et al. (2017) using a subset of samples including all replicates and a minimum of two samples per population (B.monticola N = 33, B. pratorum N = 29).For each parameter set, SNPs were called using the populations.plscript in Stacks 2.48 (Rochette et al., 2019) if they were present in at least 80% of individuals and showed a minor allele frequency above 1%.
Consistency of variant calls (SNP, allele and locus error rates) across replicates was assessed using the R-scripts provided by Mastretta-Yanes et al. (2015) to identify the best performing parameter settings.
After establishing the optimal parameter values that minimised error rates across replicate samples, these were used to build the final catalogue using eight samples from each population that showed the highest coverage, following general recommendations (Rochette et al., 2019).As additional filtering was carried out on the original Stacks output (see below), the calculation of allele and locus error rates only considered variable contigs in contrast to the original implementation by Mastretta-Yanes et al. (2015).If anything, this is expected to result in an overestimation of error rather than an underestimation.Next, de novo contigs were mapped against the Bombus terrestris genome (Bter_1.0assembly, Ensembl) using the integrate alignment tool in Stacks (Rochette et al., 2019) in order to take advantage of the linkage group and functional annotations this assembly offers.Although by aligning to the genome of a different species some information is inevitably lost, for example, where annotations are incomplete or where alignment fails in areas of high divergence, this is a necessary approach when working with non-model species (Shafer et al., 2017).Alignments of individual sample files against the catalogue, which contains all contigs resulting from the de novo assembly, were used to build an mpileup file using samtools (Li et al., 2009).The mpileup file was exported to vcf format including both variant and non-variant sites using bcftools (Li et al., 2009) and filtered in the same way as SNPs described below.
The generated vcf file was further filtered using vcftools (Danecek et al., 2011) (Table S1) et al., 2018).SNPs that showed a mean coverage across all individuals higher than double the mode average were also filtered, as these are likely to be derived from paralogous loci (Willis et al., 2017).Finally, loci significantly deviating from Hardy-Weinberg equilibrium after correction for multiple testing were also filtered, as these are likely to be affected by genotyping errors (Attia et al., 2010).Corrections for multiple testing were carried out in the R package qvalue v2.16 (Storey et al., 2015) applying a false discovery rate of 0.05 (Storey & Tibshirani, 2003).Coefficients of relatedness between all individuals were calculated in vcftools (Danecek et al., 2011) and for each pair showing a coefficient of relatedness >0.25, the individual with fewer reads was removed from all population genetic analyses (since the absence of kin is an assumption and would bias parameter estimates).

Assessment of genetic diversity, population structure and effective population size
Basic measures of population genetic diversity and pairwise F ST (Weir & Cockerham, 1984) were calculated using the R package hierfstat (Goudet, 2005).To generate confidence intervals for pairwise F ST values and inbreeding coefficients bootstrapping over loci was performed for 1000 cycles.
Population structure and admixture was assessed by principal coordinate analysis (PCoA), which is preferable to a principal component analysis when there are fewer individuals than characters (Rohlf, 1972).PCoA tests were performed in the R package ade4 v.1.7-13(Dray & Dufour, 2007) based on Euclidean distance (Gower, 1966).A Bayesian approach implemented in the programme fastStructure (Raj et al., 2014) was also implemented to assess population structure and admixture.FastStructure assumes the independence of loci and therefore only one SNP was used per RAD-tag.To find the number of components (K) that best explain the structure of the data, fastStructure was run for values of K from 1 to 10.The Python script 'chooseK.py'provided within fastStructure was then used to identify the value of K that best fitted the data.The genetic clusters identified were used in further analysis for the definition of population demes.A Mantel test (Mantel, 1967) was further conducted in Genepop to test for isolation by distance based on pairwise F ST (Raymond & Rousset, 1995).
To estimate long-term effective population size, Watterson's theta (θ W ) (Watterson, 1975) was calculated in SambaR (de Jong et al., 2021) as has been estimated for other Bombus species (Lattorff et al., 2016).This parameter is equivalent to 3N e μ in haplodiploids.
For the mutation rate μ, Liu et al.'s (2017) estimate of 3.6 Â 10 À9 (CI: 2.38 Â 10 À9 5.37 Â 10 À9 ) for a single nucleotide per generation in Bombus terrestris was used.To assess trends in the recent trajectories of the effective population size, linkage disequilibrium (LD) across different inter-locus distances can be exploited (Waples & Do, 2010).LD patterns over larger distances are informative for recent N e while LD patterns over shorter distances are informative for the more distant past (see also Nadachowska-Brzyska et al., 2022, for an informative overview).This analysis was carried out using SNeP v1.1 (Barbato et al., 2015), which is designed to estimate trends in effective population size using LD.First, Beagle 5.1 (Browning & Browning, 2007) was used to phase genotypes and then vcftools (Danecek et al., 2011) was used to calculate the correlation coefficient (r 2 ) between any two SNPs as a measure of LD used in SNeP for each population.SNeP was run using default parameters with an adjusted recombination rate of 8.7 e À008 in accordance with the recombination rate (8.7 cM/M) observed in bumblebees (Liu et al., 2017).As methods based on LD can greatly underestimate absolute values of N e in haplodiploid species (Wang et al., 2016) we emphasise that the focus of this analysis is on the relative change in N e over time.Results are presented for generations 10-100, given the uncertainties around estimates in the most recent and most distant generations (Barbato et al., 2015;Corbin et al., 2012).

Investigation of predictors of parasite richness
Relationships between observed parasite richness (number of parasite species detected), wing wear, longitude and latitude as well as species-specific habitat suitability of the sampling site (estimated using Polce et al. (2018)) and date (represented as numerical value corresponding to the day and month of sampling) were investigated using a general linear modelling approach in the stats R package framework (Chambers & Hastie, 1992).Habitat suitability estimates were derived using predictor variables such as land cover, distance from natural or semi-natural habitat and climate variables (Polce et al., 2018).Values representing the probability of presence were then extracted for each sampling site.An additional GLM was performed with observed heterozygosity as an additional potential predictor including all locations with genetic data available.In both cases, a Poisson distribution was implemented.Furthermore, relationships were investigated independently for each parasite species based on presence and absence data using a binomial distribution.Akaike's information criterion (AIC) was used to compare models with different numbers of explanatory variables (Table 2) using the step AIC function in the R-package gamlss, which implements a stepwise reduction in variables based on AIC (Rigby & Stasinopoulos, 2005).

Detection of signatures of selection F ST -outlier approach
In accordance with recommendations to combine several methods for the identification of F ST -outliers to minimise the false positive detection rate (De Mita et al., 2013), here three different approaches were adopted.BayeScan (Foll & Gaggiotti, 2008), which implements a basic regression model to differentiate between locus and population specific effects on the distribution of F ST -values, was run using default parameters.Secondly, FstHet (Flanagan & Jones, 2017), which identifies outliers based on F ST -heterozygosity relationships, was run using the betahat option for the calculation of F ST (Cockerham & Weir, 1993).This method implements a correction for sample size.
Finally, pcadapt (Luu et al., 2017), which identifies outliers as loci that deviate from patterns of population structure defined by principal component analysis, was used.Overlap between the different methods was visualised using Venn diagrams (Oliveros, 2007).Outliers that were identified by at least two approaches were investigated using gene ontology (GO) analysis.

Selective sweep detection
RAiSD (Alachiotis & Pavlidis, 2018) was used to detect signatures of selective sweeps.This programme has the advantage of combining the three most common approaches used in the detection of selective sweeps, which are based on changes in the level of polymorphism, the site frequency spectrum (SFS) and LD, into a single statistic (μ).Data were generated across sliding windows of 10 SNPs to balance information content and resolution for the detection of signals of selective sweeps (Vatsiou et al., 2016).Outliers were then identified as those above the 95th percentile of the empirical distribution.In the absence of detailed knowledge of the demography of the study populations needed to generate a simulated neutral distribution, this approach has been shown to lead to comparable results (Vatsiou et al., 2016).

Gene ontology analysis
All SNPs were annotated to specific genes if they were located within a distance of 5000 bp and therefore in high LD, which is documented to decay after 10 kb (Stolle et al., 2011) in bumblebees.GO terms for the B. terrestris genome were obtained from the Ensembl database (metazoa genes 55; Yates et al., 2022).Furthermore, orthologous sequences were identified in the Apis mellifera and Drosophila melanogaster genome using Biomart (Kinsella et al., 2011), where no GO annotation could be obtained from the B. terrestris genome.Enrichment tests were conducted using the R package topGO (Alexa & Rahnenfuhrer, 2016) using Fisher's exact test on gene count (the 'classic' algorithm, Alexa et al., 2006) and applying the 'weight' and 'elim' algorithms, which account for dependencies within the GO hierarchies (Alexa et al., 2006).A minimum node size of 10 was required to prune our hierarchy from nodes with the support of less than 10 annotated genes, which represents a frequently applied threshold (e.g., Ahrens et al., 2013;Rademacher et al., 2017).Correction for multiple testing (FDR < 5%, q-value < .05)was carried out for the 'classic' algorithm using the qvalue package in R (Storey et al., 2015), but not the 'weight' and 'elim' algorithms for which multiple testing theory does not directly apply as tests are not independent and raw p-values ≤.05 were considered as significant (Alexa & Rahnenfuhrer, 2016)   best fit to the data (Figure 3a).A further separation of Dartmoor from Long Mynd and Snowdonia becomes apparent on the second and third axes for B. monticola (Figure 2a).The presence of further substructure between these three populations is also confirmed by fastStructure analysis on this population subset (Figure 3b).For B. pratorum, the PCoA shows that variation is predominantly

Population genetic parameters
FastStructure analysis for all populations (a) and a population subset (b) for Bombus monticola and all populations for Bombus pratorum (c).explained along the first axis, which separates Antrim from all other populations sampled (Figure 2b), which is consistent with the fas-tStructure analysis identifying two clusters as the best fit to the data (Figure 3c).A separate analysis in fastStructure to assess finer population structure across all B. pratorum populations except Antrim found no support for any further population structure.
Long-term estimates of the genetic effective population sizes in B. monticola ranged from 13,444 in Dartmoor to 50,700 in Glenshee (Table 1).The long-term estimates of the genetic effective population size in B. pratorum were 43,981 for the Antrim population and 113,907 for the panmictic population across England and Wales (Table 1).Trajectories of genetic effective population size changes over the last 100 generations show a strong decrease in both species (Figure 4).This was 9-fold for the population in England and Wales in B. pratorum and 2-7 fold for B. monticola populations (Figure 4).In both species, the most recently colonised sampling site in Ireland (Antrim) exhibits the smallest population size over recent generations, and the least decrease in size (Figure 4).

Parasite prevalence
Across all sites, the prevalence of L. buchneri was significantly higher in B. pratorum than in B. monticola (Mann-Whitney-Wilcoxon test: p < .001)but no differences were found for the other parasites investigated (Figure 5).Patterns of parasite prevalence varied greatly between locations and were not correlated between B. monticola and B. pratorum (i.e., there was no pattern of parasite prevalence with geographic location, Figure S1).The number of parasites observed showed a significant relationship with habitat suitability in B. monticola, with lower suitability values being associated with higher parasite loads in both the model including all sampled locations (Table 3.1A) and the model restricted to only locations with genetic data (Table 3.1B).The 'all sample' model also indicated a positive association between age (estimated using wing wear) and parasite load.Date was also significant and showed a trend of decreasing parasite load throughout the year (Table 3.1A).In B. pratorum a significant positive effect of latitude was observed in both models (Table 3.1A, B).Significant effects of predictor variables were also identified for infections with the parasites C. bombi, V. bombi and A. bombi for both B. monticola and B. pratorum (Table 3.2-3.5).
F I G U R E 4 SNeP analysis for Bombus monticola and Bombus pratorum showing the trajectory of effective population size (Ne) change in the past.Note that Ne is plotted on a log scale and generations into the past increase along the x-axis.Results are presented for generations 10-100, given the uncertainties around estimates in the most recent and most distant generations (Barbato et al., 2015;Corbin et al., 2012).T A B L E 3 Results of the general linear model used to test relationships between: number of parasites with the predictor variables longitude, latitude, habitat suitability, wing wear and sampling date for all sampled locations (1A) and with the additional predictor of heterozygosity (1B) for all locations that included genetic information; the individual infection status with C. bombi (2), V. bombi (3), A. bombi (4) and L. buchneri (5) Venn diagram of outliers identified by each of the three methods applied for Bombus monticola (left) and Bombus pratorum (right).

Detection of signatures of selection F ST -outlier approach
Of the three methods used, Bayescan identified the lowest number of outliers in both species (Figure 6) with no outliers identified for B. pratorum.This was followed by Fsthet and then Pcadapt, which identified the largest number of outliers.For B. monticola a small proportion of outliers were consistently identified by all three methods (Figure 6).Each of these fell within 5000 bp of one of eight genes, including a transcription factor (100647190), an RNA binding protein ( 100644912), an acetyl esterase (100649264), a receptor tyrosin kinase ( 100645523) and an inorganic phosphate cotransporter (100644143) (Table 4).Another subset of outliers identified were consistently called by at least two approaches (12.5% for B. monticola and 6.8% for B. pratorum) (Figure 6) and were considered for GO analysis.
This revealed these were significantly overrepresented for the molecular function of DNA binding transcription factor activity in both species (Table 5).

Selective sweep detection
Signals of selective sweeps detected by the RaisD programme were similar across populations of the same species, but also showed some congruency across species.Of all genes within areas identified as outliers in any population, 70% and 19% were consistently found to be outliers across both populations in B. pratorum and across at least four out of five populations in B. monticola, respectively.Molecular functions and biological processes, over-represented in areas that showed signals of a selective sweep in both species, were predominantly related to neurotransmission and regulation of transcription (Table 6).

DISCUSSION
Population structure and genetic diversity Our results are largely consistent with the predicted evolutionary and population genetic effects affecting species with different demographic histories.Gene flow is much higher in the widespread and more common B. pratorum than in B. monticola, with no population structure detected across mainland UK populations.In contrast  Consistent with a larger effect of genetic drift in B. monticola, the number and proportion of polymorphic sites was on average more than three times higher in B. pratorum (Table 1).The relatively smaller proportion of polymorphic sites observed in B. monticola may be of concern given the lower adaptive genetic potential this species therefore has to respond to future environmental change.
However, estimates of heterozygosity and nucleotide diversity across polymorphic sites were comparable between both species and only the Dartmoor population showed any strong evidence of inbreeding in B. monticola.Estimates of genetic diversity are also similar to those observed in other bumblebee species assessed using a RAD-seq approach (Huml et al., 2021;Jackson et al., 2018;Lozier, 2014).

Trends in effective population size
The effective population size is formally defined as the 'the size of the ideal (Wright-Fisher) population (N) that will result in the same amount of gentic drift as in the actual population being considered' (Allendorf et al., 2013).Its relationship with the census population size (Nc) varies greatly among taxa but is generally expected to be considerably lower than the census population size in most scenarios, given that not all individuals in a population successfully breed into the next generation and that there is likely to be variance in family size in natural populations.In eusocial insects, Ne is expected to be particularly lower than Nc, since the unit of reproduction is the colony and not the individual (Chapman & Bourke, 2001).The 'threshold' of concern for Ne is a subject of much debate and the use of 'magic numbers' has been much derided (Caughley, 1994;Flather et al., 2011;Lande, 1988Lande, , 1995)).Regardless, the long-term values of Ne reported for both species are comfortably higher than the thresholds generally discussed in the literature (Lande, 1995;see also Frankham et al., 2014).What is surprising, however, is that the Ne values observed for B. pratorum are of similar orders of magnitude to B. monticola (Table 1), despite the former being generally considered as widespread and common in the UK (and Ireland).Lattorff et al.
(2016) reported much higher Ne values than those observed here for other common bumblebee species using the same approach, although as those were continental European populations this would be expected (B.terretris 2.7 Â 10 6 ; B. lapidarius 6.4 Â 10 5 , estimates that would be higher when applied to the Bombus mutation rate that has since become available (Liu et al., 2017)).We note the danger of overinterpreting raw values using single point datasets (Nadachowske-Bryska et al., 2022), and instead focus on the general pattern of decline observed in both species over recent generations (Figure 4).

Patterns of parasite prevalence
The number of infections with different parasites was best explained by the distribution of habitat suitablility in B. monticola and latitude in B. pratorum (Table 3).Habitat suitability estimates likely relate to the distribution, quality and seasonal availability of ressources.Maintenance of immunity is energetically costly (Schmid-Hempel, 2005) and pollen deprivation has been shown to reduce both the constitutive immunity (Roger et al., 2017) and the specificity of the immune response in B. terrestris (Brunner et al., 2014) although food deprivation can also lead to reduced parasite loads where these rely on host nutrition (Conroy et al., 2016;Logan et al., 2005).Consistent with our results of higher infection rates for multiple parasites in B. pratorum at higher latitudes, McArt et al. ( 2017) also identified a trend for greater prevalence of V. bombi and increased local losses in four bumblebee species in North America at higher latitudes.
Although fungicide usage showed a similar geographic pattern, and was the strongest predictor of V. bombi prevalence identified by McArt et al. (2017), a higher prevalence of pathogens at higher latitudes may be a factor contributing to the generally observed lack of bumblebees tracking changing climate conditions at northern range limits (Kerr et al., 2015).
Although inbreeding and loss of genetic diversity have been reported to represent a major factor for increased susceptibility to disease (Ekroth et al., 2019;Spielman, Brook, Briscoe, & Frankham, 2004), here we did not find any associations between parasite prevalence and genetic diversity.This contrasts with results from studies at the colony-level.For example, lower parasite loads have been observed for C. bombi and V. bombi in B. terrestris within more genetically diverse colonies (Baer & Schmid-Hempel, 1999;Huth-Schwarz et al., 2012;Liersch & Schmid-Hempel, 1998;Parsche & Lattorff, 2018).Similarly, Cameron et al. ( 2011) observed higher prevalence of V. bombi in declining populations with lower genetic diversity.Although in the absence of polyandrous mating a large proportion of colony diversity should be captured by sampling an individual worker, colony parasite loads may be significantly under-or over-estimated from just one worker, as reported in other studies (Whitehorn et al., 2014).Although we also investigated population level effects of heterozygosity on susceptibility, the number of sampled populations was too small (n = 6 and 5 for B. monticola and B. pratorum, respectively) to derive reliable conclusions.

Signatures of selection
There were some striking similarities between the two species in signatures of selective sweeps within genomic regions relating to neurotransmission (e.g., ion channel activities, synaptic transmission and sleep) and regulation of transcription (Table 6).A recent study on another species of Pyrobombus, Bombus hypnorum, also reported evidence for selection acting on regions related to ion channel activities and regulation of transcription (Huml et al., 2021).Signatures of selection on genes involved in neurotransmission have also been reported across UK populations of Bombus terrestris (Colgan et al., 2022).Sun et al. (2021) found ionotropic glutamate receptor activity, important in neurotransmission (Benton et al., 2009), to be among the top 20% of genes showing elevated dN/dS ratio across a wider range of Bombus species, which may indicate functional divergence in respect to neurotransmission in bumblebees.Such areas of high divergence among species are of particular interest, as these are likely to be involved in lineage specific adaptations (Ellegren, 2014).Ionotropic glutamate receptors were among the genes showing evidence for positive selection in both species here (Table S2).Among the biological processes showing evidence for a selective sweep, we identified regulation of cholinergic synaptic transmission (Table 6), well known to represent the primary target for neonicotinoids (Matsuda et al., 2001).
These results support a growing body of evidence for strong selection on neurological genes in bumblebees and other insects, which may be indicative of adaptive shifts in response to modern selection pressures, such as pesticides, in this group.We note, however, the need to be cautious to avoid over-or mis-interpreting genomic patterns in this context, and highlight the need for further research in this area.
A RAD-seq approach inevitably only samples a proportion of the genome.While this inevitably means that many loci are not screened this does not negate the results reported (Lowry et al., 2017;Tiffin & Ross-Ibarra, 2014).Unless a species exhibits exceeding levels of polymorphism, RAD-seq is an effective method to obtain unbiased estimates of genomic diversity across the genome (Cariou et al., 2016), and is a powerful tool for the study of population genomics and natural selection (McKinney et al., 2017).
Here, habitat suitablility strongly affected the range-restricted species, B. monticola, both in terms of the highly restricted population connectivity and parasite prevalence.This suggests that this species is particularly vulnerable to habitat degradation, and at one site (Dartmoor) we observed some potential indicators of population stress (inbreeding).The conservation actions relevant for B. monticola would seem relatively straight-forward: conservation and restoration of upland heathland habitats, already identified as a conservation priority under the EU Habitats Directive 1994 due to widespread losses over the last two centuries, particularly in the UK (Fagúndez, 2013).The picture for B. pratorum is less tractable and arguably more concerning.
Wide areas of habitat suitable for this species are available, but either that habitat is of declining quality, or other factors relating to declines are operating, or both.There are several possible and non-exclusive possibilities, for example, changes in weather patterns, long-term exposure to agrochemcials and pollutants, and decreased availability of suitable forage and nesting sites.Although the potential genomic regions under selection include targets for neonicotinoids, which may be indicative of the particular stresses that these populations are under, these are by no means exhaustive and require further study for causative links to be established.Decreases in the abundance of common species are of particular concern as these disproportionally affect ecosystem services (Winfree et al., 2015).We might take heed from the warning of Bombus sylvarum: now one of the UK's rarest species it was common across large parts of the the UK just a century ago (Sladen, 1912).Similar patterns of declines in common Bombus species have been reported in North America (Cameron et al., 2011).
Evidently, there is a presssing need for effective monitoring systems and research focus on the drivers of declines not just on rare and already threatened species, but for those species that are apparently of least conservation concern (Gaston & Fuller, 2008;Lindenmayer et al., 2011;Roycroft et al., 2021).Identifying differences in the potential of species to adaptively respond to current threats and environmental change is likely to enhance our understanding of species vulnerability to decline, and help in directing conservation efforts more effectively (Hoelzel et al., 2019;Lavergne et al., 2013;Salamin et al., 2010).

F
I G U R E 1 Sampling locations for Bombus pratorum (right) and Bombus monticola (left) with and without genetic information; species specific habitat suitability estimates from Polce et al. (2018) are also shown 2010 and September 2011 (B.monticola) and between June 2010 and in randomised batches of 12 samples including 15 technical replicates for each species (N = 324, individually tagged samples in the library).Paired-end sequencing (150 bp) was carried out on a NovaSeq6000 by BGI-Hong Kong Co. Ltd.
following the recommendations of O'Leary et al. (2018; SNPs with a minimum allele count below 3 and a minimum mean depth of coverage across all individuals below 15 were filtered).A minimum genotype depth of five reads was required to call a SNP in any individual to avoid a bias towards homozygote calls in areas of low sequencing coverage (O'Leary et al., 2018).The percentage of missing data was minimised through iteratively removing SNPs and individuals to achieve a maximum of 5% of missing data for any SNP and a maximum of 25% of missing data in any individual (O'Leary most of the variation is explained by the first three axes with four distinct clusters separated on the first axis.This is consistent with the results of fastStructure, which identified four distinct clusters as the F I G U R E 2 PCoA on Euclidian distance for Bombus monticola (a) and Bombus pratorum (b).

B
. monticola exhibited strong population structuring with low levels of gene flow.The only sampling locations where significant differentiation was not observed were Long Mynd and Snowdonia, which are less than 100 km apart and connected largely by habitat of high suitability according to the estimates ofPolce et al. (2018).
Measurements of population genetic diversity per site and species: number and proportion of polymorphic sites, number of private alleles, observed (H IS significantly larger than zero are highlighted in bold. Population pairwise F ST on the upper and confidence limits on the lower triangular of the matrix for Bombus monticola (top) and Bombus pratorum (bottom) . The background gene set used in each analysis represented only areas actually sampled (e.g., only regions within T A B L E 2 windows that qualified in the RAiSD analysis, present in the number of populations/species specified or only genes within 5000 bp of SNPs for the F ST -analysis).
outliers identified by all three methods in B. monticola Gene ontology analysis of F ST -outliers identified by at least two methods for Bombus monticola and Bombus pratorum respectively; significant p-values are highlighted in bold; results are presented where at least two evaluation metrics (elim/weight/classic) were significant Gene ontology analysis of RaisD outliers identified in each species separately and combined (in a minimum of four population clusters for B. monticola and 2 population clusters for B. pratorum); significant p-values are highlighted in bold; results are presented where at least two evaluation metrics (elim/weight/classic) were significant Downloaded from https://resjournals.onlinelibrary.wiley.com/doi/10.1111/icad.12626 by Plymouth University, Wiley Online Library on [01/02/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License T A B L E 6