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