Contrasting influence of seascape, space and marine reserves on genomic variation in multiple species

factors, space and protection status in multiple species. It stresses the need for a multi-species approach to inform marine conservation planning, opening up new perspectives at the community level.


Introduction
No-take marine protected areas (hereafter reserves) have been established around the world (Costello and Ballantine 2015) to protect marine biodiversity.These reserves are an essential tool for the conservation and management of marine resources and play a key role in mitigating the impacts of climate change on marine ecosystems (Gell and Roberts 2003, Watson et al. 2014, Roberts et al. 2017, Sala et al. 2021).Beyond individual reserves, networks of marine reserves have been advocated to tackle marine conservation and management challenges at the regional scale (Gaines et al. 2010).Although reserve design is rarely oriented towards this conservation objective (Jenkins and Stevens 2018, but see Xuereb et al. 2019), it is often assumed that reserve networks also contribute to maintain the genetic diversity of marine populations because they contribute to increase effective population size and stepping-stone dispersal across the whole network, and act as refuges for colonizing and non-migratory species (Roberts et al. 2017).The absence of a genetic diversity perspective for the design of marine reserve networks often relates to the lack of genetic data, but also to a knowledge gap on the importance of this level of diversity for conservation and management (Jenkins andStevens 2018, Andrello et al. 2022).This is unfortunate since intraspecific genetic diversity is the most fundamental level of biodiversity (May 1994), providing the potential for populations to adapt to environmental changes (Waldvogel et al. 2020, Aguirre-Liguori et al. 2021) and, in the marine realm, is directly linked to ecosystem resilience (Hughes andStachowicz 2004, Reusch et al. 2005).
Understanding the extent and drivers of genetic variation within a reserve network can contribute to identify ecologically similar and genetically connected areas that should be managed jointly (Xuereb et al. 2019), pinpoint areas with singular genetic variation (Munguía-Vega et al. 2015) or capture genetic patterns at the scale of the entire network.To date, few studies have dissected the drivers of genetic variation between marine reserves and their surroundings (Pujolar et al. 2013, Munguía-Vega et al. 2015, Sinclair-Waters et al. 2018).These studies are also relatively limited in terms of numbers of individuals and genetic markers considered (Pujolar et al. 2013, Beger et al. 2014, Munguía-Vega et al. 2015, Sahyoun et al. 2016), which reduces the power to detect weak albeit significant genetic structure in species with high gene flow (Willing et al. 2012, Benestan et al. 2015).Studies also often focus on single-species genetic patterns, which allows conservation issues that are specific to these species to be addressed (Carreras et al. 2017, Sinclair-Waters et al. 2018), but has limited power to detect overall landscape effects (Thomassen et al. 2011, Manel andHolderegger 2013).Considering multiple species allows the capture of seascape effects through congruent patterns of genetic diversity among species (Manel et al. 2020), or between genetic and species diversity (Donati et al. 2021).
To provide a better understanding of the drivers of genetic variation in the marine realm, comparative genomics approaches have been touted as the next breakthrough in this field (Gagnaire 2020).A rapidly evolving literature is demonstrating the striking parallelism of large-scale ecological and evolutionary patterns in a wide range of taxa (Stanley et al. 2018, Manel et al. 2020), providing insights into the macroecological determinants of genetic diversity (Romiguier et al. 2014, Leigh et al. 2021) and marine connectivity (Nielsen et al. 2020).Nevertheless, the role of these determinants remains to be assessed at a fine scale and within a conservation context.Furthermore, reserve networks are home to many species that display a complex mosaic of life-history traits.Consequently, reserve networks should be evaluated on the basis of multi-species rather than single-species information (Gaines et al. 2010, Magris et al. 2016, D'Aloia et al. 2017).Recent empirical studies have shown that genetic variation of co-occurring species with a spectrum of reproductive strategies may present different associations with the same environmental variables (Gamboa andWatanabe 2019, Nielsen et al. 2020).A comparative genomics approach is therefore essential for long-term and holistic management, i.e. an economically, socially and ecologically sound management that reflects the patterns and processes that characterize reserve networks from genes to seascape.
Here, we explore patterns of genomic variation in four abundant and broadly co-distributed marine species within a network of eight reserves along 950 km of coast in the northwestern Mediterranean.The Mediterranean Sea holds an exceptional biodiversity (Coll et al. 2010) and is highly threatened by both local and global human impacts (Lejeusne et al. 2010, Micheli et al. 2013).However, reserves in the Mediterranean Sea were often not only established based on conservation criteria but also using economic, political and sociological criteria.This requires a balance between provisioning proteins, local development (tourism) and protecting biodiversity (Cadoret 2021).With the recommendation of the Convention on Biological Diversity to extend the surface of reserves to 30% by 2030 (O'Leary et al. 2016), guidelines for the establishment of new reserves or the expansion of existing ones based on scientific criteria are urgently needed.Considering a network of eight reserves as a model system, we aim to evaluate the capacity of each reserve and the entire network to maintain genetic diversity.The size of the focal reserves ranges between 6 and 120 km 2 (average 52 km 2 ).We hypothesized that 1) reserves contribute to effectively preserve genetic diversity, resulting in greater genetic diversity within their boundaries, and 2) seascape variables, space and protection status influence differently the genetic make-up of the four marine species harboring different life-history traits.To address these hypotheses, we sought to 1) test for differences in genetic diversity within versus outside reserves across the whole reserve network, and 2) document the influence of seascape, space and reserve-induced protection on genetic variation of the four species.

Study area and sampling
The study focuses on four species: the white seabream (Diplodus sargus, n = 276), the striped red mullet (Mullus surmuletus, n = 311), the European spiny lobster (Palinurus elephas, n = 242) and the comber (Serranus cabrilla, n = 468) sampled in the northwestern Mediterranean Sea (Fig. 1).Individuals and GPS sampling locations were provided by artisanal fishermen between June and December 2017.We sampled a total of 139, 181, 124 and 178 geographical positions for the white seabream, the striped red mullet, the European spiny lobster and the comber, respectively.The sampling sites were separated by < 1 to > 1000 km and we sampled between 1 and 31 individuals per site.Samples are widely distributed within and between the eight reserves (Supporting information), in a region that is characterized by natural gradients in salinity, temperature and chlorophyll (Supporting information), giving an appropriate opportunity for disentangling the effect of the seascape, space and protection status on genome-wide variation.Previous analyses indicate that there is no genetic structure at the scale of the whole study area in D. sargus, M. surmuletus and P. elephas, and weak structure (F ST = 0.021) between northern and southern S. cabrilla populations (Benestan et al. 2021).Thus, the system is characterized by generally high levels of genetic connectivity.

Characterizing seascape variables and space
Sea surface temperature (SST) and salinity (SSS) were retrieved from http://marine.copernicus.eu/as the monthly minimum, maximum and average values per year and per reproductive period of each species, respectively, for the period 1987-2015.Salinity data were collected both as surface and as water-column values.Surface and benthic chlorophyll α (Chl) was retrieved from the Bio-ORACLE database (Tyberghein et al. 2012), over the same 1987-2015 period, with the same frequency (monthly minimum, maximum and average) and at a resolution of 5 arc-minutes within a 1 km radius of the sample location using the 'raster' package in R (Hijmans et al. 2013).Chlorophyll is a proxy of biological productivity that decreases from north to south and west to east (Supporting information) and is inversely related to the increase in temperature and salinity in the Mediterranean Sea (Coll et al. 2010).We performed a principal component analysis (PCA) on 24 environmental variables (Supporting information) collected per geographical position for all individuals that includes SST, SSS and Chl.We used the dudi.pca function in the R package 'ade4' (www.r-project.org) and kept the resulting first three principal component (PC) axes to calculate the environmental dissimilarity between samples.Seascape variables were then represented by the three major PCA axes in the dbRDA analysis (Results).
To address the effect of space on genetic variation, we estimated Moran's Eigenvector Maps (MEMs) using the dbmem function in the 'adespatial' R package (Dray et al. 2022).MEMs represent a spectral decomposition of the spatial relationships among the geographical positions of all samples.This decomposition generates (n − 1) eigenfunctions, which are orthogonal variables that can be used in statistical models as explanatory variables representing spatial structure at different scales across the studied area.To generate these eigenvectors, we used a distance-based MEM approach, also known as principal coordinates of neighbour matrices (PCNM).This method, however, often generates few spatial predictors (and thereby displays low statistical power) and is sensitive to irregular sampling (Bauman et al. 2018a, b).The latter is not expected to be an important issue in our case since sampling was relatively continuous over the study area.Nevertheless, to overcome the limitations of PCNM, we also tested alternative graph-based MEM approaches (i.e.Delaunay triangulation and Gabriel graph) that allow for the accounting of irregular sampling (Bauman et al. 2018a, b).The results were congruent with the MEM approach and are therefore not reported here.Euclidean distances were considered here (as opposed to in-water distances that we use below) because these are more standard in the context of MEM analysis.Only MEMs with positive eigenvalues were retained in the analysis to account for positive spatial autocorrelations generated by broad-scale processes (e.g.genetic structure at large spatial scales, phylogeographic breaks) and smaller-scale processes such as gene flow among subpopulations and genetic drift (Dray et al. 2006).A total of 12, 11, 25 and 18 positive MEMs were obtained for D. sargus, M. surmuletus, P. elephas and S. cabrilla, respectively.

Defining reserve-induced protection status and reserve uniqueness
We calculated the in-water geographic distance (i.e. the shortest distance avoiding land) between individuals and the nearest reserve boundaries based on the GPS sampling points using the 'marmap' package (www.r-project.org).To reflect the fact that reserves are expected to have protection benefits beyond their boundaries but that this effect is often spatially restricted to a few kilometers (Manel et al. 2019), including in at least two of the species considered here (Benestan et al. 2021), we defined an individual as being protected when living inside the boundaries of the reserves and its surrounding area within a 5 km buffer.This categorical variable (protected versus not protected) was then transformed into a dummy variable with a value of 0/1.The number of samples considered protected and not protected for each species and reserve are presented in the Supporting information.We also evaluated the genetic singularity, i.e. particular genetic variation that is specific to just one reserve, by coding reserve identities as dummy variables, with the values of 0/1 indicating the absence/presence of one individual within one of the eight reserves (and 5 km buffer).Seven dummy variables were created to represent the network of eight reserves (n − 1).A detailed diagram of our statistical design is provided in Fig. 2.

Building a genomic database
We used two different approaches to generate single nucleotide polymorphism (SNP) markers: restriction-site-associated Figure 2. Statistical framework developed to disentangle the influence of seascape, space and protection status on genomic variation in four species (Diplodus sargus, Mullus surmuletus, Palinurus elephas and Serranus cabrilla).The genomic dataset of each species was transformed using principal coordinates analysis (PCoA).All axes were retained for the PCoA.Geographic coordinates were used to estimate Moran's Eigenvector Maps (MEMs) with the dbmem function in the 'adespatial' package (Borcard and Legendre 2002).Only MEMs with positive eigenvalues were retained.Seascape/environmental variables were transformed using a principal component analysis (PCA) and three principal components (PCs) were kept in the dbRDA.Distance-based redundancy analysis (dbRDA) was used to associate genomic variation to three types of explanatory variables (marine reserve, seascape and space).
Page 5 of 13 DNA (RAD) sequencing for M. surmuletus and D. sargus and diversity array technology (DArT) sequencing for S. cabrilla and P. elephas.RAD libraries were prepared in-house following a protocol described in Fietz et al. (2020).Full details on SNP calling of the white seabream (Diplodus sargus; 18 576 SNPs), the European spiny lobster (Palinurus elephas, 25 230 SNPs), the striped red mullet (Mullus surmuletus; 14 312 SNPs) and the comber (Serranus cabrilla; 13 101 SNPs) are available in Benestan et al. (2021).Briefly, we selected SNPs with less than 30% missing data, a coverage between 10× and 100× and a minor allele frequency (MAF) of 0.05 to avoid biases in the genome scan (Roesti et al. 2012).Individuals with too many missing loci (> 30%) were also filtered out.For each species, SNPs were separated into putatively neutral (hereafter neutral) and outliers (Supporting information) detected by the 'pcadapt' R package (Luu et al. 2017).'pcadapt' detects SNPs that contribute disproportionately to discriminant axes of a PCA, an approach that has been successfully used in a variety of contexts (Cayuela et al. 2020, Reynes et al. 2021).We first run 'pcadapt' using a large number of PCs (K = 20) to finally retain only two PCs that best explain population structure among individuals.All SNPs were then regressed against the retained PCs and p-values associated to z-scores were assigned to all SNPs.The 1% of SNPs with lowest p-values were considered outliers.This approach does not require grouping individuals into populations, which fits well with our individual-based sampling design.

Estimating genetic diversity and pairwise relatedness between individuals
Here, we consider the observed heterozygosity within each individual as a measure of genetic diversity.These were estimated using 'vcftools' (Danecek et al. 2011) and the 'hierfstat' (Goudet and Jombart 2015) R package through the het and basics.statsfunctions, respectively.As a reduction in genetic diversity may be due to inbreeding, we calculated the relatedness among individuals using the Loiselle kinship coefficient as implemented in GenoDive (Meirmans 2020).We explored the influence of the categorical variable protected/ unprotected on genetic diversity and kinship, respectively.Since heterozygosity and kinship did not follow normality (Shapiro test, p-value < 0.001), we carried out non-parametric Wilcoxon tests to test for difference in mean genetic diversity in samples considered protected versus not protected for each species.Similarly, we tested for difference in kinship among individuals considered protected versus among individuals considered not protected for each species.

Seascape, space and reserve influence on genomic patterns
We used distance-based redundancy analysis (dbRDA) to build linear models between the pairwise Euclidean genetic distances among individuals transformed by a principal coordinates analysis (PCoA) and three types of explanatory variables: seascape, space and protection status (Fig. 2).This step allows the identification of associations between genomic variation and explanatory variables.All axes were retained for the PCoA.A PCA on linear prediction then allows visualization of the genomic environmental associations.Here, the dbRDA was conducted using the dbrda function available in the 'vegan' package on the binary genotype variation obtained from plink (Purcell et al. 2007).We first tested the global model with all positive MEM variables.Then, when the global model was significant, we applied the ordiR2step forward selection procedure of the 'vegan' package based on permutation tests considering the global p-value and adjusted R 2 before applying the selection procedure (Blanchet et al. 2008).This approach is recommended for the selection of MEMs for pattern description because it is robust to type I errors (Bauman et al. 2018b).The significance and importance of different variables in the models were defined using the anova function of the 'vegan' package with 10 000 permutations (Legendre et al. 2011).After selecting the models that were significant, we computed a partial dbRDA to disentangle the exclusive influence of each explanatory variable on the genomic regions potentially reflecting drift and migration (namely neutral SNPs) and selection (namely outliers), separately.

Seascape structure
The two first axes of the PCA carried out on the 24 environmental variables considering the 1287 individuals (all species included) explain 74.4% of the total variation in seascape.The first (namely PC1 env) and second (namely PC2 env) axes alone capture 58.3% and 16.1% of total inertia.Overall, the PCA highlights the environmental divergence between three geographical groups: Gulf of Lion (i.e.Cerbère-Banyuls, Cap de Creus), Balearic Sea (i.e.Llevant de Mallorca, Menorca and Illes Columbretes) and Alboran Sea (i.e.Cabo de Palos, Cabo de Gata Níjar and Illa de Tabarca; Fig. 3a).Illes Columbretes is spread between the Gulf of Lion and Balearic Sea groups (Fig. 3a).PC1 is related to variations in temperature and salinity (Fig. 3b), and delineates the seascape differences between the Gulf of Lion located in the north, and the Balearic Sea and Alboran Sea, located further south (Fig. 1).Chlorophyll variables contributed most to the variation observed in PC2 (Fig. 3c), which reflects the differences between the Balearic Sea and the other groups.

Genetic diversity in protected versus unprotected areas
Genetic diversity estimates in protected versus not protected areas averaged over the whole reserve network are similar for the white seabream, the striped red mullet and the European spiny lobster both for neutral and outlier loci (Wilcoxon test; p-value > 0.05; Supporting information).But it is Page 6 of 13 significantly different for S. cabrilla (Wilcoxon test; p-value < 0.001; Fig. 4; Supporting information) both for neutral and oulier loci, respectively, even after removing individuals from the Alboran Sea to correct for the weak population structure observed in this species (Benestan et al. 2021).Genetic diversity at either neutral or outlier SNPs leads to the same conclusion for all species (Supporting information).As the comparisons protected/not protected for the four species remain consistent, irrespective of the use of outlier or neutral SNPs, only the results for the neutral SNPs are presented (Fig. 4).Kinship is not significantly different among S. cabrilla individuals considered protected versus among S. cabrilla individuals considered not protected (Wilcoxon test, p-value = 0.38), indicating that lower genetic diversity inside reserves is not related to a higher degree of kinship in protected areas.

Contrasting associations with seascape, space and genetic singularity
The distance-based redundancy analysis (dbRDA) using neutral SNPs reveals genetic singularities in four reserves, namely Cabo de Palos, Cerbère-Banyuls, Cap de Creus and Illes Columbretes, which stand out as significant in the final model for one or two species (Table 1, Fig. 5).We note that, although these models are globally significant, they have a low adjusted R 2 (< 1%).The variables representing reserve-induced protection, i.e. protected/not protected category and distance to the nearest reserve (see Fig. 2 for the statistical design), are not retained in any model when all species are combined.For D. sargus, the best model includes PC1 env (temperature/ salinity) and the Cabo de Palos dummy variable (Table 1, Fig. 5a).The model displays a clear separation of individuals according to their closest reserve, and we again observe the distinction between the Gulf of Lion and Balearic reserves along the second axis of the dbRDA.We detect a greater overlap among reserves for M. surmuletus than for D. sargus, although with a minor separation of individuals belonging to the Alboran Sea (Fig. 5b).Samples from Cerbère-Banyuls and Cap de Creus tend to form distinct groups that are significant in the model (Fig. 5b).In M. surmuletus, a total of five spatial vectors are retained to build the best model; MEM1 and PC1 env are the strongest drivers (Fig. 5b), revealing the south-north dichotomy previously shown by the PCA characterizing seascape composition (Fig. 3).The best model in P. elephas includes Cerbère-Banyuls and MEM9, with a very narrow scatterplot suggesting weaker structure in this highly dispersive species (Fig. 5c).For S. cabrilla, we also recover PC2 env (chlorophyll influence) as a significant variable in the model.This axis distinguishes the Alboran Sea from the other samples (Fig. 5d).Illes Columbretes seems to harbor a unique genetic signature, but samples are unevenly distributed along the first axis of the dbRDA.The dbRDA with outlier SNPs shows no clearly delineated clusters and mostly spatial vectors are retained in the best models for all species (Table 1).All the db-axes are significant (p-values < 0.05) for all species in all the presented models (Fig. 5, Table 1).

Local genetic diversity is not higher within reserves
Little consideration is still paid to genetic diversity in reserve networks, despite a growing literature demonstrating the effect of reserve-induced protection on species size, abundance and richness (Rolim et al. 2019).Our study aims to reduce this knowledge gap by documenting local genetic diversity in protected and unprotected areas for four species, D. sargus, M. surmuletus, P. elephas and S. cabrilla.The results go against our expectation, which was that reserves may contribute to preserve local genetic diversity, even in a limited way, as for example in the European Alps (Schoville et al. 2018) or in a French Atlantic dune system (Frey et al. 2016) where genetic diversity was found to be lower in humanimpacted areas.Higher allelic richness in protected areas has also been reported in Mediterranean populations of the white seabream in two marine reserves (Pérez-Ruzafa et al. 2006).This pioneering study was conducted 15 years ago with a few microsatellite markers.Today, SNPs outperform microsatellites in estimating individual-level multilocus heterozygosity and relatedness (Lemopoulos et al. 2019), but microsatellites retain the advantage of having higher mutation rates.Our results are more in line with the study of Pujolar et al. (2013) and Sahyoun et al. (2016), who evidenced homogenous genetic diversity in the white seabream across the marine reserve of Torre Guaceto (established in 1991) and neighboring non-protected areas up to 100 km from the reserve.In our case it is important to stress that the results apply to the whole network of eight reserves, which is the scale at which the analyses were conducted, as opposed to a particular reserve.Our study also echoes the results of a study of four key Baltic Sea fish species that documented identical genetic diversity in protected and not protected areas, concluding that protected areas do not harbor higher genetic diversity (Wennerström et al. 2017).Indeed, reserves have historically been created to support the management of harvested and threatened species with the underlying goal of providing fish biomass via spillover to fisheries (Di Lorenzo et al. 2016).This management strategy may have contributed to the exclusion of some fish populations that were not harvested or threatened, limiting the role of these reserves in effectively protecting local genetic diversity of these populations.This lack of change in genetic diversity in protected versus unprotected areas for the white seabream, striped red mullet and the European spiny lobster may be, at least partially, explained by the age and size of the reserves studied.Reserve age and size play an important role in reserve efficiency (Claudet et al. 2008, Edgar et al. 2014, Manel et al. 2019) and may have contributed to limit the effect of protection on genetic diversity.Here, reserves are on average less than 25 years old, allowing for only 5-8 generations to retain and accumulate genetic variation (Supporting information).Remarkably, only 3.7% (46) of protected areas are more than 20 years old in the Mediterranean Sea.The maintenance of these reserves is crucial for conservation.Indeed, even if protection tends to maintain a high number of breeding individuals within reserves, the signature of a greater genetic diversity in the following generations is expected to build up slowly.Studies have demonstrated the stability of local genetic diversity despite population declines, for example in Atlantic salmon Salmo salar over a time frame of three to five generations (Tessier and Bernatchez 1999) and in thornback ray Raja clavata over a 40-year timeframe (Chevolot et al. 2008).Another good example is the longspine sea urchin Diadema antillarum, whose populations collapsed across most of its range due to a disease outbreak and still shows no signs of reduction in genetic diversity (Lessios 2016).Given the high levels of genetic connectivity that characterize the study species and region (Benestan et al. 2021), we suggest that more than 5-8 generations would be required to capture a reserve effect on genetic diversity.Furthermore, these reserves are mostly small (< 100 km 2 ), which limits their ability to retain and build up genetic diversity through self-recruitment.This size could also be too small to capture a higher genetic diversity and species variability than is found in unprotected areas.In other words, the mere presence of a population within a

A conservation paradox
Our results show that reserves may influence genetic diversity, but in the opposite direction to what was expected.In S. cabrilla, a species of no commercial interest, greater genetic diversity is found in unprotected areas across the whole reserve network, which poses a conservation conundrum.This pattern is observed both at neutral and outlier SNPs, indicating that it is a genome-wide phenomenon.Moreover, the fact that kinship is not higher in protected areas indicates that this effect is not due to inbreeding within reserves.The difference in genetic diversity within versus outside reserves is subtle and requires therefore high power in terms of number of markers, individuals and reserves analyzed to be detected.Nevertheless, the fact that it is highly significant suggests that it reflects a biological reality.Furthermore, sample size is consistently higher for the protected category in the four species, which if anything would create a bias in favor of higher genetic diversity in protected areas.Nevertheless, we did not observe higher diversity in protected areas in any of the four species, and on the contrary report lower genetic diversity in protected areas for S. cabrilla.Thus, if anything, our sampling bias would tend to reinforce our findings.On the other hand, the spatial extent of the sampling is broader for the non-protected category due to the small size of the reserves.However, this is not expected to affect our results for D. sargus, M. surmuletus and P. elephas, since there is no genetic structure at the scale of the whole study area in these three species (Benestan et al. 2021).In S. cabrilla, the result of lower genetic diversity in protected areas is robust to the removal of individuals from the southern populations (Alboran Sea), indicating that it is not due to the weak genetic structure between southern and northern population that is observed in this species.
A similar conservation paradox has been reported before for S. cabrilla (Louisy et al. 2012), and more broadly on community assemblages (Micheli et al. 2004, Lester et al. 2009, Takashina et al. 2012, Boulanger et al. 2021), albeit not at the genetic level.Louisy et al. (2012) observed smaller S. cabrilla individuals inside the Cerbère-Banyuls reserve, whereas body size increased strongly outside the reserve.In agreement with our results, the authors found a significant reserve effect for S. cabrilla but not for D. sargus and M. surmuletus.Serranus cabrilla is a sedentary and territorial fish (Alós et al. 2011) with high site-fidelity, while D. sargus and M. surmuletus are mobile species (Claudet et al. 2006).Serranus cabrilla may therefore be more impacted by predation within reserves, due also to its smaller size and association with the benthos, even though the three studied species are at the same trophic level (Albouy et al. 2015).Reserves are known to benefit large predatory fishes in particular, since they are highly targeted in non-protected areas, which results in higher abundance of large predators within reserves (Macpherson andRaventos 2006, Boulanger et al. 2021).We therefore hypothesize that the decrease in genetic diversity in protected areas observed in S. cabrilla is potentially due to increased predation in the reserves.Our study only constitutes a first reference point to test this hypothesis, for S. cabrilla specifically and for other ecologically similar fishes more generally.A follow-up study comparing predation rates within and outside reserves is needed to test our hypothesis.

Delimiting the influence of seascape and reserveinduced protection on genomic variation
Marine organisms exhibit a myriad of reproductive strategies, posing a real challenge for the preservation of species with contrasting life histories.We have demonstrated that heterogeneous seascapes can shape fine-scale genome-wide patterns of variation in contrasting ways across species that share a common environment.Species-specific evolutionary patterns were also evidenced by a recent seascape genomic study on three sympatric southern African marine invertebrates with distinct life histories (Nielsen et al. 2020).Here, seascape has a significant but different effect on each of these species, as expected from the seascape genomics literature (Selkoe et al. 2016, Liggins et al. 2019).For the white seabream and the striped red mullet, temperature and salinity are associated with genomic variation (Fig. 5a, b), echoing the results of a recent population genomics study on the striped red mullet in the Mediterranean Sea (Dalongeville et al. 2018).For the comber, chlorophyll, an environmental factor characterizing primary productivity, was detected as the main seascape variable (Fig. 5d).Chlorophyll is known to influence the larval stage in S. cabrilla (Álvarez et al. 2012), which may contribute to causing associations between chlorophyll and genomic variation as observed in Homarus americanus (Dorant et al. 2022).In all three species, we note that a reserve signal also tended to emerge when genomic background and seascape variables were brought together, as clearly observed in white seabream and comber on the dbRDA analysis (Fig. 5a, d).This signal is based on allelic variation among individuals, transformed into PCoA, which can be rearranged faster that heterozygosity within individuals (e.g. the Hardy-Weinberg equilibrium can be attained in a single generation), but remains nonetheless weak (adjusted R 2 < 1%).This suggests that we are only seeing the onset of the effect of protection on genomic variation.

Limitations
Adaptive genetic diversity is related to short and longterm responses to selective pressures, whereas neutral diversity is more sensitive to past population bottlenecks (Holderegger et al. 2006).Highlighting these two evolutionary processes may provide distinct conservation solutions for optimal management of marine protected areas (Hanson et al. 2020, Xuereb et al. 2021).In our study, we did not find higher differences in genetic diversity or stronger population structure when considering outlier versus putatively neutral SNPs, which goes against to the growing literature on the subject (Gagnaire et al. 2015, Diopere et al. 2018, Sandoval-Castillo et al. 2018).Nevertheless, these results are Page 10 of 13 based on a specific outlier detection method (pcadapt) that is appropriate for our individual-based approach and does not require the a priori grouping of individuals with respect to sampling location or environmental variables.Our continuous sampling scheme may have limited our ability to detect loci under selection, as few genome scan methods are available for individual-based datasets, and different outlier detection approaches can lead to substantially different results (Lotterhos and Whitlock 2014, de Villemereuil and Gaggiotti 2015, Liggins et al. 2019); further developments in outlier detection methods are required to alleviate this limitation.The outlier SNPs that we have identified may also be mostly driven by other sources of selection than the 24 environmental variables considered here.However, it is still possible to apply genotype-environment association approaches, such as redundancy analysis (RDA) (Forester et al. 2018), to detect adaptive loci.Using this approach, we detected 193 (1.9%, D. sargus), 365 (2.5%, M. surmuletus), 190 (1.4%, P. elephas) and 722 (2.8%, S. cabrilla) loci, respectively, showing evidence of selection, a few of which were also identified by pcadapt (Supporting information).These results confirm that different outlier detection methods can identify different sets of markers which only partially overlap or even do not overlap at all, as for S. cabrilla (Supporting information).Additionally, we recommend broadening the sampling scale, as in Dalongeville et al. (2018) and Boulanger et al. (2020), to achieve optimal representativeness of seascape heterogeneity and increase the power to detect local adaptation as well as long-distance dispersal (Manel et al. 2012).

Connecting conservation genomics to marine reserve design
To be efficient, reserves of a network must be representative (e.g.boundaries that are consistent with population structure), exhaustive (e.g.consider the population structure of multiple species) and irreplaceable (e.g.preserve local genetic diversity and genetic singularity).We therefore emphasize that maximizing the protection of local genetic diversity may require expanding present reserves boundaries to increase self-recruitment with the reserve boundaries.This is particularly true for the four study species that are characterized by generally high levels of genetic connectivity in the study area (Benestan et al. 2021).Conservation benefits are maximized for large reserves that ensure a high fecundity of old and large individuals (De Leo and Micheli 2015).In terms of irreplaceability, efforts should be made to preserve the reserve network in the long term, especially the genetic uniqueness of the Cabo de Palos, Cap de Creus, Illes Columbretes and Cerbère-Banyuls reserves.Effective protection measures at the multi-generational level, whose effects can be captured with genetic methods, must be sustained over time.Finally, our results reveal that conservation policies to protect the genetic variability of populations inside a network of reserves should be based on a multi-species approach because, as we demonstrate here, species can respond differently to the seascape and to the protection measures in place.

Figure 1 .
Figure 1.Map of the individual-based sampling design conducted for each species: Diplodus sargus, Palinurus elephas, Mullus surmuletus and Serranus cabrilla.Each individual is colored by its closest marine reserve.From south to north, eight marine reserves were sampled: Cabo de Gata Níjar, Cabo de Palos, Illa de Tabarca, Llevant de Mallorca, Illes Columbretes, Menorca, Cap de Creus and Cerbère-Banyuls.Individuals were sampled inside (circle) and outside (triangle) the marine reserves, considering a 5 km buffer zone around the reserve's boundaries.The number of samples (n) varies from 1 to 31 individuals per geographic coordinate.

Figure 3 .
Figure 3. Seascape structure, space and protection status of the four species using principal component analysis (PCA).(a) Scatterplots of the first and second axis of the PCA on 24 environmental variables.Each individual is colored according to reserve designation and its protection status.(b) Coefficient loadings of the top three variables contributing to the first (left panel) and second axes (right panel).The variables considered here were: sea surface temperature (SST), sea surface salinity (SSS), chlorophyll benthic (ChlB) and chlorophyll surface (ChlS).The variables that show the highest coefficient loadings (in absolute values) are colored in black (b, c).

Table 1 .
Distance-based redundancy analysis results for each genomic dataset (X: genomic variation), species (species), seascape environment (Y: seascape data) and reserve (Y: closest reserve to which the sample belongs).Seascape environment includes spatial variation captured by the Moran's Eigenvector Maps (MEMs) and environmental parameters such as temperature, salinity and chlorophyll.Significance is indicated with *p-value < 0.05, **< 0.01, *** < 0.001.Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/ecog.06127 by Instituto Espanol de Oceanografia, Wiley Online Library on [05/06/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