Strong decreases in genetic diversity despite high gene flow for a solitary bee

Many organisms live in landscapes that include fragments of natural habitat interspersed with agricultural or urban areas. In such areas, dispersal across human-dominated landscapes may be important for the maintenance of genetic diversity within habitat fragments. Understanding the extent to which dispersal occurs across fragmented habitat and how genetic diversity changes over time is particularly important for bee populations. Population sizes of many bee species are declining, which can have profound negative impacts on both agricultural pollination and the fitness of wild plants. We quantified dispersal across 100 km of mixed landscapes in southern Costa Rica for the solitary bee species Euglossa championi, and we quantified changes in genetic diversity over an eleven-year period. Our SNP-based analysis revealed high dispersal across the geographic area and no evidence that deforestation restricted dispersal. However, we found strong decreases in genetic diversity over time. Expected heterozygosity decreased by 28% over eleven years, regardless of the amount of forest that surrounded sites. Our results suggest that while high dispersal over deforested landscapes may render E. championi less susceptible to land use change than species that show restricted dispersal, the strong decrease in genetic diversity within sites over time suggests this species may warrant conservation monitoring.


Introduction
Approximately 75% of the earth's land surface has been modified by humans (Luyssaert et al. 2014).In tropical areas, for which species richness and abundance are extremely high (Basset et al. 2012), land use change can have profound negative impacts on biodiversity (Barnes et al. 2017).Land use change impacts organisms directly via reductions in resource availability and the ability to find mates (Tscharntke et al. 2005).Longer-term negative consequences of land use change can include reductions in genetic diversity due to restricted dispersal and subsequent inbreeding (Schlaepfer et al. 2018).Persistence in fragmented landscapes is therefore thought to be facilitated by the ability of individuals to successfully travel across altered habitat (Luyssaert et al. 2014), although increased dispersal may also compromise persistence due to tradeoffs between dispersal and reproduction (Baguette and Schtickzelle 2006).
Whether organisms disperse across anthropogenicallyaltered habitat depends on their mode of movement.For organisms that can fly, altered habitat may not act as a direct barrier to movement, but travel over such areas can still have negative ecological and physiological impacts on individuals (Zuluaga et al. 2022).Deforested areas are often warmer than areas of intact habitat (Vargas Zeppetello et al. 2020), which may impose heat stress on organisms (Alves de Oliveira et al. 2021).In addition, organisms may be less likely to find places to camouflage themselves in deforested areas, leading to higher predation (Dalrymple et al. 2015).Individuals may therefore be less likely to use deforested areas (Milet-Pinheiro and Schlindwein 2005), or they may evolve morphological differences that better allow them to tolerate deforested habitat (Pinto et al. 2015;Cote et al. 2017;Spaniol et al. 2020;Brito et al. 2021;Suni and Dela Cruz 2021).
Understanding effects of land use change on pollinators is particularly important because pollination is a key ecosystem service.Over 87% of angiosperms require pollination by animals (Ollerton et al. 2011).In addition, pollinators enhance seed set in over 75% of crops (Kearns et al. 1998), and wild pollinators will visit nearby agricultural areas, thus reducing the costs of food production (Garibaldi et al. 2013).The population sizes of many pollinator species are declining, with land use and climatic changes as drivers (Potts et al. 2010;Soroye et al. 2020).Pollinator decline has worrisome implications for plants, with reductions in populations of wild plant species common in areas that have experienced pollinator declines (Biesmeijer et al. 2006;Thomann et al. 2013).In addition, implications for pollination vary globally because the increasing intensity of anthropogenic land use is predicted to lead to disproportionate reductions in pollinator biodiversity in tropical areas (Millard et al. 2021).
Among pollinators, bees may be especially vulnerable to negative effects of land use changes due to genetic factors.Bees are haplodiploid, which reduces effective population sizes (N e ) by no less than 25% relative to equally-sized diploid populations.In addition, bees have single-locus sex determination whereby heterozygotes develop into females but homozygotes develop into diploid males that are usually sterile (Whiting and Whiting 1925).Diploid male production is predicted to be exacerbated by small population sizes and inbreeding (Zayed and Packer 2005), which are associated with habitat fragmentation (Schlaepfer et al. 2018).While land use changes can impact all bees, solitary species may be the most vulnerable due smaller population sizes and lower reproductive rates relative to eusocial species (Minckley et al. 1994;Séguret et al. 2016).
Here, we investigate whether genetic parameters important to population persistence have changed over a period of 11 years for the solitary bee Euglossa championi across a fragmented geographic area that spans 100 km.Our approach integrates an assessment of both spatial and temporal genetic structure, and investigates the extent to which deforested land affects genetic structure.We build on previous work on the Euglossine bee E. championi that used microsatellite markers and that was unable to detect an effect of forest or distance on dispersal (Suni and Brosi 2012;Suni et al. 2014).Estimates of gene flow are strongly influenced by the number and type of genetic markers used (Lozier 2014;Allendorf 2017;Leroy et al. 2018).Therefore, genomic studies that use thousands of loci are needed to clarify whether E. championi individuals travel widely across deforested landscapes, of whether the lack of isolation by distance or landscape was due to low resolution of the markers previously used.
Using a dataset containing thousands of SNP loci, we tested the following predictions.First, genetic structure will be higher between sites separated by deforested landscapes.Second, genetic diversity will be higher within sites that are surrounded by more forest.Third, temporal turnover in allele frequencies across years will be lower for sites in larger tracts of forest.Fourth, the amount of time elapsed between arrival of two individuals to a foraging location will be inversely correlated with their relatedness.Our analysis of spatial and temporal genetic structure together provides a framework for concomitantly investigating dispersal over a large, fragmented geographic area, as well as the extent to which individuals use small patches of habitat.

Species and sampling
Euglossa championi is found from Mexico to Ecuador at elevations from sea level to around 2000 m, is about 13 mm in length, and pollinates orchids in at least 10 genera as well as many other tropical plants (Roubik and Hanson 2004).Euglossa championi can be distinguished from other sympatric species by facial markings and a fuzzy tuft on the middle leg (Roubik and Hanson 2004).We focused our collections on E. championi in the Puntarenas Province of Costa Rica because this area is predicted to experience the most extreme deforestation of the Costa Rican provinces in future years (Stan and Sanchez-Azofeifa 2019).
We sampled a total of 381 E. championi in 2009, 2010, 2017, 2018, and 2019 within six forest fragments (Fig. 1a; Table 1).In addition to forest fragments, the sampling area is comprised of agricultural areas, towns, and pastureland.We sampled each site between the hours of 9 am and 12 pm on sunny days using the chemical baits cineole and methyl salicylate.These bates attract male bees but not females, so the samples in our dataset are males.We placed baits in forest fragments between zero and 93 m from forest edges, approximately 1.5 m off the ground on tree trunks, and netted bees as they arrived.In 2017 and 2018, we recorded the exact time that each bee arrived to each bait.We stopped sampling after no bees arrived for more than 15 min.Ethyl acetate was used to kill bees, which were then placed in vials with 95% ethanol for transfer back to the laboratory and subsequent pinning.Specimens are vouchered at the University of San Francisco, San Francisco, CA.

DNA sequencing and SNP calling
Genomic DNA was extracted from the middle leg of each specimen using DNeasy Blood and Tissue Extraction Kits (Qiagen) and was quantified using a Qubit 2.0 fluorometer (Thermo-Fisher).Double digest RADseq (ddRADseq) libraries were prepared based on Poland et al. (2012).Briefly, ddRADseq libraries were prepared using 100 ng of DNA per individual as follows.DNA was digested with the enzymes PstI and MspI (New England Biolabs).Adaptor ligation was performed using unbarcoded adaptors that were synthesized by IDT (Integrated DNA Technologies) for use with barcoded PCR primers.Ligated DNA was then cleaned with Agencourt Ampure XP beads (Beckman Coulter) and used as a template for PCR.For each sample, PCR was performed using one of 285 uniquely barcoded TrueSeq primers that were synthesized by the University of California San Francisco Center for Advanced Technology (UCSF CAT).The DNA concentration of each PCR reaction was quantified using the AccuBlue DNA Quantification Kit (Biotium), and 40 ng of each sample was pooled.The pooled DNA was then cleaned using Agencourt Ampure XP beads, and was sizeselected (300-500 bp) using a Blue Pippin Prep (Sage Science).Success in obtaining accurate target fragment  1.Panel B: Temporal genetic structure for the six sites between 2009 and 2019, calculated as F ST between samples obtained from the same site across each pair of years.Colors repre-sent the % forest surrounding each site at a radius of 24 km.For the three sites that were sampled in at least three years, lines represent regression lines depicting the relationship between the number of years between sampling events and F ST Table 1 For each site at which E. championi was sampled, the percent forest at a radius of 24 km surrounding the sampling location, the year and month sampled, the number of individuals sampled, and the expected heterozygosity (He) and nucleotide diversity (Pi) with standard errors Site names correspond to those in Fig. 1 Site Samples were demultiplexed by the UCSF CAT and quality control of the sequencing run was checked using the software FastQC v.0.11.8.The STACKS v. 1.57 pipeline (Catchen et al. 2011(Catchen et al. , 2013) ) was used to identify orthologous loci across individuals.STACKS parameter optimization was performed following (Paris et al. 2017), and the following combination was chosen: m = 3, M = 2, n = 3.We set the maximum observed heterozygosity required to process a locus to 0, as in Alonso-Garcia et al. (2021).Across years, the number of polymorphic loci recovered ranged from 6266-11,864.Pairwise F ST values, expected heterozygosity, and nucleotide diversity (π, Nei et al. 1975) for each site-year combination were calculated using the populations program of the dennovo_map.plSTACKS pipeline, and an output file was generated for the population genetic software program Genepop (Rousset 2008).Relatedness between each pair of individuals within each year was calculated by reading the Genepop file into the R (R Core Team 2021) program Adegenet (Jombart 2008), transforming it into a Genind object, and then using the R package Hierfstat (Goudet 2005) to estimate pairwise kinship coefficients between all pairs of individuals (Goudet et al. 2018).
We estimated the percent forest surrounding each sampling location and between locations using ArcGIS v.2.4 (Esri, Redlands, CA).We obtained forest cover of the study region from the Esri 2020 Land Cover dataset that corresponded to scene 17P (Karra et al. 2021).For each site we quantified the amount of forest cover within a circle of radius 24 km.We estimated the amount of forest between pairs of sampling locations in two ways.First we used ArcGIS to calculate Euclidian (straight-line) geographic distances between all possible site pairs.Euclidian distances reflect the shortest distance between sites, which may be over water.Second, we calculated "broken stick" geographic distances as in Davis et al. (2010), which are the shortest overland distances between two sites.For both types of distances we calculated the amount of forest between each pair of sites by overlaying rectangles of width 1000 m between each pair.Rectangles were centered at each pair of sites and the percent forested area was quantified within that rectangle (Figs.S1 & S2, Online Resource 1).Many sites are located nearby to the coastline.We did not clip the circular or rectangular buffers to the coastline if they extended into the Golfo Dulce or the Pacific Ocean, so water was included as deforested area.This was done to obtain a more informative estimate of the proportion of forest cover relative to other land cover types within each buffer and to reflect Euglossine bee flight path, as there is evidence that some Euglossine species have restricted dispersal over large bodies of water (da Rocha Filho et al. 2013).

Statistical analyses
To determine if gene flow is higher within forested rather than deforested areas, we performed multiple regression on distance matrices (Wang 2013) with F ST among sites as the dependent variable, and the amount of forest between pairs of sites as well as the geographic distance between pairs of sites as independent variables.We did this for each year in which samples were obtained, and for the Las Cruces site we did this for each month within 2017 since the site was visited in two different months.We estimated the overall statistical significance of the multiple comparisons using a modified false discovery rate procedure (Benjamini and Yekutieli 2001;Narum 2006).To determine if genetic diversity varied with the amount of forest surrounding a site, or over time, we performed linear mixed modeling.The dependent variable was expected heterozygosity or nucleotide diversity, percent forest and year were independent variables, and site was a random effect.
We ran two linear mixed models to determine if allele frequencies changed temporally, and to determine if the magnitude of allele frequency changes varied with the amount of forest surrounding fragments.For the first model, F ST between samples obtained at the same site but in different years was the dependent variable, and the percent of forested land at a radius of 24 km and the amount of years that had elapsed between sampling events (0-10 years) were independent variables.We also included whether the samples had been obtained in the same season (fall or spring) as an independent variable, as there was variability in whether we sampled in fall or spring across the years.Site was included as a random effect.
We examined if bees that arrived at baits were more genetically similar if they arrived closer to the same time.The rationale for this investigation was that bees might be spatially structured in kin groups, but when lured to baits the signal might disappear as all bees in a given area congregate at a centralized bait.To examine the possibility of spatial structuring based on relatedness, we conducted a Mantel test implemented using the R package Vegan (Oksanen et al. 2022).Relatedness between pairs of bees was the dependent variable, and the amount of time that elapsed between when those two bees arrived at baits was the independent variable.We did this separately for each site-day combination, and included only sites that had at least ten samples obtained on that day (three sites in 2017, and four sites in 2018, Table S1, Online Resource 1).We predicted an inverse association between relatedness and the difference in arrival time.

Results
The extent to which allele frequencies changed temporally within sites increased over time, such that the amount of time that elapsed between sampling events was positively associated with F ST (Est.= 0.0009, Chisq = 4.2, P = 0.043, Fig. 1b).The temporal change in allele frequencies was lower in sampling locations with greater amounts of forest, such that the amount of forest was negatively associated with F ST (Est.= −0.00052,Chisq = 4.8, P = 0.029, Fig. 1b).Similarity in the time at which bees arrived to sampling locations was not associated with relatedness between those bees (Table S1, Online Resource 1).
Expected heterozygosity decreased over time across sites (Est.= −0.005,Chisq = 6.8,P = 0.0092, Fig. 2).The mean expected heterozygosity fell from 0.217 in 2009 to 0.156 in 2019, representing a 28% decrease.Similarly, nucleotide diversity decreased over time (Est.= −0.005,Chisq = 6.2, P = 0.013).The mean nucleotide diversity fell from 0.225 in 2009 to 0.161 in 2019, also representing a 28% decrease.The amount of forest surrounding sites did not affect expected heterozygosity within sites (Chisq = 1.23,P = 0.27).Neither the amount of forest separating two sites nor the geographic distance between sites was a significant predictor of genetic differentiation between sites after correction for multiple comparisons (Table S2, Online Resource 1).

Discussion
Our analysis of spatial and temporal genetic structure suggests high levels of dispersal of E. championi over an area spanning more than a hundred kilometers in southern Costa Rica.There was high allelic turnover within sites, and there was no evidence of isolation by distance, whereby individuals are more genetically different the farther they are sampled apart.Furthermore, deforested areas do not seem to restrict dispersal.Similarity in the time at which individuals arrived to baits did not predict relatedness between pairs of bees.Despite high dispersal across sites, genetic diversity decreased over time within sites.
The high dispersal we found for E. championi across deforested areas and over 100 km suggests that anthropogenic alteration of land does not restrict movement.This contrasts sharply with what has been found for a number of other animals, such that habitat fragmentation is associated with lower gene flow (Tucker et al. 2018;Cosgrove et al. 2018).Flying organisms have been shown to be less susceptible to restricted movement due to habitat fragmentation than organisms that do not fly (Tucker et al. 2018), but examples of reduced dispersal are still plentiful even for flying organisms (Schtickzelle et al. 2006;Méndez et al. 2011;Stevens et al. 2018).In our case, the reasons for the high gene flow across partially deforested landscapes could include the presence of sufficient forest to provide refuge for travelling individuals.The percent of land that was forested between two pairs of sites averaged 66% (minimum 44%), which may be enough to facilitate movement.It is also possible that E. championi will fly across deforested habitat.Females of some species of Euglossine bees (genus Eulaema) will travel 2.5 km over open water (Janzen 1971), and males of some species (genus Euglossa) will leave forest to visit chemical baits (Milet-Pinheiro and Schlindwein 2005, p.).To our knowledge no studies have attempted to bait E. championi at varying distances from forest edges, but our genetic results suggest it is likely that E. championi will cross deforested habitat when foraging.
The high gene flow observed for E. championi is consistent with previous work examining the genetic structure of E. championi using microsatellites, which found a lack of genetic structure across 80 km and no effect of deforestation on gene flow (Suni and Brosi 2012;Suni et al. 2014).Some studies that have estimated genetic structure using both microsatellites and SNPs have found strong differences in the patterns observed from each marker type (Lozier 2014;Lozier et al. 2016).Differences in patterns resulting from microsatellites versus other markers can be driven by the high levels of polymorphism exhibited by microsatellites, which limit inferences that can be drawn from F ST values (Jost 2008;Meirmans and Hedrick 2011).In our case, the congruence between previous work using microsatellites and our current work using SNPs suggests that the lack of isolation by distance or landscape found previously was not due to low resolution of microsatellite markers and is instead due to a robust pattern of high gene flow across deforested areas.The high gene flow exhibited by E. championi over large and deforested geographic areas is plausible given that Euglossine bees are capable of long-distance movement (Janzen 1971).A mark-recapture study conducted in an area of intact forest observed that two male bees had traveled 45 km and 50 km away from the site at which they were marked, within a period of 17 days (Pokorny et al. 2015).Despite examples of long-range dispersal, males of several species also show high site fidelity in areas of intact habitat and in anthropogenically-altered landscapes.Mark-recapture observations conducted on Barro Colorado Island recorded the same individual bees being sighted over a period of 44 days at fragrance baits (Eltz et al. 1999).In addition, a mark-recapture study conducted in an urban area in Brazil found that 75% of individuals could be recaptured at the same trees for at least 30 days (López-Uribe et al. 2008).To our knowledge there have been no mark-recapture studies of E. championi.
Results from studies that have characterized genetic connectivity for Euglossine bees in fragmented areas suggest that species differ in gene flow across mixed landscapes.Similar to our results, microsatellite genotyping found no significant genetic differentiation between Euglossa dilemma separated by 130 km of deforested habitat (Zimmermann et al. 2011), or Euglossa imperialis separated by 226 km of partially-deforested habitat (Suni 2017).In addition, microsatellite genotyping of E. dilemma and Euglossa viridissima from sites separated by areas that were urban, agricultural, natural, or water revealed that no land types restricted gene flow across 114 km, although genetic diversity was lower on islands for E. dilemma (Soro et al. 2017).In contrast, microsatellite genotyping revealed significant genetic structuring for Eulaema bombiformis across just 14 km (Suni and Brosi 2012).In addition, microsatellite genotyping of Euglossa cordata on islands and the Brazilian mainland found that a population on the island furthest from the mainland (38 km) was genetically differentiated from those on the mainland (Boff et al. 2014).Microsatellite genotyping also found significant genetic structuring between Brazilian island and mainland populations of Euglossa stellfeldi that were 30 km apart (da Rocha Filho et al. 2013).It is possible that the differences among species in their level of genetic structure over non-forested areas can be explained by differences in the degree of specialization on particular floral resources, flight ability, the ability to camouflage outside of forest habitat, or the susceptibility of bees to increased temperature in open areas.Deforestation can select for higher dispersal (Milet-Pinheiro and Schlindwein 2005), so it is possible that the high gene flow reflects individuals traveling far distances when seeking sufficient resources or mates.
Despite high gene flow, our findings suggest that the genetic health of E. championi is decreasing over time.The 28% decrease in genome-wide levels of expected heterozygosity over a period of 11 years is striking, and could greatly decrease the adaptive potential of populations (Hoffmann et al. 2017).This loss of genetic diversity could reflect an overall decrease in the abundance of E. championi across the study area, which would remove alleles from the population.Even with no decrease in abundance, genetic drift could reduce genetic diversity within sites.Due to haplodiploidy, the effective population size of bees is lower than what would be expected from a randomly mating diploid population, so loss of heterozygosity is expected to be higher.Sampling of Euglossine bee communities that took place in 2004 in the area around the Las Cruces Biological Station revealed that E. championi was one of the most abundant species present (Brosi 2009).To our knowledge there have been no systematic assessments of the abundance of E. championi since then, so such monitoring would be worthwhile to understand if abundance is changing over time.
The extreme decreases in genetic diversity over time observed in this study are on the high end of what has been observed for other taxa.A meta-analysis that quantified changes in genetic diversity across organisms since the industrial revolution found an average of a 6% decrease in within-population genetic diversity, with the greatest loss being 28% in species living on islands (Leigh et al. 2019).That meta-analysis included one study that tracked changes in genetic diversity change for a bee species (Bombus veteranus), for which there was a 5% reduction in genetic diversity between 1895 and 1923 in Belgium (Maebe et al. 2013).Other work that investigated temporal changes in genetic diversity revealed on average a decrease in mitochondrial nucleotide diversity between 1990 and 2016 for a suite of insect species (Millette et al. 2020).Interestingly, despite the overall decrease, Millette et al. (2020) found an interaction between landscape type and year, such that there was an overall increase in genetic diversity in areas of intensive land use.In our case we found a consistent decline in genetic diversity across sites and no effect of the amount of forest surrounding sites on the way that genetic diversity changed over time.
We stress that our results do not suggest that deforestation was not a driver of the decline in genetic diversity.We found evidence of higher allelic turnover within sites that were surrounded by less forest, which would be consistent with bees spending less time foraging in these areas.We speculate that sites surrounded by less forest have fewer resources, and represent poorer habitat quality than sites surrounded by more forest.Negative effects of deforestation on genetic diversity have been found for other Euglossine bees (Boff et al. 2014) so it is worth considering the limitations of our study in detecting such a pattern should it exist.In order to understand if deforestation was a driver of the decreases in genetic diversity, it would be useful to consider our entire study area as one data point, and compare our findings to those of future studies that examine dispersal and genetic diversity within forested areas that exist over hundreds of kilometers.
Whether the decreasing genetic diversity documented here will threaten persistence of E. championi populations depends on the extent to which it results in the production of diploid males or compromises the adaptive potential of populations.High diploid male frequencies have been found for some species such as Euglossa annectans, for which more than 15% of males sampled on islands were diploid (Giangarelli et al. 2015).However, estimates of diploid male frequencies have been low for most bee species, including E. championi (Souza et al. 2010), and despite low genetic diversity and island isolation (Boff et al. 2014;Soro et al. 2017).In addition, haplodiploidy has been suggested to render bees less sensitive to genetic impoverishment than diploid species, due to deleterious alleles being purged from populations when they are exposed in haploid males (Packer and Owen 2001).Such purging of maladaptive alleles could facilitate adaptation to new ecological conditions.That said, genetic diversity can be lost quickly via drift and then take thousands of generations to recover (Nei et al. 1975), so conservation monitoring of E. championi might be warranted given ongoing habitat loss.
To our knowledge this is the first study that has quantified genetic structure for Euglossine bees using SNP loci.Our findings are consistent with previous work using microsatellite markers to assess dispersal in E. championi, and they fall in line with the body of evidence that suggests that many species in the genus Euglossa show no evidence of restricted dispersal across hundreds of kilometers.In addition, this is one of the first studies that tracks contemporary changes in genetic diversity within the same geographic area for a bee species.Most conservation genetic studies that estimate genetic diversity and its drivers quantify genetic diversity during one time period.Our monitoring of genetic diversity over more than a decade reveals that genetic diversity can decline sharply despite high gene flow.

Fig. 1
Fig. 1 Panel A: Map of sampling locations in southern Costa Rica in which E. championi was sampled from 2009 to 2019.Site names correspond to Table 1.Panel B: Temporal genetic structure for the six sites between 2009 and 2019, calculated as F ST between samples obtained from the same site across each pair of years.Colors repre-

Fig. 2
Fig. 2 Change in expected heterozygosity within sites over time for E. championi at six sites in southern Costa Rica between 2009 and 2019.Shapes represent sites