Results
The number of assembled loci with no missing data ranged from 60,812 to 117,916 (Table 1). Mean coverage (and standard deviation) per sample was 9.02 (± 2.62), 7.93 (± 0.48), 8.39 (± 0.81) for G. cuvier, C. limbatus and C. amblyrhynchos respectively. Genetic diversity was very heterogeneous, with Ɵπ varying from 0.0004 (C. melanopterus ) to 0.0021 (C. amblyrhynchos ) per site (Table 1). Tajima’s D was significantly positive for C. melanopterus (TD = 0.691, p<0.001) and C. limbatus (TD = 0.43, p<0.001), significantly negative for C. amblyrhynchos (TD = -0.23, p<0.001) and not significant forG. cuvier (TD = -0.03, p=0.059) (Table 1).
We compared the models NS, FIM, and SST (Figure 1) in the four species by means of an ABC-RF algorithm and estimated demographic parameters for the most supported model. After checking for the evolution of the out-of-bag error of the RF, model selection and parameter estimation were computed using respectively 500 and 1,000 trees in each species. We found that NS had the higher posterior probability (p =0.84) forG. cuvier , suggesting that an unstructured model fit the data better than structured models for this species (Table 1). In contrast, demographic histories of the three other species were best described by SST, with a posterior probability ranging from 0.53 to 0.88 (Table 1 and S2). The estimated median number of migrants per generationNm was 1.8 (95% CI: 0.7-3.0) for C. melanopterus , 6.6 (95% CI: 1.5-15.4) for C. limbatus, and 11.5 (95% CI: 3.0-22.0) for C. amblyrhynchos (Figure 2, Table 1). The posterior distribution of Nm strongly differed from the prior distribution and showed a clear unimodal peak with small credible intervals, and low mean square error (SME) and mean root square error (SMRE) in all three species (Figure 2, Table S3), suggesting that these estimates are highly reliable. Conversely, both Tcol and Nanc had larger SME and SMRE errors in all species (Table S3), but it was only in C. melanopterus where posterior and prior distribution could not be distinguished (Figure 2). Tcol has a clear unimodal distribution in C. amblyrhynchos but a more disperse one (and with wider credible intervals) in C. limbatus (Figure 2, Table 1).
The stairwayplots showed a nearly similar dynamic for C. amblyrhynchos and C. limbatus , characterized by a strong ancestral expansion (Figure 3). During a short time interval, theNe increased rapidly about ~2 times, reaching a value between 35,000 and 42,000 individuals. When approaching T=0, both species underwent a bottleneck but of distinct strength: Nedropped down to ~12,500 for C. amblyrhynchos but almost to 0 for C. limbatus . This is consistent with the shape of the normalized SFS, which clearly shows a stronger deficit in low frequency variants for C. limbatus compared to C. amblyrhynchos (Figure 3). Similarly to C. limbatus , C. melanopterus experienced a recent 10-fold population collapse around 20,000 years B.P. starting from a long term constant Ne . However,C. melanopterus showed no signature of ancestral expansion, consistent to the results obtained by Maisano Delser et al. (2019) usingabc-skyline method. Finally, G. cuvier displayed an ancestral expansion around 100,000 years B.P. with Ne reaching ~12,000 before dropping to ~3000 at T ~1,600 years B.P. Remarkably, the ancestral expansion retrieved by the stairwayplot (Figure 3) for both C. amblyrhynchos and C. limbatus is overlapped with the posterior distribution of Tcol estimated by the SST model (Table 1). This analogy holds too for C. melanopterus, where Tcol could not be properly estimated under the structured model (we obtained a flat posterior distribution, Figure 2) and there was no signature of ancestral expansion in the stairwayplot (Figure 3).
In order to study the shape of the SFS under the considered models and to understand how Ne changes over time assuming a panmictic population, we run sets of coalescent simulations. The first set of coalescent simulations was run under FIM and SST only (Table 2 and S4) to check if simulated data could reproduce the pattern of genetic variability (both θ estimators and TD ) observed for C. melanoptetus , C.limbatus , and C. amblyrynchos . The simulated θ values (excluding the equilibrium model) ranged between 0.001 and 0.003 per site, in line with the observed values (Table 1, 2 and S4). TD follows a U-shaped distribution for each Nmvalue as a function of Tcol , being more positive at recentTcol and at equilibrium and less positive (or negative for higherNm ) at intermediate values. Therefore, species demography withNm ~10 (and higher) and Tcol within 15k and 50k generations B.P. will have negative TD values. In contrast, species with lower Nm and very recent or very ancientTcol will have positive TD . This matches strikingly theTD observed for the three shark species and their estimated demographic parameters under SST (Table 1). We plot the normalized SFS and the stairwayplot for all scenarios presented in Table 2 (Figure 4, 5, 6, S1 and S2). First, we note that none of our scenarios, even those at equilibrium and with no variation in Nm through time, showed a normalized SFS compatible with a constant size population (Figure 4, 5, 6, S1 and S2). The normalized SFS and the reconstructedstairwayplot depend generally on the interaction betweenNm and Tcol with a dynamic strikingly similar to TD(which is indeed a summary of the SFS). For Nm=1 we observed the signature of a recent decrease in Ne for all scenarios and independently of Tcol (Figure 4). The normalized SFS showed consistently a strong deficit of low frequency variants, typical of a demographic bottleneck and in agreement with the positive TD(Figure 4 and Table 1). Furthermore, the stairwayplot could never detect the ancestral expansion for any Tcol . These results are consistent with the normalized SFS and the stairwayplot observed in C. melanopterus , which had a median Nm of ~1 and for which it was not possible to estimateTcol (Figure 2, 3 and Table 1). For growing Nm, the interplay with Tcol becomes more complex. A general result is that, once again, all scenarios were characterized by a recent decrease of Ne when looking at the stairwayplot and a deficit of singletons compared to the other low frequency classes when looking at the normalized SFS (Figure 5, 6 and S1). However, a strong signature of ancestral expansion appeared for Nm >10 andTcol between 15k and 50k generations B.P., mirroring the results of TD for which most of these scenarios displayed a negative value. Remarkably, the stairwayplot retrieved the ancestral expansion only slightly overestimating the simulated Tcol (Figure 5, 6, S1). The signature of the expansion becomes weaker for growingTcol (up to the equilibrium) even for larger Nm and the magnitude of the deficit of low frequency variants increases (Figure 5, 6, S1 and S2). These results are consistent with the stairwayplotand the normalized SFS observed in C. amblyrynchos and C. limbatus , with the former showing larger Nm and more recentTcol than the latter. This can explain why C. limbatusshowed a stronger decrease in Ne in recent times compared toC. amblyrhynchos (Figure 2, Table 1). Similar results were obtained for FIM (Figure S3).
We further tested whether a change in connectivity occurred in any of the three species best described by SST. To this end, we compared SST vs SST-CH model (Figure 1) by means of the same ABC-RF model selection framework previously adopted. The two models cannot be clearly distinguished in any of the three species since: i ) they showed similar posterior probability (~0.50); ii ) the prior error rates are large ~0.40 (Table S5);iii ) posterior distributions of Nm before and afterTch are wide and largely overlapping (Table S6); iv ) the normalized SFS closest to the observed data retrieved under the two models are very similar (Figure S4). These results can be explained either because SST and SST-CH produce a similar pattern in genetic diversity, or because we lacked enough data given our framework. To explore this issue, we ran a second set of coalescent simulations focusing on the consequences of a recent change in connectivity on the observed SFS and the reconstructed stairwayplot (Table S7 and S8). The decrease in connectivity was simulated by reducing eitherm (the migration rate per generation) or N (the effective population size of each deme) (Figure 7 and 8). As expected, we found a signature of recent population decline in all simulated scenarios, with its intensity only slightly affected by the change in Nm (Figure 7, 8, S5 and S6). However, the drop in N (Figure 8) had larger effect compared to the drop in m (Figure 7) on both the normalized SFS and the expansion time estimated by thestairwayplot . In scenarios with 100x reduction in N , thestairwayplot could not retrieve the ancestral expansion even for large long-term Nm (Figure 8). FIM-CH models displayed a behaviour similar to SST-CH models but more pronounced (Figure S7, S8, S9 and S10, Table S9 and S10). While at Tch =10 a decrease inNm slightly affected the SFS and the reconstructedstairwayplot , the consequence of the change in connectivity are more substantial at Tch =50, with a stronger deficit in singletons and a more pronounced recent decline in Ne particularly in scenarios with a 100-fold reduction of N (Figure S8, S10).