Seminatural areas act as reservoirs of genetic diversity for crop pollinators and natural enemies across Europe

Despite increasing recognition of the importance of the multiple dimensions of biodiversity, including functional or genetic diversity as well as species diversity, most conservation studies on ecosystem service‐providing insects focus on simple diversity measures such as species richness and abundance. In contrast, relatively little is known about the genetic diversity and resilience of pollinators or natural enemies of crop pests to population fragmentation and local extinction. The genetic diversity and demographic dynamics of remnant populations of beneficial insects in agricultural areas can be a useful indicator proving additional insights into their conservation status, but this is rarely evaluated. Although gene flow between agricultural and seminatural areas is key to maintaining genetic diversity, its extent and directionality remain largely unexplored. Here, we apply a pan‐European sampling protocol to quantify genetic diversity and structure and assess gene flow between agricultural and nearby seminatural landscapes in populations of two key ecosystem service‐providing insect species, the lady beetle Coccinella septempunctata, an important predator of aphids and other crop pests, and the bee pollinator Andrena flavipes. We show that A. flavipes populations are genetically structured at the European level, whereas populations of C. septempunctata experience widespread gene flow across the continent and lack any defined genetic structure. In both species, we found that there is high genetic connectivity between populations established in croplands and nearby seminatural areas and, as a consequence, they harbor similar levels of genetic diversity. Interestingly, demographic models for some regions support asymmetric gene flow from seminatural areas to nearby agricultural landscapes. Collectively, our study demonstrates how seminatural areas can serve as genetic reservoirs of both pollinators and natural enemies for nearby agricultural landscapes, acting as sources for recurrent recolonization and, potentially, contributing to enhancing ecosystem service and crop production resilience in the longer term.

or natural enemies of crop pests to population fragmentation and local extinction.The genetic diversity and demographic dynamics of remnant populations of beneficial insects in agricultural areas can be a useful indicator proving additional insights into their conservation status, but this is rarely evaluated.
Although gene flow between agricultural and seminatural areas is key to maintaining genetic diversity, its extent and directionality remain largely unexplored.Here, we apply a pan-European sampling protocol to quantify genetic diversity and structure and assess gene flow between agricultural and nearby seminatural landscapes in populations of two key ecosystem service-providing insect species, the lady beetle Coccinella septempunctata, an important predator of aphids and other crop pests, and the bee pollinator Andrena flavipes.We show that A. flavipes populations are genetically structured at the European level, whereas populations of C. septempunctata experience widespread gene flow across the continent and lack any defined genetic structure.In both species, we found that there is high genetic connectivity between populations established in croplands and nearby seminatural areas and, as a consequence, they harbor similar levels of genetic diversity.Interestingly, demographic models for some regions support asymmetric gene flow from seminatural areas to nearby agricultural landscapes.
Collectively, our study demonstrates how seminatural areas can serve as genetic reservoirs of both pollinators and natural enemies for nearby agricultural landscapes, acting as sources for recurrent recolonization and, potentially, contributing to enhancing ecosystem service and crop production resilience in the longer term.

| INTRODUCTION
With more than 80% of the world's land under the direct influence of agricultural systems, conservation strategies that include agricultural landscapes could potentially be important to buffer and connect protected areas, but also in the direct conservation of species that exist there (Garibaldi et al., 2021;Kleijn et al., 2020).Indeed, protecting those species is becoming particularly important in providing ecosystem services such as enhanced crop yield by pollinators (Klein et al., 2007) and control of pest species by natural enemies (Garratt et al., 2011).Advances have been made in transitioning to lower-intensity agriculture and increasingly adopting an ecological intensification approach (Kleijn et al., 2019), with most efforts focused on halting the loss of species.For example, it has been shown that actions such as planting flower strips can increase pollinator richness and abundance (Scheper et al., 2013), and reduced soil tilling regimes improve the survival and numbers of arthropods in agroecosystems (Thorbek & Bilde, 2004).However, in addition to restoring and maintaining species diversity, it is necessary to maintain healthy and viable populations, that is, populations with a sufficient number of individuals to maintain high genetic diversity, which is key to responding to environmental change, and therefore for the long-term persistence of natural populations (Falconer & Mackay, 1996;Frankham, 1995;Spielman et al., 2004).This calls for integrating knowledge of population genetic diversity in the management and conservation of biodiversity (Kardos et al., 2021;Willi et al., 2022).
Although it is increasingly acknowledged that it is important to consider multiple dimensions of biodiversity in conservation, most studies to date focus on simple diversity metrics at the level of species, such as richness and abundance (Cerini et al., 2023;Kardos et al., 2021;Willi et al., 2022).Although species diversity metrics have been core for evaluating the effectiveness of conservation and restoration measures, they cannot provide insights into some of the more complex mechanisms underlying biodiversity responses.Biodiversity comprises many levels of organization, from community composition and functional diversity to species metapopulation dynamics and genetic diversity (e.g., United Nations Convention on Biological Diversity; www.cbd.int).In particular, genetic diversity is an often-neglected aspect with far-reaching implications for long-term conservation (Kardos et al., 2021;e.g., Saccheri et al., 1998;Spielman et al., 2004).More genetically diverse populations may better withstand global change pressures, as these populations have a larger genetic pool from which to adapt to new environmental conditions (Frankham, 2005;Kardos et al., 2021).However, the impact of land-use changes on the genetic diversity of arthropods and the connectivity between populations from seminatural and agricultural landscapes in insect species that provide important ecosystem functions and services are still poorly known (Gonz alez et al., 2020;Keyghobadi, 2007;Schlaepfer et al., 2018;Webster et al., 2023).
Insects and other arthropods providing ecosystem services are often highly mobile organisms that disperse and forage across habitats in heterogeneous landscapes (e.g., Jaffé et al., 2016;Lecompte et al., 2016).This implies that their persistence and conservation does not only rely on local habitat quality but also on the landscape-level habitat composition and configuration, in particular, connectivity between suitable habitats (Kremen et al., 2007;Martin et al., 2019).A generally observed pattern is that croplands surrounded by a larger proportion of natural habitats tend to support relatively higher beneficial insect diversity and abundance, which facilitates ecosystem services such as pollination and pest control (Dainese et al., 2019).However, most insect conservation studies focus on a single habitat type (either agricultural or seminatural habitats) and rarely integrate different habitats within the same landscape (Marini et al., 2019).In fact, our understanding of how sourcesink metapopulation dynamics work in agricultural systems is still very limited (e.g., Ballare & Jha, 2021;Dreier et al., 2014;Jha, 2015;Ortego, Aguirre, et al., 2015).Many species are able to maintain sustainable populations in agricultural habitats, but other species may not be able to persist in simplified agricultural habitats without continuous immigration from nearby natural or seminatural areas, that is, the agricultural habitats act as population sinks and the (semi-)natural areas as population sources (e.g., Öckinger & Smith, 2007aÖckinger & Smith, , 2007b)).Either way, the benefits of having bidirectional gene flow between natural and agricultural areas can be important, even when populations would be able to persist in isolation; for example, this would promote long-term persistence of species if one area is subject to a major disturbance, as the remaining populations in neighboring habitats can repopulate and contribute to gene flow toward the disturbed area.However, in what way natural areas contribute to the genetic diversity of ecosystem service providing insect populations in nearby agricultural landscapes and vice versa has, to our knowledge, not yet been explored (Webster et al., 2023).
To improve our understanding of the importance of natural habitats for the conservation of population genetic diversity of key ecosystem service-providing insects, we used a pan-European sampling design to characterize gene flow dynamics between populations of the pest control agent Coccinella septempunctata (Coleoptera: Coccinellidae) and the bee pollinator Andrena flavipes (Hymenoptera: Andrenidae) in agricultural areas and nearby seminatural landscapes.These are two of the most common ecosystem service providers in Europe, with A. flavipes ranking 7th in a recent review of the most common bee crop pollinators in Europe (Kleijn et al., 2015) whereas C. septempunctata is a widespread predator of aphids (Jervis & Kidd, 1996;New, 1991), a ubiquitous pest of several crops.To this end, we sampled and genotyped populations of the two focal species in agricultural and nearby seminatural areas across Europe (i.e., a pairedreplicated approach), quantified their respective levels of population genetic diversity and structure, and tested alternative models of gene flow.Specifically, we first (i) tested the hypothesis that populations in seminatural areas sustain larger effective population sizes and, thus, present higher levels of genetic diversity, than those established in agricultural landscapes.We then quantified spatial patterns of genetic structure across the different European regions, and (ii) tested the hypothesis that populations of mobile beneficial insects from seminatural and agricultural landscapes regularly exchange gene flow and show limited genetic structure (e.g., Jaffé et al., 2016;Lecompte et al., 2016).Finally, we used a coalescent-based model-testing approach to evaluate alternative demographic scenarios and (iii) test the hypothesis that seminatural areas play a role as sources of gene flow for populations established in agricultural landscapes (Figure 1).Our study demonstrates how seminatural areas can serve as genetic reservoirs of both pollinators and natural enemies for nearby agricultural landscapes, acting as sources for recurrent recolonization and, potentially, contributing to enhance ecosystem services and crop production resilience in the longer term.These results, general across a network of European sites, directly appeal to practitioners and policy makers-for example, through the European Union's Common Agricultural Policy (CAP)-and emphasize the importance of implementing management actions aimed at preserving or restoring seminatural habitats for the maintenance of source populations of ecosystem service-providing insects in agricultural landscapes.

| Population sampling
In spring and summer 2021, we sampled pairs of populations of C. septempunctata and A. flavipes in agricultural landscapes (i.e., dominated by agriculturally managed land for the last 100 years) and nearby seminatural areas (i.e., dominated by seminatural habitat that was not under agricultural management at least in the last 100 years) separated by <25 km (mean: 11.9 km; range: 4.6-24.0km; Table 1; Figure 2; Figure S1).The selected species are representative of generalist ecosystem service providers within their guild and are often locally abundant and widely distributed.In seminatural areas, specimens were collected in patches of native vegetation, whereas specimens from agricultural areas were collected within crop fields and along paths and field margins.We targeted to sample 8-10 specimens per species and site area and aimed to replicate this sampling scheme in seven European countries: Estonia, Denmark, the Netherlands, Great Britain, Hungary, Switzerland, and Spain (Figure 2; Figure S1).However, low population densities in some areas often reduced our sample sizes and one of the two focal taxa could not be collected in some countries at a minimum sample size of n = 5 in agricultural and/or seminatural sites (Table 1).Consequently, our sampling scheme could be replicated in six countries for C. septempunctata and four countries for A. flavipes and final datasets included 5-10 genotyped specimens per species, country, and area (mean = 8; mode = 8; Table 1).Our pan-European sampling scheme, and the countries selected for our study, frames within the context of a multipartner project funded by the European Union's Horizon 2020 Research and Innovation Programme H2020.We spread the sampling across each site (<500 m of radius around sampling coordinates; Table 1) to avoid collecting closely related individuals (i.e., siblings or half-siblings).Specimens were collected by hand (C.septempunctata) or sweep-netting (A.flavipes) and preserved at À20 C in 1500 μL ethanol 96% until DNA extraction.We registered spatial coordinates of each sampling plot using a Global Positioning System (GPS).Further details on sampling locations are provided in Table 1.

| Genomic library preparation and processing
We extracted and purified DNA from each specimen using NucleoSpin Tissue kits (Macherey-Nagel, Düren, Germany).We processed genomic DNA into different genomic libraries using the double-digestion restrictionfragment-based procedure (ddRAD-seq) described in Peterson et al. (2012).In brief, we digested DNA with the restriction enzymes MseI and EcoRI (New England Biolabs, Ipswich, MA, USA) and ligated Illumina adaptors including unique 7-bp barcodes to the digested fragments of each individual.We pooled ligation products, sizeselected them between 350 and 450 bp with a Pippin Prep machine (Sage Science, Beverly, MA, USA), and amplified the fragments by PCR with 12 cycles using the iProofTM high-fidelity DNA polymerase (BIO-RAD, Veenendaal, The Netherlands).Single-read 201-bp sequencing was performed on an Illumina NovaSeq6000 platform.We used the different programs distributed as part of the STACKS v. 1.35 pipeline (Catchen et al., 2013) to filter and assemble our sequences into de novo loci, call genotypes, calculate genetic diversity statistics, and export input files for all downstream analyses (Methods S1).Unless otherwise indicated, all downstream analyses were performed using datasets including only the first single nucleotide polymorphism (SNP) per RAD locus (180-bp length, see Methods S1) and retaining loci with a ?(a) minimum stack depth ≥5 (m = 5) and that were represented in at least 75% of the populations (p = 9 for C. septempunctata and p = 6 for A. flavipes) and 50% of the individuals within each population (r = 0.5).

| Quantifying genetic structure
We calculated genetic differentiation between each pair of populations using the Weir and Cockerham weighted fixation index (F ST ) (Weir & Cockerham, 1984) as implemented in ARLEQUIN v. 3.5 (Excoffier & Lischer, 2010).We determined statistical significance with Fisher's exact tests after 10,000 permutations, applying a false discovery rate (FDR) adjustment (5%, q < 0.05) to control for multiple tests.We performed Mantel tests to analyze the association between genetic differentiation (F ST ) and geographical distances among populations of each focal species using ZT software with 10,000 permutations (Bonnet & van de Peer, 2002).We quantified population genetic structure and admixture using the Bayesian Markov chain Monte Carlo clustering method implemented in the program STRUCTURE v. 2.3.3 (Pritchard et al., 2000).We conducted STRUCTURE analyses hierarchically, initially analyzing data from all populations of each taxon jointly and, subsequently, running independent analyses for subsets of populations assigned to the same genetic cluster in the previous hierarchical level analysis (e.g., Gonz alez-Serna et al., 2019).To make analyses computationally tractable, we ran STRUCTURE using a random subset of 10,000 SNPs.We ran STRUCTURE analyses assuming correlated allele frequencies and admixture and without using prior population information (Hubisz et al., 2009).We conducted 15 independent runs for each value of K (from K = 1 to K = 10) to estimate the most likely number of genetic clusters with 200,000 MCMC cycles, following a burn-in step of 100,000 iterations.We retained the ten runs having the highest likelihood for each value of K and determined the number of genetic clusters that best describes our data according to log probabilities of the data (LnPr(XjK)) (Pritchard et al., 2000) and the ΔK method (Evanno et al., 2005), as implemented in STRUCTURE HARVESTER (Earl & vonHoldt, 2012).We used CLUMPP v. 1.1.2and the Greedy algorithm to align multiple runs of STRUCTURE for the same K value (Jakobsson & Rosenberg, 2007) and DISTRUCT v. 1.1 (Rosenberg, 2004) to visualize the individuals' probabilities of population membership in bar plots.Complementarily, we used the same datasets of 10,000 random SNPs to perform principal component analyses (PCA) as implemented in the R v. 4.2.1 (R Core Team, 2023) package "adegenet" (Jombart, 2008).Before running PCAs, we replaced missing data by the mean allele frequency of the corresponding locus estimated across all samples (Jombart, 2008).Results from PCAs were visualized using the R package "scatterplot3d" (Ligges & Mächler, 2003).

| Testing alternative models of gene flow
We used the simulation-based approach implemented in FASTSIMCOAL2 to evaluate alternative migration models and test the hypothesis of asymmetric gene flow between populations from agricultural landscapes and nearby seminatural areas (Excoffier et al., 2013; Figure 1a).Specifically, we tested four alternative migration models for each taxon (C.septempunctata and A. flavipes) and spatial replicate  1b).We used the site frequency spectrum (SFS) and FASTSIMCOAL2 to estimate the composite likelihood of the observed data given a specified model (Excoffier et al., 2013).For each taxon and pair of populations, we calculated a folded joint SFS considering a single SNP per locus.Because we did not include invariable sites in the SFS, we fixed the contemporary effective population size for one of the demes (θ AGR in all cases; see Figure 1b) to enable the estimation of other parameters in FASTSIMCOAL2 (Excoffier et al., 2013).The effective population size fixed in the model was calculated from the level of nucleotide diversity (π) and estimates of mutation rates per site per generation (μ).Nucleotide diversity (π) was estimated from polymorphic and nonpolymorphic loci using STACKS (Table 1).We used the mutation rate per site per generation (2.8 Â 10 À9 ) estimated for Drosophila melanogaster (Keightley et al., 2014), which is similar to the spontaneous mutation rate estimated for the butterfly Heliconius melpomene (2.9 Â 10 À9 ; Keightley et al., 2015).To remove all missing data for the calculation of the SFS, minimize errors with allele frequency estimates, and maximize the number of retained SNPs, each population group was downsampled using a Python script written by Qixin He and available on Dryad (Papadopoulou & Knowles, 2015).Note that this downsampling procedure necessarily results in the number of SNPs retained in the SFS differ across each pair of analyzed populations, as this depends on the number of available samples and the pattern and amount of missing data in each specific dataset (Table 2; e.g., Papadopoulou & Knowles, 2015).
Each model was run 100 independent times considering 100,000-250,000 simulations for the calculation of the composite likelihood, 10-40 expectation-conditional maximization (ECM) cycles, and a stopping criterion of 0.001 (Excoffier et al., 2013;Excoffier & Foll, 2011).We used an information-theoretic model selection approach based on the Akaike's information criterion (AIC) to compare the set of candidate models and identify the one that is best supported by the data (Burnham & Anderson, 2002;e.g., Ortego et al., 2021;Thome & Carstens, 2016).After the maximum likelihood was estimated for each model in every replicate, we calculated the AIC scores as detailed in Thome and Carstens (2016).AIC values for each model were rescaled (ΔAIC), calculating the difference between the AIC value of each model and the minimum AIC obtained among all competing models (i.e., the best model has ΔAIC = 0).Point estimates of the different demographic parameters for the best supported model were T A B L E 2 Alternative demographic models (detailed in Figure 1b) tested using FASTSIMCOAL2 for pairs of populations of Coccinella septempunctata and Andrena flavipes sampled in agricultural and nearby seminatural landscapes from different countries.selected from the run with the highest maximum composite likelihood.Finally, we calculated confidence intervals of parameter estimates (based on the percentile method; e.g., de Manuel et al., 2016) from 100 parametric bootstrap replicates by simulating SFS from the maximum composite likelihood estimates and re-estimating parameters each time (Excoffier et al., 2013).

| Samples and genomic data
Overall, we obtained samples of C. septempunctata and A. flavipes from six and four countries, respectively, with three countries sampling both species (Table 1).

| Quantifying genetic structure
In C. septempunctata, all pairwise F ST values were nonsignificantly different from zero (Table S1).In contrast, pairwise F ST values in A. flavipes ranged between 0 and 0.444, and several pairs of populations showed significant genetic differentiation (Table S2).Except for comparisons involving populations from the Netherlands and Switzerland, all pairs of populations of A. flavipes from different countries were significantly differentiated (Table S2).Within most countries, populations of A. flavipes from agricultural landscapes and nearby seminatural areas did not show significant genetic differentiation (Table S2).The only exception was the comparison involving the two populations from Denmark, albeit with a small F ST value (F ST = 0.03; Table S2).Simple Mantel tests revealed that genetic differentiation (F ST ) is positively correlated with geographical distances among populations in A. flavipes (r = 0.96, p < .001,n = 8) but not in C. septempunctata (r = 0.00, p = 1.000, n = 12).
STRUCTURE analyses for C. septempunctata showed that LnPr(XjK) peaked at K = 2 and decreased for higher K-values (Figure S4a).For K = 2, all individuals and populations of C. septempunctata presented a very low probability of assignment (q < 0.05) to one of the two genetic clusters (i.e., a "fictive" or "ghost" cluster sensu Guillot et al., 2005; see also Chen et al., 2007), which indicates a lack of genetic structure across analyzed populations of this taxon (Figure 3a).Clustering solutions for higher K-values (from K = 3 to K = 5) did not reveal either any genetic structure in C. septempunctata (Figure S5a).STRUCTURE analyses for A. flavipes considering all populations identified the most likely number of clusters as K = 2 according to the ΔK criterion and LnPr(XjK) reached a plateau at the same K value (Figure S4b).For K = 2, the two genetic clusters separated the populations of A. flavipes from Spain from the rest of European populations (Figure 3b).STRUC-TURE analyses at a lower hierarchical level, and excluding Spanish populations, identified the most likely number of clusters as K = 2 according to the ΔK criterion and LnPr(XjK) reached a plateau at the same K value (Figure S4c).In this case, the two genetic clusters separated the populations from Denmark from the populations from the Netherlands and Switzerland, which grouped together (Figure 3b).The new genetic clusters that arise for higher K-values (from K = 3 to K = 5) were either similarly represented in all populations of A. flavipes assigned to the same genetic cluster for K = 2 or were circumscribed to certain populations but showed very low probabilities assignment in all individuals (q < 0.15; Figure S5b,c).Principal component analyses (PCA) confirmed the results yielded by STRUCTURE for the two taxa (Figure 4).For C. septempunctata, PCA did not reveal any genetic clustering and individuals collected in different countries and sites were interspersed across the three first principal components (Figure 4a).Analyses are based on a random subset of 10,000 single nucleotide polymorphisms (SNPs).Population codes as described in Table 1.
the Netherlands and Switzerland, which in this case tended to segregate along PC2 (Figure 4c).

| Testing alternative models of gene flow
Demographic analyses in FASTSIMCOAL2 supported unidirectional gene flow from populations in seminatural areas to those in agricultural landscapes (Model D; Figure 1b) in five out of six C. septempunctata (EE, DK, NL, HU, and CH) and two out of four A. flavipes (NL and ES) pairwise comparisons (Table 2).In the rest of the comparisons (C.septempunctata in GB and A. flavipes in DK and CH), the best supported model was a scenario of unidirectional gene flow from populations in agricultural landscapes to those in nearby seminatural areas (Model C; Table 2).In the case of populations of C. septempunctata from Great Britain (GB), the second-best supported model (Model D) was statistically indistinguishable from the best supported model (Model C) (ΔAIC = 0.79), indicating that the two alternative models of unidirectional gene flow are similarly probable (Table 2).Estimates of migration rates between populations were fairly similar across all pairwise comparisons, ranging from 2.35 Â 10 À4 to 5.79 Â 10 À4 in C. septempunctata and from 2.27 Â 10 À4 to 5.67 Â 10 À4 in A. flavipes (Table 3).In both C. septempunctata and A. flavipes, estimates of gene flow from populations established in agricultural areas to those in nearby seminatural areas (i.e., Model C) were smaller than estimates obtained when the most supported model was the scenario of unidirectional gene flow in the opposite direction (i.e., Model D) (Table 3).

| DISCUSSION
Our data support contrasting patterns of genetic structure and gene flow in two key ecosystem service providing insect species, the natural enemy C. septempunctata and the pollinator A. flavipes.A panmictic structure characterizes populations of C. septempunctata across both local and regional scales.Conversely, populations of A. flavipes present a marked structure at a European scale, but a very limited or null genetic differentiation between nearby agricultural and seminatural landscapes within each country.We found no support for the hypothesis of higher levels of genetic diversity in populations from seminatural areas with respect to those established in agricultural landscapes.However, demographic model testing and coalescent-based estimates of migration rates lent partial support to the hypothesis of increased gene flow from seminatural to agricultural landscapes (Figure 1a).For C. septempunctata, models considering anisotropic gene flow in the predicted direction received comparatively higher statistical support than the rest of the scenarios across most landscape replicates (i.e., countries), while the evidence on this respect for A. flavipes was limited to 50% of the regions.Collectively, these results emphasize the occurrence of gene flow from seminatural to cultivated areas and confirm the importance of preserving seminatural habitats at the landscape level for the maintenance of source populations of ecosystem service-providing insects (Garibaldi et al., 2021).
4.1 | Genetic diversity and structure at regional and landscape scales Analyses of genetic differentiation (F ST ) and structure (STRUC-TURE and PCAs) revealed marked differences in the spatial distribution of genomic variation in the two focal taxa (Figures 3 and 4).Populations of C. septempunctata were panmictic and presented no evidence for genetic differentiation at both European and landscape scales, a phenomenon likely explained by the long-distance migratory behavior linked to seasonal aggregations characterizing many coccinellids (Hagen, 1962).Long-distance dispersal in this taxon is exemplified by unrestricted gene flow between British and mainland European populations, which indicates the limited effects of seawater as a barrier to dispersal (i.e., The English Channel is 34 km at its narrowest).These results are in line with the findings of a previous microsatellitebased study covering the entire distribution of the species (Lecompte et al., 2016).This study revealed a lack of genetic structure, very limited differentiation (F ST < 0.05), and a weak pattern of isolation-by-distance across populations sampled from Portugal to India.Only eastern Asia (China, Japan) and some North African populations (Algeria) were assigned to different genetic clusters and presented a significant genetic differentiation (F ST = 0.10-0.30),which suggests either long-term isolation of these peripheral populations (e.g., in different glacial refugia) or the evolution of reproductive isolation and an incipient process of cryptic speciation (Lecompte et al., 2016).In contrast, analyses for A. flavipes revealed a marked genetic structure and differentiation (F ST = 0.10-0.44)at a European scale, with more peripheral populations from Jutland (DK) and Iberian (ES) peninsulas forming distinct genetic clusters with limited gene flow between them as well as with the rest of analyzed populations (Figures 3 and 4).Likewise, although central European populations from NL and CH (separated by >400 km) clustered together in STRUCTURE analyses (Figure 2; Figure S5) and showed nonsignificant genetic differentiation (F ST = 0; Table S2), they tended to segregate in PCAs (Figure 4).These results reveal more limited gene flow in A. flavipes than in C. septempunctata at a continental scale, which has likely resulted in the development of genetic structure across populations from most studied regions.In contrast, populations from agricultural and nearby seminatural landscapes within countries presented weak genetic differentiation in A. flavipes (F ST < 0.03).Only the two sites within Denmark showed significant differentiation, albeit with a very low estimate (F ST = 0.03).Remarkably, populations of A. flavipes from Denmark showed the lowest levels of genetic diversity compared to the rest of European populations, which suggests that demographic stochasticity, founder effects, genetic drift, and low effective population sizes characterizing leading-edge peripheral populations could have contributed to fine spatial-scale (<5 km) genetic structure between populations from agricultural and nearby seminatural landscapes within Jutland Peninsula (Hampe & Petit, 2005;Lira-Noriega & Manthey, 2014) Collectively, these results are in line with the findings of spatially explicit landscape genetic studies suggesting a limited impact of habitat fragmentation and structure on the distribution of spatial patterns of genetic variation in different pollinators (Barbosa et al., 2022;Exeler et al., 2010;Jackson et al., 2018;Jha & Kremen, 2013).Aligned with inferences from analyses of genetic structure, testing of alternative demographic models strongly rejected the scenario of strict isolation (i.e., total lack of gene flow; Model A) between populations from agricultural and natural areas in the two focal taxa (Table 2).Connectivity between nearby areas of different land use is expected to contribute to prevent genetic drift and loss of genetic diversity through time, which can explain similar levels of genetic variation across populations from natural and agricultural landscapes (Table 1).

| Anisotropic gene flow between agricultural and seminatural landscapes
Asymmetric gene flow has been previously documented linked to passive dispersal in the direction of prevailing winds or ocean currents (e.g., Bertola et al., 2020;Kling & Ackerly, 2021) and associated to differences in the production of propagules in the context of central-peripheral and source-sink metapopulation dynamics (Kirkpatrick & Barton, 1997;e.g., Hauser et al., 2019;Holliday et al., 2012).Our model-based approach provided mixed support to the hypothesis of directional gene flow from populations in seminatural areas to those sustained in nearby agricultural landscapes (Figure 1).Coalescent-based model testing in the natural enemy C. septempunctata largely supported unidirectional gene flow in the predicted direction (Model D), suggesting that seminatural areas genetically subsidize conspecific populations established in nearby croplands.These results are in line with previous capture-mark-recapture and observational studies pointing to the role of seminatural landscapes as population sources for beneficial insects (Öckinger & Smith, 2007a(Öckinger & Smith, , 2007b)).Although analyses in A. flavipes consistently supported models fitting asymmetric migration rates over those considering symmetric gene exchange, the scenario of unidirectional gene flow from seminatural to agricultural landscapes was the most supported one in only half of our spatial replicates.Remarkably, the bidirectional migration model received low statistical support in all spatial replicates for both C. septempunctata and A. flavipes, suggesting that asymmetric gene flow and source-sink dynamics govern the demography of the focal taxa in natural-agricultural mosaic landscapes (Hanski, 1998).Under this scenario, source populations with a net growth supply individuals to sink populations, reducing the chances of extinction of small populations and contributing to increasing the stability and maintaining the levels of genetic diversity of the whole metapopulation (Howe et al., 1991;e.g., Ortego et al., 2008;Saccheri et al., 1998).It must be noted, however, that source-sink dynamics in highly anthropized landscapes are only expected to be established among generalist taxa with relatively high dispersal ability (e.g., Keller et al., 2013;Ortego, García-Navas, et al., 2015), such as the two widespread ecosystem service-providing species analyzed in our study (Jervis & Kidd, 1996;Kleijn et al., 2015;New, 1991).Rare species with more limited dispersal or increased sensitivity to population subdivision (e.g., pools of regionally threatened species) will likely experience dispersal depression, population isolation, reduced levels of genetic diversity, and increased chances of local extinction under the severe fragmentation of natural habitats that characterize most agricultural landscapes (e.g., Jauker et al., 2009;Ortego, Aguirre, et al., 2015;Schtickzelle et al., 2006).

| CONCLUSIONS
Collectively, our results support high genetic connectivity between populations established in agricultural and nearby seminatural landscapes, indicating that (i) the latter may serve-at least in some cases-as reservoirs of both natural enemies and pollinators for surrounding agricultural areas, (ii) act as sources for recurrent recolonization, and (iii) potentially contribute to enhance ecosystem services and crop production.This highlights the importance of conserving complex agricultural systems which embed seminatural areas within productive regions, as those can boost gene flow exchange and contribute to demographic reinforcement of populations of ecosystem service-providing insects and increase the chances that they persist and maintain their important functional role through time.
Schematic illustrating the hypothesis of unidirectional gene flow from seminatural to nearby agricultural landscapes for populations of the pest predator Coccinella septempunctata (shown in panel a) and the pollinator Andrena flavipes.Panel (b) shows the four alternative gene flow models tested using FASTSIMCOAL2.Parameters include mutation-scaled ancestral effective population sizes (θ ANC ), contemporary effective population sizes for populations from agricultural (θ AGR ) and seminatural (θ NAT ) landscapes, timing of population split (T DIV ), and symmetric (m) and unidirectional (m AN and m NA ) migration rates.
(i.e., country), including a null-model considering no gene flow between populations from agricultural and seminatural areas (Model A), a model of symmetric gene flow (Model B), and two models of unidirectional gene flow, one fitting unidirectional gene flow from agricultural to seminatural areas (Model C) and another exclusively considering unidirectional gene flow in the opposite direction (Model D) (Figure Map showing the location of study sites in seven European countries, each including an agricultural (brown dots) and a nearby seminatural (green triangles) plot.Geographical coordinates are given in Table1and satellite images for each sampling site are presented in FigureS1.CH, Switzerland; DK, Denmark; EE, Estonia; ES, Spain; GB, Great Britain; HU, Hungary; NL, Netherlands.

F
For A. flavipes, PCA separated Spanish populations from the rest of European populations along the first principal component (PC1) (Figure 4b,c).The second principal component (PC2) separated the populations from Denmark from those sampled in the Netherlands and Switzerland, which also tended to segregate along the third principal component (PC3) (Figure 4b).Analyses for A. flavipes excluding Spanish populations provided similar results, with PC1 separating populations from Denmark from those sampled in Results of genetic assignments based on the Bayesian method implemented in the program STRUCTURE for populations of (a) Coccinella septempunctata and (b) Andrena flavipes sampled in agricultural (agr.) and nearby seminatural (nat.)landscapes from different countries.STRUCTURE analyses for A. flavipes were performed for all populations jointly (top) and excluding populations from Spain (bottom).Analyses are based on a random subset of 10,000 single nucleotide polymorphisms (SNPs).Each individual is represented by a vertical bar, which is partitioned into K colored segments showing the individual's probability of belonging to the cluster with that color.Thin vertical black lines separate individuals from different populations.Population codes as described in Table 1..)NL (nat.)DK (nat.)EE (nat.)HU (nat.)CH (nat.)ES (nat.)I G U R E 4 Principal component analyses (PCA) of genetic variation for populations of (a) Coccinella septempunctata and (b, c) Andrena flavipes sampled in agricultural (agr.; dots) and nearby seminatural (nat.; triangles) landscapes from different countries.Principal component analyses for A. flavipes were performed for (b) all populations jointly and (c) excluding populations from Spain.
Parameters inferred from coalescent simulations with FASTSIMCOAL2 under the most supported demographic model (Model C or Model D, as indicated in Table2) for pairs of populations of Coccinella septempunctata and Andrena flavipes sampled in agricultural and nearby seminatural landscapes from different countries.