# Ecological and biogeographic features shaped the complex evolutionary history of an iconic apex predator (Galeocerdo cuvier)

Pierre Lesturgie, Hugo Lainé, Arnaud Suwalski, Pascaline Chifflet-Belle, Pierpaolo Maisano Delser, Eric Clua, Sébastien Jaquemet, Hélène Magalon, Stefano Mona

► To cite this version:

Pierre Lesturgie, Hugo Lainé, Arnaud Suwalski, Pascaline Chifflet-Belle, Pierpaolo Maisano Delser, et al.. Ecological and biogeographic features shaped the complex evolutionary history of an iconic apex predator (Galeocerdo cuvier). BMC Ecology and Evolution, 2022, 22 (1), pp.147. 10.1186/s12862-022-02100-y. http://hal-03923924

HAL Id: hal-03923924

https://hal.science/hal-03923924

Submitted on 23 Feb 2023

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>
RESEARCH

Open Access

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

Check for
updates

# Ecological and biogeographic features shaped the complex evolutionary history of an iconic apex predator (*Galeocerdo cuvier*)

Pierre Lesturgie<fnref n="1" />, Hugo Lainé<fnref n="1" />, Arnaud Suwalski<fnref n="1" />,², Pascaline Chifflet-Belle<fnref n="1" />,², Pierpaolo Maisano Delser³, Eric Clua²,⁴, Sébastien Jaquemet⁵, Hélène Magalon⁴,⁵† and Stefano Mona<fnref n="1" />,²,⁴*†

**Abstract**

**Background:** The tiger shark (*Galeocerdo cuvier*) is a large iconic marine predator inhabiting worldwide tropical and subtropical waters. So far, only mitochondrial markers and microsatellites studies have investigated its worldwide historical demography with inconclusive outcomes. Here, we assessed for the first time the genomic variability of tiger shark based on RAD-seq data for 50 individuals from five sampling sites in the Indo-Pacific (IP) and one in the Atlantic Ocean (AO) to decipher the extent of the species' global connectivity and its demographic history.

**Results:** Clustering algorithms (PCA and NMF), *F*<sub>ST</sub> and an approximate Bayesian computation framework revealed the presence of two clusters corresponding to the two oceanic basins. By modelling the two-dimensional site frequency spectrum, we tested alternative isolation/migration scenarios between these two identified populations. We found the highest support for a divergence time between the two ocean basins of ~ 193,000 years before present (B.P) and an ongoing but limited asymmetric migration ~ 176 times larger from the IP to the AO (*Nm* ~ 3.9) than vice versa (*Nm* ~ 0.02).

**Conclusions:** The two oceanic regions are isolated by a strong barrier to dispersal more permeable from the IP to the AO through the Agulhas leakage. We finally emphasized contrasting recent demographic histories for the two regions, with the IP characterized by a recent bottleneck around 2000 years B.P. and the AO by an expansion starting 6000 years B.P. The large differentiation between the two oceanic regions and the absence of population structure within each ocean basin highlight the need for two large management units and call for future conservation programs at the oceanic rather than local scale, particularly in the Indo-Pacific where the population is declining.

**Keywords:** Agulhas leakage, Coalescent modelling, Demographic history, Population genomics, RAD-seq, Tiger shark

**Background**

Predation plays a fundamental role in the top-down regulation of ecosystem dynamics, with apex predators being key actors in promoting species diversity [1]. However, in marine ecosystems, many predatory species have declined across their ranges [2]. Efforts to develop conservation programs need be tailored to the appropriately scaled units of managements for the species under investigation [3]. Recent advances in DNA sequencing technologies allow the characterization of thousands of independent loci giving the power to assess the

<footnote>†Hélène Magalon and Stefano Mona jointly supervised this work</footnote>

<footnote>*Correspondence: stefano.mona@mnhn.fr</footnote>

<footnote n="1">¹ Institut de Systématique, Evolution, Biodiversité (ISYEB), EPHE-PSL, Université PSL, MNHN, CNRS, SU, UA, Paris, France
Full list of author information is available at the end of the article</footnote>
genetic diversity of any target model or non-model species, which can inform management policies. However, genetic diversity assessment should be complemented by the reconstruction of species connectivity and historical demography to better establish conservation priorities. For instance, understanding how populations are spatially connected as well as the divergence time between lineages is essential to decipher at which geographic scale a species should be managed. Reconstructing the evolutionary history of a species is often a complex task that requires an educated choice of the most likely model of evolution, often selected among a reduced selection of biologically meaningful models. Unfortunately, selection of an inappropriate model can yield misleading estimates of critically important parameters as more data are collected. This has important implications for conservation genetic applications that rely on accurate estimates of genetic diversity and changes in effective population size through time [4–7].

The tiger shark (*Galeocerdo cuvier*, Péron & Lesueur, 1822) is a large and iconic apex marine predator, that is considered “Near Threatened” by the International Union for Conservation of Nature (IUCN). Though the tiger shark is a predominantly coastal species, its distribution includes tropical and subtropical waters worldwide [8]. The species is heavily impacted by fisheries [9] and shark control programs in the Indo-Pacific [10]. Indirect estimates have suggested an annual number of tiger shark catches between 50,000 and 300,000 individuals [11], raising conservation concerns. Though not directly endangered by global climate change, the species is likely to extend its habitat range poleward as a consequence of the rising in annual sea surface temperatures [12], which may increase the potential for greater trans-oceanic movements [13]. Despite being found predominantly along the coast, tiger sharks spend considerable time in pelagic waters and telemetry studies have shown that they can cross oceanic expanses [14–16], but no evidence of contemporary migration between the Indo-Pacific and the Atlantic Ocean has yet been found.

Knowledge of these ecological traits is important to devise meaningful evolutionary models, but it is not sufficient. Large marine predators with continuous distributions can be intuitively considered as capable of high dispersal due to the absence of clear physical barriers [17]. Nevertheless, there are both examples of panmictic species, such as the blue shark *Prionace glauca* [18] or the mako shark *Isurus oxyrinchus* [19], and examples of species structured according to ocean basins such as the Galapagos shark *Carcharhinus galapagensis* and the dusky shark *Carcharhinus obscurus* [20]. In the tiger shark, there has been contrasting evidences about the degree of population structure, the

extent of genetic diversity and particularly about the historical demography [21–27]. All studies recognize the existence of a clear separation between the Indo-Pacific (IP) and the Atlantic Ocean (AO), with genome-wide data supporting low to no population structure within each basin [27]. However, it remains unclear whether the two basins hold two allopatric species as originally proposed by Naylor et. al (2012) or two divergent lineages as proposed by Bernard et. al (2016). An accurate characterization of divergence and migration is additionally still lacking: Bernard et al. (2016) provided a divergence time computed on mtDNA using a non-equilibrium model implemented in *MDIV* [28] but no confidence interval could be determined. At the same time, migration rates were estimated using the equilibrium model implemented in *MIGRATE* [29]. Yet a global analysis estimating simultaneously all parameters is warranted. Furthermore, an in-depth analysis of the historical demography in both regions is still lacking, with only [24] supporting a recent decrease in effective population size in two sampling sites from the IP. Using the wealth of data provided by RAD-seq, we sequenced a total of 50 sharks from six sites in the IP and one in the AO (Fig. 1) in order to shed light on the complex evolutionary history of the tiger shark. We first investigate the extent of genetic diversity, the level of population structure and historical demography in all sampling sites, and finally tested alternative evolutionary scenarios to model the divergence and migration between IP and AO by fitting the observed two-dimensional site frequency spectrum (2D-SFS) with coalescent simulations. These analyses are necessary not only to reconstruct the evolutionary history of the tiger shark but also to better inform conservation strategies.

## Results

### RAD-seq sequencing

The average number of reads retained per individual after the quality filtering and demultiplexing step was 4,011,430 (±1,314,894 s.d.). After a first round of de novo assembly and filtering using STACKS v2.5, the depth of coverage was low with a mean of 12.65 (±6.36 s.d.), which motivated the use of the genotype-free allele frequency estimation pipeline implemented in *ANGSD* [30] rather than the direct call. The final number of loci (variable and fixed) was highly variable between sampling locations (Table 1) ranging from 16,953 to 118,591 in Brazil (BRA) and New Caledonia (NCA) sampling sites respectively, with a number of SNPs with no missing data following a similar trend (from 5868 to 26,075 for BRA and NCA, respectively).
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Map of the sampling sites. From west to east: Brazil (BRA: n = 7), Reunion Island (RUN; n = 15), North Coast of Australia (AUS<sub>Ni</sub>; n = 8), Coral Sea (COR; n = 5), East Coast of Australia (AUS<sub>E</sub>; n = 7) and New Caledonia (NCA; n = 8)</figcaption></figure>

<table><caption>**Table 1** Sample size (n), mean pairwise difference (θ<sub>π</sub>), Watterson theta (θ<sub>w</sub>), Tajima's D (TD), and total number of loci (monomorphic included) (n<sub>loci</sub>) and SNPs (n<sub>SNP</sub>) without missing data for all sampling sites (ranged from west to east)</caption><thead><tr><th></th><th>n</th><th>θ<sub>π</sub> (10<sup>-3</sup>)</th><th>θ<sub>w</sub> (10<sup>-3</sup>)</th><th>TD<sup>a</sup></th><th>n<sub>loci</sub></th><th>n<sub>SNP</sub></th></tr></thead><tbody><tr><td>BRA</td><td>7</td><td>0.97</td><td>1.16</td><td><strong>0.14</strong></td><td>16,953</td><td>5868</td></tr><tr><td>RUN</td><td>15</td><td>0.59</td><td>0.73</td><td><strong>-0.19</strong></td><td>71,214</td><td>19,971</td></tr><tr><td>AUS<sub>N</sub></td><td>8</td><td>0.76</td><td>0.97</td><td><strong>-0.19</strong></td><td>38,420</td><td>11,627</td></tr><tr><td>COR</td><td>5</td><td>0.58</td><td>0.72</td><td>0.04</td><td>97,736</td><td>18,407</td></tr><tr><td>AUS<sub>E</sub></td><td>7</td><td>0.63</td><td>0.77</td><td>0.02</td><td>49,380</td><td>11,385</td></tr><tr><td>NCA</td><td>8</td><td>0.57</td><td>0.7</td><td>-0.03</td><td>118,591</td><td>26,075</td></tr></tbody></table>

AUS<sub>E</sub> East Coast of Australia, AUS<sub>N</sub> North Coast of Australia, BRA Brazil, COR Coral Sea, NCA New Caledonia, RUN Reunion Island

<sup>a</sup> Tajima's D values in bold are significantly different from 0 (P < 0.001)

## Population structure

Population structure was investigated using datasets allowing up to 20% of missing data per SNP. Thus, after filtering, the remaining number of SNP was 24,454 for the Principal Component Analysis (PCA) and the non-negative matrix factorization (nmf) inference, and ranged from 8785 to 15,977 per population pair for the F<sub>ST</sub> computation. Pairwise F<sub>ST</sub> highlighted a moderate differentiation between Indo-Pacific (IP) and Atlantic Ocean (AO) sampling sites, with values ranging from 0.117 to 0.129 and systematically significant (P ≤ 0.001, Fig. 2A and Additional file 1: Table S1). Conversely, the average F<sub>ST</sub> between IP sites was 0.023 (ranging from 0.018 to 0.029) and not statistically significant for the majority of pairwise comparisons (Additional file 1: Table S1). The Mantel test, computed between IP sampling sites only (given the evidence of a clear genetic discontinuity between AO and IP), showed no correlation between genetic and geographic distances (r = 0.005, P = 0.62, Additional

file 1: Fig. S1). Clustering analyses were consistent with the observed pattern of differentiation. First, the *nmf* algorithm selected *K* = 2 ancestral populations corresponding to IP and AO (Fig. 2). Individuals from IP had a probability ancestry to cluster 1 ranging from 92.6 to 100% whereas individuals from AO showed a probability ancestry to cluster 2 ranging from 84.6 to 100%. Average ancestry of cluster 2 in IP individuals was only 0.7% while average ancestry of cluster 1 in AO individuals was 4.4%. Second, the PCA clearly segregated AO from IP individuals, with 38.71% of the total variance explained by the first axis (Fig. 3A). Individuals from Reunion Island (RUN), the IP site closest to the AO, did not show more proximity to the AO in the PCA or a higher contribution from cluster 2 than other IP individuals, nor did they show a lower pairwise F<sub>ST</sub> with the AO than the other IP sites. (Figs. 2, 3A and Additional file 1: Fig. S2A). When computed on IP individuals only, the PCA identified a single cluster (Fig. 3B and Additional file 1: Fig. S2B) and
<figure><img src="image_4.png" /></figure>

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

the *nmf* did not show any meaningful geographic clusters with *K*=2 (Additional file 1: Fig. S3). We further applied an Approximate Bayesian Computation (ABC) framework using a 500 trees random forest for all sampling sites after checking for the evolution of the out-of-bag error rate. This coalescent framework allows to detect genetic structure using a single sampling location by testing whether its gene genealogy yield signatures of a Stepping Stone (SS), Finite Island (FIM) or Non-Structured (NS) model (Additional file 1: Fig. S4). The model

selection (Additional file 1: Table S2) highlighted NS as the most supported model with a posterior probability ranging from 0.48 to 0.89 in the IP sampling sites and of 0.62 for BRA (Additional file 1: Table S3).

## Genetic diversity and variation of *N*<sub>e</sub>

Genetic diversity values were very similar among sampling sites, with BRA being slightly more variable than the IP counterpart (Table 1). Tajima's *D* (*TD*) was significantly positive for BRA (*TD*=0.137; *P* ≤ 0.001), while
significantly negative ($P \le 0.001$) for Northern Australia (AUS$_N$) and Reunion Island (RUN) and not significantly different from zero for the other populations (Table 1). Except for the AUS$_N$ population, the reconstructions of the effective size ($N_e$) through time by the STAIRWAYPLOT were very similar among IP locations: an ancestral expansion occurred between ~100,000 and ~200,000 years before present (BP). bringing the median $N_e$ to ~10,000 followed by a very recent bottleneck ~2000 to ~4000 years BP (Fig. 4). The STAIRWAYPLOT for AUS$_N$ displayed a different signal, with an ancestral $N_e$ median value similar to the one retrieved in the other sampling sites (~10,000) followed by a strong and recent expansion that raised the modern $N_e$ to ~35,000, contrasting with the recent decrease observed for the other IP

sampling sites. The demographic history reconstructed for BRA was slightly more complex with the ancestral $N_e$ of ~12,000 first decreasing to ~9000 at ~40,000 years BP and then increasing (between ~4000 and 6000 years BP) to a modern $N_e$ of ~20,000 (Fig. 4).

## Population divergence and migration rate estimation

In the light of the absence of population structure signatures within the IP and the AO and the genetic distinctness between them, we tested five Isolation-Migration (IM) models to determine the divergence time and the migration pattern between the two oceanic regions (Fig. 5). The likelihood distribution computed over 100 replicates was similar for models IM-bsc and IM-full,

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Variation of the effective population size ($N_e$) through time and its 75% confidence interval estimated by the STAIRWAYPLOT for the AO (panel A) and IP (panel B) sampling sites. AUS$_E$ East Coast of Australia, AUS$_N$ North Coast of Australia, BRA Brazil, COR Coral Sea, NCA New Caledonia, RUN Reunion Island</figcaption></figure>

<figure><img src="image_7.png" /></figure>
but the AIC values supported IM-bsc as the model with the highest probability (Additional file 1: Fig. S5). The distribution of the likelihood evaluated in each model under the Maximum Likelihood (ML) parameters proved that they can be distinguished based on the available data (Additional file 1: Fig. S5). The two oceanic regions appeared connected, though the migration rate is limited and strongly asymmetric, being ~176 times higher from IP to AO (Nm ~ 3.9) than vice versa (Nm ~ 0.02) (Table 2 and Additional file 1: Fig. S6). Going backward in time the populations from the two regions merged ~193,000 years BP (90% CI: [77,000; 355,000]). The ancestral population size was almost half of those estimated in both IP and AO derived populations (Table 2 and Additional file 1: Fig. S6), indicating an ancestral expansion, consistent with the observed STAIRWAYPLOT dynamics for IP populations. An increase in effective size was estimated in the AO starting ~45,000 years BP (90% CI: [14,000; 53,000]), bringing the effective size from 1,406 (90% CI: [1,178; 4,781]) to 16,810 (90% CI: [12,885; 40,036]) which is consistent with the observed expansion in the STAIRWAYPLOT of BRA. We observed a decrease in Ne in the IP from 49,000 to 17,000, but the timing was poorly estimated. Moreover, ancient and modern Ne in the IP showed largely overlapping confidence intervals (Table 2).

Discussion

To shed light on the evolutionary history of the tiger shark, we sequenced thousands of loci in 50 individuals following a double digest RAD-seq protocol. We handled low coverage issues by applying an appropriate framework based on genotype free likelihood estimation of allele frequencies [30]. The first result is the unambiguous

presence of two highly divergent genetic clusters, corresponding to the Indo-Pacific (IP) region and the Atlantic Ocean (AO), and the signature of very weak population structure within each of them. Despite the large panel of SNPs that could potentially detect fine spatial structure compared to previous work based on microsatellites, we did not find any barrier to gene flow within the IP, but rather a signature consistent with a large panmictic population or a meta-population characterized by very large amount of gene flow (Figs. 2, 3, Additional file 1: Figs. S1, S2, S3). This strongly confirms previous findings [23, 24, 27] and contradicts the conclusions of [22], who found significant evidence for population structure. Using a larger amount of genomic information, we provide evidence for the presence of a single mating population in the IP based on the following observations: (1) the PCA and *nmf* analyses display one single cluster in the IP; (2) the *F*<sub>ST</sub> values computed between sampling sites did not exceed 0.029 (as a comparison, an average value of ~0.124 was found between IP and AO) with no signature of isolation by distance (Additional file 1: Fig. S1); (3) all IP sharks showed a similar amount of genetic contribution from the AO cluster in the *nmf* analysis (Fig. 2). These results are consistent with one of two explanations: either (1) tiger sharks randomly mate within IP or (2) the number of migrants exchanged each generation between sampling sites is so large to erase any signature of genetic structure. The absence of multiple sampling sites in AO prevented us to perform similar analyses in AO. To test the presence of a single panmictic population, we therefore followed an ABC strategy based on coalescence simulations to reconstruct the evolutionary history of BRA population and assess whether the patterns of genetic variation within the BRA samples are better described by

<table><caption>**Table 2** Maximum Likelihood (ML), 90% confidence interval (5% lower bound and 95% upper bound) and search bounds for the parameters estimated by FASTSIMCOAL under the IM-bsc model</caption><thead><tr><th>Parameter</th><th>ML</th><th>5% Lower bound</th><th>95% Upper bound</th><th>Parameter bounds</th></tr></thead><tbody><tr><td>N<sub>anc</sub></td><td>9274</td><td>1310</td><td>80,372</td><td>U:{100;100,000}</td></tr><tr><td>N<sub>modIP</sub></td><td>17,591</td><td>16,643</td><td>31,915</td><td>U:{100;100,000}</td></tr><tr><td>N<sub>modAO</sub></td><td>16,810</td><td>12,885</td><td>40,036</td><td>U:{100;100,000}</td></tr><tr><td>T<sub>SI</sub> a</td><td>153,090</td><td>27,460</td><td>769,750</td><td>U:{100;1,000,000}</td></tr><tr><td>T<sub>SAO</sub> a</td><td>45,980</td><td>14,030</td><td>53,090</td><td>U:{100;1,000,000}</td></tr><tr><td>N<sub>ancIP</sub></td><td>48,823</td><td>3103</td><td>107,816</td><td>U:{100;100,000}</td></tr><tr><td>N<sub>ancAO</sub></td><td>1406</td><td>1178</td><td>4781</td><td>U:{100;100,000}</td></tr><tr><td>T<sub>DIV</sub> a</td><td>193,850</td><td>77,110</td><td>355,910</td><td>U:{100;1,000,000}</td></tr><tr><td>N<sub>AO/IP</sub> <sup>b</sup></td><td>0.022</td><td>0.008</td><td>0.124</td><td>U:{0;50}</td></tr><tr><td>N<sub>IP/AO</sub> <sup>b</sup></td><td>3.873</td><td>3.171</td><td>8.919</td><td>U:{0;50}</td></tr></tbody></table>

<footnote><sup>a</sup> Times are expressed in years using a mutation rate of 1.93 × 10⁻⁸ per generation per site and a generation time of 10 years</footnote>

<footnote><sup>b</sup> Number of migrants per generation are expressed in forward</footnote>

<footnote><sup>c</sup> uniform distribution</footnote>
unstructured or meta-population models. This approach has been successfully applied in both empirical and simulation-based works [7, 31, 32] and represents an alternative (or complementary, when possible) method to infer the presence of population structure. As found for the IP, the NS (No Structure) model had the most support within the AO (Additional file 1: Table S3). Even though more sampling sites and SNPs would likely have provided tighter estimates, this result is most consistent with a single mating population in the AO, despite our single sampling site does not allow us to infer the geographic extent of that population. We note that population structure has previously been reported within the AO based on mitochondrial markers [25, 26]. However, results in the current study are consistent with the genome wide study of [27] which suggested low to no population structure in the AO. The differences between the inferred mitochondrial structure and genomic DNA signals could be due to sex-biased dispersal, as female philopatry has been proposed for some shark species [33–36]. However, inconsistencies between mitochondrial and autosomal data are widespread in nature and it has been suggested to causively interpret the mitochondrial variation in the light of the demographic history of a species, as other evolutionary forces such selection may act as confounding factors [37]. Furthermore, the result found herein (i.e., low to no population structure in each of the two oceanic regions) is consistent with the fact that tiger sharks have been documented to move large distances across oceanic basins. However, it does not explain why a large predator that is capable of covering distances of several thousands of kilometers [13–15] could not erase the genetic differentiation between the two regions.

One possible explanation was proposed by [21], who suggested the presence of two allopatric subspecies in IP and AO. This hypothesis was later refuted by [22], who still agreed on a long-term genetic isolation between the two oceanic regions but proposed some genetic exchanges. By harnessing the power of RAD-seq genome wide data, our coalescent modelling could not only disentangle the two hypotheses, but also provide quantitative estimates of the tempo and mode of divergence between the two populations. Comparing five Isolation/Migration (IM) scenarios, we found that the most supported model included a divergence around ~193,000 years BP (90% CI: [77,000; 356,000]) between the two regions (Table 2), which nevertheless remained in contact since then through a very limited (3.9 individuals per generation; 90% CI: [3.2; 8.9]) and asymmetric gene flow ~176 times higher from IP to AO than the opposite direction. The low number of migrants *Nm* exchanged each generation (Table 2) and the asymmetric exchange are consistent with the clustering results and the *F*<sub>ST</sub> values (Fig. 2),

which differentiated the two regions and clearly highlighted a higher, but weak, genetic contribution from IP to AO than vice versa (Fig. 2). These results support the idea that populations from the AO and from the IP represent two lineages [22], rather than two allopatric species [21]. A permeable barrier to gene flow between the two oceanic regions is therefore responsible for the observed pattern of divergence and asymmetric migration. The presence of this barrier can be explained by the ecology of the tiger shark and by the environmental conditions governing the Indian-Atlantic water exchange, the so-called Agulhas leakage. As a tropical to sub-tropical species, tiger sharks prefer warm water and they show the peak of swimming activities at ~22 °C [12] so that their movement from the Atlantic to the Indo-Pacific is hampered by the upwelling of cold water off South-Western Africa (the Benguela current). However, the AO receives warm water from the IP through the Agulhas leakage [38], which can account for the asymmetric migration reported, consistent with the pattern observed in other tropical sharks, bony fishes and turtles [39–42]. The Agulhas leakage has not been constant through time, with an increasing intensity in the Holocene, preceded by a period of stasis and a strong peak in the late Pleistocene around 130,000 years BP [43, 44]. This variation in Agulhas leakage intensity could have influenced the relation between the two basins. However, we could not distinguish pulses of migrations in our data since the model IM-bsc was preferred to those accounting for variation in migration rate through time. Model selection procedure robustly selected the IM-bsc model (Additional file 1: Fig. S5), which is neither the most nor the least parameters rich, supporting the idea that our results are not an artefact of incorrect modelling (see also below). In the future there will still be space to improve our estimates: more data will help refining the confidence interval of each parameter and ameliorating the calibration of the molecular clock.

We found divergent demographic histories in the two oceanic regions examined (Fig. 4). First, we note that since both regions are most likely described by non-structured models and the migration rates between them are very low, it is possible to directly interpret the results of unstructured model such as the STAIRWAYPLOT [45]. It is important to stress this point, because population structure, if not accounted for, can generate signatures that erroneously look like changes in effective population size [32, 46, 47]. All the Indo-Pacific populations (except AUSN) underwent a recent bottleneck between ~2000 and 4000 years BP, which was robust to the pooling of all IP sampling points but AUSN (but with larger incertitude in recent times due to the lower number of SNPs, see Additional file 1: Fig. S7). Despite the caveat regarding
the interpretation of *N*<sub>e</sub> variation in recent times due to the inclusion of singletons, this is barely consistent with what was recently proposed by [24] based on 25 microsatellites combined with mitochondrial DNA. Here we refined their estimates and better characterized the intensity of the bottleneck with a non-parametric approach (the STAIRWAYPLOT) exploring a large parameter space with genome wide data. The inferred demographic history of the AO is strikingly different from that of the IP. The estimated *N*<sub>e</sub> was ~20,000 following a population expansion occurring between ~ 4000 and ~ 6000 years BP (Fig. 4). These values are consistent with the estimates retrieved by the IM-bsc model (Table 2). The strong signature of population expansion recovered implies that the tiger shark is profiting from recent environmental changes in AO, in contrast to the IP population. Consistently, [26, 48] found a recent demographic trend suggesting an expansion rather than a decrease in AO, while the recent demographic trends in the IP appear to have been the result of intense pressure from fisheries and shark-control programs [9–11]. More investigations are needed to determine the origin of the difference between the two oceans. The applications of methods based on linkage disequilibrium applied to whole genome data will help detect more recent events [49], which will be important for planning conservation strategies. Given the very low genetic structure within IP, we would expect AUS<sub>N</sub> to have the same demographic history than the other sampling sites in the IP. We cannot exclude a scenario where AUS<sub>N</sub> represents an isolated population experiencing its own demographic history that separated from the rest of the IP too recently to accumulate divergence. If confirmed, this would highlight the presence of independent lineages in the IP, with important consequence for conservation programs. However, more data is needed to shed light on this topic, both in terms of individuals and in terms of genomic coverage: ultimately only whole genome sequencing will give the opportunity to confidently resolve this issue.

## Conclusions

Reconstructing the evolutionary history of a species relies on the application of a realistic demographic model, which is mostly unknown in empirical studies. Here, we investigated the evolutionary history of the tiger shark and found that it is characterized by an asymmetrical migration between the AO and the IP and a signature of random mating within each region. These findings let us model each oceanic region as a single unstructured population and evaluate competing demographic scenarios to investigate their divergence time and migration rates. The two regions are separated by the Benguela barrier, but our estimates strongly suggest that the Agulhas

leakage allows an overwhelmingly asymmetric migration between them, by far stronger from IP to AO than in the opposite direction. While we confirmed that the tiger shark is likely undergoing a reduction in *Ne* in the Indo-Pacific, we show that it probably underwent a strong expansion in the Atlantic Ocean. Even if a better calibration of the molecular clock and full genome analyses would still be needed to confirm our results, our findings support the existence of two management units. This implies that local conservation or shark control programs will have very limited impact on the dynamics of the species, which needs to be managed at the ocean basin level, demanding considerable communication efforts among different countries and coordination as suggested for other megafaunal organisms [50].

## Material and methods

### Sampling, library preparation and sequencing

A total of 50 tiger shark individuals (*Galeocerdo cuvier*) from both the Indo-Pacific (IP) and the Atlantic Ocean (AO) were sampled off (from west to east) Brazil (BRA), Reunion Island (RUN), North Coast of Australia (AUS<sub>N</sub>; North Territory), East Coast of Australia (AUS<sub>E</sub>; Sunshine Coast), Coral Sea (COR) and New Caledonia (NCA). Sharks were grouped into six populations based on their sampling site (Fig. 1; Table 1). Total genomic DNA was extracted from muscle tissue or fin clips preserved in 96% ethanol using QIAGEN DNeasy Blood and Tissue kit (Qiagen, Hilden, Germany) according to the manufacturer's protocols. Double-digest restriction-associated DNA (ddRAD) libraries were prepared following [51] using EcoRI and MspI restriction enzymes, a 400-bp size selection, and a combination of two indexes and 24 barcodes to pool 48 individuals per lane. The genomic libraries obtained were sequenced with a HiSeq 2500 Illumina sequencer (single-end, 125 bp).

### RAD-seq de novo assembly

Raw reads were first demultiplexed and quality filtered through the *process_radtags.pl* pipeline in *STACKS v.2.5* [52]. In the absence of a reference genome of *G. cuvier* or of closely related species, RAD-seq loci (125 bp sequences) were de novo assembled under the *denovo_map.pl* pipeline in *STACKS*. Preliminary results based on parameters *m*=3 (minimum read depth to create a stack), *M*=3 (number of mismatches allowed between loci within individuals), and *n*=3 (number of mismatches allowed between loci within catalogue) found an average depth of ~10x (see Results). Despite the absence of a clear cut-off indicating an acceptable coverage value above which genotype calling may be considered reliable, simulation results suggest that ~10x may produce inconsistent calling under different algorithms [53]. To prevent
with custom script from the SFS. Significance of *TD* was evaluated after 1000 coalescent simulations of a constant population model with size θ<sub>w</sub>. We reconstructed the variation in the effective population size (*N*<sub>e</sub>) through time by running the STAIRWAYPLOT v.0.2 software [45] with singletons, where the composite likelihood is evaluated as the difference between the observed (folded) SFS and its expectation under a specific demographic history.

Population divergence and migration rate estimation

Based on the results of the previous analyses, we devised
five alternative Isolation/Migration (IM) models of diver-
gence between IP and AO regions using the composite
likelihood method implemented in FASTSIMCOAL [59].
We presented in Fig. 3 the model richest in parameters,
the remaining four representing simplified versions
nested within it. Hereafter, a brief description of the five
models going from the most complex to the simplest:
(a) **IM-full**: the two ocean regions with their respective
modern effective population sizes, N<sub>modIP</sub> and N<sub>modAO</sub>
diverged at T<sub>div</sub> from an ancestral population of effec-
tive population size N<sub>anc</sub>. Due to the STAIRWAYPLOT
results, we allowed the two modern effective population
sizes N<sub>modIP</sub> and N<sub>modAO</sub> to change to N<sub>ancIP</sub> and N<sub>ancAO</sub>
following an exponential dynamic in T<sub>SIP</sub> and T<sub>SAO</sub> years
respectively. Migration is defined by two time periods:
m<sub>1</sub> representing the migration rate occurring between
time 0 until T<sub>mig</sub> and m<sub>2</sub> between time T<sub>mig</sub> until T<sub>div</sub>. The
migration matrix in each time period is asymmetric: for
instance, m<sub>1AO/IP</sub> represents the forward migration rate
from AO to IP and m<sub>1IP/AO</sub> from IP to AO. In summary,
the model is defined by thirteen parameters: five effective
population sizes, four migration rates and four historical
events; (b) **IM-anc**: same as IM-full with ancestral migra-
tion only between T<sub>mig</sub> and T<sub>div</sub>. The model is defined
by eleven parameters, the two m<sub>1</sub> migration rates being
removed; (c) **IM-rec**: same as IM-anc, but with recent
migration only occurring between time 0 until T<sub>mig</sub>, keep-
ing only the two m<sub>1</sub> migration rates; (d) **IM-bsc**: the clas-
sic model where migration is constant from time 0 until
T<sub>div</sub> (i.e. m<sub>1</sub> = m<sub>2</sub> = m [28]). We modelled the variation
in effective size of the two regions similarly to the other
models, for a total of ten parameters; (e) **IM-div**: a pure
divergence model with no migration. This is defined by
eight parameters: the five effective population sizes and
three historical events (T<sub>div</sub>, T<sub>SIP</sub>, T<sub>SAO</sub>). The analyses
are based on the folded 2D-SFS computed by ANGSD
between six individuals from Brazil (representing the
AO) and six from the Indo-Pacific (IP). This sample size
was chosen to obtain a balanced design and to maxim-
ise the number of SNPs shared among the two ocean
basins. Similarly, in each basin, we selected the individu-
als presenting the smaller proportion of missing data to

further increase the number of joint SNPs. To maximize
the observed 2D-SFS we applied the following options in
FASTSIMCOAL: -N 300,000 (number of coalescent simula-
tions), -L 40 (number of expectation-maximization (EM)
cycles), and -C 10 (minimum observed SFS entry count
considered for parameter estimation). For all model
parameters we used wide search ranges with uniform
distributions (Table 2). We ran each model 100 times in
order to determine the maximum likelihood parameters
and to perform model selection using the Akaike’s infor-
mation criterion comparing the best run of each model
[59]. To check the robustness of the model selection pro-
cedure and to take into account the presence of linked
sites in our dataset, we further examined the likelihood
distribution obtained based on 100 expected 2D-SFS sim-
ulated under the parameters estimated in the best run of
each model, each approximated with 10⁶ coalescent sim-
ulations. This procedure is needed to take into account
the variance in the likelihood estimation given our dataset: if the distributions obtained by the various models do
overlap, the difference in the estimated likelihoods of our
models is not significant [63]. Finally, we determined the
confidence interval of the parameter estimated under the
best run of our best model by parametric bootstrapping.
The 2D-SFS was bootstrapped 100 times using FASTSIM-
COAL and each of these datasets was analysed under the
same conditions as the original data (one hundred inde-
pendent runs for each dataset). Calibrating the molecular
clock is crucial to obtain accurate estimates of demo-
graphic parameters and historical events, but it is chal-
lenging when fossil records and/or orthologous loci from
an outgroup are lacking. Here, all demographic infer-
ences were performed using the RAD-seq mutation rate
of μ = 1.93 × 10⁻⁸ per site and per generation previously
used for the tiger shark [7], and the generation time was
set to 10 years [24, 64].

**Supplementary Information**

The online version contains supplementary material available at https://doi.org/10.1186/s12862-022-02100-y.

**Additional file 1: Figure S1.** Isolation by distance (IBD) plot within the Indo-Pacific. Pairwise genetic distances ($F_{ST}$/(1-$F_{ST}$)) are plotted against geographic distances between Indo-Pacific sampling sites. **Figure S2.** Principal Component Analysis (PCA) computed with: (A) all individuals (n = 50) and (B) Indo-Pacific individuals only (n = 43). The axes represented in both panels are the first and the third component. **Figure S3.** Ancestry proportions retrieved using the *nmf* algorithm with K = 2 ancestral populations for Indo-Pacific samples performed with *PCAANGSD*. **Figure S4.** Evolutionary scenarios used to investigate the population structure of the Atlantic Ocean based on data from Brazil population through an Approximate Bayesian Computation (ABC) framework. NS (No Structure) is an unstructured model where the modern effective size ($N_{mod}$) instantaneously changes to $N_{anc}$ at time shift $T_s$ generations. FIM (Finite Island Meta-population) represents a finite island meta-population model with 100 demes that have been instantaneously colonised $T_{col}$ generations ago, from an ancestral population of size $N_{anc}$. Demes are allowed to
exchange migrants with any other. SS (Stepping-Stone) is similar to FIM but the migrants are only exchanged between the four nearest neighbours in a two-dimensional grid. **Figure S5.** Akaike Information Criterion (AIC) values for the five isolation/migration models and the associated ranking on the x-axis. Boxplots represent the likelihood distribution of the data evaluated under the best parameter estimates for each of the five models (presented in Fig. 2) after 100 replicates. The models are presented from the richest in parameters (IM-full, 13 parameters) to the poorest (IM-div, 8 parameters). **Figure S6.** Maximum likelihood for the parameter estimated by fastsimcoal under model IM-bsc, representing two populations from each ocean basin with an effective size that changed $T_{SIP}$ and $T_{SAO}$ years ago from a modern effective size ($N_{modIP}$ and $N_{modAO}$) to an ancestral effective size ($N_{ancIP}$ and $N_{ancAO}$). The two populations are connected by an asymmetrical number of migrants constant from 0 to $T_{div}$ ($Nm_{IP→AO}$ and $Nm_{AO→IP}$) and diverged $T_{div}$ years ago from an ancestral population of size $N_{anc}$.

**Acknowledgements**

We thank Serge Planes for providing DNA samples of all specimens but Reunion Island, Gavin Naylor for critical reading of the manuscript, and Thibaut Caley for discussion above the Agulhas current. We thank Andrea Benazzo for help with ANGSD. We are grateful to the Genotoul bioinformatics platform Toulouse Midi-Pyrenees (Bioinfo Genotoul; http://bioinfo.genotoul.fr) for providing computing resources. All data and analyses from Reunion Island were produced with the help of A. Pirog and financially supported by DEAL Réunion under the project ECoReCo-Run.

**Author contributions**

S.M. and H.M. conceived the project. S.M., PL and H.L. analysed the data with input from P.M.D. E.C., H.M. and S.J. provided reagents and samples. A.S. and P.C.B performed the molecular lab work. PL., S.M. and H.M. wrote the manuscript with input from all the others. All authors read and approved the final manuscript.

**Funding**

This work was supported by two ATM grants (2016 and 2017) from the Muséum National d'Histoire Naturelle to Stefano Mona.

**Availability of data and materials**

Fastq sequence files are available from the GenBank at the National Center for Biotechnology Information short-read archive database (BioProject ID: PRJNA887936).

**Declarations**

**Ethics approval and consent to participate**

This research (including experimental protocols) was conducted in accordance with Permit No. 6024-4916/DENV/SMer (New Caledonia), Permit No. QS2010 GS065 and Permit No. 143005 (Australia), and Permit No. TREL21187155/S41 (Reunion). Samples from Brazil were provided by Jennifer Schults from Hawaii, in accordance with local regulations. All methods were carried out in accordance with relevant guidelines and regulations (certificate of compliance ABSCH-IRCC-FR-259518-1) and in accordance with ARRIVE guidelines (https://arriveguidelines.org).

**Consent for publication**

Not applicable.

**Competing interests**

The authors declare no conflict of interest.

**Author details**

Received: 8 May 2022 Accepted: 9 December 2022

Published online: 16 December 2022

**References**

1. Terborgh JW. Toward a trophic theory of species diversity. Proc Natl Acad Sci U S A. 2015;112(37):11415-22.

2. Myers RA, Worm B. Rapid worldwide depletion of predatory fish communities. Nature. 2003;423(6937):280-3.

3. Palsbøll PJ, Bérubé M, Allendorf FW. Identification of management units using population genetic data. Trends Ecol Evol. 2007;22(1):11-6.

4. Chikhi L, Sousa VC, Luisi P, Goossens B, Beaumont MA. The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes. Genetics. 2010;186(3):983-95.

5. Mona S, Ray N, Arenas M, Excoffier L. Genetic consequences of habitat fragmentation during a range expansion. Heredity (Edinb). 2014;112(3):291-9.

6. Mazet O, Rodríguez W, Chikhi L. Demographic inference using genetic data from a single individual: separating population size variation from population structure. Theor Popul Biol. 2015;104:46-58.

7. Lesturgie P, Planes S, Mona S. Coalescence times, life history traits and conservation concerns: an example from four coastal shark species from the Indo-Pacific. Mol Ecol Resour. 2022;22(2):554-66. https://doi.org/10.1111/1755-0998.13487.

8. Compagno LJV. FAO species catalogue, Vol. 4: Sharks of the world: an annotated and illustrated catalogue of shark species known to date. Part 2—Carcharhiniformes. In: FAO Fisheries synopsis. Rome, Italy; 1984. p. 503-12.

9. Temple AJ, Kiszka JJ, Stead SM, Wambiji N, Brito A, Poonian CNS, et al. Marine megafauna interactions with small-scale fisheries in the southwestern Indian Ocean: a review of status and challenges for research and management. Rev Fish Biol Fish. 2018;28(1):89-115.

10. Sumpton W, Taylor S, Gribble N, McPherson G, Ham T. Gear selectivity of large-mesh nets and drumlines used to catch sharks in the Queensland Shark Control Program. Afr J Mar Sci. 2011;33(1):37-43.

11. Clarke SC, McAllister MK, Milner-Gulland EJ, Kirkwood GP, Michielsens CGJ, Agnew DJ, et al. Global estimates of shark catches using trade records from commercial markets. Ecol Lett. 2006;9(10):1115-26. https://doi.org/10.1111/j.1461-0248.2006.00968.x.

12. Payne NL, Meyer CG, Smith JA, Houghton JDR, Barnett A, Holmes BJ, et al. Combining abundance and performance data reveals how temperature regulates coastal occurrences and activity of a roaming apex predator. Glob Chang Biol. 2018;24(5):1884-93.

13. Holland KN, Anderson JM, Coffey DM, Holmes BJ, Meyer CG, Royer MA. A perspective on future tiger shark research. Front Mar Sci. 2019;6(FEB):1-7.

14. Werry JM, Planes S, Berumen ML, Lee KA, Braun CD, Clua E. Reef-fidelity and migration of tiger sharks, Galeocerdo cuvier, across the coral sea. PLoS ONE. 2014;9(1): e83249.

15. Lea JSE, Wetherbee BM, Queiroz N, Burnie N, Aming C, Sousa LL, et al. Repeated, long-distance migrations by a philopatric predator targeting highly contrasting ecosystems. Sci Rep. 2015;5(1):11202.

16. Meyer CG, Anderson JM, Coffey DM, Hutchinson MR, Royer MA, Holland KN. Habitat geography around Hawaii's oceanic islands influences tiger shark (Galeocerdo cuvier) spatial behaviour and shark bite risk at ocean recreation sites. Sci Rep. 2018;8(1):4945.

17. Palumbi SR. Genetic divergence, reproductive isolation, and marine speciation. Annu Rev Ecol Syst. 1994;25:547-72.

18. Bailleul D, Mackenzie A, Sacchi O, Poisson F, Bierne N, Arnaud-Haond S. Large-scale genetic panmixia in the blue shark (Prionace glauca): a single worldwide population, or a genetic lag-time effect of the "grey zone" of differentiation? Evol Appl. 2018;11(5):614-30.

19. Corrigan S, Lowther AD, Beheregaray LB, Bruce BD, Cliff G, Duffy CA, et al. Population connectivity of the Highly Migratory Shortfin Mako (Isurus oxyrinchus Rafinesque 1810) and implications for management in the Southern Hemisphere. Front Ecol Evol. 2018;6(NOV):1-15.

20. Corrigan S, MaisanoDelser P, Eddy C, Duffy C, Yang L, Li C, et al. Historical introgression drives pervasive mitochondrial admixture between two species of pelagic sharks. Mol Phylogenet Evol. 2017;110(December):122-6.

<footnote n="1">¹Institut de Systématique, Evolution, Biodiversité (ISYEB), EPHE-PSL, Université PSL, MNHN, CNRS, SU, UA, Paris, France. ²EPHE, PSL Research University, Paris, France. ³Department of Zoology, University of Cambridge, Cambridge, UK.</footnote>

<footnote n="4">⁴Laboratoire d'Excellence CORAIL, Papetoai, French Polynesia. ⁵UMR ENTROPIE (Université de La Réunion/IRD/CNRS), Université de La Réunion, Saint Denis, France.</footnote>
21. Naylor GJP, Caira JN, Jensen K, Rosana KAM, White WT, Last PR. A DNA sequence based approach to the identification of shark and ray species and its implications for global elasmobranch diversity and parasitology. *Bull Am Mus Nat Hist*. 2012;367:1–262.

22. Bernard AM, Feldheim KA, Heithaus MR, Wintner SP, Wetherbee BM, Shivji MS. Global population genetic dynamics of a highly migratory, apex predator shark. *Mol Ecol*. 2016;25(21):5312–29.

23. Holmes BJ, Williams SM, Otway NM, Nielsen EE, Maher SL, Bennett MB, et al. Population structure and connectivity of tiger sharks (*Galeocerdo cuvier*) across the Indo-Pacific Ocean basin. *R Soc Open Sci*. 2017;4(7):170309.

24. Pirog A, Jaquemet S, Ravigné V, Cliff G, Clua E, Holmes BJ, et al. Genetic population structure and demography of an apex predator, the tiger shark *Galeocerdo cuvier*. *Ecol Evol*. 2019;9(10):5551–71. https://doi.org/10.1002/eece3.5111.

25. Carmo CB, Ferrette BLS, Camargo SM, Roxo FF, Coelho R, Garla RC, et al. A new map of the tiger shark (*Galeocerdo cuvier*) genetic population structure in the western Atlantic Ocean: hypothesis of an equatorial convergence centre. *Aquat Conserv*. 2019;29(5):760–72.

26. Andrade FRS, Afonso AS, Hazin FHV, Mendonça FF, Torres RA. Population genetics reveals global and regional history of the apex predator *Galeocerdo cuvier* (carcharhiniformes) with comments on mitigating shark attacks in north-eastern brazil. *Mar Ecol*. 2021;42(2):1–16. https://doi.org/10.1111/maec.12640.

27. Bernard AM, Finnegan KA, PavinskiBitar P, Stanhope MJ, Shivji MS. Genomic assessment of global population structure in a highly migratory and habitat versatile apex predator, the tiger shark (*Galeocerdo cuvier*). *J Hered*. 2021;112(6):497–507.

28. Nielsen R, Wakeley J. Distinguishing migration from isolation: a Markov chain Monte Carlo approach. *Genetics*. 2001;158(2):885–96.

29. Beerli P, Felsenstein J. Maximum likelihood estimation of a migration matrix and effective population sizes in *n* subpopulations by using a coalescent approach. *Proc Natl Acad Sci U S A.* 2001;98(8):4563–8.

30. Korneliussen TS, Albrechtsen A, Nielsen R. ANGSD: analysis of next generation sequencing data. *BMC Bioinform*. 2014;15(1):1–13.

31. Peter BM, Wegmann D, Excoffier L. Distinguishing between population bottleneck and population subdivision by a Bayesian model choice procedure. *Mol Ecol*. 2010;19(21):4648–60.

32. MaisanoDelser P, Corrigan S, Duckett D, Suwalski A, Veuille M, Planes S, et al. Demographic inferences after a range expansion can be biased: the test case of the blacktip reef shark (*Carcharhinus melanopterus*). *Heredity* (Edinb). 2019;122(6):759–69. https://doi.org/10.1038/s41437-018-0164-0.

33. Pardini AT, Jones CS, Noble LR, Kreiser B, Malcolm H, Bruce BD, et al. Sex-biased dispersal of great white sharks. *Nature*. 2001;412(6843):139–40.

34. Keeney DB, Heupel MR, Hueter RE, Heist EJ. Microsatellite and mitochondrial DNA analyses of the genetic structure of blacktip shark (*Carcharhinus limbatus*) nurseries in the northwestern Atlantic, Gulf of Mexico, and Caribbean Sea. *Mol Ecol*. 2005;14(7):1911–23.

35. Tillett BJ, Meekan MG, Field IC, Thorburn DC, Ovenden JR. Evidence for reproductive philopatry in the bull shark *Carcharhinus leucas*. *J Fish Biol*. 2012;80(6):2140–58.

36. Mourier J, Planes S. Direct genetic evidence for reproductive philopatry and associated fine-scale migrations in female blacktip reef sharks (*Carcharhinus melanopterus*) in French Polynesia. *Mol Ecol*. 2013;22(1):201–14.

37. Ballard JWO, Whitlock MC. The incomplete natural history of mitochondria. *Mol Ecol*. 2004;13(4):729–44.

38. Beal LM, De RWPM, Biastoch A, Zahn R, Wcrp S, Working I. On the role of the Agulhas system in ocean circulation and climate. *Nature*. 2011;472(7344):429–36.

39. Gaither MR, Bowen BW, Rocha LA, Briggs JC. Fishes that rule the world: circumtropical distributions revisited. *Fish Fish*. 2016;17(3):664–79.

40. Maduna SN, Rossouw C, da Silva C, Soekoe M, Bester-van der Merwe AE. Species identification and comparative population genetics of four coastal houndsharks based on novel NGS-mined microsatellites. *Ecol Evol*. 2017;7(5):1462–86.

41. Reid BN, Naro-Maciel E, Hahn AT, FitzSimmons NN, Gehara M. Geography best explains global patterns of genetic diversity and postglacial co-expansion in marine turtles. *Mol Ecol*. 2019;28(14):3358–70.

42. van der Zee JP, Christianen MJA, Bérubé M, Nava M, Schut K, Humber F, et al. The population genomic structure of green turtles (*Chelonia mydas*) suggests a warm-water corridor for tropical marine fauna between the

Atlantic and Indian oceans during the last interglacial. *Heredity* (Edinb). 2021;127(6):510–21.

43. Caley T, Giraudeau J, Malaizé B, Rossignol L, Pierre C. Agulhas leakage as a key process in the modes of Quaternary climate changes. *Proc Natl Acad Sci U S A.* 2012;109(18):6835–9.

44. Caley T, Peeters FJC, Biastoch A, Rossignol L, van Sebille E, Durgadoo J, et al. Quantitative estimate of the paleo-Agulhas leakage. *Geophys Res Lett*. 2014;41(4):1238–46.

45. Liu X, Fu YX. *Airway Plot 2: demographic history inference with folded SNP frequency spectra*. *Genome Biol*. 2020;21(1):1–9.

46. Chikhi L, Rodríguez W, Grusea S, Santos P, Boitard S, Mazet O. The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: insights into demographic inference and model choice. *Heredity* (Edinb). 2018;120(1):13–24.

47. Mazet O, Rodríguez W, Grusea S, Boitard S, Chikhi L. On the importance of being structured: instantaneous coalescence rates and human evolution—lessons for ancestral population size inference? *Heredity* (Edinb). 2016;116(4):362–71.

48. Peterson CD, Belcher CN, Bethea DM, Driggers WB, Frazier BS, Latour RJ. Preliminary recovery of coastal sharks in the south-east United States. *Fish Fish*. 2017;18(5):845–59.

49. Kerdoncuff E, Lambert A, Achaz G. Testing for population decline using maximal linkage disequilibrium blocks. *Theor Popul Biol*. 2020;134:171–81.

50. Barkley AN, Gollock M, Samoilys M, Llewellyn F, Shivji M, Wetherbee B, et al. Complex transboundary movements of marine megafauna in the Western Indian Ocean. *Anim Conserv*. 2019;22(5):420–31.

51. Peterson BK, Weber JN, Kay EH, Fisher HS, Hoekstra HE. Double digest RADseq: an inexpensive method for de novo snp discovery and genotyping in model and non-model species. *PLoS ONE*. 2012;7(5):e37135.

52. Rochette NC, Rivera-Colón AG, Catchen JM. Stacks 2: analytical methods for paired-end sequencing improve RADseq-based population genomics. *Mol Ecol*. 2019;28(21):4737–54.

53. Fountain ED, Pauli JN, Reid BN, Palsball PJ, Peery MZ. Finding the right coverage: the impact of coverage and sequence quality on single nucleotide polymorphism genotyping error rates. *Mol Ecol Res*. 2016;16(4):966–78.

54. Khimoun A, Doums C, Molet M, Kaufmann B, Peronnet R, Eyer PA, et al. Urbanization without isolation: the absence of genetic structure among cities and forests in the tiny acorn ant *Temnothorax nylanderi*. *Biol Lett*. 2020;16(1):20190741. https://doi.org/10.1098/rslb.2019.0741

55. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. *Bioinformatics*. 2009;25(14):1754–60.

56. Meisner J, Albrechtsen A. Inferring population structure and admixture proportions in low-depth NGS data. *Genetics*. 2018;210(2):719–31.

57. Nielsen R, Korneliussen T, Albrechtsen A, Li Y, Wang J. SNP calling, genotype calling, and sample allele frequency estimation from new-generation sequencing data. *PLoS ONE*. 2012;7(7):e37558.

58. Mantel N. The detection of disease clustering and a generalized regression approach. *Cancer Res*. 1967;27(2):209–20.

59. Excoffier L, Dupanloup I, Huerta-Sánchez E, Sousa VC, Foll M. Robust demographic inference from genomic and SNP data. *PLoS Genet*. 2013;9(10):e1003905. https://doi.org/10.1371/journal.pgen.1003905.

60. Pudlo P, Marin JMM, Estoup A, Cornuet JMM, Gautier M, Robert CP. Reliable ABC model choice via random forests. *Bioinformatics*. 2016;32(6):859–66.

61. Watterson GAA. On the number of segregating sites in genetical models without recombination. *Theor Popul Biol*. 1975;7(2):256–76.

62. Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. *Genetics*. 1989;123(3):585–95.

63. Meier JI, Sousa VC, Marques DA, Selz OM, Wagner CE, Excoffier L, et al. Demographic modelling with whole-genome data reveals parallel origin of similar *Pundamilia cichlid* species after hybridization. *Mol Ecol*. 2017;26(1):123–41. https://doi.org/10.1111/mec.13838.

64. Cortés E. Incorporating uncertainty into demographic modeling: application to shark populations and their conservation. *Conserv Biol*. 2002;16(4):1048–62.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.