Introduction

The Coral Triangle (sensu Carpenter et al. 2011) and the Great Barrier Reef (GBR) are among the most species-rich marine regions in the world (McWilliam et al. 2018; Mouillot et al. 2014; Renema et al. 2008). The species assemblages of such biodiversity hotspots have been compared to less diverse regions to infer how processes like continental drift and climatic fluctuations have generated and maintained biodiversity (Keith et al. 2013; Leprieur et al. 2016; Pellissier et al. 2014; Renema et al. 2008). Yet, biodiversity begins at the population level: patterns of genetic diversity within and amongst populations are linked to their demographic histories (Fang et al. 2021, 2020), which are largely shaped by past geological and climatic events (Kotlík et al. 2018; Wei et al. 2020). The geographic distribution of genetic diversity—i.e., the level of biodiversity at which evolution takes place—can therefore shed light on the generation and maintenance of diverse populations in biodiversity hotspots.

Part of the Arafura Shelf connecting Australia and New Guinea (Fig. 1a) most recently re-submerged only ~7000 years ago to form the Torres Strait (Woodroffe et al. 2000). This reconnected the GBR (in the Pacific Ocean) with the western Coral Triangle and north-western Australia (in the Indian Ocean) for the first time in over 100,000 years. Before this reconnection (during previous Pleistocene climate fluctuations), western Australia and the Coral Triangle were refugia for coral reef-associated species (Cowman and Bellwood 2013; Greenstein and Pandolfi 2008; Pellissier et al. 2014). On the other side of the Torres Strait—following millennia of changing sea levels and intense sedimentation—the GBR only began growing into its current form and size around 10,000 years ago (Carter and Johnson 1986; Marshall and Davies 1984; Webster et al. 2018). Therefore, the most recent submergence of the Torres Strait could have aided the recolonization of the GBR by marine fauna from western Australian and Coral Triangle refugia and increased the genetic diversity of species with populations that persisted in both regions.

Fig. 1: Bathymetric map of our Indo-Pacific study region with Carcharhinus amblyrhynchos sampling location information.
figure 1

Sampling locations are indicated by red circles, and the number of specimens collected at each location is indicated in parentheses. Grey shading denotes land area above mean sea level; light blue denotes area no deeper than 200 m; darker blue shading denotes area deeper than 200 m. Inset (a) is a bathymetric profile of southern Indonesia, New Guinea, and northern Australia during the last glacial maximum. Area that is currently below sea level but was exposed (i.e., land) during the last glacial maximum is coloured green. Inset (b) shows the wider Indian and Pacific Oceans and inset (c) shows the New Caledonia area in more detail.

Full size image

The role played by the cyclical emergence and submergence of the Torres Strait in shaping genetic diversity of populations on either side of the strait remains unclear. Even when dispersal-related life history traits are accounted for across a variety of organisms, the level of genetic differentiation amongst populations of species present on opposite sides of the Torres Strait varies widely (Gaither and Rocha 2013; Keyse et al. 2018; Liggins et al. 2016; Mirams et al. 2011). With up to six-metre differences in the sea level on each side of the Torres Strait due to the different tidal regimes in adjacent basins (Wolanski et al. 1988), the corresponding eddies and currents can be both strong and highly variable. Furthermore, net water movement through the Torres Strait both switches directions and varies threefold in strength seasonally (Margvelashvili et al. 2008; Saint-Cast and Condie 2006). All available genetic and hydrological evidence therefore suggests that dispersal through the strait is highly stochastic.

Large marine predators such as sharks can undertake active migrations. Therefore, despite the irregular oceanographic conditions, sharks can likely move through the Torres Strait region more consistently than organisms with planktonic larval stages. Catch and acoustic tagging studies have demonstrated that several species of coastal and coral reef associated sharks have specific habitat requirements such that they tend to disperse mainly along continental shelves (Chin et al. 2012; Dudgeon et al. 2013; Espinoza et al. 2015a, 2015b, 2014). Accordingly, genetic connectivity amongst coastal and coral reef-associated shark species in this region is facilitated by shallow shelfal areas (like the Torres Strait), while hindered by deep oceanic expanses (Dudgeon et al. 2009; Momigliano et al. 2017; Ovenden et al. 2009; Vignaud et al. 2014). When the shallow Torres Strait became land, however, shark populations on both sides of Australia may have diverged, although this hypothesis has not yet been explicitly tested.

Many previous coastal and reef shark studies have relied mainly on mtDNA and microsatellites (Boissin et al. 2019; Dudgeon et al. 2009; Ovenden et al. 2009; Portnoy et al. 2016; Vignaud et al. 2014) which often may not provide the necessary statistical power to resolve detailed evolutionary history (although see Portnoy et al. (2016) for a counterexample). While these studies may provide valuable insights into the genetic diversity and connectedness of reef shark populations, recent studies have highlighted the improved resolution of genomic sequencing (e.g., Kraft et al. 2020). Recent genomic studies of the blacktip reef shark (Carcharhinus melanopterus) demonstrated how studies with fewer markers may miss the signature of an ancient range expansion (Maisano Delser et al. 2018, 2016). Further study of populations using genome-wide markers has become a priority in the management and conservation of reef sharks (Domingues et al. 2018; Dudgeon et al. 2012; Heupel et al. 2019), and the more intensive evolutionary analyses they allow could provide insight into the origin of species as well as their current distribution and patterns of genetic diversity.

In this study, we used genome-wide single nucleotide polymorphism (SNP) data from populations of the endangered grey reef shark (Carcharhinus amblyrhynchos) (Simpfendorfer et al. 2020) to shed light on both their historical demographics and contemporary effective population sizes (Ne). In particular, we aimed to reconstruct the demographic history of this species, as well as determine its centre of origin and the potential legacy that the Torres Strait land bridge has had on patterns of genetic diversity. If C. amblyrhynchos has undergone a range expansion similar to the one inferred for C. melanopterus (Maisano Delser et al. 2018, 2016), demographic history models could detect an initial expansion across multiple lineages of the species, bottlenecks (decreases in Ne) due to founder effects, and possibly a spatial gradient of genetic diversity. If other major demographic events have occurred since the inception of the species—such as populations on both sides of the Torres Strait intermittently exchanging migrants—then demographic history models may infer more recent changes in Ne and migration rates (and possibly erase genetic signatures of earlier events). Estimates of coalescent Ne from demographic history models were compared with estimates of contemporary Ne from a common conservation genetics method to assess both the ancient and recent past of this endangered species.

Materials and methods

Sampling and genotyping

The genetic data in this study come from 513 C. amblyrhynchos samples collected across 17 locations in the Indian and Pacific Oceans (Fig. 1b, c). The capture and tissue sampling procedures for all individuals have been described in previous studies (Boussarie 2019; Momigliano et al. 2017, 2015a). The sequencing and SNP discovery was performed by Diversity Arrays Technology Pty. Ltd. (Canberra, Australia) using the standard DArTseq protocol (Sansaloni et al. 2011). After calling and filtering, the dataset contained 7952 loci with a total of 14,935 SNPs across ~551,040 base pairs (bp). For more details about the bioinformatics and filtering, see the Supplementary Information.

Genetic diversity

Species range expansions are expected to leave transient spatial patterns of genetic diversity. More recently colonized areas tend to have lower genetic diversity than the region of origin due to the sequential bottlenecks that characterise range expansion (Harpending and Rogers 2000; Mona 2017). This transient spatial pattern of decreasing genetic diversity is expected to persist until migration-drift equilibrium is reached. A simple test to identify a centre of origin was proposed by Ramachandran et al. (2005): genetic diversity is expected to be negatively correlated with distance from a putative centre of origin, with the most likely centre of origin being the location with the strongest negative correlation coefficient. Combining all samples from a given location, we computed the number of segregating sites (S), pairwise differences between sequences (π) and Tajima’s D directly from its site frequency spectrum (SFS) using the software package moments (Jouganous et al. 2017). We report values of π per site (i.e., divided by the total number of sequenced bases; Supplementary Information). A geographical lattice of 22,880 points (every 0.5 degrees between 70.25°E–179.75°E and 25.75°N–25.75°S) was constructed and then filtered down to include only points located in the ocean shallower than 200 metres since C. amblyrhynchos is very rarely observed away from continental shelves (with the exception of coral reefs in isolated seamounts) (Last and Stevens 2009). For each of the included 2028 points, the shortest distance by sea to each sampling location was calculated using the r packages ‘gDistance’ (van Etten 2018) and ‘marmap’ (Pante and Simon-Bouhet 2013). The distance from the point to each of these sites was then rank-correlated using Spearman’s rho (ρ) with the π at those sites using r version 3.6 (R Core Team 2019). A second, finer-scale run was performed in a narrower geographic area considering only 15 sampling locations (omitting the two oceanic archipelagos in the Indian Ocean—Cocos (Keeling) and Chagos) where an initial lattice of 125,584 points (every 0.125 degrees between 91.625°E–174.875°E and 20.875°N–25.875°S) was filtered the same way to produce 7521 points for the correlation analysis. Maps and figures were produced using r as described in the Supplementary Information.

Demographic history models

Studies on the demographic history of populations are often constrained by the necessity of using predefined models to test alternative demographic scenarios. Instead of using pre-defined models to detect up to one or two changes in the size of individual populations (e.g., Maisano Delser et al. 2016, 2018), we used stairway plot v2 (Liu and Fu 2015), a model-flexible method that calculates the composite likelihood of the observed SFS under several possible size change scenarios across the entire history of a population.

stairway plot infers the most likely timing and magnitude of changes in Ne at various epochs, with the transition between epochs corresponding to coalescent events in the population. stairway plot initially calculates the population-scaled mutation rate θ, from which Ne can be derived if μ (the germline mutation rate) is known (θ = 4•Ne•μ). Previous studies on the congener C. melanopterus used a uniform prior distribution of 8.05–8.54 • 10–9 for μ in an Approximate Bayesian Computation analyses, which they estimated using samples of Carcharhinus galapagensis from opposite sides of the Central American Isthmus (Maisano Delser et al. 2018, 2016). Here, we used the midpoint of this range (8.295 • 10–9). After producing a randomly-thinned, minor allele (folded) SFS for each population (Supplementary Information), stairway plot was run using the recommended (nseq − 2) ÷ 4 break points (where nseq is the number of haplotypes sampled) and 200 replicate runs for each population (except Cocos (Keeling) which had too few sequenced specimens to be considered in this analysis).

We further evaluated the joint demographic history of the Misool and North GBR populations using the method implemented in moments. This program is largely based on dadi (Gutenkunst et al. 2009), but uses direct computation of the frequency spectrum without solving the diffusion system. After creating the two-dimensional site frequency spectrum (2D-SFS) of the Misool and North GBR populations (Supplementary Information), we tested five types of gene flow models: the strict isolation model (SI), the isolation with migration model (IM), the two-epoch model (2EP), the secondary contact model (SC) and the ancient migration model (AM) (Sousa and Hey 2013). Combined with each type of gene flow scenario, we also modelled different plausible historical population growth scenarios, including ancestral expansions (AE) and bottlenecks (B). Visual representations of these models can be found in Fig. S1 and Table S1 contains a list of all model names and their parameters. The SI model assumes a split at time T1, with no migration. The IM model includes continuous, asymmetric migration following the population split at T1, while the 2EP model allows for migration rates to change at time T2. The SC model has a period of strict isolation starting at the time of split (T1) followed by gene flow starting at time T2, while the AM model assumes migration at the initial stages of divergence but strict contemporary isolation. If gene flow has not occurred between Misool and the North GBR since they were colonized (or stopped shortly after colonization), then the SI and AM models should fit the data well. On the other hand, if there is ongoing gene flow between these two regions, the IM, 2EP, or SC models should fit the data better. If this gene flow has been constant since both regions were colonized, the IM model should not fit the data any worse than the SC or 2EP models. If gene flow was interrupted or diminished by the presence of the Torres Strait land bridge, the SC or 2EP models should fit the data best. Since overlooked demographic events can in some cases bias demographic inference (Momigliano et al. 2021), and since both stairway plot and Tajima’s D values suggest the possibility of an ancestral population expansion, we also considered an additional version of each main model that included a population expansion at time TAE prior to the split (AE models). Since stairway plot also suggested a possible recent expansion in the North GBR (potentially corresponding with its regrowth following the last glacial maximum), we tested variations of each model that included exponential growth following a bottleneck starting at the time of split T1 (B models). We further tested variations of each model that included both ancestral expansion and exponential growth in the North GBR. For each basic scenario, we therefore tested four variants (e.g., for the IM model, we tested the following variants: IM, IMB, IMAE, IMAEB).

Model parameters were optimized in 10 independent optimization routines, following the pipeline developed by Portik et al. (2017) for dadi, adapted to moments by Momigliano et al. (2021), and further modified by Le Moan et al. (2021). Each independent routine was composed of four rounds of subsequent parameter optimization. In the first round a Broyden–Fletcher–Goldfarb–Shanno algorithm optimization (function “optimize.log”; max-iter = 10) was run for 30 sets of threefold randomly-perturbed parameters. The parameter values with the highest composite likelihoods were then chosen as starting parameters from a second round of optimization, this time run for 20 sets of two-fold randomly-perturbed parameters and increasing max-iter to 20. The same procedure applied to round two was used in round three and four, this time using two-fold randomly-perturbed parameters and one-fold randomly perturbed parameters. Convergence was checked using the best replicate of each of the 10 independent optimizations. Model choice was performed as described in Rougeux et al. (2017) and Momigliano et al. (2021). We calculated the Akaike information criterion (AIC) and ΔAIC for each model from the replicate with the highest likelihood. Then we calculated Akaike weight of evidence (wAIC) as outlined in Rougeux et al. (2017). Finally, we checked the fit of the best scoring models by plotting the residuals between the empirical 2D-SFS and the 2D-SFS produced by the models.

Contemporary effective population size

Contemporary Ne was estimated for each sampled site using the ldne linkage disequilibrium method (Waples and Do 2008) implemented in NeEstimator 2.1 (Do et al. 2014). This is a commonly used method to estimate recent Ne in population genetic studies for conservation purposes (e.g., Reid-Anderson et al. 2019). Simulation studies have shown that it can reliably detect changes in Ne from restriction-site associated sequencing data (Nunziata and Weisrock 2018). The mean linkage disequilibrium (r2) between all pairs of (putatively) unlinked loci in n sampled individuals from a given cohort is compared to the linkage disequilibrium that would be expected from sampling n individuals from an idealized (infinite) population (in which there is no actual linkage disequilibrium other than what is introduced by finite sampling bias). The difference between the observed and expected linkage disequilibrium indicates the level of genetic drift occurring in the sampled cohort, which is inversely proportional to its Ne. An estimate of infinite Ne is possible if the observed linkage disequilibrium between pairs of loci in the sampled cohort is indistinguishable from the expected linkage disequilibrium in an equivalent sample size from an idealized population. This method of estimating Ne is therefore most useful for small population sizes (tens to hundreds), and the upper 95% confidence interval limit on estimates can be infinite for estimates in the hundreds or low thousands (Lonsinger et al. 2018; Moran et al. 2019).

An adjustment of this method’s final estimate of Ne should be applied when a large number of loci are used. This depends on the size of the species’ genome, which has not yet been determined for C. amblyrhynchos. However, diploid chromosome numbers of 74–86 have been estimated for other species in the family Carcharhinidae (Gregory 2018), so an intermediate value of 40 haploid chromosomes was used for this correction. The recommended pcrit (minor allele frequency cut-off) of 0.05 was used to avoid upwardly biasing estimates of Ne (Waples et al. 2016). Since confounding factors such as overlapping generations (Waples et al. 2014) cannot be controlled for in this dataset (see Reid-Anderson et al. (2019) for additional considerations), the estimates provided here should not be taken as valid, standalone assessments of Ne, but rather as calculations used to compare the inferred Ne of each sampling location relative to each other and to the other methods of estimating Ne in this study.

Results

Genetic diversity

Overall, there were only slight differences in nucleotide diversity between sampling locations. While Entrecasteaux and several other New Caledonian sampling locations had the highest number of segregating sites (Table S2), Scott Reef and the three locations closest to it—Misool, Rowley Shoals, and Ningaloo Reef—all had the highest nucleotide diversity with π > 0.00147 (Table S3). The Indian Ocean archipelago locations (Chagos and Cocos (Keeling)) had both the lowest number of segregating sites and lowest nucleotide diversity, while Matthew had the next lowest nucleotide diversity. Tajima’s D was negative at all locations.

In the diversity-distance correlation analysis, the most likely centres of origin were located slightly east of Misool in the Cenderawasih Bay (Fig. 2). The ρ of −0.87 for these locations is indicative of a strong negative correlation in genetic diversity with distance away from the Coral Triangle, particularly towards the east. The second run without Indian Ocean archipelago sampling locations confirmed this: a steep gradient in ρ was observed through the Torres Strait (Fig. S2), and the correlation between nucleotide diversity and degrees East (ρ = −0.839) was only slightly weaker than the strongest correlation found at any point (ρ = −0.857).

Fig. 2: Spatial distribution of the Spearman’s ρ (rho) coefficient between genetic diversity and distance by sea from all 17 Carcharhinus amblyrhynchos sampling locations (black crosses) for 2028 points throughout the Indo-Pacific region shallower than 200 m.
figure 2

The lowest values of rho (indicating the most likely centre of origin) were −0.87 and located east of Misool in the Cenderawasih Bay of Western New Guinea (white X’s).

Full size image

Demographic history models

A strong signal of expansion 40,000–90,000 generations ago was inferred by stairway plot for all locations except Chagos (Fig. 3). Immediately after this ancestral expansion, the Ne of the three locations in western Australia were inferred to be at least 20% larger than all others. More recent expansions were inferred at the Misool and North GBR populations. The steady decrease after the initial expansion inferred for Matthew is consistent with what has previously been interpreted as a founder effect in studies on a closely-related reef shark C. melanopterus (Maisano Delser et al. 2018). Less extreme, more recently-initiated declines were inferred for Walpole, Grand Astrolabe, and Ningaloo Reef.

Fig. 3: Carcharhinus amblyrhynchos demographic histories inferred by stairway plot.
figure 3

Sampled locations are grouped into three different regions, and those with trends distinct from others in their region are labelled. Cocos (Keeling) has been removed due to insufficient sample size. The dashed vertical lines at 40,000 and 90,000 generations ago on the x-axis span the period during which an expansion was inferred across nearly all sampled locations.

Full size image

Two-population models that included contemporary migration (IM, SC and 2EP models) always had lower AIC scores (better fit) than models with no contemporary gene flow (Fig. 4). The SI models performed worst, and they are not represented in Fig. 4 since their very high AIC would obscure the differences between the rest of the models (Table S1). In general, models that included growth in the North GBR performed better than models that did not, and models that included an ancestral expansion performed better than models that assumed an ancestral population at drift-equilibrium at the time of split between the North GBR and Misool (Fig. 4). The model with the lowest AIC was SCB (secondary contact coinciding with exponential growth in the GBR), however Akaike weights for the SCB, IMAEB, SCAE, and SCAEB models were fairly similar.

Fig. 4: Modelling results from moments for the Misool and North GBR populations of Carcharhinus amblyrhynchos.
figure 4

SI = strict isolation, IM = isolation with migration, twoep = two periods with different levels of migration, SC = secondary contact, AM = ancestral migration (see Fig. S2 for visual representations of relevant models). a AIC from the best 3 optimizations for all tested models, with the exclusion of SI gene flow (GF) scenarios (as SI models had such high AIC that their inclusion in the graph would have obscured differences among better fitting models). The x-axis groups models by whether or not they include an ancestral expansion (AE), exponential growth in the North GBR following a bottleneck (B), both (AEB), or neither (basic); (b) Weight of evidence (wAIC) for each model, in decreasing order on the x-axis. Four models had wAIC above 0.1: SCB, IMAEB, SCAE and SCAEB. See Table S1 for more detailed results and all parameter estimates for each model.

Full size image

The four best models (i.e., those with the lowest AIC values: SCB, IMAEB, SCAE, and SCAEB) all show a good fit to the data, with small and normally-distributed residuals (Fig. 5). However, they each represent a fairly different demographic history. The second, third, and fourth best models (respectively IMAEB, SCAE and SCAEB) inferred an ancestral population expansion between 79,000 and 87,000 generations ago, in line with the results from stairway plot. These three models (IMAEB, SCAE and SCAEB) indicate the Misool and North GBR groups diverged between 12,200 and 13,500 generations ago, and the SC models indicate that secondary contact was established again 1400–2400 generations ago. The SCB model suggests a much older divergence time of ~145,000 generations ago—even older than the oldest coalescence event inferred within any single population by stairway plot (~120,000 generations ago).

Fig. 5: Graphical representation of the empirical 2D-SFS for Carcharhinus amblyrhynchos in Misool and the North GBR, as well as the inferred parameters and residuals from the four best models in moments.
figure 5

Time parameters are given as thousands of generations ago (kga).

Full size image

Contemporary effective population size

Effective population sizes were estimated to be largest in the most northern and western locations sampled in New Caledonia (Fig. 6). Chesterfield had the largest Ne at 5529, followed by adjacent sampling sites Northern Lagoon and Entrecasteaux at 4974 and 4401, respectively (Table S3). The Ne at each of these three sites was more than six times larger than all others east of the Torres Strait, and even the lower limit of their jackknife 95% confidence intervals overlapped only with the Ne estimated for Ningaloo Reef (3885), which was twice as large as any other location west of the Torres Strait (Fig. 6, Table S3). The only other locations with Ne > 1000 were Scott Reef and Misool which, along with Ningaloo Reef, were the three locations with the highest nucleotide diversity (Fig. 6, Table S3). Matthew had the second lowest Ne at 161 and was the only location for which a finite upper 95% confidence interval limit was estimated (547), while the South GBR had the lowest Ne at only 70 (Fig. 6, Table S3).

Fig. 6: Effective population size estimates at each Carcharhinus amblyrhynchos sampling location using the linkage-disequilibrium (ldne) method implemented in NeEstimator.
figure 6

Locations are listed in increasing order of longitude on x-axis and Ne is represented on the y-axis (on a log2 scale). Shaded area indicates jackknife 95% confidence intervals (all upper limits were infinite except for Matthew).

Full size image

Discussion

This genomic analysis of 513 C. amblyrhynchos individuals from 17 locations indicates that a range expansion and glacial cycles may have affected the distribution of diversity within the species. However, these ancient events do not appear to have had much of an effect on modern populations, as their estimated contemporary effective population sizes do not correspond with levels of genetic diversity or coalescent history.

Grey reef shark genetic diversity decreased with distance away from the Coral Triangle and north-western Australia. A similar genetic diversity gradient was observed in C. melanopterus (Maisano Delser et al. 2018). While these gradients are typically inferred to indicate a species’ centre of origin, in this case, they could also indicate a centre of overlap or accumulation of diversity from the adjacent Indian and Pacific Oceans. Maisano Delser et al. (2018) found further support for the centre of origin hypothesis in C. melanopterus by calculating directionality indices (Peter and Slatkin 2013), demonstrating that significantly higher proportions of shared derived alleles were found outside of north-western Australia. This region could be both a centre of origin and a centre of accumulation for C. amblyrhynchos as these hypotheses are not mutually exclusive (Bowen et al. 2013; Briggs and Bowen 2013; Cowman and Bellwood 2013; Gaither and Rocha 2013). Further evidence that C. amblyrhynchos originated from the Coral Triangle comes from the detection of an initial expansion across almost all sampled locations 40,000–90,000 generations ago by stairway plot, followed by larger post-expansion Ne west of the Torres Strait (in the Coral Triangle and north-western Australia). Tajima’s D was negative across all locations, which also supports the hypothesis of an expansion. However, this last result should be interpreted with caution as our sample sizes are not homogeneous and Tajima’s D can be affected by sample size (Subramanian 2016). The early and steady decline inferred by stairway plot for Matthew matches the founder effects observed in many Pacific island populations of C. melanopterus using a similar ABC-skyline approach (Maisano Delser et al. 2018). Indeed, our modelling results are consistent with previous studies of a C. melanopterus range expansion: the C. melanopterus expansion inferred to have occurred around 59,000 generations ago (Maisano Delser et al. 2016) is in the middle of the expansion range inferred in this study, and the highest Nm (the product of Ne and migration rate) values through time were also found in populations from northern and western Australia (Maisano Delser et al. 2018). The two closest relatives of C. melanopterus (Carcharhinus cautus and Carcharhinus fitzroyensis) are geographically restricted to Western Australia, Queensland, and New Guinea (Lyle 1987). Accordingly, Maisano Delser et al. (2018) suggest that the cautus-melanopterus-fitzroyensis clade may have originated in this region. The results of our study indicate that a third, more distantly related species of the same genus shares a similar demographic history. This is consistent with species assemblage studies showing that north-western Australia and the Coral Triangle were refugia for species-level diversity during the Pleistocene (Greenstein and Pandolfi 2008; Ludt and Rocha 2015; Pellissier et al. 2014), and that the Coral Triangle is a centre of speciation (Bowen et al. 2013; Briggs 1999; Cowman and Bellwood 2013).

While the geographic distribution of genetic diversity in C. amblyrhynchos is similar to that of C. melanopterus (strongly correlated with distance away from the Coral Triangle), the differences between locations were less extreme in C. amblyrhynchos. Maisano Delser et al. (2018) found genetic diversity differences of greater than 310% between continental and island C. melanopterus populations, and as much as a 38% difference between continental Australian populations. Furthermore, putative founder effects (localized bottlenecks) were detected across almost all C. melanopterus populations outside of Australia (Maisano Delser et al. 2018). Only a few bottlenecks were observed in this study, while the greatest disparity in nucleotide diversity between locations was 208% overall and only 8.7% between continental Australian locations. Different sequencing and associated bioinformatic approaches (for restriction site-associated sequencing data here versus target gene capture in Maisano Delser et al. (2016, 2018)) may make our results less comparable. However, the more uniform genetic diversity in C. amblyrhynchos relative to C. melanopterus could also be explained by the two species’ differing dispersal potential. Migration increases local genetic diversity and decreases structuring across the range of a species, leading towards migration-drift equilibrium during a range expansion (Mona et al. 2014). Although both C. amblyrhynchos and C. melanopterus generally exhibit site fidelity (Bonnin et al. 2021; Chin et al. 2013; Espinoza et al. 2015a, 2015b; Mourier et al. 2012; Papastamatiou et al. 2009; Speed et al. 2011), C. amblyrhynchos has a larger average dispersal distance (Dwyer et al. 2020), activity space (Speed et al. 2016), and utilizes offshore resources more frequently (McCauley et al. 2012; Papastamatiou et al. 2018). Furthermore, only C. amblyrhynchos has been observed to undertake repeated migrations of several hundred kilometres (Bonnin et al. 2019)—sometimes more than 900 kilometres (White et al. 2017)—across oceanic expanses.

In single population models, inferred changes in Ne may also be due to changes in migration, which, as discussed above, brings additional diversity into the population via gene flow. While Maisano Delser et al. (2016, 2018) accounted for this by using metapopulation structure in their single population models, we followed up single population modelling with joint population modelling to explicitly consider different historical gene flow scenarios between locations on each side of the Torres Strait. Three of the best four models were of isolation followed by secondary contact between Misool and the North GBR, implying not only that migrant exchange between these two locations is ongoing, but also that they have previously been isolated from each other. However, the inferred ages of demographic changes were quite different between models. The SCB model had the lowest AIC, yet yielded a very old divergence time of ~143,000 generations ago. This seems unrealistic given that it is ten times older than the divergence times inferred in the next three best models and is even older than all coalescent events inferred by stairway plot. Furthermore, it was the only model of the best four to not include an ancestral population expansion (as indicated by stairway plot and negative Tajima’s D values at all locations). Momigliano et al. (2021) showed that when true demographic scenarios include both an ancestral expansion and subsequent growth in a newly diverged population, models including only new growth (like the SCB model) can erroneously produce a very good fit and misestimate parameters. Thus, given the multiple lines of evidence in this study that suggest an ancestral expansion has occurred, we believe it is likely that the SCB model is providing an oversimplified scenario.

To compare the timing of geological events (in years) to the timing of demographic events inferred by coalescent models (in generations, when a mutation rate is given), an estimated generation time for the sampled organism is required. Maisano Delser et al. (2016, 2018) used seven years as the generation time for C. melanopterus. If the same value was used as the generation time in our best two-population models (excluding SCB), the divergence of the Misool and North GBR populations would be inferred to have occurred during the last glacial period, and the SC models would suggest that secondary contact occurred around when the Torres Strait resubmerged 7000 years ago (Woodroffe et al. 2000). However, this figure of 7 years used to calibrate previous models was actually female maturation time, not generation time. Generation time also factors in the timing of subsequent pups throughout the female’s lifespan and can effectively be considered the average difference in age of a parent and all of its offspring. Therefore, the generation time of iteroparous organisms (that reproduce more than once during their lifespan) is somewhere between their age of sexual maturity and death (typically closer to the latter age as older, larger individuals are generally more fecund). In C. amblyrhynchos, female maturity was observed after a minimum of 11 years, and as fecundity increased with age, generation time was estimated to be 16.4 years (Robbins 2006). Using this generation time to calibrate our models, all events are inferred to have occurred before the last glacial maximum (more than 21,000 years ago). Note that because of uncertainties in estimating generation times and mutation rates, these dates should be regarded as approximations.

Regarding the grey reef shark populations on either side of the Torres Strait, we can confidently draw two conclusions. The first is that there is likely ongoing gene flow between them (as indicated by the poor performance of the ancient migration and strict isolation models). The second is that, as shown by stairway plot, the sizes of these two populations, the level of gene flow between them, or both, have increased since the initial expansion of the species. Four of the top five joint history models and stairway plot indicate that the North GBR population expanded since its divergence from the Misool population. This may correspond with the expansion of the GBR habitat following the last glacial maximum (Carter and Johnson 1986; Marshall and Davies 1984; Webster et al. 2018). However, the Ne increase in Misool inferred by stairway plot instead seems more likely to be due to resumed gene flow through the Torres Strait. If it were instead due to an expansion in habitat availability similar to that of the GBR (or increased gene flow from a different sampled region), we would have expected to observe similar increases in Ne in western Australian populations around the same time (given their similar bathymetry and history as refugia). The possibility of migrant exchange between the Coral Triangle and the Coral Sea via northern New Guinea and the Bismarck Sea could confound these modelling results. However, the most consistent genetic break identified amongst a variety of organisms with larval stages in the Coral Triangle was in northern New Guinea: southeast of Halmahera Eddy near the Cenderawasih Bay (Carpenter et al. 2011). This implies that connectivity through this region is low for organisms with larval dispersal, yet a similar analysis of organisms without larval dispersal such as sharks has not been completed. To fill this knowledge gap, further sampling in northern New Guinea, the Bismarck Sea, and the Solomon Islands could reveal unknown lineages with different demographic histories, as well as how gene flow differs between larval and non-larval organisms near barriers such as the Halmahara Eddy and the Torres Strait.

The large estimates and infinite 95% confidence interval upper bounds for most linkage-disequilibrium estimates of Ne imply that most populations are large, and that their sizes cannot be confidently estimated or statistically distinguished from one another. Despite this, the largest estimated population sizes were upwards of 60 times larger than the smallest, arguably warranting further consideration. Despite the higher genetic diversity and ancient Ne inferred via coalescent modelling for locations west of the Torres Strait, linkage-disequilibrium estimates of Ne show that the largest populations are currently found around remote reefs and atolls in the north-western part of the New Caledonian archipelago. These locations, due to their isolation from centres of human population, are also likely among the least anthropogenically-impacted worldwide (Cinner et al. 2018; Januchowski‐Hartley et al. 2020; Juhel et al. 2019, 2018; MacNeil et al. 2020) and are fully protected by strict no-entry rules within the Coral Sea Nature Park. This supports previous findings that reef shark abundance and diversity is higher in isolated coral reef ecosystems (Juhel et al. 2019, 2018; MacNeil et al. 2020; Sandin et al. 2008). On the other end of the spectrum, the Ne estimated for the South GBR is considerably smaller than all other estimates in this study. This is likely the result of several factors, possibly linked to the habitat structure from which these individuals were sampled (Momigliano et al. 2017, 2015a). Reef habitat is less contiguous towards the southern end of the GBR with greater gaps between reefs, unlike further north where most reefs are highly connected and likely perceived as contiguous habitat by sharks across a much wider range (Espinoza et al. 2015a; Momigliano et al. 2015b). This site is also at the extreme southern end of the C. amblyrhynchos distribution on the GBR, further limiting potential population size. Habitat area is likely a major limiting factor in the oceanic islands of Chagos, Cocos (Keeling), Walpole, and Matthew as well—which in addition to having low diversity and inferred Ne in all analyses performed here—are also furthest from the inferred point of origin for C. amblyrhynchos.

Most alarmingly, studies using multiple assessment methods have detected reef shark population declines on the GBR (Hisano et al. 2011; Robbins et al. 2006), with catch per unit effort of reef sharks in this region having decreased by 82% between 1962 and 2016 (Roff et al. 2018). This figure is similar to the much lower abundance and diversity of reef sharks in southern New Caledonia (near the capital city Nouméa) relative to remote north-western New Caledonia (Juhel et al. 2019, 2018). These apparent localized declines (detected using the linkage disequilibrium method) have not resulted in a loss of genetic diversity though as π in the south GBR and Southern Lagoon is greater than the π observed in Northern Lagoon (which is estimated to have a contemporary Ne of nearly 5000 in this study). This is consistent with recent findings that human-induced extinction can exceed the pace at which losses of genetic diversity are detectable—even with genome-wide sequencing (Roycroft et al. 2021). Furthermore, our coalescent analyses (stairway plot) also did not detect any recent declines, inferring the history of the south GBR and Southern Lagoon populations to be similar to those in north-western New Caledonia (with the largest contemporary Ne observed in this study). Simulation studies have shown that the linkage disequilibrium method of estimating Ne detects declines (in populations of 250–1000) two times sooner than coalescent methods (i.e., following only half as many generations since the initiation of a decline (Nunziata and Weisrock 2018)). Thus, it seems likely that the reef shark population declines previously documented in this area (Hisano et al. 2011; Robbins et al. 2006; Roff et al. 2018) are not yet detectable using coalescent methods, but that future sampling efforts may detect these declines. Overall, these contrasting results show that the various methods employed in this study, which are routinely used to infer population declines and range expansions, are not interchangeable. Researchers and managers looking to assess the current status as well as the histories of populations should therefore consider the appropriate methods to achieve their goals.

Conclusions

While reef shark populations are declining worldwide, the scientific community has called for increased use of genomic tools to assess their status (Domingues et al. 2018; Dudgeon et al. 2012; Heupel et al. 2019). By analysing the diversity and Ne of C. amblyrhynchos across 17 locations using genome-wide SNP data, we identified the Coral Triangle and north-western Australia as the most likely region from which this species originated. Although joint demographic history modelling of populations on both sides of the Torres Strait yielded inconclusive results, single population history models indicate that the intermittent land barrier here may have obstructed gene flow. We also reinforce the conclusions of previous studies, this time using genetic methods, that the largest current reef shark populations are found in isolated, protected coral reef ecosystems (Juhel et al. 2019, 2018; MacNeil et al. 2020; Robbins et al. 2006; Sandin et al. 2008). This study demonstrates which methods can be used to disentangle long term Ne from contemporary population sizes and illustrates how genetic diversity in marine biodiversity hotspots is shaped by ancient demographic processes.