Seascape Genetics of a Globally Distributed, Highly Mobile Marine Mammal: The Short-Beaked Common Dolphin (Genus Delphinus)

Identifying which factors shape the distribution of intraspecific genetic diversity is central in evolutionary and conservation biology. In the marine realm, the absence of obvious barriers to dispersal can make this task more difficult. Nevertheless, recent studies have provided valuable insights into which factors may be shaping genetic structure in the world's oceans. These studies were, however, generally conducted on marine organisms with larval dispersal. Here, using a seascape genetics approach, we show that marine productivity and sea surface temperature are correlated with genetic structure in a highly mobile, widely distributed marine mammal species, the short-beaked common dolphin. Isolation by distance also appears to influence population divergence over larger geographical scales (i.e. across different ocean basins). We suggest that the relationship between environmental variables and population structure may be caused by prey behaviour, which is believed to determine common dolphins' movement patterns and preferred associations with certain oceanographic conditions. Our study highlights the role of oceanography in shaping genetic structure of a highly mobile and widely distributed top marine predator. Thus, seascape genetic studies can potentially track the biological effects of ongoing climate-change at oceanographic interfaces and also inform marine reserve design in relation to the distribution and genetic connectivity of charismatic and ecologically important megafauna.


Introduction
Identifying environmental conditions underlying the division of species into smaller units is central for understanding ecological and evolutionary processes and for the conservation management of biodiversity. In highly mobile species that are distributed across continuous environments with few barriers to dispersal, it is expected that persistent gene flow will stifle genetic differentiation and speciation. Nevertheless, there is growing recognition that gene flow can be limited even in the absence of geographical barriers, both in terrestrial and aquatic environments [1,2]. A detailed knowledge of how landscape characteristics structure populations has therefore become an important focus of molecular ecological research [3], leading to the emerging field of landscape genetics [3,4]. This multidisciplinary approach aims to complement genetic data with lines of evidence from other areas such as spatial statistics and landscape ecology in order to understand the effects of the landscape on the spatial distribution of genetic diversity [3,5,6]. Although extensively applied in terrestrial systems, this approach has been used less frequently in the marine environment [4]; but see [7,8].
The study of connectivity in marine systems can be challenging due to the absence of obvious barriers to dispersal and generally large population sizes of marine organisms that often resist genetic divergence, leading to low statistical power to detect population structure [8,9]. Therefore, the use of an integrative approach such as the one used in landscape genetics (or 'seascape genetics' when applied to the marine environment) has provided valuable insights into which factors may be shaping genetic structure in the world's oceans [7,10]. Biogeographic barriers and environmental variables such as ocean currents, upwelling, variation in sea surface temperature and salinity are some of the factors that have been proposed to explain genetic diversity and structure in marine organisms [9,10,11]. However, most of these studies have been conducted in organisms with larval dispersal. In active marine dispersers such as sharks and dolphins, where dispersal potential is dependent upon individual vagility, the interplay of environmental features and genetic structure has remained largely untested (but see [12]). Although differences in salinity, temperature and productivity levels have been suggested to explain genetic discontinuities in dolphins [13,14,15,16], a direct relationship between such oceanographic features and genetic structure has only been recently evaluated for two coastal dolphin species with limited distribution: the franciscana (Pontoporia blainvillei) [12] and the humpback dolphin (Sousa chinensis) [17]. These authors found that heterogeneity in chlorophyll concentration, water turbidity and temperature likely influenced the occurrence of genetically distinct populations of these species along the coast of Argentina and in the Western Indian Ocean, respectively.
In this study we use as model a highly mobile, widely distributed cetacean species belonging to the genus Delphinus, the short-beaked common dolphin. Common dolphins occur in all oceans from tropical to temperate waters. Two species and four subspecies are currently recognized: the short-beaked common dolphin, Delphinus delphis Linnaeus, 1758, distributed in continental shelf and pelagic waters of the Atlantic and Pacific Oceans; the long-beaked common dolphin, Delphinus capensis Gray, 1828, distributed in nearshore tropical and temperate waters of the Pacific and southern Atlantic waters; D. d. ponticus Barabash, 1935, restricted to the Black sea; and D. c. tropicalis van Bree, 1971, restricted to the Indian Ocean [18]. However, due to discordance between morphological and genetic characters, the phylogenetic relationships and taxonomy within the genus, particularly in regard to the specific status of the long-beaked form, are still under debate (Amaral et al. unpublished data;[19]).
Short-beaked common dolphins are known to occur in large groups of dozens to hundreds of individuals. Although their social structure is still poorly understood, individuals seem to group irrespective of genetic relationships, with possible gender and age segregation [20]. However, there is a gap in knowledge if these findings are representative for common dolphins in other geographic regions. The movements of common dolphins are thought to be largely determined by those of their potential prey (e.g. [21]) and their diet varies between locations and seasons [21,22]. Nonetheless, they generally depend on small, mesopelagic shoaling fishes such as scombroids and clupeoids, and squids [21,22]. It has been suggested that short-beaked common dolphins often prefer specific water masses [15,23,24] and in the Eastern Tropical Pacific they occur preferentially in upwelling-modified waters [23].
Genetic studies conducted so far have shown significant genetic differentiation among populations inhabiting different oceans and different coasts of the Atlantic Ocean [19,25]. However, within each side of the Atlantic Ocean, no genetic structure has been detected, suggesting a lack of strong dispersal barriers in these areas [25,26]. Within the Pacific Ocean, results from regional   [15,27]). However, a direct evaluation of the influence of oceanographic variables on the genetic structure of this species has never been carried out. Our aim is to assess the relative influence of key oceanographic variables on population subdivision of short-beaked common dolphins at a range of medium to large spatial scales, including within ocean basins and across oceans. To achieve this aim we have sampled populations inhabiting the Atlantic, Pacific and Indian Oceans and used remote sensing data under a seascape genetics approach. The global distribution, high mobility, and putatively close association of short-beaked common dolphins with water masses, makes them an excellent model species to test for interactions between variation in environmental factors and genetic structure, contributing towards an understanding of ecological processes affecting population connectivity in the sea.

Ethics Statement
This study was conducted according to relevant national and international guidelines. No ethics approval was considered necessary because the animals were not handled directly. Permissions for collecting samples were obtained separately in countries where it was required (Macquarie University Animal Ethics Committee, Australia; Southwest Fisheries Science Center Ethics Advisory Committee, USA; Institute for Nature Conservation and Biodiversity, Portugal; and Department of Conservation, New Zealand). CITES permits numbers used to export/import samples were: 07US168545/9,  (Table 1). All tissue samples were obtained from either stranded animals (103 samples) or from skin biopsies (178 samples) collected from free-ranging dolphins. Tissues were stored either in ethanol or in 20% DMSO/saturated NaCl.

DNA extraction and microsatellite genotyping
Genomic DNA was isolated from skin or muscle using a standard proteinase K digestion and two phenol-chloroform and one chlorofom-isoamyl extractions followed by ethanol precipitation [28] for samples originated from stranded animals or, alternatively, using a salting-out protocol [29] for samples originated from biopsies. DNA quality and concentration was verified using Thermo Scientifc NanoDrop 1000 Spectrophotometer (Thermo Fisher Scientific Inc.). Samples from NEPAC and NWATL were provided as DNA by the Southwest Fisheries Science Center, Marine Mammal and Turtle Research Sample Collection (SWFSC-NOAA, La Jolla, CA).
All samples were genotyped at 14 polymorphic microsatellite loci: 7 tetranucleotide (Tur4_80, Tur4_87, Tur4_92, Tur4_105, Tur4_141, Tur4_142; [30] and Dde59 [31] and 7 dinucleotide (Dde66, Dde70; [31]), KW2, KW12 [32], EV1 [33], MK6 and MK8 [34]. The forward primer for each primer pair was labelled with a M13 tag [35]. Fluorescent dyes were also labelled with the M13 tag. Amplification reactions contained 50-100 ng DNA, 16 GoTaqH reaction buffer (Promega), 2.5 mM MgCl 2 , 0.2 mM dNTPs, 0.1 mM of each primer and 1 U GoTaqH Taq DNA polymerase (Promega). The thermal cycler profile for the tetranucleotide loci and Dde66 and Dde70 consisted of initial denaturation at 94uC for 3 min followed by a touchdown profile for 5 cycles with the annealing temperature starting at 63uC and decreasing 2uC per cycle, followed by 30 cycles with an annealing temperature of 53uC, and a final extension step at 72uC for 10 min. The tetranucleotide loci were amplified in multiplex after optimization. For the remaining dinucleotide loci, conditions followed the original publications. All reactions included both positive and negative controls. Following amplification, samples were mixed with an internal size standard (LIZ 500) and run on an ABI 3130 Genetic Analyzer. The GeneMapper v.4.1 software (Applied Biosystems, CA) was used for sizing of allele fragments.

Data analysis
Genetic diversity. The program Micro-checker v.2.2.3 [36] was used to check for the presence of genotyping errors such as scoring errors due to stuttering, large allele dropout or evidence for null alleles. Departures from Hardy-Weinberg Equilibrium were tested for each population using the Fisher exact test in Genepop v.4.0 [37]. Genepop was also used to test for linkage disequilibrium between loci. Samples were grouped into 7 putative populations according to their geographical origin as described above. Genetic diversity measures such as mean number of alleles per locus and observed (H O ) and expected (H E ) heterozygosities were calculated in Arlequin v.3.5.1 [38] and allelic richness (A R ) calculated using FSTAT v.2.9.3 [39].
Genetic differentiation. Three different measures of population differentiation were used: the fixation index F ST , estimated using FSTAT [39]; the analogous R ST , estimated using Genepop v.4.0 [37]; and the statistic Jost's D [40], estimated using SMOGD v.1.2.5 [41]. The latter has been shown to provide a more accurate measure of differentiation when using highly polymorphic microsatellite loci [40]. Additionally, we tested for a mutation effect on genetic structure by randomly reassigning allele sizes while keeping allele identity the same [42]. The test was conducted in SPAGEDI v.1.3 through 10,000 permutations. R ST values significantly larger then F ST values indicate that mutation, in addition to drift and gene flow, has contributed to frequency differences among samples, which in some cases can be interpreted as phylogeographic signal [42].
In order to visualize relationships among putative populations based on genetic variation, we performed a principal component analysis (PCA) on a table of standardised allele frequencies using the adegenet and ade4 packages in R [43]. In addition, we performed an analysis of nonmetric multidimensional scaling (MDS, [44]) on each of the genetic distance matrices using the PRIMER computer package [45].
An analysis of molecular variance, AMOVA [46] was conducted in Arlequin to assess population structure. Different hierarchical levels were tested, considering differences occurring between populations in different oceans and within the same ocean basin.
A Bayesian approach to identify the number of populations (K) present in the dataset was implemented in the program STRUCTURE v.2.3.3 [47,48]. The admixture and the correlated allele frequencies models were implemented since we expect that allele frequencies in the different populations are likely to be similar due to migration or shared ancestry. Sampling locations were used as prior to help detect population structure [49]. Ten independent runs of K between 1 and 8 were run with 400 000 ''burn in'' and 4 million MCMC replicates. The maximum loglikelihood values from all runs corresponding to each given K were checked for consistency and averaged. The K with the highest averaged maximum log-likelihood was considered the most likely number of clusters that better explains our dataset. CLUMMP v.1.1.2 [50] was used to summarize parameters across 10 runs and distruct v.1.1 [51] was used to produce the corresponding graphical output.
Isolation by distance. Isolation by distance (IBD) was evaluated using a Mantel test implemented in the program IBDWS v.3.16 [52]. Genetic distance matrices given by F ST / (12F ST ) were regressed against the logarithm of geographical distances following a two-dimensional model [53]. R ST and Jost's D values were also used. Geographic distances were measured in Google Earth by using set points and measuring either straight-line distance across oceans, or the shortest geographical distance along continental margins. The set points were chosen so as to represent the middle point of the area of distribution where the samples were collected.
Environmental predictors of genetic structure. Three different oceanographic variables were used as predictors of the observed genetic differences between short-beaked common dolphin populations. These were night-time sea surface temperature (SST, uC), chlorophyll concentration (CHL, mg/m 3 ) and water turbidity measured as diffuse attenuation coefficient at 490 nm (KD490, m 21 ). These variables, here obtained from remote sensing data, have been previously related to habitat heterogeneity [54] and associated with genetic differences in other dolphin species [17]. Furthermore, the oceanographic variables chosen have a wide geographic coverage through remote sensing, making them ideal for a global approach. Seven oceanic regions, corresponding to the sampling areas for short-beaked common dolphins, were used for the extraction of these oceanographic variables to assess association with patterns of genetic differentiation. Polygons were defined considering the possible range of common dolphins within that oceanic region, with the last side being the coastline. For NWATL the area was defined between 46uN, 38uN and 57uW; for CEATL between 34uN, 32uN and 16uW; for NEATL between 60uN, 35uN and 0u; for NEPAC between 45uN, 25uN and 108uW; for SWPAC_NZ between 32uS, 44uS and 180uW; for SWPAC_AUS between 26uS, 44uS and 156uE; and for SEIND between 31uS, 37uS and 140uE. In order to account for possible influence of area choice in the final results, areas restricted to where samples from free-ranging animals originally came from or from published distributional data were considered and re-analysed. Since no differences were found in the final results, only analyses including the areas defined above are presented, which account for a possible wider ranging distribution of common dolphins. Monthly averaged data of the three variables, with a 4 km spatial resolution was obtained from Ocean Color Web (http://oceancolor.gsfc.nasa.gov/) for the period from July 2002 to October 2010 and processed using MATLAB software (www. mathworks.com). Data collected during this time period provide a characterization of the oceanographic features for each region and are robust to inter-annual oscillations (Supplementary Material, Figure S1). Data analysis included the construction of temperature, chlorophyll and turbidity maps for each region, where each pixel of the map corresponds to the eight-year average value for a 4 km grid. These maps were visually inspected to detect geographical areas of environmental heterogeneity. Monthly averages for each oceanic region were then statistically analysed using a paired t-test to detect differences among those regions. Total averages for the 8 year-period for each factor and each sampled region were subsequently used to examine environmental and genetic associations (details below). Environmental distances were calculated as pairwise differences in mean temperature, chlorophyll and turbidity between regions. Pairwise F ST , R ST and Jost's D were used as genetic distances.
All analyses were carried out at different spatial scales: at a large scale, all oceans included; each ocean considered in separate, i.e. all populations within the Atlantic and all populations within the Pacific Ocean and the population in the Southeast Indian Ocean; and at a medium scale, the North and Central Atlantic populations (hereinafter referred to as North Atlantic) and the South Pacific and Southeast Indian Ocean populations (hereinafter referred to as South Indo-Pacific).
Seascape genetics. Associations between genetic and environmental factors were examined using a hierarchical Bayesian method implemented in GESTE [55], which estimates individual F ST values for each local population and then relates them to environmental factors via a generalized linear model. Here we used 10 pilot runs of 1,000 iterations to obtain the parameters of the proposal distribution used by the MCMC, and an additional burn-in of 5610 6 iterations with a thinning interval of 20. The model with the highest posterior probability is the one that best explains the data [55].
Additionally, we used the BIOENV procedure of [56] as implemented in PRIMER v.5 [45] and as described in [57] to examine which predictor variable would provide the best model to explain the population genetic structure observed in the data. This procedure calculates the value of Spearman's rank correlation coefficient (r) between a genetic distance matrix (response matrix) with a distance matrix calculated as the Euclidean distance among one or more predictor variables. It then calculates the value of r using every possible combination of predictor variables until it finds the ''best fit'', corresponding to the combination of predictor variables whose Euclidean distance matrix yields the highest value of r [56]. We used three different response matrices corresponding to F ST , R ST and Jost's D distance matrices to identify the best one, two or three-variable fits.
Mantel tests [58] were also used to test for correlations between the pairwise genetic and environmental distances. Partial Mantel tests were used to control the effect of geographical distances in these potential correlations. These tests were performed using the package vegan in R.

Genetic Diversity
In total 281 short-beaked common dolphin samples were genotyped at 14 microsatellite loci (Table 1). Results from Micro-  Checker and the Fisher exact test suggested deviations from Hardy-Weinberg equilibrium (HWE) in 4 loci. Two of these (Tur91 and Tur80) showed deviations in only one population each and were therefore included in subsequent analyses, whereas the other two (Tur141 and Dde66) showed deviations in 4 and 2 populations, respectively. These deviations are due to a deficit of heterozygotes (significant F IS values, Table 1). To test whether results would be affected by the inclusion of these two loci, estimates of genetic variability and differentiation were carried out with and without them. Since no major differences in results were observed (data not shown), all 14 loci were used in subsequent analyses. These deviations are likely not related with the fact that some samples originated from strandings and others from biopsies. In fact, it has been recently shown that no apparent differences occur when testing population structure in common dolphins using samples originated from carcasses or from free-ranging dolphins [59].
Levels of genetic diversity, given by mean number of alleles, allelic richness and expected and observed heterozygosities were high for most populations (Table 1). Significant F IS values were obtained for populations from NE Pacific and SW Pacific Australia and New Zealand, which can be due to the presence of population sub-structure (i.e. Wahlund effect). In fact, this is known to be the case for common dolphins inhabiting those regions ( [15,27]; Stockin et al. unpublished).

Genetic differentiation
Pairwise F ST and R ST comparisons showed significant levels of differentiation among all putative populations (Table 2), although the extent of that differentiation differed for each index. Jost's D values tended to be higher than F ST and R ST values. R ST also tended to be higher than F ST . Since R ST is based on allele size, the differences observed indicate that mutation, in addition to drift or gene flow may be affecting the differentiation between these populations. This result was confirmed using SPAGEDI. The overall R ST value was significantly higher than the overall F ST value (P = 0.042).
Taken as a whole, the fixation indices showed high levels of differentiation between short-beaked populations inhabiting different ocean basins. The SEIND and NEPAC populations showed the highest levels of differentiation when compared with all other short-beaked populations. Contrasting to the inter-ocean basin differentiation, lower levels of differentiation were observed between short-beaked populations inhabiting the same ocean basins.
The first two principal components of the PCA analysis explained 84.35% of the variance in allele frequencies among putative populations (Figure 2). The first principal component shows a clear separation between populations inhabiting the Indo-Pacific and the Atlantic Oceans. The second principal component further shows some structure within the Indo-Pacific region, with the SEIND and NEPAC populations appearing separated from the SWPAC_AUS and SWPAC_NZ populations.
Non metric MDS analyses using the three different genetic indices also show a clear separation from populations inhabiting the Atlantic, the Pacific and Indian oceans, with the exception of the analysis using R ST , which grouped the NEPAC population with Atlantic ones (Figure 3). The analyses using F ST and Jost's D show a closer proximity among the short-beaked populations inhabiting the North Atlantic, and also of the populations inhabiting the Pacific Ocean.
Results obtained in STRUCTURE using the correlated allele frequency model resulted in a peak of maximum ln P(K) at K = 3 ( Figure 4, Supplementary Table S2). These clusters correspond to populations inhabiting the three ocean basins: the Atlantic (including the NEATL, NWATL and CEATL populations), the Pacific (including the NEPAC, SWPAC_AUS and SWPAC_NZ populations) and the Indian Ocean including the SEIND population ( Figure 4).
The AMOVA analysis showed that the highest levels of differentiation were obtained when populations were divided by eastern versus western regions within ocean basins (F CT = 0.03425, P,0.0001) ( Table 3).

Isolation by distance
The relationship between geographic and genetic distance was only observed when populations inhabiting all oceans were considered in the analysis and when F ST and Jost's D values were used (Table 4). This relationship was not detected when R ST values were used, nor when finer spatial scales were considered.  Data on sea surface temperature (SST), chlorophyll concentration (CHL) and water turbidity (KD490) was gathered for the seven oceanic regions where short-beaked common dolphins were sampled: NEATL, CEATL, NWATL, NEPAC, SWPAC_AUS, SWPAC_NZ and SEIND ( Figure 5). Paired t-tests showed significant differences in the 8 year average values of SST between most regions with exception of the comparison between NEATL and NWATL, between NEPAC and SWPAC (both AUS and NZ), and between NEPAC and SEIND, where differences were not statistically significant (P,0.01, see Supplementary Material, Table S1). In the SST maps, all regions are heterogeneous, having regions of colder and warmer waters ( Figure 5). Nevertheless, NEATL and NWATL regions are dominated by colder waters when compared with other regions, which are dominated by warmer waters, such as SWPAC_AUS and SWPAC_NZ. Significant differences were not detected in mean CHL values between NEPAC and SWPAC (both AUS and NZ) and between NEPAC and SEIND, as well as among SEIND, SWPAC_AUS and SWPAC_NZ. All other comparisons were significant. Despite this, in the CHL maps, clear differences can be seen among the regions located in the Pacific Ocean. Chlorophyll concentrations are higher in the NEPAC region closer to the coast when compared to the SWPAC_AUS and SWPAC_NZ regions. Regarding turbidity mean values, these were only not significant in the comparisons among SWPA-C_AUS, SWPAC_NZ and SEIND (Table S1). Patterns seen in the maps are similar to the ones obtained for the CHL maps ( Figure 5).

Seascape genetics
Hierarchical Bayesian analyses implemented in GESTE identified the model including the constant as the best one in all spatial scales considered ( Table 5). The second best model for all analyses was the one including KD490, though the third and fourth models (including CHL and SST) all had very similar posterior probability values. Higher posterior probabilities were obtained when medium spatial scales were analysed. Positive signals of the regression coefficients were obtained for the association between CHL and genetic differentiation in the Pacific Ocean and South Indo-Pacific Ocean populations, and for the association between KD490 and genetic differentiation in the Pacific Ocean populations (Table 5). Regarding SST, positive signals of the regression coefficients were obtained for all populations across all oceans, for the North Atlantic populations, and for the South Indo-Pacific populations (Table 5). Therefore, genetic isolation of populations within the Pacific Ocean increases with differences in CHL and KD490 among regions, whereas genetic isolation of populations within the Atlantic Ocean increases with differences in SST among regions. In the South Indo-Pacific region, both CHL and SST increase genetic isolation among populations. The percentage of variation that remained to be explained (indicated by sigma values) was however moderate ( Table 5).
The BIOENV procedure found strong positive correlations between oceanographic predictors and genetic differentiation for the analyses conducted at medium spatial scales (Table 6). For the populations within the Atlantic Ocean and within the South Indo-Pacific, CHL and KD490 showed stronger correlation with genetic distance. For the larger spatial scales considered (across all oceans and within the Pacific Ocean), a strong negative correlation between CHL and KD490 with rank genetic distance was found ( Table 6).
Mantel tests and Partial Mantel tests between genetic and environmental distances were not statistically significant for any comparison, even considering different spatial scales (results not shown). Failures of these tests to detect relationships between genetic and environmental data have been previously described [60,61] and could explain the unsuccessful use with our datasets.

Discussion
We used a seascape approach to investigate the interaction between a set of oceanographic variables and population structure in a highly mobile, widely distributed top marine predator, the short-beaked common dolphin. We show that sea surface temperature, chlorophyll concentration and water turbidity seem to be important factors in explaining the observed patterns of genetic structure in these dolphins, more than geographical distance alone, particularly when medium spatial scales were considered.

Genetic structure
The overall global pattern of genetic structure obtained here supports previous studies [19]: higher levels of differentiation were obtained across large geographical scales, between different ocean basins, and lower levels were obtained when medium geographical scales were considered, within the same ocean basin. While results from STRUCTURE showed a clear differentiation between ocean basins, the AMOVA analysis resulted in higher F CT estimates for partitioning of short-beaked populations among regions within each ocean basin. The low levels of divergence found between populations inhabiting the same ocean basin may have affected the power of the program STRUCTURE to detect such differentiation, even using recently developed algorithms that account for weak differentiation [49]. Nonetheless, the PCA and the NMDS plots also indicate some level of differentiation within ocean basins, which seems to be stronger among the Pacific Ocean populations. Multivariate analysis does not require strong assumptions about the underlying genetic model, such as Hardy-Weinberg equilibrium or the absence of linkage disequilibrium [43]. The high levels of differentiation found for the SEIND population (southern Australia) were surprising given the comparatively shorter distance separating this population from the Southwest Pacific populations (off New South Wales, southeastern Australia), even considering that the region where the SEIND population was sampled (off South Australia) falls into a different biogeographic region (see [62] to the one of the SWPAC_AUS population. Such high differentiation was also reported by [27] when comparing individuals from this region to individuals from southeastern Tasmania (Southwest Pacific) -in that case oceanographic features affecting the distribution of target prey were suggested to be the likely explanation for the genetic differentiation found. Our study corroborates this previous finding (see below).  Table 5. Posterior probabilities of the four most probable models for the GESTE analysis of environmental associations with genetic structure (population specific F ST ) of short-beaked common dolphins.

Isolation by distance
A pattern of isolation by distance was only observed when large spatial scales were considered, indicating that the stronger genetic differentiation observed in short-beaked common dolphins from different oceans may be an effect of geographic distance. Isolation by distance has been reported for other cetacean species, such as in the harbour porpoise [63] and in bottlenose dolphins [64]. Conversely, when medium geographic scales were considered (i.e. within each ocean basin), no isolation by distance effect was detected, and genetic differentiation could be explained by oceanographic variables. This pattern has also been described for common dolphins at small geographical scales, along the eastern Australian coast [15], for bottlenose dolphins in South Australia where a temperature and salinity front coincides with the boundary between two distinct genetic populations [13], and for pilot whales, where ecological factors, such as SST, were more important in explaining genetic structure than geographic separation [14]. In franciscana and humpback dolphins, environmental factors were also more important in explaining genetic structure than distance at small geographical scales [12,17].  Oceanographic predictors All oceanographic variables tested, CHL, KD490 and SST, showed an association with population genetic structure in shortbeaked common dolphins. These associations were strongest at the medium spatial scales considered. In the Pacific Ocean, CHL and KD490 were the environmental predictors that were most strongly associated with increased genetic isolation in short-beaked common dolphins. Conversely, in the Atlantic Ocean, SST was the strongest predictor associated with population divergence. Although no significant statistical differences in the 8-year average values of CHL and KD490 were detected among regions in the Pacific Ocean, a visual inspection of the regional maps shows heterogeneity in these variables among regions ( Figure 5). Heterogeneity in SST, CHL and KD490 is also seen among Atlantic Ocean regions, although our results suggest that only SST seems to explain genetic differentiation of short-beaked common dolphins in this area. Marine productivity and SST are important variables for habitat occupancy and dispersal in cetaceans [65,66] and have been shown to influence population structure in Franciscana [12] and in humpback dolphins [17]. Here, we suggest that they are also important drivers of population structure in common dolphins. A direct causality is however difficult to establish. For example, it has been suggested that ecological factors such as prey behaviour rather than inherent sensitivity to environmental factors, could account for the relationship between SST and population structure in pilot whales [14,66,67]. Similarly, differences in prey distribution and abundance between regions rather than SST differences themselves are suggested to account for genetic differentiation of bottlenose dolphins in South Australia [13] and short-beaked common dolphins in southern [27] and southeastern Australia [15]. We suggest that a similar process may account for the patterns obtained in this study. Since dolphins feed high in the food chain, a statistical association with oceanographic variables that do not directly affect the individuals, but rather affect their prey, is expected to be weak [23]. This could also explain the fact that analyses performed in GESTE did not result in a single best-chosen model and that the percentage of variability that remained to be explained in the data was moderate.
Chlorophyll concentration, water turbidity and SST are routinely used to map ocean primary productivity (e.g. [68]). Due to the bottom-up processes that control marine ecosystems [69], these variables have been related to prey distribution and abundance, and to the occurrence of top marine predators (e.g. [70,71]). Distribution and abundance of prey has been suggested as the main factor dictating seasonal migrations in several species of delphinids, including short-beaked common dolphin (e.g. [21]). Moreover, short-beaked common dolphins feed primarily on small mesopelagic schooling fish such as sardines and anchovies [21,22]. These fishes are filter feeders and occur in association with nutrient rich waters (e.g. [72]), and could explain the dolphins' preference for certain oceanographic conditions. We further suggest that a behavioural mechanism such as specialization for local resources could also explain the patterns observed. Resource specialization is a common mechanism driving population structure in delphinds [73]. Moreover, dietary segregation is known to occur in short-beaked common dolphins. In the Bay of Biscay, Northeast Atlantic Ocean, common dolphins inhabiting neritic and oceanic waters feed on different prey species [74]. Feeding specialization leading to local adaptation has also been suggested as driving speciation of the short and long-beak forms [19] and as important triggers for the process of population divergence and speciation in the genera Tursiops and Stenella [75,76]. Perhaps the best studied example within delphinids are killer whales (Orcinus orca), where resource partitioning and foraging specializa-tions of sympatric populations occurring in the North Pacific have lead to the evolution of distinct lineages [77]. Short-beaked common dolphins could therefore be locally adapted to the existent prey species and only move within certain regions following prey migration. Seasonal migrations are known to occur in the Northeast Pacific [78] and Southwest Indian Ocean [79]. Further investigation is however required to support this hypothesis.
There are also other factors that may account for population divergence in common dolphins that were not assessed in this study. Fine-scale oceanic processes, for example, have recently been suggested to affect connectivity in common dolphins [15]. A proper assessment of its direct relationship with genetic structure requires knowledge on hydrodynamic modelling and will certainly be the aim of forthcoming studies. Demographic and historical processes can also contribute to population structure and should also be integrated in future analyses.

Implications for conservation and management
The results presented here are of particular importance for marine conservation management and design of marine protected areas (MPA). MPAs are usually designed to protect coastal regions that are either important habitats, as part of the marine ecosystem, or biodiversity hotspots [80]. Marine predators are often used as indicators for MPA design, because their protection aids in protecting the more complex environments they use [81,82,83]. Although several studies have described the distribution and occurrence of cetacean species in relation to different habitat variables (e.g. [84,85,86]), only a few have found a direct correlation between oceanographic variables and population structure [12,17]. In this study, by showing how marine productivity correlate with population structure in short-beaked common dolphins, we highlight the importance of using seascape genetic studies to inform MPA design in relation to distribution and genetic connectivity of charismatic and ecologically important megafauna. Furthermore, we highlight how such an approach can track the biological effects of ongoing climate-change and prevent the loss of top marine predators [87]. Table S1 Mean pairwise difference between average values of a) sea surface temperature (SST), b) chlorophyll concentration (CHL) and c) water turbidity (KD490) obtained for each oceanographic region where short-beaked common dolphins were sampled for this study, with significant values of paired t-tests indicated in bold. (XLS)

Table S2
Individual runs for the Bayesian analysis implemented in the program STRUCTURE with a burn-in phase of 4610 5 and 4610 6 MCMC replicates. The log-likelihood of the data (LnP(D)) for each run and an average across 10 runs for each K are shown. The K with the highest averaged maximum log-likelihood was considered the most likely number of clusters that better explains our dataset (in bold). (XLS)