Sardines at a junction: Seascape genomics reveals ecological and oceanographic drivers of variation in the NW Mediterranean Sea

By evaluating genetic variation across the entire genome, one can address existing questions in a novel way while raising new ones. The latter includes how different local environments influence adaptive and neutral genomic variation within and among populations, providing insights into local adaptation of natural populations and their responses to global change. Here, under a seascape genomic approach, ddRAD data of 4609 single nucleotide polymorphisms (SNPs) from 398 sardines (Sardina pilchardus) collected in 11 Mediterranean and one Atlantic site were generated. These were used along with oceanographic and ecological information to detect signals of adaptive divergence with gene flow across environmental gradients. The studied sardines constitute two clusters (FST = 0.07), a pattern attributed to outlier loci, highlighting putative local adaptation. The trend in the number of days with sea surface temperature above 19°C, a critical threshold for successful sardine spawning, was crucial at all levels of population structuring with implications on the species' key biological processes. Outliers link candidate SNPs to the region's environmental heterogeneity. Our findings provide evidence for a dynamic equilibrium in which population structure is maintained by physical and ecological factors under the opposing influences of migration and selection. This dynamic in a natural system warrants continuous monitoring under a seascape genomic approach that might benefit from a temporal and more detailed spatial dimension. Our results may contribute to complementary studies aimed at providing deeper insights into the mechanistic processes underlying population structuring. Those are key to understanding and predicting future changes and responses of this highly exploited species in the face of climate change.


| INTRODUC TI ON
Identifying genomic signatures of selection provides insights not only into the local adaptation of natural populations but also into their adaptive responses to global change (Carroll et al., 2014;Hauser & Carvalho, 2008;Heino et al., 2015). This can be better achieved through population genomic studies, by separating the effects of neutral processes (drift-migration) from those of selection, where frequency shifts point to rapid changes in environmental conditions (Hoffmann & Willi, 2008).
Modern genomic tools provide the means to assess the population structure of species at high resolution and gain a detailed view of their status (Greenbaum et al., 2019), supplying new means for examining the genetic basis of local adaptation (Price et al., 2020).
Combining state-of-the-art genomic approaches with spatially explicit information on the main habitat features allows linking genotypes to environmental variables, within the field of seascape genomics (Riginos et al., 2016). This enhances our understanding of the processes affecting inhabitant populations, including adaptive responses to environmental fluctuations reflecting either environmental heterogeneity or climate change (Selkoe et al., 2010), and providing valuable clues on how these populations cope with environmental change, in the face of ocean warming (Liggins et al., 2019).
The Mediterranean Sea has experienced a long history of overexploitation with one of the world's highest percentage of stocks fished at unsustainable levels (FAO, 2020) and anthropogenic perturbations (Coll et al., 2012). It is strongly affected by pollution (García-Rivera et al., 2017;UNEP/MAP, 2015). In addition to the human-induced loss of habitat and biodiversity through direct exploitation (Coll et al., 2015), indirect effects such as climate change and Lessepsian migration/invasive species, pose additional threats that will probably increase in severity in the future Ramírez et al., 2017).
A shift towards sustainable fisheries in the Mediterranean Sea is more urgent than ever and requires robust fisheries and stock assessments, which means that biological and management units are aligned (Reiss et al., 2009;Vasilakopoulos et al., 2014). This is intended to preserve stocks with adaptive diversity that can contribute to the long-term maintenance of productive fisheries and ecosystems (Antoniou & Magoulas, 2014).
Because of their commercial value and ecological role within marine communities, small pelagic fish are important species in the Mediterranean Sea FAO, 2018;Palomera et al., 2007). Overall, declines of commercially important small pelagic fish populations, particularly of the European pilchard (Sardina pilchardus, (Walbaum, 1792); Clupeidae; Clupeiformes, hereafter sardine), have been observed in the Mediterranean Sea mostly attributed to the combined effects of climate impacts and fishing pressure Pennino et al., 2020;Peristeraki et al., 2019;Ramírez et al., 2018Ramírez et al., , 2021Saraux et al., 2019).
Highly mobile sardine, similarly to other clupeids, display localized batch spawning behaviour and migratory patterns that may result in restricted gene flow between populations and subsequent divergence (Bacha et al., 2014). This has already been demonstrated within the Mediterranean and to a lesser extent in the eastern Atlantic, mainly reflecting the levels of environmental heterogeneity of the two regions (Caballero-Huertas et al., 2022, and references therein;Ganias et al., 2004). The prevailing conditions are distinct in the two areas in terms of productivity, circulation patterns, wind events, water temperature, and biogeochemical and biological components but also topographic features that greatly impact the species' spatiotemporal variability in spawning (Ganias, 2009;Lima et al., 2022). Sardine spawning is temperature-dependent and conducted under a narrow range of temperatures (preferred range 12-14°C, also occurring up to 19°C) along their distribution range (Palomera et al., 2007), taking advantage of shelf upwelling systems with cold water and enhanced productivity, and under the influence of favourable winds that minimize offshore transport (Abdelouahab et al., 2021;Santos, Brito, et al., 2006;Santos, Ré, et al., 2006;Stratoudakis et al., 2007). Spawning is presumably triggered by decreasing water temperature in the autumn (<19°C) and completed as temperature increases again in the spring (Fernández-Corredor et al., 2021;Palomera et al., 2007;Stratoudakis et al., 2007). Multiple lines of evidence suggest a poleward migration and spawning habitat e Innovación, Grant/Award Number: RYC2020-030078-I; European Regional Development Fund Handling Editor: Paul A. Hohenlohe processes. Outliers link candidate SNPs to the region's environmental heterogeneity.
Our findings provide evidence for a dynamic equilibrium in which population structure is maintained by physical and ecological factors under the opposing influences of migration and selection. This dynamic in a natural system warrants continuous monitoring under a seascape genomic approach that might benefit from a temporal and more detailed spatial dimension. Our results may contribute to complementary studies aimed at providing deeper insights into the mechanistic processes underlying population structuring. Those are key to understanding and predicting future changes and responses of this highly exploited species in the face of climate change.

K E Y W O R D S
climate change, ddRAD, outliers, population genomics, SNPs shift as a response to rising temperatures (Alheit et al., 2012;Lima et al., 2022;Palomera et al., 2007). Genetic studies conducted on sardines in the Mediterranean employing different types of markers have mostly reported shallow phylogeographical structure, low genetic differentiation, and weak or absent population structure, with signs of late Pleistocene expansion (Kasapidis, 2014, and references therein). Significant genetic structure was only revealed in one particular locus and a microsatellite probably linked to it. The pattern observed coincided with a genetic cline tightly connected to the local thermal and hydrodynamic characteristics of the Atlantic-Mediterranean junction (Chlaida et al., 2009;Kasapidis et al., 2012;Laurent et al., 2007;Patarnello et al., 2007).
Within the Mediterranean, the Almeria-Oran oceanographic front (AOF) is a sharp semipermanent boundary between surface water masses formed by the convergence of the Atlantic Ocean and the Mediterranean Sea. It is a key oceanographic feature, situated at the Northeastern Alboran Sea that separates the Mediterranean waters from northeastward inflow of Atlantic waters, reaching the Balearic Islands and forming a second density front in the northern part of the archipelago (Galarza et al., 2009). Its hydrological characteristics are closer to those of the Northeastern Atlantic than to the Western Mediterranean (Bacha et al., 2014), indicative of enhanced gene flow between the Alboran Sea and the Atlantic (Northeast and Moroccan Atlantic Ocean, Baibai et al., 2012). The AOF constitutes a barrier to sardine dispersal that could be held responsible for founder effects, followed by genetic drift with selection acting on the local scale driving adaptations mostly related to minimum sea surface temperature (Baltazar-Soares et al., 2021). The AOF as a barrier that separates the Atlantic and the Mediterranean is commonly observed in other fish species, such as Sardinella aurita, Mugil cephalus and Xiphias gladius (Borsa et al., 1997); Merluccius merluccius (Lundy et al., 1999); three sparid species, Lithognatus mormyrus, Spondyliosoma cantharus and Dentex dentex (Bargelloni et al., 2003); horse mackerel (Trachurus trachurus) populations (Abaunza et al., 2008) and Dicentrarchus labrax (Naciri et al., 1999).
Despite the wealth of previous studies, the role of multiple stressors in shaping the stock structure of sardines within the Mediterranean Sea remains elusive. This prohibits inferences of local adaptive variation and the potential of populations to respond to environmental changes, with implications for the conservation and sustainable management of this declining fish stock. Spatial heterogeneity is expected to generate adaptive divergence by imposing challenges to the organisms inhabiting the area, further enhanced by changing environmental conditions (Grummer et al., 2019;Sandoval-Castillo & Beheregaray, 2020). The latter, and especially changes in sea surface temperature and the occurrence of winter upwelling events, have been held responsible for the recent decline of sardine recruitment (Santos et al., 2018) which in conjunction with spawning and larval behaviour constitute the main drivers of recruitment variability (Lowerre-Barbieri et al., 2017, and references therein).
Considering all the above, and given the fragmented information on sardine biology, ecology and evolution (Coll & Bellido, 2019), with this study we aim to provide a better understanding of the processes affecting sardine populations and to define environmental variables with high impact in shaping the species genomic patterns. Populations may show evident signatures of local adaptation regardless of gene flow (Diopere et al., 2018;Moody et al., 2015).
Furthermore, adaptive variation is expected to be maintained in cases where the selection of locally adapted alleles is stronger than the influx of nonlocally adapted alleles (Bal et al., 2021). Given the above, we conducted an exploratory analysis aimed at finding the most relevant variables to the adaptive potential of sardine as well as to its susceptibility to climate change. By analysing genome-wide single nucleotide polymorphisms (SNPs) within a seascape genomics approach, we evaluated the influence of climatological ocean conditions and their trends, as well as landscape characteristics, biological components and exploitation based on different gear types, on the spatial patterns of population genetic structure across the environmentally heterogeneous region of the Central-Western Mediterranean and adjacent Atlantic waters. This set of variables includes well-established species-specific drivers and variables never tested before. Our approach enabled the assessment of the species population structure as well as genome-environment associations by disentangling the patterns of adaptive and neutral genetic variation. We further attempted to explain whether observed patterns in size and body condition, that is the north/south gradient in size and body condition observed in the Western Mediterranean Bachiller et al., 2020;Brosset et al., 2017;Coll & Bellido, 2019), reflect random genetic processes and demography, or selection imposed by environmental forcing (both the biotic and abiotic attributes of the areas under study), or their combined effect.
The methodology and findings of this study can contribute to better understand processes shaping population differentiation and to evaluate the adaptive potential of species in response to environmental change. It is also anticipated to shed light on how the interplay between natural selection and gene flow affects local adaptation in this broadcast spawner, especially at times of increasing anthropogenic influence on the natural world (Flanagan et al., 2018).  Table S1).

| Sampling, DNA extraction, library construction and genotyping
Genomic DNA was extracted from tissue samples and libraries were prepared following the double-digest restriction siteassociated DNA (ddRAD) sequencing approach of Kess et al. (2016).
Sequenced reads were analysed using stacks version 2.4 pipeline (Catchen et al., 2013), to quality control the reads, identify the genomic loci sequenced, genotype each individual and conduct basic population genetic analyses against the genome of sardine (Louro et al., 2019). A detailed description on sampling, library construction and genotyping is provided in the Supporting Information.

| Population genetic analyses
Population genetic structure was assessed by both Bayesian and multivariate ordination methods. First, a model-based clustering was performed with structure version 2.3.4 (Pritchard et al., 2000) under the correlated allele frequency model allowing admixture. In cases where K > 1, samples allocated to clusters with high membership coefficient (q ≥ 0.9) were further analysed with structure for each cluster and with the same settings as above, in an attempt to examine whether further substructuring occurs.
Second, a discriminant analysis of principal components (DAPC) was performed with the adegenet version 2.1.1 package (Jombart, 2008), in the R environment (R Core Team, 2020) as a means to infer population subdivision of the samples under study, with an independent of population genetics model method (for a detailed description see Appendix S1).
Genetic diversity between the 12 sampling localities and the identified clusters was then compared in terms of observed (H O ) and expected heterozygosity (H E ), as a means of genetic differentiation using the stacks population program. Furthermore, levels of differentiation among localities and clusters were assessed by the fixation index F ST , with statistical significance evaluated through 10,000 permutations with the software stratag R package (Archer et al., 2017), and p values were adjusted for multiple testing through false discovery rate (FDR) correction (Benjamini & Hochberg, 1995).
A detailed description of the procedure followed is given in the Supporting Information.

| The data sets
Four types of genomic data sets were analysed to detect whether the different levels of population structure and the observed contrast in physical condition are driven by distinct environmental drivers. The first one included all studied samples ("all samples," i.e. a data set that displayed population structure, see Results). The second and third data sets corresponded to the population clusters discovered in the studied area, that is the "Atlantic" [ATL] and "Mediterranean" [MED] clusters, respectively (Section 3). Finally, the fourth data set contained samples from northern and southern sampling sites of the Western Mediterranean Sea ("northern vs. southern sites," i.e., GSA07a, GSA07b, GSA06a vs. GSA06c) that displayed remarkable differences in their physical condition, namely fish length and weight. This difference had a north/south gradient with larger and heavier sardines found in the southernmost areas (Table S1; Figure 1). Although those differences might be attributed to the sampling strategy itself (both time and gear employed for sampling), as well as the inclusion of six immature samples at the F I G U R E 1 Geographical distribution of sampling sites. Colours in pies indicate the proportion of males, females and immatures captured, while pie size is proportional to mean length of sardines. Bathymetry of the studied area is also indicated with a gradient of blue colours with darker colours for deeper seas.

GSA05
Gulf of Cadiz northern sites, this seems less likely given that sampling was conducted during the same oceanographic surveys, and the consistency of such differences with recent studies reporting a similar trend in maximum size and body condition of sardine in the area Bachiller et al., 2020;Brosset et al., 2017 (Guillot et al., 2014). To minimize the occurrence of false positives, loci with a log 10 Bayes factor (logBF) of ≥3 were interpreted as having an outstanding statistical dependence with a certain environmental variable and therefore likely to belong to a genomic region under selection.
Environmental variables were selected considering evolutionary history, climate variability or change, and fishing activity as the dominant processes affecting small pelagic fish. Variables were selected to include all previously reported variables and their trends with significant impact on sardines, such as temperature (Garrido et al., 2017;Gordó-Vilaseca et al., 2021), salinity (Santos, Ré, et al., 2006) and hydrodynamics (for a review see Caballero-Huertas et al., 2022), as well as variables related to the biology of the species whose impact has never previously been assessed. Given that the impacts of climate change are a rather complex issue that refers to more than one parameter (e.g., an increase in temperature), we tested whether there are any environmental variables related to topography, hydrodynamics, and biochemical and biological components at different layers (i.e., surface, water column, bottom), different types of exploitation tools and food availability acting as selective agents in sardines of the Northwestern Mediterranean (Table 1; Table S2, see Supplementary Information for a detailed description). These were conducted under an exploratory framework to pinpoint significant factors shaping sardine genomic patterns that have previously been neglected, with subsequent implications for future directions of global change.
Outlier detection was also conducted by a principal component analysis (PCA)-based approach on individual genotype data using the pcadapt r package (Luu et al., 2017) and by redundancy analysis (RDA) implemented in the r package vegan, which was employed as a multivariate method to detect environmental variables relevant for sardine as well as outlier loci (Oksanen et al., 2012). The RDA was conducted in two data sets with different levels of population structure. Multiple data sets were analysed given the not well-known performance of RDA in systems with increased levels of population structure or metapopulation dynamics .
The two data sets included "all samples" and MED while for the remaining two the analysis was not possible due to the low number of sites that resulted in collinearity of all the environmental variables.
Environmental covariates were selected to minimize collinearity among them using variance inflation factors (VIF < 10). Significant constrained axes were identified using 999 permutations of the response data and a p-value threshold of .05. We identified candidate adaptive loci as SNPs loading ±3 SD from the mean loading of these significant RDA axes. We then identified the covariate most strongly correlated with each candidate SNP (i.e., highest correlation coefficient), to group candidates by potential driving environmental variables. Detailed information for this section is provided in the Supporting Information.

| Neutral and putatively adaptive genomic variation
To minimize the detection of false positives, SNPs that were selected by all methods (ginland, pcadapt, RDA) were considered as outliers.
Furthermore, the highest number of overlapping SNPs among pairs of methods was also considered indicative of the signal in our data and provided better insight into the acquired results. Similarly, a putative neutral SNP set included only SNPs that were not highlighted as outliers by any of the methods employed. This yielded two different sets of loci, namely putatively under selection and putatively neutral. The two sets of loci were analysed with structure version 2.3.4 (Pritchard et al., 2000), using the same settings as above, in order to define whether different signals of population structure could be revealed.

| Outlier loci functional annotation
The two groups of outlier loci (i.e., selected by all methods and by the two methods with the highest number of overlapping loci) were further analysed to identify which of these loci were located within known genic regions and, for those found in intergenic regions, which is the closest gene using the genome annotation files (gff) from the sardine genome portal (https://bioin forma tics.psb.ugent.be/gdb/ Spil/). The protein sequences of the genes recovered were used in a blastp similarity search (e-value threshold 10 −8 ) against the swissprot database. To identify potential biological functions involved in the studied populations' response to environmental variables, a functional annotation analysis was conducted (for details see Supporting Information).
All computationally intensive analyses of the study were conducted at the IMBBC HPC facility (Zafeiropoulos et al., 2021).

| ddRAD sequencing and data analysis
We sequenced one billion paired reads from 398 sardine individu-

| Assessing the population genetic structure
According to structure, the optimal number of clusters that best fit the data was K = 2 (with all ten replicates providing identical The layer source is also provided. Further detailed information can also be found in Table S2. To detect adaptations of the two differentiated gene pools, as a strict measure of defining as pure clusters as possible, seascape genomic analyses were conducted based on structure's allocation to clusters (i.e., on the ATL and MED clusters), without taking into consideration the admixed samples.

| Seascape genomic analyses
The overlapping outliers among all three methods were five, with RDA analysis conducted in only two out of the four data sets. The highest number of overlapping outlier loci was that between ginland and pcadapt that reached 196 and included the five loci that In blue, the cluster that mostly contains samples from the Atlantic (similar to ATL cluster of structure analyses) and in red, the cluster that mostly contains samples from the Mediterranean (similar to MED cluster of structure analyses).
overlapped among all employed methods. ginland and RDA shared six loci while pcadapt and RDA had 28 loci in common ( Figure S1).
Across data sets with all levels of population structure (i.e., "all samples" data set that displayed population structure, and ATL and MED data sets at the lower level of population structure), there were no outlier loci in common when ginland was used. By contrast, when pcadapt was used, a core of 33 loci was shared among the three data sets. The two data sets analysed with RDA shared 14 loci ( Figure S1).
In the "all samples" data set, ginland detected 262 outlier loci with logBF >3, indicative of strong evidence of selection (Kass & Raftery, 1995). These outlier loci were associated with 39 environmental variables with none of them related to topography or fishing effort. Three groups of loci exhibited shared patterns of association across multiple environmental variables ( Figure 4)  sites were positively related to nitrate, plan curvature and photosynthetically active radiation (parmean) respectively. Furthermore, a distinction of northern (GSA07a, GSA07b, GSA06a) versus southern sites (GDS06c) was evident with a positive relationship with the number of days with SST < 12°C (sst_tr12) and with maximum annual SST (sst_max_sl) respectively (Figure 8a,b).  Table S2. In the ATL data set, ginland detected 37 outlier loci with logBF >3 in 24 environmental variables. In this data set, ginland highlighted variables related to topography covarying with SNP data. A group of loci exhibited shared patterns of association across the majority of environmental variables, including those related to topography ( Figure 6). pcadapt analysis detected 957 loci for the same data set, of which five were common with ginland's list of outlier loci.
The "northern vs. southern" data set yielded 45 outlier loci with logBF >3 in ginland analysis, involving seven environmental variables. One variable (i.e., sst_max_sl) displayed correlation with most of the highlighted SNPs (23 SNPs). The next most important variables in terms of the number of highlighted loci were cum_impact †, sst_tr19 and tot_impact (Figure 7). One SNP was found to covary with clim_impact †. pcadapt for the same data set detected 218 outliers, three of which in common with ginland.
Common outliers among employed approaches and data sets are depicted in Figure S1. A more detailed description of the results is provided in the Supporting Information.
Overall, from all examined data sets and all employed methods (i.e., ginland, RDA and pcadapt), 1607 SNPs were highlighted as potential outliers. Almost all environmental variables used in this study correlated with at least one SNP and up to 156 (see Table S4). Four variables that did not correlate with any of the SNPs in any of the data sets were: bathymetry, biogeo01, biogeo13 and Berline's re- Population structure analyses on ginland and pcadapt common outliers (n = 196, including the 5 common SNPs by all methods) estimated K = 2 clusters and allocated the samples to the MED and ATL clusters as with the complete loci data set. By contrast, F I G U R E 8 Triplots of the "all samples" data set for RDA axes 1 and 2 (a), and RDA axes 1 and 3 (b). Triplots of the MED data set for RDA axes 1 and 2 (c), and RDA axes 1 and 3 (d). The dark grey cloud of points at the centre of each plot represents the SNPs. The coloured points represent individual sardines with coding by sampling site. Blue vectors represent environmental predictors (see Table 1 for abbreviations).
when using only the neutral set of loci (i.e., 3002 loci [deducting 1607 outliers from the total 4609 loci]), no signal of population structure was detected.

| Functional annotation
Out of the five outlier loci in common through all approaches, four were found in scaffolds with genes in their vicinity: Heterogeneous nuclear ribonucleoprotein A1 (hnRNP A1), G protein-activated inward rectifier potassium channel 2 (GIRK-2), Myosin regulatory light chain 2, ventricular/cardiac muscle isoform (MLC-2) and Transmembrane protein 178B. Functional annotation analyses of the 196 outlier loci highlighted in any of the data sets by ginland and pcadapt approaches revealed that 156 are located within scaffolds that contain predicted genes. Few of them (i.e., 12) were found within genic regions while for the rest the closest gene was identified (top blast hit per gene is given in Table S5). Ten loci were located within the same gene with the rest being quite distant (Table S5). Gene Ontology (GO) terms were retrieved for 140 of these genes ( Figure S2).
See also Appendix S1 for a detailed description of the Results.

| Understanding the population structure of sardines
The present study indicates the existence of two clearly differentiated gene pools of sardines that meet in the Western Mediterranean (i.e., ATL and MED clusters). The results of population structure analysis tools (i.e., structure and DAPC) were consistent despite their methodological differences, indicative of the strong signal in our data. The two clusters are geographically well defined, in agreement with the geographical distribution of the two sardine subspecies recognized by morphological and molecular characters (Atarhouch et al., 2007;Fonseca et al., 2021;Parrish et al., 1989). Admixed individuals are found throughout the Mediterranean sampling sites, with the AOF having a transition role between the two regions, in agreement with previous studies (e.g., Atarhouch et al., 2007). Our results suggest that the AOF acts as the actual barrier separating the Atlantic and the Mediterranean sardines allowing for some interaction as reported in previous studies (Santos et al., 2018, and references therein). This pattern coincides with that observed in other fish species (e.g. Abaunza et al., 2008;Bargelloni et al., 2003;Borsa et al., 1997;Duranton et al., 2018Duranton et al., , 2020Naciri et al., 1999;Patarnello et al., 2007). Within the Mediterranean, it was not possible to detect any substructure corresponding to that well described for various species of the Tunisian-Siculo barrier (Duranton et al., 2018;López-Márquez et al., 2021). This might be due to the inclusion of only one site  Figure S1).

| Seascape genomics provides insights into the main environmental conditions affecting the evolution of sardines
Water temperature and derived environmental products indicating changing temperature conditions were highlighted as extremely important variables driving natural selection and local adaptation of sardines in the study area. The trend in the number of days with SST above 19°C was the sole environmental variable highlighted in all data sets. This is indicative of the extreme conditions experienced by the species within the past few years (as in European sea bass, Anastasiadi et al., 2021), impacting its reproduction (Palomera et al., 2007) and spawning success potential (Gordó-Vilaseca et al., 2021). This variable is critical for the reproduction of sardines in the Western Mediterranean Sea, where spawning starts when water temperature falls below 19°C (Palomera et al., 2007) and salinity is over a threshold (~25 PSU, Lima et al., 2022). Furthermore, temperature could be a proxy of other variables, namely of coastal upwelling strength (e.g., Zelle et al., 2004), and therefore larval dispersal (Gordó-Vilaseca et al., 2021), as well as recruitment success (Garrido et al., 2017;Gordó-Vilaseca et al., 2021;Solari et al., 2010). The impact of SST on sardine abundance has been demonstrated in California as well,
debunking previous hypotheses of competitive exclusion with anchovies as the driving factor in population size for both species (Sugihara et al., 2012).
Apart from SST, nutrient availability, sea currents and temperature velocity of change were important variables for the differentiation and local adaptation of Atlantic and Mediterranean clusters.
This was evident given that those variables were the most important in all but one data set (northern vs. southern sites data sets).
These variables are expected to mostly affect early life stages, which are particularly vulnerable to environmental variability, and thus the strength of recruitment. In sardines, such environmental forcing acting upon the egg and larval stages has already been demonstrated (Garrido et al., 2017), where recruitment strength is affected by differences in the prevailing environmental conditions in the studied areas. In the main recruitment hotspots of the species in Atlantic-Iberian waters, high recruitment years have been associated with high chlorophyll-a (Chla) concentration (as an indicator of food availability) and low SST (reflecting the optimal range for larval development) during periods of sardine larvae development (Garrido et al., 2017;Gordó-Vilaseca et al., 2021). Temperature and food availability also affect the intensity of reproduction and quality of the eggs produced (Garrido & van der Lingen, 2014).
The approaches employed identify different numbers of outlier loci with few of them in common. This is frequently observed (Gagnaire et al., 2015) and

| The genomic profile of the outlier loci
Functional annotation analyses of the outlier-associated gene set revealed that the genes are linked to key biological processes associated with environmental pressures ( Figure S2). It is worth noting that the retrieved top GO categories (e.g., cellular process, metabolism, biological regulation and response to stimulus) have been reported in other studies that measure transcriptomic responses under experimentally induced heat stress (e.g., Huang et al., 2018Huang et al., , 2020Li et al., 2021). Those findings reflect a similar functional profile of the mechanisms selected during sardines' adaptation to sea warming.
Reproduction and sex determination are the main processes greatly impacted by temperature in fish (Vandeputte et al., 2020). Among others, our gene data set included the proteins Gametogenetin-binding protein 2, AT-rich interactive domaincontaining protein 4B, Myosin-7 and Thioredoxin reductase 3, all involved in the process of spermatogenesis, and Cilia-and flagella-associated protein 52 involved in sperm mobility, reflecting a possible impact of the significant environmental parameters of our analysis to sardine reproduction. Another key process that seems to be highly impacted is muscle formation which has been described as being affected by temperature (Balbuena-Pecino et al., 2019;Georgakopoulou et al., 2010)

| Sardines from the Mediterranean and the Atlantic respond differently to environmental pressures
Τwo data sets (i.e., all samples and ATL cluster data sets) displayed shared patterns of association of individual SNPs with the majority of environmental variables. Loci showing multiple environmental associations were found in the same loci cluster using hierarchical clustering, with these clusters representing proxies for functional modules (Kess et al., 2020;Lotterhos et al., 2018). Such modular genetic architectures are favoured when adapting to climate on a landscape scale (Kliebenstein, 2011) where complex spatial and temporal environments with both abiotic and biotic challenges occur at different spatial scales (Le Nagard et al., 2011), as well as when multiple traits are under a combination of directional (among populations) and stabilizing (within populations) selection (Griswold, 2006;Le Corre & Kremer, 2003;Wagner & Altenberg, 1996). Such modularity can be beneficial not only for facilitating evolvability in response to a variable environment (Edwards & Weinig, 2011), but also by allow- ing the system to function as a set of independent groups of genes, allowing a species to better sample its potential phenotypic and genotypic diversity (Leroi, 2000). Sardines undoubtedly face spatially and temporarily complex environments within and between the two areas (i.e., Mediterranean Sea and Atlantic Ocean). Whether sardines possess phenotypic and genomic characteristics that could be indicative of directional and stabilizing selection among and within populations, respectively, demands more data and testing. Current data and analyses favour a hypothesis of among-population directional selection in sardines (i.e., as indicated by the 43 SNPs, which appear monomorphic in the Atlantic cluster, and with high frequency of the fixed allele in the Mediterranean cluster). However, given that our sampling scheme includes only one site from the Atlantic, no firm conclusion can be drawn. areas of river runoffs such as GSA07 (Palomera et al., 2007).
The Atlantic cluster seems to be mostly affected by the seascape and especially plan curvature, as indicated by both ginland and RDA analyses, probably reflecting adaptations to its habitat which is different from that within the Mediterranean. An extensive body of water such as the Eastern Atlantic displays fewer physical barriers associated with coastal relief (Santos et al., 2018), and with bottom topography that impacts its oceanographic elements such as currents, gyres, upwelling and oceanic fronts (Castelao & Luo, 2018). A significant relationship between the extent of the potential sardine habitat, larval food availability (defined by depth, SST and Chla) and sardine recruitment has already been found in the Atlantic coasts off Morocco (Machu et al., 2009). It was suggested that years of low recruitment occurred when the optimal spawning period in terms of temperature and salinity did not coincide with periods of high food availability. However, further investigations over wider geographical areas in the Atlantic Ocean are necessary to draw firm conclusions on the main physical, biological or environmental drivers of the genetic variability in the Atlantic cluster.
Sardines commonly occur in highly productive areas of enhanced food availability. In the Atlantic Ocean, these major feeding areas are largely associated with upwelling offshore systems (Checkley et al., 2017). In general, the offshore occurrence of sardines in the Atlantic Ocean has been held responsible for their larger size with respect to other small pelagic fishes (i.e., anchovies), but also for their longer migrations and lower levels of population structure (Checkley et al., 2017). In contrast, marine productivity and highly productive marine areas in the Mediterranean Sea are largely driven by river discharges (e.g., Feuilloley et al., 2020). Sardines stay closer to the coastline and feed in areas directly influenced by river runoff (Costalago et al., 2015). Based on these lines of evidence, we hypothesize that the larger size of the individuals of the ATL cluster might be due to its adaptation to wander in offshore areas. An indication of the distinct selective pressures acting on the different hierarchical clustering levels is provided by almost zero overlap in the lists of loci across the different data sets (i.e., zero overlap in ginland, 33 loci in pcadapt and 14 loci in RDA across the data sets, Figure S1).

| Is fishing effort driving the adaptation of sardines?
Covariation of SNPs with variables reflecting fishing effort was observed in the cumulative impacts variable, that is cum_impact † (Ramírez et al., 2018, where climate impacts on sardine biomass and spawning are combined with fishing effort footprint (ginland analysis), as well as in fishing effort (fishing_ef, when RDA was employed). Such findings agree with multiple lines of evidence where fishing may alter fluctuations in sardine stocks, but may not be the primary source of change (Checkley et al., 2017;Saraux et al., 2019) and the fact that environmental parameters (ICES, 2015) and in particular temperature and food are essential factors for sardines in general. Unarguably, fishing effort is affecting sardine stocks and might become detrimental to their sustainability (Ramírez et al., 2018

| Is the north/south gradient solely affected by the environment or does it have a heritable basis?
In RDA, the distinction of northern from southern sites was evident in terms of the number of days with SST < 12°C (sst_tr12) and the maximum annual SST (sst_max_sl), respectively (Figure 8a,b). According to Van Beveren et al. (2014), the recent decrease in the abundance of small pelagic fishes in the Gulf of Lions was related to slower growth and the disappearance of older individuals rather than to a decrease in recruitment strength or overfishing. The authors proposed that the current decline in sardine biomass could be due to qualitative and/or quantitative modifications in planktonic production (i.e., a bottom-up control) or mass mortalities of adults due to an epidemic disease (Van Beveren et al., 2014). This is further supported by the findings of recent studies where the decrease in body condition was not related to selective processes (Saraux et al., 2019), with spatial-seasonal factors at local scales being pivotal (Lloret-Lloret et al., 2022). It has also been observed that the spatial differences in the feeding ecology of adult sardines might affect their lipid composition, and also the fatty acid reserves transferred to their progeny (Garrido et al., 2007(Garrido et al., , 2008, with direct implications for their condition. Such phenotypic differences are not always in line with observed genetic divergence, as in the case of sardines off Moroccan coasts (Atarhouch et al., 2006(Atarhouch et al., , 2007, which were also attributed to responses to local environmental variables.
Unfortunately, the quality of the food was not directly included in our analysis due to the lack of suitable variables and no conclusions on this can be drawn.
Overall, our observed responses to local environments comply with the environmental conditions prevailing in those areas and their deterioration . According to the marine regionalization of Zhao and Costello (2019), there are extreme changes in PAR, waves and wind in the northern sites (i.e., GSA07a, GSA07b and GSA06a) during winter (December-February) but high PAR during summer (June-August) at the southern sites with consequent effects on productivity and temperature, especially during the reproductive period. The above, coupled with rapid changes in temperature (Ramírez et al., 2018, could be held responsible for such differences in physical conditions between the north and the south that might be reinforced by periods of lower connectivity with areas subject to milder pressures of extreme environmental variability and where healthier fish are found, such as the GSA06c sampling site (as in Hidalgo et al., 2019). These intensively fluctuating conditions constrain species either to tolerate such variation, adjust their activity, growth and reproduction in response to seasons, or evolve into distinct biogeographical assemblages (Zhao & Costello, 2019).

| Conclusions
Insights into the main environmental parameters that drive local adaptation and adaptive responses to changes were gained under a seascape genomics approach separating the signals of neutral and putative adaptive variation. The results obtained are not only relevant for the biology of the species, but also for understanding natural, as well as human-induced, evolutionary processes. They can assist planning of an adapted fisheries management in this area by predicting future trajectories of biodiversity and setting conservation priorities. Despite admixture between the two groups (i.e., MED and ATL), their discrete management will maintain locally adapted groups, enhancing their sustainability. Our findings provide evidence for a dynamic equilibrium in which population structure is maintained by physical and/or biological (including predation and behaviour) factors under the opposing influences of migration and selection. Drastic changes may occur in such dynamic systems that require continuous monitoring under a seascape genomics approach. Including a temporal as well as a more detailed spatial dimension will be beneficial. This study can offer valuable information for sustainable species management, providing important clues on the definition and delineation of stocks and their potential to adapt.
Furthermore, it may assist the prediction of the evolutionary responses to ongoing and future climate change scenarios. Finally, our results may further contribute to studies aimed at providing deeper insights into the mechanistic processes underlying population structuring of sardines; this is a key step for understanding and predicting future changes and responses by this highly exploited species in the face of climate change. and C.S.T. conceived and supervised the study and wrote the paper.

AUTH O R CO NTR I B UTI O N S
All authors approved the final version of the paper.

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw reads have been deposited to European Nucleotide Archive (ENA) under the accession id PRJEB45992. A VCF file with all filtered and unlinked SNPs used in the analyses of this study is provided as an additional file.