Population structure, connectivity, and demographic history of an apex marine predator, the bull shark Carcharhinus leucas

Abstract Knowledge of population structure, connectivity, and effective population size remains limited for many marine apex predators, including the bull shark Carcharhinus leucas. This large‐bodied coastal shark is distributed worldwide in warm temperate and tropical waters, and uses estuaries and rivers as nurseries. As an apex predator, the bull shark likely plays a vital ecological role within marine food webs, but is at risk due to inshore habitat degradation and various fishing pressures. We investigated the bull shark's global population structure and demographic history by analyzing the genetic diversity of 370 individuals from 11 different locations using 25 microsatellite loci and three mitochondrial genes (CR, nd4, and cytb). Both types of markers revealed clustering between sharks from the Western Atlantic and those from the Western Pacific and the Western Indian Ocean, with no contemporary gene flow. Microsatellite data suggested low differentiation between the Western Indian Ocean and the Western Pacific, but substantial differentiation was found using mitochondrial DNA. Integrating information from both types of markers and using Bayesian computation with a random forest procedure (ABC‐RF), this discordance was found to be due to a complete lack of contemporary gene flow. High genetic connectivity was found both within the Western Indian Ocean and within the Western Pacific. In conclusion, these results suggest important structuring of bull shark populations globally with important gene flow occurring along coastlines, highlighting the need for management and conservation plans on regional scales rather than oceanic basin scale.


| INTRODUC TI ON
Delineating populations and their connectivity by gene flow (i.e., effective dispersal) is of primary importance for the conservation and management of endangered and/or exploited species (Begg, Friedland, & Pearce, 1999;Moritz, 1994;Palsbøll, Bérubé, & Allendorf, 2007). In marine species, genetic analyses allow for stocks to be defined, species exploitation status to be assessed, and genetic diversity underlying recruitment potential and species adaptability to be preserved (Begg et al., 1999;Hilborn, Quinn, Schindler, & Rogers, 2003;Palumbi, 2003). Once genetically distinct groups (i.e., populations) that may be managed independently are identified, estimating abundance and the number of individuals effectively exchanged among populations is needed to assess population viability and resilience (Frankham, 2010;Schwartz, Luikart, & Waples, 2007).
Studies of population structure and connectivity are challenging because commonly used direct approaches (mark-recapture, satellite, and acoustic tracking) are often difficult to use for pelagic marine species. This difficulty leads to small sample sizes (Grothues, 2009) and an underestimation of individual movements (Ng, Able, & Grothues, 2007;Thorsteinsson, 2002). Therefore, indirect methods based on the conceptual framework of population genetics have been increasingly used to address ecological and evolutionary questions in such species. First, genetic methods allow the assessment of population structure resulting from evolutionary forces shaping allele frequencies within and among populations (mutation, genetic drift, migration, and selection; Wright, 1931). At neutral loci, while gene flow homogenizes allele frequencies and limits population differentiation, genetic drift promotes population differentiation by randomly fixing alleles (Hartl & Clark, 1997). Second, genetic methods can provide estimates of the effective population size (N e ; Wright, 1931). This parameter represents the size of an idealized Wright-Fisher population affected by genetic drift at the same rate per generation found in the population of interest. Combined with the mutation rate (µ), N e provides an estimate of population genetic diversity (4N e µ for the diploid autosomal part of the genome and N e µ for the haploid mitochondrial genome). N e is also related to the number of breeders per generation (Waples, Antao, & Luikart, 2014) and has been shown to correlate with a population's ability to adapt to environmental changes (Hare et al., 2011). N e has thus been increasingly used in conservation and management to estimate the health status of a population and its ability to recover when depleted (Frankham, Briscoe, & Ballou, 2010).

Marine species characterized by large populations commonly
show weak genetic structuring at neutral loci. In large populations, even a low dispersal rate can lead to weak population genetic structure because the number of migrants is not negligible. Also, genetic drift is limited in these species due to their large population sizes (Bailleul et al., 2018;Gagnaire et al., 2015;Palumbi, 1992). Weak genetic structuring may result from the existence of large isolated populations or, conversely, the existence of one large panmictic population. Identifying which situation is driving population structure can be challenging, but recent developments are providing the necessary analytical resolution. Incorporation of migration into simulation models, combined with new approximate Bayesian computation algorithm relying on random forest (i.e., ABC-RF), allows comparisons and selection of alternative demographic models that best fit the observed dataset (Pudlo et al., 2016;Raynal et al., 2017). ABC-RF provides estimates of the posterior probability of the selected model and the parameters of interests, such as migration rates between populations and effective population size (Pudlo et al., 2016;Raynal et al., 2017). For both model choice and parameter estimates, ABC-RF is more accurate and requires a smaller number of simulated datasets than previous ABC methods (Fraimout et al., 2017;Pudlo et al., 2016;Raynal et al., 2017).
Many large sharks face considerable exploitation, and populations have declined globally in recent decades (Dulvy et al., 2014).
The bull shark Carcharhinus leucas is caught in recreational, subsistence, and targeted commercial fisheries, as well as bycatch computation with a random forest procedure (ABC-RF), this discordance was found to be due to a complete lack of contemporary gene flow. High genetic connectivity was found both within the Western Indian Ocean and within the Western Pacific.
In conclusion, these results suggest important structuring of bull shark populations globally with important gene flow occurring along coastlines, highlighting the need for management and conservation plans on regional scales rather than oceanic basin scale.
Estimates of long-term effective population size of bull sharks vary among studies and locations, but are likely in the order of 100,000 individuals (Karl et al., 2011;Testerman, 2014), which represents potentially greater genetic diversity than other shark species, for which estimates of N e are in the order of 10,000-50,000 individuals (e.g., basking sharks, Hoelzel et al., 2006; sicklefin lemon shark, Schultz et al., 2008). This may suggest that (a) bull shark populations are not severely depleted and/or (b) that fishery pressures are too recent to be detected through genetic analyses (Karl et al., 2011;Testerman, 2014).
To date, few studies have investigated bull shark genetic structure and have relied either on (a) extensive sampling on a global scale using only nuclear markers (Testerman, 2014), or (b) a locally intensive sampling (either restricted to the Atlantic or Northern Australia), using relatively few nuclear and mitochondrial markers (3-5 microsatellites along with 1 or 2 mitochondrial genes; Karl et al., 2011;Tillett et al., 2012). Thus, improving our understanding of bull shark population structuring and connectivity across ocean basins is needed. Combining the information from two types of molecular markers (25 microsatellite loci and three mitochondrial genes [CR, nd4, and cytb]), we analyzed the genetic variation in 370 bull sharks from 11 locations in the Western Indian Ocean, the Western Pacific, and the Western Atlantic, including both continental coasts and oceanic islands ( Figure 1). By including new locations and increasing the number of markers presenting different modes of evolution, our objective was to combine classical population genetic analyses with coalescent-based approximate Bayesian computation approaches (Beaumont, 2010;Csilléry, Blum, Gaggiotti, & François, 2010). This was performed to delineate bull shark populations and assess their demographic history and connectivity, using model selection analyses to refine the evolutionary history of this species. Specifically, the objectives were to: 1. Expand our understanding of the genetic structure previously documented by Testerman (2014) in order to delineate genetic clusters at different scales (e.g., within vs. among ocean basins) that should be managed separately; 2. Decipher whether contemporary migration occurs among defined clusters; and 3. Estimate the effective population sizes of these clusters.

| Sampling
Tissue samples were collected in the Western Indian Ocean (WIO), the Western Pacific (WP), and the Western Atlantic (WA; Figure 1).
In the WIO, samples came from continental coasts and oceanic is-

| Laboratory procedures
Genomic DNA was extracted using Qiagen DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) following manufacturer instructions.
PCRs were performed in a total volume of 25 µl: 1× of MasterMix (Applied Biosystems), 0.3 µM of forward and reverse primers, and 1.6 ng/µl of genomic DNA. The thermocycling program for CR contained an initial denaturing step at 94°C for 5 min, 35 cycles × (94°C for 30 s, 56°C for 30 s, 72°C for 1 min 30 s), and a final extension step at 72°C for 5 min. For cytb, the same program was used, except that the annealing temperature was set to 53°C. For nd4, the annealing temperature was 50°C and the elongation step was 45 s. Amplicons were sequenced directly with primers used for PCR on a capillary sequencer ABI 3730XL (Applied Biosystems) by Genoscreen.

| Genetic diversity analysis
Among the individuals from Madagascar, 12 out of 25 samples were kept for data analyses, because the remaining samples extracted from teeth exhibited high amounts of missing data (more than 50%) due to low-quality DNA.  (Weir & Cockerham, 1984) were estimated using Fstat v.2.9.3.2 (Goudet, 1995). Departure from Hardy-Weinberg equilibrium (HWE) at each microsatellite locus was tested using 5,000 permutations in Fstat v.2.9.3.2 (Goudet, 1995). The mean allelic richness A r and the mean private allelic richness A rp were calculated using a rarefaction method, as implemented in hp-rare v.1.0 (Kalinowski, 2005). This method accounts for differences in sample size by standardizing A r and A rp values across sampled locations by resampling the lowest number of genotypes available (i.e., 12 haploid gene copies or six diploid genotypes in Rodrigues Island) in each location.
Detection of partitioning schemes within the concatenated sequence CR-nd4-cytb and of substitution models was performed using partitionFinDer v.2.1.1 (Guindon et al., 2010;Lanfear, Frandsen, Wright, Senfeld, & Calcott, 2017). We used Beast v.1.8.4 (Drummond, Suchard, Xie, & Rambaut, 2012) to reconstruct phylogenetic relationships and infer divergence times on the mitochondrial concatenated sequence CR-nd4-cytb. Bayesian Markov chain Monte Carlo (MCMC) analyses were performed assuming a HKY85 + I model of substitution as the latter was shown to best fit the data. The rate of variation among sites was modeled with a discrete gamma distribution with four rate categories. We assumed an uncorrelated lognormal relaxed clock to account for rate variation among lineages. To minimize prior assumptions about demographic history, we adopted an extended Bayesian skyline plot (EBSP) approach in order to integrate data over different demographic histories. Trees were calibrated using two methods. First, an analysis was performed adding a sequence of S. lewini (mitochondrion available in GenBank; accession number JX827259), and the tree was calibrated using the divergence date between Carcharhinus and Sphyrna genera, 38 millions years ago (Mya), estimated from fossil data (Maisey, 1984). Second, the tree was calibrated using the closure of the Isthmus of Panama as the divergence time of bull shark populations from the Western Atlantic and the Indo-Pacific, 3.1-3.5 Mya (Coates, Collins, Aubry, & Berggren, 2004;Coates et al., 1992). For each analysis, a normal prior distribution was set for the calibrated node (mean ± SD: 38 ± 7 and 3.5 ± 0.4, respectively). Evolutionary model parameters were then estimated, with samples drawn from the posterior every 10 5 MCMC steps over a total of 10 8 steps from five independent runs. The first 10 7 steps were discarded as burn-in.

| Population genetic structure
Two complementary clustering methods were used to investigate population structure in the bull shark. First, Bayesian clustering analyses were performed using structure v.2.3.4 (Falush, Stephens, & Pritchard, 2003;Pritchard, Stephens, & Donnelly, 2000). For any given number of clusters (K) between 1 and 10, individual assignment probabilities to each cluster were determined so as to minimize departures from HWE within clusters and maximize LD between them. Two analyses were performed, with and without the LOCPRIOR model, which uses prior sampling location information in the Bayesian clustering to detect genetic population structure that might be weaker (Hubisz, Falush, Stephens, & Pritchard, 2009).
Conditions were set to 10 6 chain length after a burn-in of 5 × 10 5 , and 10 chains were run for each K assuming correlated allele frequencies and the admixture model. For a given K, distinct modes were identified, and for each mode and each individual, the assignment probabilities to each cluster were averaged using CluMpak (Kopelman, Mayzel, Jakobsson, Rosenberg, & Mayrose, 2015). Second, a discriminant analysis of principal components (DAPC Jombart, Devillard, & Balloux, 2010), which does not rely on HWE or LD contrary to structure, was performed to check whether similar clustering patterns were identified. This method summarizes the genetic variation of the microsatellite allele frequencies using a principal component analysis as a prior step to a discriminant analysis and defines clusters such as to minimize variations within them and maximize differentiation between them. DAPC was applied using the adegenet package (Jombart, 2008) for R (R Core Team, 2017). Methods traditionally used to detect the most likely number of clusters according to the analysis performed (Structure and DAPC) might provide different outputs for the same dataset. To cope with these inconsistencies, we chose to consider the highest number of clusters and the individual assignments that were retrieved by both analyses. Moreover, in a hierarchical approach, these analyses were repeated on each cluster found separately. Commonly, using Structure and DAPC, when the finest level of structuring is reached, adding a supplementary cluster leads to inconclusive assignments with individuals assigned to several clusters in the same proportions.
Analyses of molecular variance (AMOVAs; Cockerham, 1969Cockerham, , 1973 were performed to estimate the genetic variation due to the partitioning in clusters (identified with the TCS haplotype network for the mitochondrial data and with structure and DAPC for microsatellite data), the variation within clusters among sampling locations, and the variation within sampling locations. AMOVAs were performed with arlequin v.3.5.1.2 (Excoffier & Lischer, 2010), and significance of fixation indices was tested using a nonparametric approach with 10,000 permutations (Excoffier, Smouse, & Quattro, 1992).

| Neutrality tests
To test for departures from a constant population size (Ramos-Onsins & Rozas, 2000), the summary statistics Tajima's D (Tajima, 1989) and Fu's F S (Fu, 1997) were estimated from the concatenated mitochondrial dataset with Arlequin v.3.5.1.2 (Excoffier & Lischer, 2010), with significance tested implementing 10 5 simulated samples. It was therefore not included in the ABC-RF analysis. Pooling individuals from different sampling locations, even with nonsignificant pairwise differentiation values, may bias results (Lombaert et al., 2014). Hence, the two regions were represented by the sampling location with the highest number of individuals, that is Reunion Island (RUN) for WIO and Eastern Australia (Clarence River, AUS1) for WP. Four demographic scenarios were built (Figure 2), all of them starting with an ancestral population from which both observed populations diverged. Scenarios then differed as to the occurrence of migration during divergence. Scenario 1 assumed constant recurrent migration from the split to present. In Scenario 2, the split was followed by a period of recurrent migration, itself followed by a period of isolation. In Scenario 3, populations diverged in isolation.

Demographic scenarios
Finally, Scenario 4 assumed that populations first went through a period of isolation before engaging in a period with recurrent migration. In all scenarios, recurrent migration was bidirectional but not necessarily symmetric.

Model choice
For each scenario, we simulated 200,000 microsatellite and mitochondrial datasets using FastsiMcoal (Laval & Excoffier, 2004). To account for both types of markers having different sample sizes, we applied a two-step procedure (bash scripts available upon request).
Microsatellite datasets were first simulated with parameters drawn in the prior distributions described in Appendix S2 (Tables S2.1 and   S2.2). The mitochondrial datasets were subsequently simulated using the same historical parameters (divergence times, starting time, and ending time of the migration period) as for microsatellites, but with different sample sizes and, importantly, different demographic and genetic parameters (effective sizes, migration rates, and mutation rates). We thus estimated different migration rates and effective sizes for microsatellite and mtDNA. Because of the lack of knowledge on effective sizes and historical divergence of bull shark populations, broad parameter ranges were chosen. Simulated datasets were described using 19 summary statistics (Appendix S2) related to the genetic polymorphism of both types of loci using ArlsuMstat (Excoffier & Lischer, 2010). For both markers, we computed the mean number of alleles over loci K and the mean of Nei's gene diversity H for each population and the pairwise F ST between populations.
For microsatellite markers only, the mean over loci of the modified Garza-Williamson index MGW were computed for each population and the mean delta mu-square δµ 2 (square difference in mean microsatellite allele length between pairs of populations) between F I G U R E 2 Graphical representations of the four scenarios depicting possible divergence histories for each pair of bull shark populations: FLO-RUN, FLO-AUS1, and RUN-AUS1. The time was measured backward in generations before present. In black, is represented the ancestral population of effective population size N anc ; in dark gray, population 1 of effective population size N 1 and in light gray, population 2 of effective population size N 2 . Double arrows represent bidirectional migration events. t 2 , time of divergence; t 1 , start and end of the isolation period for Scenario 2 and Scenario 4, respectively

Parameter estimations
Parameter values were subsequently inferred using ABC random forests as developed by Raynal et al. (2017), using 100,000 datasets simulated under the best scenario. To test the performance of the method in estimating parameters, we used 1,000 pseudo-observed datasets on which the estimation procedure was applied to measure the precision of the estimation procedure. From these values, the 95% confidence interval (CI) and the normalized mean square error NMSE were computed. Parameter inference analyses were replicated two times to ensure consistency of ABC-RF analyses.

| Genetic diversity analysis
Null alleles were detected for several loci in several sampling loca-

| Genetic clustering
Structure clustering analysis suggested that the genetic structure is best explained by two clusters. For the microsatellite dataset without the LOCPRIOR model, a clear clustering was observed at

| Genetic differentiation
AMOVAs were conducted with the previously obtained clusters

| Neutrality tests
Considering the concatenated mitochondrial dataset, no evidence of any historical population expansions or contractions was found with tests of selective neutrality, either by considering all locations separately or by grouping them in the clusters identified (all Tajima's D and Fu's F S not significantly different from zero; Appendix S8).

| Bayesian analyses using both microsatellite and mtDNA data
The PCAs on the space of the summary statistics and the analysis of the distribution of each summary statistics revealed that all scenarios could produce simulated datasets mirroring the observed dataset. On the PCAs of the summary statistics, the point representing the observed dataset fell within the cloud of points representing the simulated ones (Appendix S9). Also, most often the observed summary statistics fell well within the distribution obtained from the simulations (Appendix S10).
In all 10 replicates, Scenario 3 had the highest percentage of votes with 38.98% ± 0.97, while Scenario 1 the lowest with 5.73% ± 0.73 (Table 3). Performing Tukey's post hoc tests, we confirmed that Scenario 3 had a significantly higher percentage of votes compared to all others (all p < .001), while no significant differences were found between Scenario 2 and Scenario 4 (p = .15).
Parameter values were thus estimated using data simulated under Scenario 3 only.
Using 1,000 pseudo-observed datasets, we found that for effective population sizes, the estimation procedure had very low bias and good precision over the whole prior range with low NMSE values, ranging from 0.02 to 0.03 for contemporary populations and 0.14 to 0.16 for the ancestral unsampled population (Table 4 and Appendix S11). Using these estimations, effective population sizes in number of genes estimated from the microsatellite data ranged from RUN. Other parameters were less well resolved (Appendix S12), and values will not be interpreted.

| D ISCUSS I ON
Using a combination of markers following different models of evolution and appropriate inference methods may help reach a better understanding of genetic structure and connectivity. Here, hierarchical sampling (inter-and intra-ocean basins) and the use of both mtDNA sequences and microsatellite markers allowed us to test for the existence of migration between populations and to estimate effective

| An ancient divergence between the Atlantic and the Western Indian Ocean/Western Pacific
Both mtDNA sequences and microsatellite markers showed high differentiation, suggesting a complete absence of gene flow between the Western Atlantic and the Western Indian Ocean/Western Pacific since their divergence. This result is congruent with previous research on bull shark using microsatellites, which identified three isolated genetic clusters, one in Indo-Australia, one in Fiji, and one in the Atlantic Ocean (Testerman, 2014). This possibly long-dating genetic divergence may have enabled the emergence of biological differences between Atlantic Ocean bull shark populations on one side, and those of the Indian and Pacific Oceans (Indian/Pacific Oceans) on the other. In the Indian/Pacific Oceans, individuals are larger, both in terms of maximum length  and size at maturity (Cliff & Dudley, 1991) than those from the Gulf of Mexico (Branstetter & Stiles, 1987;Cruz-Martinez, Chiappa-Carrara, & Arenas-Fuentes, 2005).
Divergence times were inferred based on a molecular clock estimate and should thus be regarded as qualitative indicators, rather than precise values. The use of the divergence between Sphyrna and Carcharhinus genera, or of the Isthmus of Panama closure as the divergence date between the Atlantic and the Indian/Pacific bull shark populations, yielded mutation rates similar to those observed in other shark species using several different fossil records (Duncan et al., 2006;Gubili et al., 2014;Karl, Castro, & Garla, 2012;Schultz et al., 2008). Using two different calibration dates, we estimated the ter rivers and lakes (Heupel & Simpfendorfer, 2008;Thorson, 1976).
The lack of samples from the Eastern Pacific did not allow us to test this hypothesis or the presence of any relationships between animals from the Eastern and the Western Pacific. Yet, populations from these two regions might be genetically structured because of the East Pacific Barrier, in place since 65 Mya (Grigg & Hey, 1992). This biogeographical barrier is characterized by depths over 5,000 m over a wide oceanic distance (~7,000 km), limiting longitudinal dispersal across the Pacific Ocean (Briggs, 1995). Nevertheless, some gene flow among these three regions may have occurred after the formation of the East Pacific Barrier, via the southern tip of Africa, before the formation of the Benguela Upwelling System.
The Benguela Upwelling System may be more constraining than the closure of the Isthmus of Panama for the bull shark, which is more sensitive to cold temperatures than species for which some gene flow after the formation of this current has been highlighted (e.g., tiger shark Galeocerdo cuvier [Bernard et al., 2016], dusky shark Carcharhinus obscurus [Benavides et al., 2011], or scalloped hammerhead shark Sphyrna lewini [Duncan et al., 2006]). Bull sharks remain in warmer waters, favoring temperatures of 24-26°C (Smoothey et al., 2016), and found less frequently in waters less than 18°C (Brunnschweiler et al., 2010;Lea et al., 2015;Matich & Heithaus, 2012 Note: N e is expressed in number of genes; N e (sat), effective population size estimated using microsatellite data; N e (seq), effective population size estimated using mtDNA; OOB-MSE, out-of-bag mean square error; NMSE, normalized mean square error; NMAE, normalized mean absolute error; CI, 95% confidence interval.
disrupted the migratory behavior of bull sharks and led to the divergence of the Atlantic and Indian Ocean populations. Additional samples from the Eastern Pacific and the Southern Atlantic (both Eastern and Western) are needed to further investigate the worldwide phylogeography of the bull shark.

| Negligible gene flow between Western Indian Ocean and Western Pacific
We observed a high differentiation at mtDNA sequences and a low differentiation at microsatellite markers between the Western Indian and the Western Pacific Oceans, a finding that had not been identified in previous studies. For example, Testerman (2014) Bernard et al., 2016;Karl et al., 2011;Pardini et al., 2001;Portnoy et al., 2015). But sex-biased dispersal is not the only possible cause of a higher differentiation in mtDNA sequences as compared to microsatellite markers. In addition to their difference in modes of evolution, nonpanmictic mating systems may affect differentially the levels of differentiation at both types of markers. ABC random forest procedure, which is regarded as one of the most precise Bayesian methods to identify demographic histories (Fraimout et al., 2017;Pudlo et al., 2016;Raynal et al., 2017), offers a mean to formally test for the evolutionary forces underlying genetic population structure, including migration regimes. To do so and to account for sex-biased dispersal, we independently estimated migration rates and effective sizes for both types of markers. Analyses revealed that the scenario with no gene flow between the Western Indian Ocean and the Western Pacific populations since their isolation best explained the observed data. Indeed, while scenarios with migration were designed to allow sex-biased dispersal, they were chosen significantly less to explain the observed data than the scenario with no migration over 10 independent replicate analyses. This may reflect either an absence of gene flow or dispersal events that are rare enough not to be detected. For populations of large sizes (N e > 10 3 ), rare effective dispersal events may be sufficient to homogenize allelic frequencies, leading to F ST estimates nonsignificantly different from zero (in the order of 10 −3 ) while maintaining high mitochondrial differentiation (Hauser & Carvalho, 2008;Mariani & Bekkevold, 2014).
To increase juvenile survival, females may exhibit high fidelity to their breeding areas and nurseries, which are typically good foraging areas and offer protection from large predators (Branstetter, 1990;Castro, 1993;Heupel, Carlson, & Simpfendorfer, 2007;Springer, 1967). These breeding sites are sometimes the same as the natal places of females, as these latter represent suitable habitats for parturition (Heupel et al., 2007;Hueter et al., 2005).
Female philopatry to nursery areas has notably been demonstrated in the lemon shark Negaprion brevirostris in the Bahamas by reconstructing parental genotypes (microsatellites) through sampling juveniles in specific nurseries over several decades: Some females returned to the nursery to give birth, sometimes 14 to 17 years after being born (Feldheim et al., 2013). In contrast, males may exhibit roaming behaviors and undertake migration, possibly to avoid inbreeding depression, and demographic and environmental stochasticity, especially in polygynous systems (Henry, Coulon, & Travis, 2016), as may occur for the bull shark (A. Pirog, personal communication). It is thus possible that female philopatry to nursery areas also occurs in the bull shark as hypothesized in the Western Atlantic (Karl et al., 2011) and in Australia (Tillett et al., 2012). Furthermore, no direct evidence of bull sharks moving between the Western Indian Ocean and the Western Pacific has been documented using satellite tracking or conventional tagging.
While it may be due to relatively small sample sizes, it may also illustrate the absence, or at least extremely low occurrence, of bull shark migration across the Indian and Pacific Oceans. Pleistocene, each lasting approximately 100,000 years, followed by shorter interglacial periods of about 10,000 years (Dawson, 1992;Martinson et al., 1987), fluctuations in sea levels were as great as 100 m during this time period (Shackleton, 1987). These within Australian waters (Dudgeon, Broderick, & Ovenden, 2009;Ovenden, Kashiwagi, Broderick, Giles, & Salini, 2009). It is possible that the deep waters of the Timor Trench (2,000-3,000 m) and the strong Indonesian through-flow current along the Makassar and Lombok Straits induced the genetic subdivisions observed between Indonesian and Australian waters (Dudgeon et al., 2012(Dudgeon et al., , 2009Ovenden et al., 2009), and thus limits gene flow between the Indian and Pacific Oceans.  . As such, long-distance migration of adult bull sharks may genetically link ecosystems within these regions. Each movement study also highlights the fidelity of bull sharks to specific sites at discrete times, as shown in Reunion Island , in New Caledonia (Werry & Clua, 2013), in Australia , in Fiji, and in the Bahamas (Brunnschweiler & Baensch, 2011;Brunnschweiler et al., 2010

| Effective population sizes
Changes in population size were not detected with neutrality tests performed with mitochondrial data. Estimates of effective population sizes (N e ) from the microsatellite dataset were ca. 3,000-4,000 individuals for the Western Indian Ocean and the Western Pacific.
We obtained much lower estimates from the mitochondrial dataset, with biparentally inherited ones (autosomal markers), and a deviation from that expectation may reflect sex-biased dispersal patterns, social organization, or specific mating systems. Chesser and Baker (1996) showed that in panmictic populations and in systems with single paternity, the effective size of maternally and paternally inherited markers was one-half of that of biparentally inherited markers and that social structure, sex-biased dispersal, or different mating systems usually lower the effective size of autosomal markers while lowering or uppering maternally and paternally inherited markers.
The bull shark has recently been shown to be a polyandrous species (Pirog et al., 2015). Sugg and Chesser (1994) showed that multiple paternity increases the effective sizes of diploid genes. However, because all the progeny will receive the maternally inherited genes from the female regardless the sire, multiple paternity should not affect the dynamics of the genes that are transmitted by the female (Chesser & Baker, 1996). Estimates of N e inferred using mtDNA may thus be more accurate than those estimated using microsatellite data.
Using the mismatch distribution of the mitochondrial control region, Tillett et al. (2012)  proportional to the number of loci (Felsenstein, 2006;Pluzhnikov & Donnelly, 1996). It may also be due to the scale of the region studied, as our estimates were obtained using samples from one locality to represent an entire region. Our mitochondrial estimates were nevertheless lower than those previously inferred.
Estimates of effective population size using genetic markers are increasingly used for fisheries stock assessments (Ovenden et al., 2016). It has been postulated that an N e of at least 500 individuals is needed for a population to adapt to environmental changes  although others estimate that at least 5,000 breeding individuals may be required (Lande, 1995). Avoiding deleterious allele accumulation may require an N e above 1,000 individuals Palstra & Ruzzante, 2008) and inbreeding depression may occur if N e falls below 50 individuals (Frankham et al., 2010). Our estimates (3,000-4,000 with microsatellite data; 380-1,800 with mtDNA) are nearly in the same range as the basking shark Cetorhinus maximus (i.e., 8,200;Hoelzel et al., 2006), but lower than estimates for the lemon shark N. brevirostris (26,000 to 52,000 in the Atlantic) and the sicklefin lemon shark Negaprion acutidens (34,000 to 52,000 in the Western Pacific; Schultz et al., 2008), and much lower than for the tope shark, Galeorhinus galeus (138,000;Chabot & Allen, 2009). All of these species are considered either (a) globally Vulnerable on the IUCN Red List or (b) subjected to a loss of genetic diversity due to a bottleneck (e.g., basking sharks). This may therefore be the case for bull sharks, especially if taking into account mtDNA estimates, and populations may even be depleted. Obtaining more precise population estimates requires greater knowledge of the reproductive biology of the bull shark, notably the number of individuals that successfully reproduce in a generation (or reproductive cycle), the age at maturity, and the mating system (Ovenden et al., 2016).

| CON CLUS ION
Here, we highlight several key findings about the global population structure of bull sharks that will inform management and conservation issues: 1. The genetic isolation between bull shark populations from the Western Atlantic and from the Western Indian Ocean/Western Pacific implies that the Western Atlantic populations should be managed separately.
2. Low gene flow, and maybe even complete isolation, has also been evidenced between bull shark populations from the Western Indian Ocean and the Western Pacific, despite a low nuclear differentiation. It implies that these populations should also be man-