Skip to main content

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

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), FST 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.

Peer Review reports

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 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,5,6,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,15,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,22,23,24,25,26,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.

Fig. 1
figure 1

Map of the sampling sites. From west to east: Brazil (BRA: n = 7), Reunion Island (RUN; n = 15), North Coast of Australia (AUSN; n = 8), Coral Sea (COR; n = 5), East Coast of Australia (AUSE; n = 7) and New Caledonia (NCA; n = 8)

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 v.2.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).

Table 1 Sample size (n), mean pairwise difference (θπ), Watterson theta (θw), Tajima’s D (TD), and total number of loci (monomorphic included) (nloci) and SNPs (nSNP) without missing data for all sampling sites (ranged from west to east)

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 FST computation. Pairwise FST 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 FST 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 FST 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 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).

Fig. 2
figure 2

Heat map representing the pairwise Reynold’s FST values between sampling sites (A) and ancestry proportions retrieved using the nmf algorithm with K = 2 ancestral populations (B). Both analyses were performed with PCAngsd. Values in the upper triangle of the heat map are the pairwise FST values and significance is displayed on the lower triangle: non-significant (NS) or p < 0.001 (*)

Fig. 3
figure 3

Principal Component Analysis (PCA) computed with: (A) all individuals (n = 50) and (B) Indo-Pacific individuals only (n = 43)

Genetic diversity and variation of N e

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 ≤ 0.001) for Northern Australia (AUSN) and Reunion Island (RUN) and not significantly different from zero for the other populations (Table 1). Except for the AUSN population, the reconstructions of the effective size (Ne) 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 Ne to ~ 10,000 followed by a very recent bottleneck ~ 2000 to ~ 4000 years BP (Fig. 4). The stairwayplot for AUSN displayed a different signal, with an ancestral Ne 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 Ne 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 Ne of ~ 12,000 first decreasing to ~ 9000 at ~ 40,000 years BP and then increasing (between ~ 4000 and 6000 years BP) to a modern Ne of ~ 20,000 (Fig. 4).

Fig. 4
figure 4

Variation of the effective population size (Ne) through time and its 75% confidence interval estimated by the stairwayplot for the AO (panel A) and IP (panel B) sampling sites. AUSE East Coast of Australia, AUSN North Coast of Australia, BRA Brazil, COR Coral Sea, NCA New Caledonia, RUN Reunion Island

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, 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).

Fig. 5
figure 5

Model IM-full, the most parameter-rich model (13 parameters) representing two populations from each ocean basin with an effective size that changed \({T}_{{s}_{IP}}\) and \({T}_{{s}_{AO}}\) years ago from a modern effective size (\({N}_{{mod}_{IP}}\) and \({N}_{{mod}_{AO}}\)) to an ancestral effective size (\({N}_{{anc}_{IP}}\) and \({N}_{{anc}_{AO}}\)). The two populations are connected by an asymmetrical migration rate allowed to change \({T}_{mig}\) years ago (respectively from \({m}_{{1}_{IP/AO}}\) and \({m}_{{1}_{AO/IP}}\) to \({m}_{{2}_{IP/AO}}\) and \({m}_{{2}_{AO/IP}}\)) and diverged \({T}_{div}\) years ago from an ancestral population of size \({N}_{anc}\). The remaining four models are nested within IM-full, having less migration rate parameters: IM-anc is similar to IM-full but only the ancestral migration occurs (i.e., between \({T}_{mig}\) and \({T}_{div}\)); IM-rec is similar to IM-full but only the recent migration occurs (i.e., between 0 and \({T}_{mig}\)); IM-bsc considers the migration constant from 0 to \({T}_{div}\); and IM-div is a strict divergence model with no migration

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

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 FST 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 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,34,35,36]. However, inconsistencies between mitochondrial and autosomal data are widespread in nature and it has been suggested to cautiously 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,14,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 FST 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,40,41,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 Ne 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 Ne 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,10,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 AUSN to have the same demographic history than the other sampling sites in the IP. We cannot exclude a scenario where AUSN 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 (AUSN; North Territory), East Cost of Australia (AUSE; 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 this, we used a genotype free estimation of allele frequencies implemented in the software angsd v.0.923 (Analysis of Next Generation Sequencing Data; [30]), which has been proven to be a more efficient method for low to medium coverage next-generation sequencing (NGS) data than SNPs calling algorithms [30]. We describe below the bioinformatics steps required to apply angsd to RAD-seq data from a non-model organism and the downstream population genetic analyses applied to the filtered datasets.

Assembly pipeline and filtering

Angsd requires a reference sequence to work, which we were lacking. To circumvent this issue, we followed the approach described in [54] by creating an artificial reference sequence from loci previously assembled by stacks under the parameter m = 3, N = 3, M = 3 (based on the results of Mona, Bertorelle, Benazzo, in preparation). To this end, we concatenated the consensus sequences of each locus spaced by a stretch of Ns and then map reads back from individual fastq files using the bwa-mem algorithm with default parameters [55]. Using custom bash scripts coupled with angsd, we then discarded: (i) sites with coverage < 3x (-minIndDepth = 3, corresponding to m in the first assembly performed by stacks) and/or of low quality (based on the per base alignment score, -baq = 1 flag); (ii) low quality bases and poorly aligned reads (-minQ and -minMapQ and -C flags with default values); (iii) SNPs present in the last 5 bp of each locus and SNPs genotyped as heterozygous in 80% or more of the individuals; (iv) loci with more than 5 SNPs that might be the result of paralog RAD loci alignment on the reference. Specific filters were further added according to the downstream analyses performed.

Population structure

A single reference sequence was created for all populations and we retained sites shared by at least 80% of the samples. The PCA was computed with PCAngsd v.0.97 based on genotype likelihood [56]. Admixture was then investigated by running the non-negative matrix factorization algorithm (nmf) implemented in PCAngsd which is based on the same covariance matrix inferred for the PCA. The number of ancestral populations (K) was automatically chosen by PCAngsd to be e + 1, where e is the optimal number of significant principal components depicting population structure, resulting from the Velicier’s minimum average partial test run on the covariance matrix. The sparseness regularization parameter α (used to reduce the noise in low depth NGS data) that best fitted the data was tested between 0 and 100 and it was chosen by comparing the resulting likelihood following [56]. We generated pairwise site allele frequency likelihood files and then computed FST with the realSFS program in angsd [57] using SNPs with a minor allele frequency ≥ 0.05 (-minMaf flag). The significance of each pairwise FST comparison was evaluated with 1000 permutations by randomly allocating individuals to one of the two populations. We finally tested isolation by distance (IBD) using a Mantel test [58] and plotted the relationship between genetic vs. geographic distances.

We applied an approximate Bayesian computation (ABC) approach similar to previous studies [7, 31, 32] in all sampling sites to further investigate the presence of population structure. This approach is particularly helpful in the Atlantic Ocean (AO) where only one locality was sampled (Brazil; BRA population). Briefly, we designed three demographic models (Additional file 1: Fig. S4): (1) NS (No Structure) which represents a panmictic population where \({N}_{mod}\), the modern effective size instantly changes to \({N}_{anc}\), the ancestral effective size, at \({T}_{s}\)(time shift) generations; (2) FIM (Finite Island Meta-population) which represents a finite island meta-population model composed of 100 demes exchanging symmetrically \(Nm\) migrants per generation with each other. All demes were instantaneously colonised, \({T}_{col}\) generations ago, from an ancestral population of size \({N}_{anc}\); (3) SS (Stepping-Stone) which represents a stepping-stone model where the 100 demes are arranged in a two-dimensional grid and where migration is only allowed symmetrically in both directions between the four nearest neighbouring demes. We performed 50,000 coalescent simulations under each model using fastsimcoal v.2.6.0.3 [59] extracting parameters from the prior distributions displayed in Additional file 1: Table S1. Model selection was evaluated by the random forest classification method implemented in the abcRF package in R [60]. We used the SFS, θπ and TD as summary statistics and further added the first two axes of the Linear Discriminant Analysis in the dataset as suggested by [60] to increase the classification method accuracy. The number of trees was chosen by checking the evolution of the out-of-bag error.

Genetic diversity and effective population size variation

We created one reference sequence per population in order to maximise the number of loci assembled. We filtered the sites with missing data by setting the -minInd flag in angsd to the total number of individuals in each population. The filtered dataset was then used to generate a site allele frequency likelihood (saf) file, where genotype likelihoods were computed using the SAMtools method (-GL = 1 flag). The folded site frequency spectrum (SFS) was directly computed from the filtered saf datasets through the realSFS program [57]. Nucleotide diversity (θπ), Watterson’s theta based on segregating sites (θw;[61]) and Tajima’s D (TD; [62]) were computed with custom script from the SFS. Significance of TD was evaluated after 1000 coalescent simulations of a constant population model with size θw. We reconstructed the variation in the effective population size (Ne) 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 divergence 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}_{mo{d}_{IP}}\) and \({N}_{mo{d}_{AO}}\), diverged at \({T}_{div}\) from an ancestral population of effective population size\({N}_{anc}\). Due to the stairwayplot results, we allowed the two modern effective population sizes \({N}_{mo{d}_{IP}}\) and \({N}_{mo{d}_{AO}}\) to change to \({N}_{an{c}_{IP}}\) and \({N}_{an{c}_{AO}}\) following an exponential dynamic in \({T}_{{s}_{IP}}\) and \({T}_{{s}_{AO}}\) years respectively. Migration is defined by two time periods: \({m}_{1}\) representing the migration rate occurring between time 0 until \({T}_{mig}\) and \({m}_{2}\) between time \({T}_{mig}\) until\({T}_{div}\). The migration matrix in each time period is asymmetric: for instance, \({m}_{{1}_{AO/IP}}\) represents the forward migration rate from AO to IP and \({m}_{{1}_{IP/AO}}\) 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 migration only between \({T}_{mig}\) and\({T}_{div}\). The model is defined by eleven parameters, the two \({m}_{1}\) migration rates being removed; (c) IM-rec: same as IM-anc, but with recent migration only occurring between time 0 until\({T}_{mig}\), keeping only the two \({m}_{1}\) migration rates; (d) IM-bsc: the classic model where migration is constant from time 0 until \({T}_{div}\) (i.e. \({\mathrm{m}}_{1}={\mathrm{m}}_{2}=\mathrm{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}_{div},{T}_{{s}_{IP}}\), \({T}_{{s}_{AO}}\)). 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 maximise the number of SNPs shared among the two ocean basins. Similarly, in each basin, we selected the individuals 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 simulations), -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 information criterion comparing the best run of each model [59]. To check the robustness of the model selection procedure 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 simulated under the parameters estimated in the best run of each model, each approximated with 106 coalescent simulations. 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 fastsimcoal and each of these datasets was analysed under the same conditions as the original data (one hundred independent runs for each dataset). Calibrating the molecular clock is crucial to obtain accurate estimates of demographic parameters and historical events, but it is challenging when fossil records and/or orthologous loci from an outgroup are lacking. Here, all demographic inferences were performed using the RAD-seq mutation rate of μ = 1.93 × 10–8 per site and per generation previously used for the tiger shark [7], and the generation time was set to 10 years [24, 64].

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).

References

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

    Article  CAS  Google Scholar 

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

    Article  CAS  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  CAS  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    CAS  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

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

    Article  Google Scholar 

  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.

    Article  CAS  Google Scholar 

  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.

    Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  CAS  Google Scholar 

  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.

    Article  Google Scholar 

  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/ece3.5111.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

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

    Article  CAS  Google Scholar 

  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.

    Article  CAS  Google Scholar 

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

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  CAS  Google Scholar 

  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.

    Article  CAS  Google Scholar 

  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.

    Article  CAS  Google Scholar 

  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.

    Article  Google Scholar 

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

    Article  Google Scholar 

  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.

    Article  CAS  Google Scholar 

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

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  CAS  Google Scholar 

  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.

    Article  Google Scholar 

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

    Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  CAS  Google Scholar 

  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.

    Article  Google Scholar 

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

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  CAS  Google Scholar 

  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.

    Article  CAS  Google Scholar 

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

    Article  CAS  Google Scholar 

  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/rsbl.2019.0741

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

    Article  CAS  Google Scholar 

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

    Article  Google Scholar 

  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.

    Article  CAS  Google Scholar 

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

    CAS  Google Scholar 

  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.

    Article  CAS  Google Scholar 

  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.

    Article  CAS  Google Scholar 

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

    Article  CAS  Google Scholar 

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

    Article  CAS  Google Scholar 

  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.

    Article  CAS  Google Scholar 

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

    Article  Google Scholar 

Download references

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.

Funding

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

Author information

Authors and Affiliations

Authors

Contributions

S.M. and H.M. conceived the project. S.M., P.L 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. P.L., S.M. and H.M. wrote the manuscript with input from all the others. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Stefano Mona.

Ethics 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. TREL2118715S/541 (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.

Additional information

Publisher's Note

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

Supplementary Information

Additional file 1: Figure S1.

Isolation by distance (IBD) plot within the Indo-Pacific. Pairwise genetic distances (FST/(1-FST)) 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 PCAngsd. 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 (\({\mathrm{N}}_{\mathrm{mod}}\)) instantaneously changes to \({\mathrm{N}}_{\mathrm{anc}}\), at time shift \({\mathrm{T}}_{\mathrm{s}}\) generations. FIM (Finite Island Meta-population) represents a finite island meta-population model with 100 demes that have been instantaneously colonised \({\mathrm{T}}_{\mathrm{col}}\) generations ago, from an ancestral population of size \({\mathrm{N}}_{\mathrm{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 \({\mathrm{T}}_{{\mathrm{s}}_{\mathrm{IP}}}\) and \({\mathrm{T}}_{{\mathrm{s}}_{\mathrm{AO}}}\) years ago from a modern effective size (\({\mathrm{N}}_{{\mathrm{mod}}_{\mathrm{IP}}}\) and \({\mathrm{N}}_{{\mathrm{mod}}_{\mathrm{AO}}}\)) to an ancestral effective size (\({\mathrm{N}}_{{\mathrm{anc}}_{\mathrm{IP}}}\) and \({\mathrm{N}}_{{\mathrm{anc}}_{\mathrm{AO}}}\)). The two populations are connected by an asymmetrical number of migrants constant from 0 to \({\mathrm{T}}_{\mathrm{div}}\) (\({\mathrm{Nm}}_{\mathrm{IP}\to \mathrm{AO}}\) and \({\mathrm{Nm}}_{\mathrm{AO}\to \mathrm{IP}}\)) and diverged \({\mathrm{T}}_{\mathrm{div}}\) years ago from an ancestral population of size \({\mathrm{N}}_{\mathrm{anc}}\).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Lesturgie, P., Lainé, H., Suwalski, A. et al. Ecological and biogeographic features shaped the complex evolutionary history of an iconic apex predator (Galeocerdo cuvier). BMC Ecol Evo 22, 147 (2022). https://doi.org/10.1186/s12862-022-02100-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12862-022-02100-y

Keywords

  • Agulhas leakage
  • Coalescent modelling
  • Demographic history
  • Population genomics
  • RAD-seq
  • Tiger shark