# Climate differently influences the genomic patterns of two sympatric marine fish species

Emilie Boulanger, Laura Benestan, Pierre-Edouard Guérin, Alicia

Dalongeville, David Mouillot, Stéphanie Manel

► To cite this version:

Emilie Boulanger, Laura Benestan, Pierre-Edouard Guérin, Alicia Dalongeville, David Mouillot, et al.. Climate differently influences the genomic patterns of two sympatric marine fish species. Journal of Animal Ecology, 2022, Understanding climate change response in the age of genomics, 91 (6), pp.1180-1195. 10.1111/1365-2656.13623. hal-03589729

HAL Id: hal-03589729

https://hal.umontpellier.fr/hal-03589729v1

Submitted on 17 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Climate differently influences the genomic patterns of two sympatric marine fish species

Boulanger Emilie <fnref n="1" />,<fnref n="2" />,*, Benestan Laura <fnref n="1" />, Guerin Pierre-edouard <fnref n="1" />, Dalongeville Alicia <fnref n="2" />, Mouillot David <fnref n="2" />,<fnref n="3" />, Manel Stéphanie <fnref n="1" />,*

## Abstract :

1-Climate influences population genetic variation in marine species. Capturing these impacts remains challenging for marine fishes which disperse over large geographic scales spanning steep environmental gradients. It requires the extensive spatial sampling of individuals or populations, representative of seascape heterogeneity, combined with a set of highly informative molecular markers capable of revealing climatic-associated genetic variations.

2- We explored how space, dispersal and environment shape the genomic patterns of two sympatric fish species in the Mediterranean Sea, which ranks among the oceanic basins most affected by climate change and human pressure. We hypothesized that the population structure and climate-associated genomic signatures of selection would be stronger in the less mobile species, as restricted gene flow tends to facilitate the fixation of locally adapted alleles.

3- In order to test our hypothesis, we genotyped two species with contrasting dispersal abilities: the white seabream (Diplodus sargus) and the striped red mullet (Mullus surmuletus). We collected 823 individuals and used genotyping by sequencing (GBS) to detect 8,206 Single Nucleotides Polymorphisms (SNPs) for the seabream and 2,794 for the mullet. For each species, we identified highly differentiated genomic regions (i.e. outliers) and disentangled the relative contribution of space, dispersal and environmental variables (climate, marine primary productivity) on the outliers' genetic structure to test the prevalence of gene flow and local adaptation.

4- We observed contrasting patterns of gene flow and adaptive genetic variation between the two species. The seabream showed a distinct Alboran sea population and panmixia across the Mediterranean Sea. The mullet revealed additional differentiation within the Mediterranean Sea that was significantly correlated to summer and winter temperatures, as well as marine primary productivity. Functional annotation of the climate-associated outlier SNPs then identified candidate genes involved in heat tolerance that could be examined to further predict species' responses to climate change.

5- Our results illustrate the key steps of a comparative seascape genomics study aiming to unravel the evolutionary processes at play in marine species, in order to better anticipate their response to climate

<footnote n="1">¹ CEFE Univ Montpellier CNRS EPHE-PSL University IRD Montpellier France</footnote>

<footnote n="2">² MARBEC Univ Montpellier CNRS Ifremer, IRD Montpellier France</footnote>

<footnote n="3">³ Institut Universitaire de France Paris France</footnote>

<footnote>* Corresponding author : Emilie Boulanger, email address : emilie.boulanger@hotmail.com ; Stéphanie Manel, email address : stephanie.manel@ephe.psl.eu</footnote>
change. Defining population adaptation capacities and environmental niches can then serve to incorporate evolutionary processes into species conservation planning.

**Keywords**: climate change, comparative seascape genomics, gene flow, local adaptation, Mediterranean Sea
INTRODUCTION

The spatio-temporal variation of environmental variables (e.g. temperature) influences microevolutionary processes such as gene flow, drift and selection, thus contributing to shaping species' genomic patterns. Steep climatic gradients can act as environmental barriers reducing gene flow between nearby localities (e.g. Stanley et al. 2018), while local environmental conditions can act as strong selective pressures inducing local adaptation (Manel and Holderegger 2013; Dayan 2018). Climate can also influence genetic diversity by promoting or impeding gene flow and the demographic events necessary to the colonization of new environmental niches and habitats (Mittelbach et al. 2007; Manel et al. 2020). Climate change is fundamentally pushing species to avoid extinction by either adapting to new conditions within their current range (i.e. local adaptation), or by moving (i.e. dispersal, gene flow) to suitable habitats (Moritz and Agudo 2013; Gienapp 2020). Understanding how space, dispersal and environmental variables shape genetic variation is therefore crucial to identify the climate-adaptive potential of populations, assess their vulnerability to climate change, and delineate relevant ecological and evolutionary conservation units (Pauls et al. 2013; Gagnaire et al. 2015; Benestan 2019; Capblancq et al. 2020).

Landscape and seascape genomic methods offer a relevant framework for decoding gene flow and local adaptation patterns along with their drivers (Balkenhol et al. 2017; Liggins et al. 2019). These genome-wide marker approaches make it possible to identify highly differentiated regions of the genome, i.e. outlier loci, which improve the detection of genetic structure in species with high gene flow (D'Aloia et al. 2020), and help identify signatures of selection (Lewontin and Krakauer 1973). Associating outliers and potentially adaptive loci to environmental heterogeneity can then provide information on the population's local adaptation (Schoville et al. 2012; Rellstab et al. 2015; Capblancq et al. 2020).

Because the environment is inherently spatially structured, the patterns of isolation-by-distance and isolation-by-colonization, resulting from genetic (e.g. gene flow vs. drift) and demographic processes respectively, can coincide with environmental gradients (Bierne et al. 2011; Orsini et al. 2013), and confound patterns of isolation-by-environment driven by differential selection (Wang and Bradburd 2014). Partitioning the variance explained by spatial and environmental predictors across a wide geographical and environmental range makes it possible to disentangle their relative contributions to genetic variation (Peres-Neto et al. 2006; Vandamme et al. 2014; Dayan 2018). The detected outlier loci associated to environmental conditions can then
be mapped onto the species' annotated genomes when available, and used to detect candidate
genes and interpret their function (Manel et al. 2016). Finally, when possible, introducing these
candidate genes into species distribution models that predict the influence of climate change
makes it possible to take species- and population-specific eco-evolutionary responses into account
(Fitzpatrick and Keller 2015; Razgour et al. 2019; Waldvogel et al. 2020; Aguirre-Liguori et al.
2021; Nielsen et al. 2021).

Marine organisms occupy an open ocean with few physical boundaries (Pascual et al. 2017; Benestan et al. 2021). The strong dispersal potential of mobile species, either in their larval or adult stages (or both) (Almany et al. 2017; Manel et al. 2019) often translates into high gene flow and low genetic structure (Selkoe et al. 2016; Hoey and Pinsky 2018; Vandamme et al. 2021). Despite this weak genetic structure, temperature and salinity emerge as important drivers of marine population structure (Benestan et al. 2016; Selkoe et al. 2016; Stanley et al. 2018; Cayuela et al. 2020). Comparative genomic approaches can then be used to decipher the influence of species' traits on their adaptive responses (Nielsen et al. 2020; Torrado et al. 2020; Nielsen et al. 2021) but such datasets are rarely available for entire biogeographic regions or oceanic basins (Gagnaire 2020; Leigh et al. 2021).

Here, we explored how environmental gradients across the Mediterranean Sea can
influence the genomic patterns of two fish species while accounting for spatial and larval dispersal
patterns. The Mediterranean Sea is a unique hotspot of biodiversity (Coll et al. 2010) but is
projected to experience higher than average impacts of climate change under increasing
temperatures, salinity, and heat waves (Marbà et al. 2015; Darmaraki et al. 2019). Linear trends
showing increases in temperature and salinity have been recorded in the Mediterranean Sea since
1900 (Vargas-Yáñez et al. 2010; Borghini et al. 2014) with consequences for fish communities
(Givan et al. 2018). Local adaptation of fish populations in the Mediterranean Sea seems to be
closely linked to the thermal (from north to south) and saline (from west to east) gradients, as
previously reported for the striped red mullet (Dalongeville, Benestan, et al. 2018) and two
wrasses (Torrado et al. 2020). Deciphering the effect of thermal and saline gradients on
Mediterranean fish populations is the first step towards understanding the eco-evolutionary
processes at play, and anticipating the effects of climate change. More broadly, monitoring how
species cope with climate change in this semi-enclosed sea could serve as an example of what
might happen globally (Lejeusne et al. 2010).
We focus on two socio-economically important species, the striped red mullet *Mullus surmuletus*, and the white seabream *Diplodus sargus*. Both species share similar distribution ranges across the Mediterranean Sea and Eastern Atlantic and are heavily targeted by small-scale fisheries (Claudet et al. 2010). Although they share similar larval duration periods (PLD, Macpherson & Raventos, 2006), these species have contrasting movement patterns at the adult stage. The seabream engages in seasonal adult migration, reported up to 90 km (Abecasis et al. 2009), while the mullet's movement as an adult is assumed to be more restricted (Macpherson and Raventos 2006). Here we used individual-based sampling on an exceptionally large spatial scale and a genotyping by sequencing (GBS) approach providing a high number of genome-wide markers. By applying an original comparative seascape genomics approach, our objectives were to (i) describe the spatial distribution of genetic diversity, (ii) define outlier SNPs to improve the detection of population structure, (iii) disentangle the drivers influencing population genetic structure and (iv) identify candidate genes potentially involved in local adaptation to environmental conditions affected by climate change. We hypothesize that population structure and environmental adaptive signatures would be lower in the more mobile species, i.e. the seabream, as gene flow tends to counteract natural selection.

# MATERIAL & METHODS

## Study species and sampling effort

The white seabream and striped red mullet are coastal fishes that spawn from March till June/July (Reñones et al. 1995; Morato et al. 2003). Both produce pelagic larvae that remain in the water column for approximately 20 to 35 days (Macpherson and Raventos 2006). The seabream inhabits rocky reefs and seagrass beds, whereas the mullet mostly occupies sandy and rocky habitats (Froese and Pauly 2019).

We analysed 526 mullet and 297 seabream individuals collected in 2014 from 64 and 59 sites, respectively (Dalongeville, Andrello, et al. 2018) (Figure 1). The samples covered a large area, roughly 15 degrees latitude, 40 degrees longitude and over 3,800 km across the entire Mediterranean Sea covering all seven marine ecoregions (Figure 1). Fin clips were taken from small-scale fishery catches which required no permit or ethical approval, and conserved in 96% ethanol prior to being stored at 4°C.
## Genetic data and SNP calling

The SNP datasets were produced using a Genotyping by Sequencing (GBS) approach, in order to genotype a large number of individuals for a large number of SNPs covering the whole genome. Genomic DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen). The equipment was cleaned with 100% alcohol between each sample to avoid contamination. After quality and concentration assessment, GBS libraries for the mullet and seabream were constructed using restriction enzyme ApeKI (Elshire et al. 2011), and Pst1/Bfa1, respectively. Libraries of 95 individuals per lane were prepared and sequenced at the Institute of Genomic Diversity at Cornell University using the Illumina HiSeq 2500 (100 bp, single-end reads). More details can be found in Dalongeville et al. (2018).

Raw GBS sequences were filtered and trimmed based on quality (Supplementary Methods). Demultiplexing was performed using the process_rad-tags program in STACKS v.1.48 (Catchen et al. 2011). Sequences were mapped on the published reference genomes (Fietz et al. 2020: genome size = 785 and 613 Mbp, contig N50 = 1101 and 384 kbp for seabream and mullet respectively) with BWA2 (Vasimuddin et al. 2019). Variant calling was performed using the FreeBayes software (Garrison and Marth 2012) as implemented in the dDocent pipeline (Puritz et al. 2014). Additional filters were further applied with VCFtools (Danecek et al. 2011) and vcflib (https://github.com/vcflib/vcflib#vcflib) following the method set out in O'Leary et al. (2018). In short, individuals with more than 50% missing data were removed, as well as Single Nucleotide Polymorphisms (SNPs) with a genotype call rate below 95%. The SNPs were further filtered based on allele balance, mapping quality ratio and quality/depth ratio (Table S1).

## Detecting outlier SNPs

We applied PCAapt to detect outlier loci that were more differentiated than under a neutral model (Luu et al. 2017). PCAapt is currently the only individual-based method available with no prior hypotheses on which environmental factors drive differentiation (Liggins et al. 2019). PCAapt conducts a Principal Component Analysis (PCA) on the individual genotypes and retains only the principal components that best explain genetic structure across individuals. All SNPs are then regressed against the retained ordination axes and outlier SNPs are selected based on their significant correlation with these axes. P-values for each SNP were obtained using the R
function *pcadapt()*, and false discovery rates (FDR) were estimated using the function *qvalue()*
(Storey et al. 2019). A conservative FDR cut-off ($\alpha$ = 0.005) was applied for outlier detection.

Genetic diversity was assessed based on the full dataset. Genetic structure was investigated separately for the outlier and non-outlier SNPs. For the non-outlier dataset, only SNPs in Hardy-Weinberg Equilibrium (HWE) were kept (Beaumont and Nichols 1996).

Genetic diversity

We estimated global genetic diversity as the expected heterozygosity ($H_e$) for all individuals with the function *Hs()* (R package *adegenet*, Jombart 2008). Because our sampling design is not stratified in populations, we further investigated the spatial variation of $H_e$ using a grid-based approach modified from Larranaga et al. (2017) and Thomas et al. (2012). Individuals were replicated to nearby cells at a resolution of approximately 80 km, and $H_e$ was calculated for each cell (Supplementary Methods). The resulting grid contained on average 56 ± 33 individuals per cell for the mullet, and 37 ± 25 individuals for the seabream. $H_e$ was calculated on bootstrapped subsamples of 20 individuals per cell to avoid sampling bias. We further calculated the observed individual heterozygosity ($H_o$) as the proportion of heterozygous loci per individual using VCFtools, and compared values between ecoregions and species.

Genetic structure

We inferred population structure using the Discriminant Analysis of Principal Components
(DAPC) with *adegenet*. DAPC maximizes genetic differences between groups while minimizing
variation within groups (Jombart et al. 2010). We did not use the sampling sites as prior
information, but selected the optimal number of clusters using the Bayesian Information Criterion
(BIC) with the function *find.clusters()*. We calculated the pairwise $F_{ST}$ between all DAPC clusters
following Weir and Cockerham (1984) using the function *stamppFst* from the *StAMPP* package
(Pembleton et al. 2013). We also calculated the effective population size ($N_e$) for each of the
DAPC clusters inferred with non-outlier SNPs using $N_e$Estimator v2.1 under a random mating
model (Do et al. 2014). $N_e$ was calculated on non-outlier SNPs to avoid potential biases from loci
under selection (Waples 2006; Candy et al. 2015) and a minimum allele frequency of 0.02
(Waples and Do 2010).
# Characterizing the seascape features influencing population structure

## Spatial variables

Marine geographic distances were calculated as the minimum distances constrained by water (Supplementary Methods). Pairwise geographic distances were then converted to distance-based Moran's eigenvector maps (dbMEM), which are uncorrelated (orthogonal) vectors that can be used as site-based variables in regression analyses (Dray et al. 2006; Legendre and Legendre 2012). The dbMEM were computed from a principal coordinate analysis (PCoA) on the truncated distance matrix. The new eigenvectors with the highest eigenvalues describe broad-scale spatial structures and were retained for subsequent analyses (Supplementary Methods, Figure S1a-b). The dbMEMs were calculated using the function *dbmem()* from the *adespatial* package with the default truncation threshold (Dray et al. 2020).

## Larval dispersal ability

Larval dispersal was quantified using a biophysical model which estimates the probability of larval connection between each pair of sites (Andrello et al. 2015). We used a pelagic larval duration of 30 days and released larvae every 3 days from 1 May until 28 May, which corresponds to the spawning season of the study species. The same larval dispersal estimates were used for both species considering their respective sampled sites (Reñones et al. 1995; Morato et al. 2003). Pairwise dispersal probabilities were then converted to Asymmetric Eigenvector Maps (AEM). The AEM are orthogonal vectors similar to dbMEM but account for the directionality of larval dispersal (Blanchet et al. 2008). Probabilities were translated into a site-by-edge binary matrix where two sites are connected if the dispersal probability is > 0 (Supplementary Methods). Orthogonal AEMs were then computed from the binary matrix using the function *aem()* from the *adespatial* package. Similarly to dbMEM, the new AEM with the highest eigenvalues describe broad-scale dispersal structures and were retained for subsequent analyses (Figure S1c).

## Environmental variables

Monthly mean sea surface salinity, temperature and chlorophyll a (included as a proxy of marine primary productivity) values around each site were retrieved from the Bio-ORACLE database for the period 1987-2015 (Tyberghein et al. 2012). Salinity and chlorophyll a values were averaged over the months and years, whereas temperature values were averaged over the years by season (i.e. spring, summer, autumn and winter) using the *raster* package (Hijmans 2020). To
avoid problems with collinearity, we only kept the winter and summer mean temperature to represent extreme values (i.e. the coldest and warmest).

## Disentangling drivers of genetic variation

We used a distance-based redundancy analysis (dbRDA) on Euclidean genetic distances estimated between individuals to characterize the relative contribution of spatial (dbMEM), dispersal (AEM vectors) and environmental (salinity, chlorophyll a, winter and summer temperature) variables (Legendre and Legendre 2012) to genetic variation. Starting with all the selected variables, we first undertook variable selection using the *vegan's* `ordiR2step()` function (Oksanen et al. 2016) to retain the most informative explanatory variables and minimize variance-inflation. In the second step, variables with the largest variance inflation factor (VIF) were sequentially removed until all VIF values were < 3 and all variables were sufficiently uncorrelated (Figure S2). The selected set of explanatory variables was then tested in a final dbRDA model. We also conducted multiple partial dbRDA to test the effect of each environmental variable after accounting for the effect of spatial and dispersal variables. We then performed variance partitioning to quantify the separate and shared proportion of the genetic variation explained by the spatial (SPACE), dispersal (DISP), and environmental (ENV) variables using the *vegan's* `varpart()` function. Model significance, as well as the significance of each dbRDA axis were tested using ANOVA-like permutation tests with 9999 permutations as implemented in the *vegan's* `anova.cca()` function (Legendre et al. 2011).

## Functional annotation & climatic variance

We retrieved the genomic position of the outlier SNPs from the annotated seabream and mullet genomes (Fietz et al. 2020). For SNPs positioned in annotated coding regions, we extracted the flanking DNA sequences centered on the variant position (200 bp length) and aligned them against the Swiss-Prot protein sequence database using BLAST+ (Camacho et al. 2009). We only retained sequences ranking first or second, with a minimum identity of 70% and an E-value below $10^{-6}$. We then assessed whether SNP variations were synonymous or not using ANVAGE (https://github.com/Grelot/anvage). Non-synonymous variants were kept to increase the possibility of selected SNPs having functional implications (Rellstab et al. 2015). Finally, we tested whether the individual genotypes of these candidate SNPs correlated to the environmental variables identified by the dbRDA.
# RESULTS

## Successful SNP calling

We genotyped a total of 8,206 filtered and informative SNPs from 297 white seabreams (2.74 % missing data) and 2,794 SNPs from 467 striped red mullets (2.83 % missing data) (Table S1). In these datasets, we detected 413 (5.0%) and 291 (10.5%) outlier SNPs for the seabream and the mullet, respectively. The remaining non-outlier SNPs, after filtering for HWE, contained 7,655 SNPs for the seabream and 2,462 SNPs for the mullet.

## Genetic diversity patterns

The expected heterozygosity ($H_e$) calculated on all individuals was 0.26 and 0.24 for the seabream and the mullet, respectively. Spatial interpolation of $H_e$ revealed that the genetic diversity of seabream was mostly stable across the Mediterranean Sea with peaks in the Western Mediterranean and Aegean Sea (Figure 2a). The genetic diversity of mullet was highest in the Alboran Sea and gradually decreased eastward (Figure 2b). Individual observed heterozygosity ($H_o$) confirmed these patterns (Figure 2c) and was significantly different between species, ecoregions, and their interaction (two-factor ANOVA, all p-values < 0.001, Table S2).

## Higher resolution for population structure with outlier SNPs

The DAPC on non-outlier SNPs highlighted almost complete panmixia for both species. Nevertheless, we found a slight difference between them as the lowest BIC value was $K = 2$ for the white seabream and $K = 1$ for the striped mullet. For the white seabream, $K = 2$ indicated that most individuals belonged to one cluster (Ne = 4556, Table S3) (Figure S3a) while 27 individuals mostly from the south of the Alboran Sea formed the second genetic cluster (Ne = 255, Table S3) (pairwise $F_{ST} = 0.02$). For the mullet, we found one single cluster, except for a few scattered individuals when forcing $K = 2$ (pairwise $F_{ST} = 0.006$) (Ne = 9,489 and 607 respectively, Table S3) (Figure S3b).

The genetic structure observed from seabream outliers revealed three clusters, two of which were similar to those obtained with non-outlier SNPs (Figure 3a) but with stronger differentiation (pairwise $F_{ST} = 0.14 – 0.58$, Table S4). For the mullet, the use of outlier SNPs increased our resolution and we uncovered a stronger genetic structure (Figure 3b). The BIC
suggested the optimal numbers of genetic clusters as between K = 3 (BIC = 1,260) and K = 7 (BIC = 1,250) (Figure S5). All pairwise FST values were higher than with non-outlier SNPs (0.05 - 0.53, Table S4). DAPC for K=3 showed that the first and second clusters represented a more gradual genetic differentiation (in orange and yellow respectively; Figure 3b). These were comprised of individuals from the Western and Eastern Mediterranean underlining a progressive genetic shift between the two basins. The third cluster (in blue; Figure 3b) was more strongly differentiated and comprised individuals from the Alboran Sea, and also from the Gulf of Lion and the Northern Adriatic Sea. Increasing the number of clusters mostly subdivided the Western Mediterranean basin. When selecting K=7, individuals from the initial third cluster (when K=3) were separated into two clusters: one comprising of individuals from the Alboran Sea (in dark blue; Figure 3b), and the second one of individuals from the Gulf of Lion and the Northern Adriatic Sea (in light blue; Figure 3b).

## Climate-associated genetic variation

Based on their eigenvalues, the first seven and five dbMEM were selected for the seabream and mullet, respectively, while the first six AEM were selected for both species. The distance-based redundancy analysis on non-outlier SNPs produced significant (p-value = 0.0001) but low fits for both species ($R^2_{adj}$ = 0.006 for the seabream and 0.004 for the mullet) and is therefore not discussed in any further detail. The dbRDAs conducted on outlier SNPs were globally significant (p-value = 0.001) and explain 14.5% and 8% of genetic variation in the seabream and mullet respectively (Figure 4). Only the first dbRDA axis was significant in both cases, accounting for 15% and 7.9% of the overall (unadjusted) outlier based genetic variation.

For the seabream, salinity, summer and winter temperature, along with three geographic-distance vectors (dbMEM2, dbMEM3 and dbMEM6) and three larval-dispersal vectors (AEM1, AEM3 and AEM5) were retained (Figure S5). The environmental partition alone explained 19% of the dbRDA variance (Figure 5a), while salinity and summer temperature were significant predictors after accounting for the variance explained by spatial and dispersal variables in additional partial dbRDAs (p-value < 0.05; Table 1).

For the mullet, chlorophyll a, summer - and winter temperature along with one geographic-distance vector (dbMEM5) and five larval-dispersal vectors (AEM1, AEM3, AEM4, AEM5 and AEM6) were retained (Figure S6). The environment explained 18% of the dbRDA
variance (Figure 5b), while summer and winter temperatures and chlorophyll a were significant
predictors after accounting for the variance explained by spatial and dispersal variables in partial
dbRDAs (p-value < 0.05; Table 1).

Candidate genes

For the white seabream, nine of the 413 outlier SNPs (2%) were located on annotated
coding regions and returned known protein-encoding genes. Six were non-synonymous variants,
five of which were significantly related to summer temperature (Table S6, Figure S7). Among
these, we found genes encoding for the Zinc finger MYND domain-containing protein 10
(summer temperature ~ genotype p-value < 0.01, Figure 6a), and the Pumilio homolog 2 (p-value
< 0.001). For the mullet, 19 out of the 291 outlier SNPs (7%) were positioned on annotated coding
regions and returned known protein-encoding genes. Eight were non-synonymous variants, seven
of which were significantly related to summer temperature and four to winter temperature (Table
S6, Figure S8). Among these, we found genes encoding for tubulin-specific chaperone D (p-value
< 0.001, Figure 6b) and protein arginine methyltransferase (p-value < 0.001).

DISCUSSION

Decrypting climate-associated genetic variation is essential to accurately predict the impact
of climate change on living organisms (Bay et al. 2017; Bay et al. 2018; Razgour et al. 2019).
Here, we assessed the extent of gene flow and adaptive genetic variation in two fish species with
different movement strategies, the white seabream *Diplodus sargus* and the striped red mullet
*Mullus surmuletus*, in order to unpack the role of gene flow and climate on local adaptation. We
found that both species are influenced by the same biogeographic barrier to gene flow and that
their respective structures were driven by a combination of spatial, dispersal and environmental
variables. Functional annotation of outlier SNPs identified candidate genes associated with
temperature suggesting local adaption to temperature for both species. The structure based on
outlier SNPs is stronger for mullet in the Western Mediterranean basin suggesting greater local
adaptation for the less mobile species. Overall, our results illustrate important advances in
seascape genomics and spatial ecology which may well facilitate climate change vulnerability
assessment.
Neutral processes structuring genetic variation and diversity

The spatial dbMEM and larval dispersal AEM together explained 35% and 60% of the modeled genetic variation in seabream and mullet, respectively. This indicates that a combination of gene flow between connected populations and drift within disconnected populations, e.g. in the Alboran Sea, partly drives the genetic structure. As population sizes are mostly large, both in the fishes studied (Table S3) and in marine species in general, we expected genetic drift to have a small effect and gene flow a much stronger effect (Allendorf et al. 2010).

In our study, the Almeria-Oran Front (AOF) divides the Mediterranean Sea into two seabream populations: the south of the Alboran Sea and the rest of the Mediterranean Sea (Figure 3). The AOF is formed by the influx of colder and fresher Atlantic water into the Mediterranean Sea through the Strait of Gibraltar, which influences the oceanographic dynamics and creates a near-permanent oceanic front (Patarnello et al. 2007). The influx of fresher water simultaneously creates an abrupt change in temperature and salinity, which gradually increases eastward (Figure 1c). The AOF acts as a barrier to gene flow and biogeographic break for numerous species (Patarnello et al. 2007; Pascual et al. 2017) including the mullet (Galarza et al. 2009). However, it was not thought to differentiate seabream populations until now (Bargelloni et al. 2005; González-Wangüemert et al. 2010). This break was only detected for the mullet by using outlier SNPs. In contrast with the seabream there is a lack of samples from the south of the Alboran Sea (Figure 1), which probably weakens the inference of a well-defined Alboran mullet population.

The AOF forms a semi-permeable barrier which allows two genetic pools to mix in the contact zone: the Mediterranean pool and the Atlantic pool (Carreras et al. 2020). This is reflected in the spatial distribution of genetic diversity which peaks in the Western Mediterranean for the seabream (Figure 2a). In the absence of strong biogeographic barriers in the main Mediterranean basin, the seabream shows high levels of gene flow impeding further differentiation and homogenizing genetic diversity. This lack of spatial variation in genetic diversity across the Mediterranean Sea is supported by a published study of 31 Mediterranean fishes, including our two species, using microsatellite and mitochondrial markers (Dalongeville et al. 2016). For the mullet, on the other hand, genetic diversity peaks in the Alboran Sea and gradually decreases eastward (Figure 2). These patterns coincide with the genetic structure which is divided into several clusters in the Western Mediterranean, but consists of one main cluster in the Eastern Mediterranean (Figure 3). This dichotomy could be explained by the introgression of Atlantic
haplotypes gradually decreasing as the distance from the Atlantic-Mediterranean contact zone
increases (Duranton et al. 2019).

The influence of climate on genetic variation and potential local adaptation

The influx of fresher Atlantic water which generates the AOF, simultaneously creates an abrupt change in both temperature and salinity. In addition to neutral processes, local adaptation of the two seabream populations on both sides of the front could maintain or even promote their differentiation, as the environmental variables explain up to 19% of the modeled genetic variance in the seabream (Figure 5a). These results contribute to the mounting evidence from studies of numerous different species, that the Atlantic-Mediterranean break is enabled by differential selection and local adaptation to contrasting temperature and salinity conditions (Milano et al. 2014; Pascual et al. 2016; Carreras et al. 2020). The lack of further differentiation and local adaptation in seabream in the main Mediterranean basin concurs with our initial hypothesis, as the seabream is capable of exceptional long-distance dispersal (Abecasis et al. 2009; Aspillaga et al. 2016).

The mullet, on the other hand, is less mobile and shows finer genetic differentiation within
the main Mediterranean basin. The spatial dbMEM, larval AEM, temperature, salinity and
productivity patterns were insufficiently correlated to quantify their relative and shared importance
(Figure S1). Taken together, the environmental variables explained 18% of the modeled genetic
variation in the mullet not explained by spatial and dispersal structures (Figure 5). The
combination of elevated productivity and cold winter temperatures correlated with the
differentiation and clustering of geographically distant individuals in the Alboran Sea (DAPC
cluster 3), in the Gulf of Lion (cluster 7) and the Northern Adriatic Sea (cluster 7) (Figure 4b).
This differentiation breaks the gradual pattern expected under isolation-by-distance and their
genetic similarity is most probably due to a shared local adaptation. A very similar pattern of local
adaptation was previously identified for the ocellated wrasse, a reef-associated fish dispersing with
a larval stage of 7 to 13 days (Torrado et al. 2020). Individuals sampled from the Alboran Sea and
the Gulf of Lion were grouped together by redundancy analysis and identified as locally adapted
to temperature and productivity (Torrado et al. 2020). Productivity, through food availability, and
temperature affect the development, survival and settlement of larvae (Robitzch et al. 2016;
Raventos et al. 2021), acting as environmental selective pressures and prompting local adaptation in fishes (Cayuela et al. 2020; Torrado et al. 2020).

Identifying the genomic signatures of selection associated with temperature from outliers can lead to the discovery of candidate genes that may be involved in species' adaptive response to climate change. Furthermore, these genes confirm the adaptive nature of the associated outlier SNPs (Manel et al. 2010). The tubulin specific chaperone D gene, detected in mullet as having potentially locally-adapted alleles, is associated with the stress-response to varying temperatures and salinities in fishes (Whitehead et al. 2011; Avarre et al. 2014) and mussels (Tomanek et al. 2012). Homozygote genotypes A/A at the SNP located in the chaperone D gene are found in individuals living in environments with the lowest winter temperatures, whereas homozygote G/G genotypes are more prevalent at higher temperatures (Figure 6b). Zinc finger proteins, for which we found two associated SNPs in seabream, are involved in heat and salinity stress response mechanisms in plants (e.g. Sakamoto et al. 2004; Droll et al. 2013) and the thermal tolerance of fish (Yu et al. 2018). Homozygote G/G genotypes are most prevalent at higher summer temperatures, whereas T/T homozygotes mostly occur at lower temperatures (Figure 6a). Finally, combining genomic information with phenotypic and fitness data would allow for the functional interpretation of potential genotype-environment associations and further unravel adaptive responses to a changing seascape (Gagnaire and Gaggiotti 2016; Waldvogel et al. 2020).

## From theory to practice: conservation, climate and genomics

Marine Protected Areas (MPAs) are considered to be a major conservation tool for maintaining and restoring biodiversity and ecosystem functions (Halpern et al. 2010; Edgar et al. 2014). Networks of MPAs are established to promote connectivity between MPAs and ensure the persistence of populations within the network (Gaines et al. 2010). Networks of MPAs then help maintain gene flow and allow the spread of advantageous alleles (Xuereb et al. 2019). Genomic studies provide crucial insights for designing efficient MPA networks, as they provide information on gene flow and the (mal)adaptation of populations to their changing environment. Spatial conservation planning can make use of genomic metrics (e.g. adaptive genetic diversity) to incorporate species' evolutionary potential into the design of MPA networks, by prioritizing the protection of populations adapted to specific environmental conditions (Nielsen et al. 2017; Hanson et al. 2020; Xuereb et al. 2020). These locally-adapted populations can act as sources of
adaptive alleles for the meta-population (Selmoni et al. 2020), thus promoting species' adaptation. For example, genotypes adapted to high temperatures can spread from the Levantine Sea to the western Mediterranean basin via the movement of individuals, allowing western populations to adapt to the expected temperature increase. They can also be used for assisted migration, by transplanting locally adapted individuals in order to increase the species' adaptive potential under climate change (Bay et al. 2017). As climate change causes rapid environmental shifts, adaptive management, under which management strategies can be modified based on new knowledge acquired and assessment of the effectiveness of previous actions (Katsanevakis et al. 2017), has been proposed as the most efficient solution to preserve marine biodiversity in the long term (Frazão Santos et al. 2020; Rilov et al. 2020). Genomic data will thus be required to monitor the effectiveness of conservation strategies in maintaining species' adaptability over time, and can also be incorporated into predictive models to forecast species' adaptive responses to predicted environmental variations induced by climate change (Razgour et al. 2019).

## Acknowledgements

We are grateful to all the fishermen and colleagues involved in the sampling and molecular analyses. We thank Marco Andrello for contributing his codes for the spatial genetic diversity analyses. This research was funded by the 'Total foundation' through the SEACONNECT project, the 'ANR' through the RESERVEBENEFIT project and the SEAMONT project, PSL Environment through the project RESMOD, and EB was partially funded by the Region Occitanie. The bioinformatic analyses made use of the Montpellier Bioinformatics Biodiversity (MBB) platform, supported by the ANR program 'Investissements d'avenir' (ANR-10-LABX-04-01).

## Data Availability Statement

All the codes required to reproduce the bio-informatic processing and statistical analyses are available at https://gitlab.mbb.univ-montp2.fr/seaconnect (Boulanger et al. 2021a). The SNP data files have been uploaded to Dryad: https://doi.org/10.5061/dryad ffbg79cvs (Boulanger et al. 2021b).

## Author Contributions
SM, AD, LB, EB and DM conceived the study and designed methodology; AD collected the data; PEG conducted the bioinformatic analyses; EB conducted the statistical analyses; EB, LB and SM wrote the first draft. All authors made critical contributions to the manuscript and approved the final manuscript for publication.

## Conflict of interest

The authors declare that they have no conflicts of interest.

## References

Abecasis D, Bentes L, Erzini K. 2009. Home range, residency and movements of Diplodus sargus and Diplodus vulgaris in a coastal lagoon: Connectivity between nursery and adult habitats. Estuar Coast Shelf Sci. 85(4):525-529. doi:10.1016/j.ecss.2009.09.001.

Aguirre-Liguori JA, Ramírez-Barahona S, Gaut BS. 2021. The evolutionary genomics of species' responses to climate change. Nat Ecol Evol. 5(10):1350-1360. doi:10.1038/s41559-021-01526-9.

Allendorf FW, Hohenlohe PA, Luikart G. 2010. Genomics and the future of conservation genetics. Nat Rev Genet. 11(10):697-709. doi:10.1038/nrg2844.

Almany GR, Planes S, Thorrold SR, Berumen ML, Bode M, Saenz-Agudelo P, Bonin MC, Frisch AJ, Harrison HB, Messmer V, et al. 2017. Larval fish dispersal in a coral-reef seascape. Nat Ecol Evol. 1(6):1-7. doi:10.1038/s41559-017-0148.

Andrello M, Jacobi MN, Manel S, Thuiller W, Mouillot D. 2015. Extending networks of protected areas to optimize connectivity and population growth rate. Ecography (Cop). 38(3):273-282. doi:10.1111/ecog.00975.

Aspillaga E, Bartumeus F, Linares C, Starr RM, López-Sanz À, Díaz D, Zabala M, Hereu B. 2016. Ordinary and extraordinary movement behaviour of small resident fish within a mediterranean marine protected area. PLoS One. 11(7):1-19. doi:10.1371/journal.pone.0159813.

Avarre JC, Guinand B, Dugué R, Cosson J, Legendre M, Panfili J, Durand JD. 2014. Plasticity of gene expression according to salinity in the testis of broodstock and F1 black-chinned tilapia, Sarotherodon melanotheron heudelotii. PeerJ. 2014(12):1-20. doi:10.7717/peerj.702.

Balkenhol N, Dudaniec RY, Krutovsky K V, Johnson JS, Cairns DM, Segelbacher G, Selkoe KA, von der Heyden S, Wang IJ, Selmoni O, et al. 2017. Landscape Genomics:
Understanding Relationships Between Environmental Heterogeneity and Genomic Characteristics of Populations. In: Rajora OP, editor. Population Genomics: Concepts, Approaches and Applications. Cham: Springer International Publishing. p. 261–322.

Bargelloni L, Alarcon JA, Alvarez MC, Penzo E, Magoulas A, Palma J, Patarnello T. 2005. The Atlantic-Mediterranean transition: Discordant genetic patterns in two seabream species, Diplodus puntazzo (Cetti) and Diplodus sargus (L.). Mol Phylogenet Evol. 36(3):523–535. doi:10.1016/j.ympev.2005.04.017.

Bay RA, Harrigan RJ, Underwood V Le, Gibbs HL, Smith TB, Ruegg K. 2018. Genomic signals of selection predict climate-driven population declines in a migratory bird. Science (80- ). 359(6371):83–86. doi:10.1126/science.aan4380.

Bay RA, Rose NH, Logan CA, Palumbi SR. 2017. Genomic models predict successful coral adaptation if future ocean warming rates are reduced. Sci Adv. 3(11):1–10. doi:10.1126/sciadv.1701413.

Beaumont MA, Nichols RA. 1996. Evaluating loci for use in the genetic analysis of population structure. Proc R Soc London Ser B Biol Sci. 263(1377):1619–1626. doi:10.1098/rspb.1996.0237.

Benestan L. 2019. Population Genomics Applied to Fishery Management and Conservation. In: Oleksiak MF, Rajora OP, editors. Population Genomics: Marine Organisms. Springer, Cham. p. 399–421.

Benestan L, Quinn BK, Maaroufi H, Laporte M, Clark FK, Greenwood SJ, Rochette R, Bernatchez L. 2016. Seascape genomics provides evidence for thermal adaptation and current-mediated population structure in American lobster (Homarus americanus). Mol Ecol. 25(20):5073–5092. doi:10.1111/mec.13811.

Bierne N, Welch J, Loire E, Bonhomme F, David P. 2011. The coupling hypothesis: Why genome scans may fail to map local adaptation genes. Mol Ecol. 20(10):2044–2072. doi:10.1111/j.1365-294X.2011.05080.x.

Blanchet FG, Legendre P, Borcard D. 2008. Modelling directional spatial processes in ecological data. Ecol Modell. 215(4):325–336. doi:10.1016/j.ecolmodel.2008.04.001.

Borghini M, Bryden H, Schroeder K, Sparnocchia S, Vetrano A. 2014. The Mediterranean is becoming saltier. Ocean Sci. 10(4):693–700. doi:10.5194/os-10-693-2014.

Boulanger E, Benestan L, Guerin P-E, Dalongeville A, Mouillot D, Manel S. 2021a. Codes from: Climate differently influences the genomic patterns of two sympatric marine fish species,
https://gitlab.mbb.univ-montp2.fr/seaconnect.

Boulanger E, Benestan L, Guerin P-E, Dalongeville A, Mouillot D, Manel S. 2021b. Data from: Climate differently influences the genomic patterns of two sympatric marine fish species, Dryad Digital Repository, https://doi.org/10.5061/dryadFFFg79cvs.

Boulanger E, Dalongeville A, Andrello M, Mouillot D, Manel S. 2020. Spatial graphs highlight how multi-generational dispersal shapes landscape genetic patterns. Ecography (Cop). 43(8):1167-1179. doi:10.1111/ecog.05024.

Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. 2009. BLAST+: architecture and applications. BMC Bioinformatics. 10(1):421. doi:10.1186/1471-2105-10-421.

Candy JR, Campbell NR, Grinnell MH, Beacham TD, Larson WA, Narum SR. 2015. Population differentiation determined from putative neutral and divergent adaptive genetic markers in Eulachon (Thaleichthys pacificus, Osmeridae), an anadromous Pacific smelt. Mol Ecol Resour. 15(6):1421-1434. doi:10.1111/1755-0998.12400.

Capblancq T, Fitzpatrick MC, Bay RA, Exposito-Alonso M, Keller SR. 2020. Genomic Prediction of (Mal)Adaptation across Current and Future Climatic Landscapes. Annu Rev Ecol Evol Syst. 51:245-269. doi:10.1146/annurev-ecolsys-020720-042553.

Carreras C, García-Cisneros A, Wangensteen OS, Ordóñez V, Palacín C, Pascual M, Turon X. 2020. East is East and West is West: Population genomics and hierarchical analyses reveal genetic structure and adaptation footprints in the keystone species Paracentrotus lividus (Echinoidea). Divers Distrib. 26(3):382-398. doi:10.1111/ddi.13016.

Catchen JM, Amores A, Hohenlohe P, Cresko W, Postlethwait JH. 2011. Stacks: Building and Genotyping Loci De Novo From Short-Read Sequences. G3 Genes, Genomes, Genet. 1(3):171-182. doi:10.1534/g3.111.000240.

Cayuela H, Rougemont Q, Laporte M, Mérot C, Normandeau E, Dorant Y, Tørresen OK, Hoff SNK, Jentoft S, Sirois P, et al. 2020. Shared ancestral polymorphisms and chromosomal rearrangements as potential drivers of local adaptation in a marine fish. Mol Ecol. 29(13):2379-2398. doi:10.1111/mec.15499.

Claudet J, Osenberg CW, Domenici P, Badalamenti F, Milazzo M, Falcón JM, Bertocci I, Benedetti-Cecchi L, García-Charton J-A, Goñi R, et al. 2010. Marine reserves: Fish life history and ecological traits matter. Ecol Appl. 20(3):830-839. doi:10.1890/08-2131.1.

Coll M, Piroddi C, Steenbeek J, Kaschner K, Lasram FBR, Aguzzi J, Ballesteros E,
Bianchi CN, Corbera J, Dailianis T, et al. 2010. The biodiversity of the Mediterranean Sea: Estimates, patterns, and threats. PLoS One. 5(8). doi:10.1371/journal.pone.0011842.

D'Aloia CC, Andrés JA, Bogdanowicz SM, McCune AR, Harrison RG, Buston PM. 2020. Unraveling hierarchical genetic structure in a marine metapopulation: A comparison of three high-throughput genotyping approaches. Mol Ecol. 29(12):2189-2203. doi:10.1111/mec.15405.

Dalongeville A, Andrello M, Mouillot D, Albouy C, Manel S. 2016. Ecological traits shape genetic diversity patterns across the Mediterranean Sea: A quantitative review on fishes. J Biogeogr. 43(4):845-857. doi:10.1111/jbi.12669.

Dalongeville A, Andrello M, Mouillot D, Lobreaux S, Fortin M-J, Lasram F, Belmaker J, Rocklin D, Manel S. 2018. Geographic isolation and larval dispersal shape seascape genetic patterns differently according to spatial scale. Evol Appl. 11(8):1437-1447. doi:10.1111/eva.12638.

Dalongeville A, Benestan L, Mouillot D, Lobreaux S, Manel S. 2018. Combining six genome scan methods to detect candidate genes to salinity in the Mediterranean striped red mullet (Mullus surmuletus). BMC Genomics. 19(1):1-13. doi:10.1186/s12864-018-4579-z.

Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, et al. 2011. The variant call format and VCFtools. Bioinformatics. 27(15):2156-2158. doi:10.1093/bioinformatics/btr330.

Darmaraki S, Somot S, Sevault F, Nabat P, Cabos Narvaez WD, Cavicchia L, Djurdjevic V, Li L, Sannino G, Sein D V. 2019. Future evolution of Marine Heatwaves in the Mediterranean Sea. Clim Dyn. 53(3-4):1371-1392. doi:10.1007/s00382-019-04661-z.

Dayan DI. 2018. Clinal Adaptation in the Marine Environment. In: Oleksiak MF, Rajora OP, editors. Population Genomics: Marine Organisms. Springer. p. 221-247.

Do C, Waples RS, Peel D, Macbeth GM, Tillett BJ, Ovenden JR. 2014. NeEstimator v2: Re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Mol Ecol Resour. 14(1):209-214. doi:10.1111/1755-0998.12157.

Dray S, Bauman D, Blanchet G, Borcard D, Clappe S, Guenard G, Jombart T, Larocque G, Legendre P, Madi N, et al. 2020. adespatial: Multivariate Multiscale Spatial Analysis. R Packag version 03-8.

Dray S, Legendre P, Peres-Neto PR. 2006. Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). Ecol Modell. 196(3-4):483-493. doi:10.1016/j.ecolmodel.2006.02.015.
Duranton M, Bonhomme F, Gagnaire P-A. 2019. The spatial scale of dispersal revealed by admixture tracts. Evol Appl. 12(9):1743-1756. doi:10.1111/eva.12829.

Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buckler ES, Mitchell SE. 2011. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS One. 6(5):1-10. doi:10.1371/journal.pone.0019379.

Fietz K, Trofimenko E, Guerin P-E, Arnal V, Torres-Oliva M, Lobréaux S, Pérez-Rúzafa A, Manel S, Puebla O. 2020. New genomic resources for three exploited Mediterranean fishes. Genomics. 112(6):4297-4303. doi:10.1016/j.ygeno.2020.06.041.

Fitzpatrick MC, Keller SR. 2015. Ecological genomics meets community-level modelling of biodiversity: Mapping the genomic landscape of current and future environmental adaptation. Ecol Lett. 18(1):1-16. doi:10.1111/ele.12376.

Frazão Santos C, Agardy T, Andrade F, Calado H, Crowder LB, Ehler CN, García-Morales S, Gissi E, Halpern BS, Orbach MK, et al. 2020. Integrating climate change in ocean planning. Nat Sustain. 3(7):505-516. doi:10.1038/s41893-020-0513-x.

Froese R, Pauly D. 2019. FishBase.

Gagnaire P-A. 2020. Comparative genomics approach to evolutionary process connectivity. Evol Appl. 13(6):1320-1334. doi:10.1111/eva.12978.

Gagnaire P-A, Broquet T, Aurelle D, Viard F, Souissi A, Bonhomme F, Arnaud-Haond S, Bierne N. 2015. Using neutral, selected, and hitchhiker loci to assess connectivity of marine populations in the genomic era. Evol Appl. 8(8):769-786. doi:10.1111/eva.12288.

Gagnaire P-A, Gaggiotti OE. 2016. Detecting polygenic selection in marine populations by combining population genomics and quantitative genetics approaches. Curr Zool. 62(6):603-616. doi:10.1093/cz/zow088.

Gaines SD, White C, Carr MH, Palumbi SR. 2010. Designing marine reserve networks for both conservation and fisheries management. Proc Natl Acad Sci U S A. 107(43):18286-193. doi:10.1073/pnas.0906473107.

Galarza JA, Turner GF, Macpherson E, Rico C. 2009. Patterns of genetic differentiation between two co-occurring demersal species: the red mullet (Mullus barbatus) and the striped red mullet (Mullus surmuletus). Can J Fish Aquat Sci. 66(9):1478-1490. doi:10.1139/F09-098.

Garrison E, Marth G. 2012. Haplotype-based variant detection from short-read sequencing. arXiv:1207 [q-bioGN].

Gienapp P. 2020. Opinion: Is gene mapping in wild populations useful for understanding
and predicting adaptation to global change? Glob Chang Biol. 26(5):2737-2749.
doi:10.1111/gcb.15058.

Givan O, Edelist D, Sonin O, Belmaker J. 2018. Thermal affinity as the dominant factor changing Mediterranean fish abundances. Glob Chang Biol. 24(1):e80-e89.
doi:10.1111/gcb.13835.

González-Wangüemert M, Cánovas F, Pérez-Ruzafa A, Marcos C, Alexandrino P. 2010. Connectivity patterns inferred from the genetic structure of white seabream (Diplodus sargus L.). J Exp Mar Bio Ecol. 383(1):23-31. doi:10.1016/j.jembe.2009.10.010.

Hanson JO, Marques A, Veríssimo A, Camacho-Sanchez M, Velo-Antón G, Martínez-Solano Í, Carvalho SB. 2020. Conservation planning for adaptive and neutral evolutionary processes. J Appl Ecol. 57(11):2159-2169. doi:10.1111/1365-2664.13718.

Hijmans RJ. 2020. raster: Geographic Data Analysis and Modeling.

Hoey JA, Pinsky ML. 2018. Genomic signatures of environmental selection despite near-panmixia in summer flounder. Evol Appl. 11(9):1732-1747. doi:10.1111/eva.12676.

Jombart T. 2008. adegenet: a R package for the multivariate analysis of genetic markers.
Bioinforamtics. 24:1403-1405. doi:10.1093/bioinformatics/btn129.

Jombart T, Devillard S, Balloux F. 2010. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 11(1):94.
doi:10.1186/1471-2156-11-94.

Katsanevakis S, Mackelworth P, Coll M, Fraschetti S, Mačić V, Giakoumi S, Jones P, Levin N, Albano P, Badalamenti F, et al. 2017. Advancing marine conservation in European and contiguous seas with the MarCons Action. Res Ideas Outcomes. 3:e11884.
doi:10.3897/rio.3.e11884.

Larranaga N, Albertazzi FJ, Fontecha G, Palmieri M, Rainer H, van Zonneveld M, Hormaza JI. 2017. A Mesoamerican origin of cherimoya (Annona cherimola Mill.): Implications for the conservation of plant genetic resources. Mol Ecol. 26(16):4116-4130.
doi:10.1111/mec.14157.

Legendre P, Legendre L. 2012. Numerical Ecology. 3rd ed. Amsterdam: Elsevier Science BV.

Legendre P, Oksanen J, ter Braak CJF. 2011. Testing the significance of canonical axes in redundancy analysis. Methods Ecol Evol. 2(3):269-277. doi:10.1111/j.2041-210X.2010.00078.x.

Leigh DM, van Rees CB, Millette KL, Breed MF, Schmidt C, Bertola LD, Hand BK,
Hunter ME, Jensen EL, Kershaw F, et al. 2021. Opportunities and challenges of macrogenetic studies. Nat Rev Genet. 0123456789. doi:10.1038/s41576-021-00394-0.

Lejeusne C, Chevaldonné P, Pergent-Martini C, Boudouresque CF, Pérez T. 2010. Climate change effects on a miniature ocean: the highly diverse, highly impacted Mediterranean Sea. Trends Ecol Evol. 25(4):250-260. doi:10.1016/j.tree.2009.10.009.

Lewontin RC, Krakauer J. 1973. Distribution of gene frequency as a test of the theory of the selective neutrality of polymorphisms. Genetics. 74(1):175-195.
doi:10.1093/genetics/74.1.175.

Liggins L, Treml EA, Riginos C. 2019. Seascape Genomics: Contextualizing Adaptive and Neutral Genomic Variation in the Ocean Environment. In: Population Genomics: Marine Organisms. p. 171-218.

Luu K, Bazin E, Blum MGB. 2017. pcadapt: an R package to perform genome scans for selection based on principal component analysis. Mol Ecol Resour. 17(1):67-77.
doi:10.1111/1755-0998.12592.

Macpherson E, Raventos N. 2006. Relationship between pelagic larval duration and geographic distribution of Mediterranean littoral fishes. Mar Ecol Prog Ser. 327(Planes 2002):257-265. doi:10.3354/meps327257.

Manel S, Guerin P, Mouillot D, Blanchet S, Velez L, Albouy C, Pellissier L. 2020. Global determinants of freshwater and marine fish genetic diversity. Nat Commun. 11(1):692.
doi:10.1038/s41467-020-14409-7.

Manel S, Holderegger R. 2013. Ten years of landscape genetics. Trends Ecol Evol. 28(10):614-621. doi:10.1016/j.tree.2013.05.012.

Manel S, Joost S, Epperson BK, Holderegger R, Storfer A, Rosenberg MS, Scribner KT, Bonin A, Fortin MJ. 2010. Perspectives on the use of landscape genetics to detect genetic adaptive variation in the field. Mol Ecol. 19(17):3760-3772. doi:10.1111/j.1365-294X.2010.04717.x.

Manel S, Loiseau N, Andrello M, Fietz K, Goñi R, Forcada A, Lenfant P, Kininmonth S, Marcos C, Marques V, et al. 2019. Long-Distance Benefits of Marine Reserves: Myth or Reality? Trends Ecol Evol. 34(4):342-354. doi:10.1016/j.tree.2019.01.002.

Manel S, Perrier C, Pratlong M, Abi-Rached L, Paganini J, Pontarotti P, Aurelle D. 2016. Genomic resources and their influence on the detection of the signal of positive selection in genome scans. Mol Ecol. 25(1):170-184. doi:10.1111/mec.13468.

Marbà N, Jordà G, Agustí S, Girard C, Duarte CM. 2015. Footprints of climate change on
Mediterranean Sea biota. Front Mar Sci. 2(AUG):1-11. doi:10.3389/fmars.2015.00056.

Milano I, Babbucci M, Cariani A, Atanassova M, Bekkevold D, Carvalho GR, Espiñeira M, Fiorentino F, Garofalo G, Geffen AJ, et al. 2014. Outlier SNP markers reveal fine-scale genetic structuring across European hake populations (Merluccius merluccius). Mol Ecol. 23(1):118-135. doi:10.1111/mec.12568.

Mittelbach GG, Schemske DW, Cornell H V., Allen AP, Brown JM, Bush MB, Harrison SP, Hurlbert AH, Knowlton N, Lessios HA, et al. 2007. Evolution and the latitudinal diversity gradient: Speciation, extinction and biogeography. Ecol Lett. 10(4):315-331. doi:10.1111/j.1461-0248.2007.01020.x.

Morato T, Afonso P, Lourinho P, Nash RDM, Santos RS. 2003. Reproductive biology and
recruitment of the white sea bream in the Azores. J Fish Biol. 63(1):59-72. doi:10.1046/j.1095-8649.2003.00129.x.

Moritz C, Agudo R. 2013. The future of species under climate change: Resilience or decline? Science (80-). 341(6145):504-508. doi:10.1126/science.1237190.

Nielsen ES, Beger M, Henriques R, Selkoe KA, Heyden S. 2017. Multispecies genetic objectives in spatial conservation planning. Conserv Biol. 31(4):872-882. doi:10.1111/cobi.12875.

Nielsen ES, Henriques R, Beger M, von der Heyden S. 2021. Distinct interspecific and intraspecific vulnerability of coastal species to global change. Glob Chang Biol.(October 2020):1-17. doi:10.1111/gcb.15651.

Nielsen ES, Henriques R, Beger M, Toonen RJ, Von Der Heyden S. 2020. Multi-model seascape genomics identifies distinct environmental drivers of selection among sympatric marine species. BMC Evol Biol. 20(1):1-17. doi:10.1186/s12862-020-01679-4.

O'Leary SJ, Puritz JB, Willis SC, Hollenbeck CM, Portnoy DS. 2018. These aren't the loci you'e looking for: Principles of effective SNP filtering for molecular ecologists. Mol Ecol. 27(16):3193-3206. doi:10.1111/mec.14792.

Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O'Hare RB, Simpson GL, Solymos P, et al. 2016. vegan: Community Ecology Package. R Packag version 24-1.

Orsini L, Vanoverbeke J, Swillen I, Mergeay J, De Meester L. 2013. Drivers of population genetic differentiation in the wild: isolation by dispersal limitation, isolation by adaptation and isolation by colonization. Mol Ecol. 22(24):5983-5999. doi:10.1111/mec.12561.
Pascual M, Palero F, García-Merchán VH, Macpherson E, Robainas-Barcia A, Mestres F, Roda T, Abelló P. 2016. Temporal and spatial genetic differentiation in the crab Liocarcinus depurator across the Atlantic-Mediterranean transition. Sci Rep. 6(July):1–10.
doi:10.1038/srep29892.

Pascual M, Rives B, Schunter C, Macpherson E. 2017. Impact of life history traits on gene flow: A multispecies systematic review across oceanographic barriers in the Mediterranean Sea. Chiang T-Y, editor. PLoS One. 12(5):e0176419. doi:10.1371/journal.pone.0176419.

Patarnello T, Volckaert FAMJ, Castilho R. 2007. Pillars of Hercules: Is the Atlantic-Mediterranean transition a phylogeographical break? Mol Ecol. 16(21):4426–4444.
doi:10.1111/j.1365-294X.2007.03477.x.

Pauls SU, Nowak C, Bálint M, Pfenninger M. 2013. The impact of global climate change on genetic diversity within populations and species. Mol Ecol. 22(4):925–946.
doi:10.1111/mec.12152.

Pembleton LW, Cogan NOI, Forster JW. 2013. StAMPP: an R package for calculation of genetic differentiation and structure of mixed-ploidy level populations. Mol Ecol Resour.
13(5):946–952. doi:10.1111/1755-0998.12129.

Peres-Neto PR, Legendre P, Dray S, Borcard D. 2006. Variation partitioning of species data matrices: Estimation and comparison of fractions. Ecology. 87(10):2614–2625.
doi:10.1890/0012-9658(2006)87[2614:VPOSDM]2.0.CO;2.

Puritz JB, Hollenbeck CM, Gold JR. 2014. dDocent : a RADseq, variant-calling pipeline designed for population genomics of non-model organisms. PeerJ. 2(1):e431.
doi:10.7717/peerj.431.

Raventos N, Torrado H, Arthur R, Alcoverro T, Macpherson E. 2021. Temperature reduces fish dispersal as larvae grow faster to their settlement size. Eizaguirre C, editor. J Anim Ecol. 88080:1365-2656.13435. doi:10.1111/1365-2656.13435.

Razgour O, Forester B, Taggart JB, Bekaert M, Juste J, Ibáñez C, Puechmaille SJ, Novella-Fernandez R, Alberdi A, Manel S. 2019. Considering adaptive genetic variation in climate change vulnerability assessment reduces species range loss projections. Proc Natl Acad Sci U S A.
116(21):10418–10423. doi:10.1073/pnas.1820663116.

Rellstab C, Gugerli F, Eckert AJ, Hancock AM, Holderegger R. 2015. A practical guide to environmental association analysis in landscape genomics. Mol Ecol. 24(17):4348–4370.
doi:10.1111/mec.13322.
Reñones O, Massuti E, Morales-Nin B. 1995. Life history of the red mullet Mullus surmuletus from the bottom-trawl fishery of the Island of Majorca (north-west Mediterranean). Mar Biol. 123:411-419.

Rilov G, Fraschetti S, Gissi E, Pipitone C, Badalamenti F, Tamburello L, Menini E, Goriup P, Mazaris AD, Garrabou J, et al. 2020. A fast-moving target: achieving marine conservation goals under shifting climate and policies. Ecol Appl. 30(1):1-14. doi:10.1002/eap.2009.

Robitzch VSN, Lozano-Cortés D, Kandler NM, Salas E, Berumen ML. 2016. Productivity and sea surface temperature are correlated with the pelagic larval duration of damselfishes in the Red Sea. Mar Pollut Bull. 105(2):566-574. doi:10.1016/j.marpolbul.2015.11.045.

Schoville SD, Bonin A, François O, Lobreaux S, Melodelima C, Manel S. 2012. Adaptive Genetic Variation on the Landscape: Methods and Cases. Annu Rev Ecol Evol Syst. 43(1):23-43. doi:10.1146/annurev-ecolsys-110411-160248.

Selkoe KA, D'Aloia CC, Crandall ED, Iacchei M, Liggins L, Puritz JB, Von Der Heyden S, Toonen RJ. 2016. A decade of seascape genetics: Contributions to basic and applied marine connectivity. Mar Ecol Prog Ser. 554:1-19. doi:10.3354/meps11792.

Selmoni O, Rochat E, Lecellier G, Berteaux-Lecellier V, Joost S. 2020. Seascape genomics as a new tool to empower coral reef conservation strategies: An example on north-western Pacific Acropora digitifera. Evol Appl. 13(8):1923-1938. doi:10.1111/eva.12944.

Stanley RRE, DiBacco C, Lowen B, Beiko RG, Jeffery NW, Van Wyngaarden M, Bentzen P, Brickman D, Benestan L, Bernatchez L, et al. 2018. A climate-associated multispecies cryptic cline in the northwest Atlantic. Sci Adv. 4(3):eaaq0929. doi:10.1126/sciadv.aaq0929.

Storey JD, Bass AJ, Dabney A, Robinson D. 2019. qvalue: Q-value estimation for false discovery rate control. R Packag version 2180.

Thomas E, van Zonneveld M, Loo J, Hodgkin T, Galluzzi G, van Etten J. 2012. Present Spatial Diversity Patterns of Theobroma cacao L. in the Neotropics Reflect Genetic Differentiation in Pleistocene Refugia Followed by Human-Influenced Dispersal. PLoS One. 7(10):e47676. doi:10.1371/journal.pone.0047676.

Tomanek L, Zuzow MJ, Hitt L, Serafini L, Valenzuela JJ. 2012. Proteomics of hyposaline stress in blue mussel congeners (genus mytilus): Implications for biogeographic range limits in response to Climate change. J Exp Biol. 215(22):3905-3916. doi:10.1242/jeb.076448.

Torrado H, Carreras C, Raventos N, Macpherson E, Pascual M. 2020. Individual-based population genomics reveal different drivers of adaptation in sympatric fish. Sci Rep. 10(1):12683.
doi:10.1038/s41598-020-69160-2.

Tyberghein L, Verbruggen H, Pauly K, Troupin C, Mineur F, De Clerck O. 2012. Bio-ORACLE: A global environmental dataset for marine species distribution modelling. Glob Ecol Biogeogr. 21(2):272-281. doi:10.1111/j.1466-8238.2011.00656.x.

Vandamme S, Raeymaekers JAM, Maes GE, Cottenie K, Calboli FCF, Diopere E, Volckaert FAM. 2021. Reconciling seascape genetics and fisheries science in three codistributed flatfishes. Evol Appl. 14(2):536-552. doi:10.1111/eva.13139.

Vandamme SG, Maes GE, Raeymaekers JAM, Cottenie K, Imsland AK, Hellemans B, Lacroix G, Mac Aoidh E, Martinsohn JT, Martínez P, et al. 2014. Regional environmental pressure influences population differentiation in turbot (Scophthalmus maximus). Mol Ecol. 23(3):618-636. doi:10.1111/mec.12628.

Vargas-Yáñez M, Moya F, García-Martínez MC, Tel E, Zunino P, Plaza F, Salat J, Pascual J, López-Jurado JL, Serra M. 2010. Climate change in the Western Mediterranean Sea 1900-2008. J Mar Syst. 82(3):171-176. doi:10.1016/j.jmarsys.2010.04.013.

Vasimuddin M, Sanchit M, Heng L, Srinivas A. 2019. Efficient architecture-aware acceleration of BWA-MEM for multicore systems. Proc - 2019 IEEE 33rd Int Parallel Distrib Process Symp IPDPS 2019.:314-324. doi:10.1109/IPDPS.2019.00041.

Waldvogel A, Feldmeyer B, Rolshausen G, Exposito-Alonso M, Rellstab C, Kofler R, Mock T, Schmid K, Schmitt I, Bataillon T, et al. 2020. Evolutionary genomics can improve prediction of species' responses to climate change. Evol Lett. 4(1):4-18. doi:10.1002/evl3.154.

Wang IJ, Bradburd GS. 2014. Isolation by environment. Mol Ecol. 23(23):5649-5662. doi:10.1111/mec.12938.

Waples RS. 2006. A bias correction for estimates of effective population size based on linkage disequilibrium at unlinked gene loci. Conserv Genet. 7(2):167-184. doi:10.1007/s10592-005-9100-y.

Waples RS, Do C. 2010. Linkage disequilibrium estimates of contemporary Ne using highly variable genetic markers: A largely untapped resource for applied conservation and evolution. Evol Appl. 3(3):244-262. doi:10.1111/j.1752-4571.2009.00104.x.

Weir BS, Cockerham CC. 1984. Estimating F-statistics for the analysis of population structure. Evolution (N Y). 38(6):1358-1370. doi:10.1088/0022-3719/13/31/023.

Whitehead A, Galvez F, Zhang S, Williams LM, Oleksiak MF. 2011. Functional genomics of physiological plasticity and local adaptation in killifish. J Hered. 102(5):499-511.
doi:10.1093/jhered/esq077.

Xuereb A, D'Aloia CC, Andrello M, Bernatchez L, Fortin M. 2020. Incorporating putatively neutral and adaptive genomic data into marine conservation planning. Conserv Biol. 0(0):cobi.13609. doi:10.1111/cobi.13609.

Xuereb A, D'Aloia CC, Daigle RM, Andrello M, Dalongeville A, Manel S, Mouillot D, Guichard F, Côté IM, Curtis JMR, et al. 2019. Marine Conservation and Marine Protected Areas. In: Population Genomics: Marine Organisms. p. 423–446.
<figure><img src="image_2.png" /></figure>

FIGURES

<caption>**Figure 1.** a) Sampling sites for white seabream (*Diplodus sargus*) and striped red mullet (*Mullus surmuletus*). Black dots indicate sites where both species were sampled, whereas silhouettes indicate sites where only the respective species was sampled. Mediterranean ecoregions are shown in different shades. AOF denotes the Almeria-Oran Front. b-c) Variation in winter temperature (b) and salinity (c) levels at sampled sites across ecoregions. The letters</caption>
indicate grouping after a Kruskal-Wallis test, where different letters indicate significantly different means.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Spatial distribution of genetic diversity. The maps show the spatial interpolation of expected heterozygosity ($H_e$) of (a) seabream and (b) mullet, calculated per cell based on imputed individuals. (c) Individual heterozygosity ($H_o$) confirms the spatial variation of the seabream (dots) and mullet (triangles) genetic diversity across ecoregions (mean ± standard error). The letters indicate grouping of ecoregions after a Kruskal-Wallis test, where different letters indicate significantly different means.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Discriminant analysis of principal components with outlier loci of (a) white seabream and (b) striped red mullet. Individuals are ordered by their longitudinal position within each ecoregion (ALB = Alboran Sea, WMED = Western Mediterranean, ADR = Adriatic Sea, ION = Ionian Sea, TUN = Tunisian Plateau, AEG = Aegean Sea, LEV = Levantine Sea).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Distance-based redundancy analysis for outlier loci of (a) white seabream and (b) striped red mullet. Arrows represent environmental variables that drive the observed population structure. MEM vectors are distance-based Moran Eigenvector Maps representing geographical isolation at different spatial scales while Asymmetric Eigenvector Maps (AEM) represent isolation by larval dispersal. Dots represent individuals while colours correspond to their clustering obtained with the DAPC in Figure 3.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Relative importance of the spatial (SPACE; dbMEMs), dispersal (DISP; AEMs) and environmental (ENV; salinity, chlorophyll a, summer – and winter temperature) variables in explaining outlier genetic variation (i.e. Euclidean genetic distances estimated between individuals) of (a) white seabream and (b) striped red mullet. The surface of the different sections of the Euler diagram is proportional to the R² of each partition.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Example of two SNPs located on candidate genes, (a) zmynd10 in the seabream</figcaption></figure>

and (b) TBCD in the mullet, and the prevalence of alleles along temperature gradients.
<table><caption>**Table 1.** Partial dbRDA results for white seabream and striped red mullet using outlier SNPs. Variables are selected through an ordiR2step procedure, and their significance and model fit are assessed in separate partial dbRDAs. \* indicates significance levels from ANOVA with 9999 permutations (\*P < 0.05, \*\*P < 0.01, \*\*\*P < 0.001).</caption><thead><tr><th>Species</th><th>Variable</th><th>P model</th><th>R<sup>2</sup> adj model</th></tr></thead><tbody><tr><td rowspan="3">White seabream<br><i>Diplodus sargus</i></td><td>Salinity</td><td>0.039 *</td><td>0.002</td></tr><tr><td>Summer temperature</td><td>1e-04 ***</td><td>0.014</td></tr><tr><td>Winter temperature</td><td>0.054</td><td>0.001</td></tr><tr><td rowspan="2">Striped red mullet<br><i>Mullus surmuletus</i></td><td>Chlorophyll a</td><td>1e-04 ***</td><td>0.01</td></tr><tr><td>Summer temperature</td><td>1e-04 ***</td><td>0.007</td></tr></tbody></table>