A globally threatened shark, Carcharias taurus, shows no population decline in South Africa

Knowledge about the demographic histories of natural populations helps to evaluate their conservation status, and potential impacts of natural and anthropogenic pressures. In particular, estimates of effective population size obtained through molecular data can provide useful information to guide management decisions for vulnerable populations. The spotted ragged-tooth shark, Carcharias taurus (also known as the sandtiger or grey nurse shark), is widely distributed in warm-temperate and subtropical waters, but has suffered severe population declines across much of its range as a result of overexploitation. Here, we used multilocus genotype data to investigate the demographic history of the South African C. taurus population. Using approximate Bayesian computation and likelihood-based importance sampling, we found that the population underwent a historical range expansion that may have been linked to climatic changes during the late Pleistocene. There was no evidence for a recent anthropogenic decline. Together with census data suggesting a stable population, these results support the idea that fishing pressure and other threats have so far not been detrimental to the local C. taurus population. The results reported here indicate that South Africa could possibly harbour the last remaining, relatively pristine population of this widespread but vulnerable top predator.

The size of a population is a determining factor in how evolutionary and demographic processes affect its longterm survival. In conservation biology, small populations are of concern as they are typically more susceptible to demographic, environmental and stochastic genetic events [1][2][3] . The loss of genetic diversity and accumulation of deleterious alleles can lead to a decreased ability to adapt to changing environmental conditions, ultimately impairing long-term persistence and resulting in an elevated extinction risk 4,5 . The concept of effective population size (N e ) was originally introduced by Wright 6 and can be interpreted as the size of an ideal population that experiences the same amount of genetic drift as the population under consideration. 'Ideal' in this context refers to a constant population of equal sex ratio, in which all individuals mate randomly, and are equally likely to produce offspring 6 . In wild populations, several factors, such as skewed sex ratio and non-random variance in reproductive success and mating, influence the amount of genetic diversity lost, and N e is often reduced in comparison to actual, or census population size (N c ) [7][8][9] .
Several studies have attempted to establish the ratio of N e and N c , with the expectation that the estimation of one parameter would allow the inference of the other 10 . Although this is an active area of research, ratios remain uncertain, as they are influenced by a variety of factors which interact in reducing N e . Reported average N e /N c ratios across a range of species from different taxa range from approximately 0.1 to 0.5 [11][12][13] . N e /N c ratios in teleost fish with high fecundity, high juvenile mortality and sweepstake recruitment (type III survivorship), can be as low as 10 -514, 15 (but see Waples et al. 2018 16 ). In contrast, in elasmobranchs (sharks, rays and skates), N e closely approximates N c 17,18 . This is in line with theoretical predictions based on their life history traits such as longevity, late maturity and low fecundity rate 19  www.nature.com/scientificreports/ as shown by Chevolot et al. 2006 20 , where N e /N c ratios for the thornback ray Raja clavata were estimated to be very low (between 9 × 10 -5 and 1.8 × 10 -3 ). Effective population size N e can be estimated from genetic and demographic data, or a combination thereof 16,21,22 . As it is usually difficult to collect enough data on demographic parameters in wild populations, genetic approaches have become widely popular, and are particularly useful when studying rare or elusive organisms, such as elasmobranchs 23 . In light of current levels of global fisheries exploitation, estimates of contemporary N e and demographic changes through time can provide valuable information to guide effective management and conservation of vulnerable elasmobranch species and populations [24][25][26] . As a coastal species, the ragged-tooth shark Carcharias taurus is not targeted by large off-shore fishing industries, but it is nonetheless threatened by small-scale fishing operations such as artisanal and recreational fisheries 27,28 . Population declines have been observed over much of the species' distribution, and the species is listed as Critically Endangered by IUCN Red List criteria in South America 29 and East Australia 30 , and is thought to be locally extinct in the Mediterranean 31 . A large number of sharks has been caught in the beach meshing program along the Australian East Coast, which is believed to be a significant driver in the dramatic decline experienced by C. taurus in that region 32,33 . In South Africa, the KwaZulu-Natal beach meshing program has also been identified as a having a negative impact on this species, due to its very low intrinsic rate of population increase 34 . Despite the potential impact of these threats, no decline in catch rate or body size of C. taurus was detected in the competitive shore-fishery 35 or beach meshing catches 34 . Population trends estimated from mark-recapture data also remained stable over a 20-year study period 36 . In the same study, it was estimated that an average of 6,800 juvenile and 16,700 adult C. taurus inhabit the South African coast. However, this modelling approach was applied to a shark species for the first time, and some uncertainties remain around the parameters of the model 36 . To gain a better understanding of whether population trends inferred from catch data present an accurate picture of the health of this species, we used genetic data to explore the demographic history of the South African C. taurus population. In particular, we aimed to confirm that the population has not yet suffered any drastic declines from anthropogenic pressures. This represents the first genetic assessment of its kind for this population and may help to direct future management decisions.

Methods
Sampling and data generation. Genetic data for Carcharias taurus were generated using opportunistically collected samples originating from both juvenile and adult sharks. Research fishing trips were conducted at various nursery sites along the South African South Coast from 2003 to 2009. Teams of 2-4 research anglers used conventional fishing tackle to catch C. taurus from the shore. Following capture, sharks were restrained in a specially constructed plastic stretcher while taking size measurements and fin clip samples from the first dorsal fin using a sterile scalpel. To minimize stress, sampling and release were completed in less than five minutes. Samples were stored in 80% ethanol until processing. Additional tissue samples from the East Coast were provided by the KwaZulu-Natal Sharks Board beach meshing program. More detail on the net installations, service and operations of the program can be found in Dudley & Simpfendorfer 2006 34 . Capture, sampling and release of sharks was conducted following guidelines endorsed by the Port Elizabeth Museums animal ethics committee, and experimental protocols applied in this study were approved by the Ethics Committee of the Faculty of Science, University of Johannesburg. DNA extraction and amplification of 12 microsatellite markers was conducted following protocols described in Klein et al. 2019 37 for a total of 189 samples. MICROCHECKER version 2.2.3 38 was used to examine genotypes for stuttering, null alleles and large allele drop-out. Linkage disequilibrium (LD) between pairs of loci and significant deviations from Hardy-Weinberg-Equilibrium (HWE) were tested in GENEPOP version 4 39 . For both tests, the dememorization number was set to 1,000, the number of batches to 500 and the number of iterations per batch to 10,000. P values were corrected for multiple comparisons with the B-Y method 40 using the p.adjust function in R 41 . Because inferences of past population size changes can be strongly affected by genetic structure 42,43 , it was first confirmed that the species is genetically homogeneous in South African waters. For this purpose, population pairwise F ST and D est values were computed in GENALEX version 6.5 44 . P values were generated under 1,000 permutations and adjusted using B-Y correction. Further, to visualize genetic distances among individuals, a principal coordinate analysis (PCoA) was conducted in GENALEX. Subsequent analyses were aimed at exploring the demographic history of C. taurus in South Africa using two complementary methods that evaluate temporal changes in N e .

Approximate Bayesian computation. Approximate Bayesian computation (ABC) as implemented in
DIYABC version 2.0 45 was used to compare several pre-defined models of historical change in N e . ABC is a Bayesian approach in which coalescent-based simulations are compared to observed data using summary statistics to identify models that fit the observed data most closely. Comparison of different demographic models was preceded by exploratory analyses to establish suitable priors. To this end, a scenario with a single, genetically homogeneous population was simulated (1 × 10 6 data sets) to determine posterior distributions of N e and mutation rate (µ). A full description of this analysis is provided in the Supplementary Information ('DIYABC: preliminary analysis and determination of suitable priors').
Following this preliminary analysis, three competing demographic scenarios were constructed: a stable population (scenario 1), a population expansion (scenario 2), and a decline in population size (scenario 3) (Fig. 1). Priors for N e and µ were refined based on posterior values resulting from the previous analysis. To explore the approximate time when changes in population size occurred, we performed two analyses using the same three models, but with different priors for the timing of population size change (t). First, a wide, uniform prior was chosen to encompass recent events such as the arrival of European settlers and the start of commercial fisheries in South Africa, as well as natural historical events that preceded the onset of European settlement (analysis 1). www.nature.com/scientificreports/ To ensure that the genetic signal of a recent, anthropogenic population decline was not masked by more significant changes in population size earlier in the history of the population, the same analysis was repeated with the prior for t restricted to the onset of commercial fisheries until present (analysis 2). All prior values are listed in Supplementary Table S3. To convert time in generations to years, an average generation time of 9-10 years was assumed 46 . A data set consisting of 3 × 10 6 iterations was simulated (1 × 10 6 simulations per scenario). Summary statistics calculated included mean number of alleles across loci, mean gene diversity across loci, mean allele size variance across loci, and the mean ratio between the number of alleles and the range in allele size (Garza-Williamson index 47 ). To identify the most likely demographic scenario, relative posterior probabilities with 95% confidence intervals (CI) were estimated using a logarithmic method as described by Cornuet  The extra composite parameter N ratio = θ/θ anc was also estimated, which characterizes the strength of the population size change, with values < 1 indicating a contraction and values > 1 an expansion. Hence, significant departures from the expectations of a stable population were evident when 95% CI of N ratio did not encompass 1 43 . To ensure that enough points were sampled near the top of the likelihood surface, 13 successive iterations were run with 500 parameter points and 200 simulated ancestral histories per point. Scaled population size (θ, θ anc ) and time (D) parameters were converted to effective population sizes (N e and N anc ) and time in years (T) using a microsatellite mutation rate of µ = 4.23 × 10 -4 as estimated by DIYABC (see results), which is consistent with rates previously used in sharks 52,53 . As in the ABC analysis, a generation time of 9-10 years was assumed.

Results
None of the microsatellite markers displayed stuttering, null alleles or large allele drop-out. No significant LD was detected between pairs of loci. One marker (Ct243) deviated significantly from HWE when all samples were analysed jointly (P value after BY correction: 0.026), but no deviations were observed when testing by sampling location. Population pairwise F ST and D est ranged from 0.005 to 0.020 and from 0 to 0.048, respectively (Supplementary Table S1) but none of the comparisons were significant after correction for multiple comparisons. The PCoA further indicated genetic homogeneity, with the first two axes only explaining 4.34% and 3.63% of the variation in the data, and overlap of groups from different sampling sites being evident ( Supplementary  Fig. S1). These results indicate that Carcharias taurus in South Africa comprises a single population and justify the inclusion of all available data in the subsequent demographic analyses.

Approximate Bayesian computation.
In the exploratory ABC analysis step, suitable prior ranges for the model comparison were successfully established (see Supplementary Table S2 and Fig. S2 for results). For Analysis 1, which included a wide prior for the time at which the size of the population changed (t), scenario 2 (population expansion), had the highest posterior probability which was close to the maximum possible value of 1.0 with very narrow 95% CIs (Fig. 2a). When t was restricted to the past ~ 400 years (analysis 2), Scenario 1 (constant population size) had the highest support when estimating posterior probabilities (Fig. 2b). In both analyses, scenario 3 (population decline) received no support. Because simulations with t restricted to the   Table S4, Fig. S3, Fig. S4).

Discussion
Estimates of contemporary effective population size may assist in guiding management decisions for exploited and vulnerable elasmobranch populations, as they provide information on the size of the breeding population and its genetic health 18,24 . Here, two methods used to reconstruct the demographic history of Carcharias taurus in South Africa provided no evidence for recent, anthropogenic population decline, and instead, they both identified a historical population expansion with parameter estimates that were largely congruent between the two methods. The scenario recovered here may undoubtedly oversimplify the true demographic history of this population but, as advocated by Cabrera and Palsbøll 2017 54 , it is reasonable to focus on a few simple, contrasting demographic scenarios that capture key demographic events. Based on the ABC analyses, the population expansion occurred ~ 32,960 years ago, where a historic N e of only 1,100 individuals expanded to a current effective size of approximately 15,900 individuals. Results obtained through IS are concordant with a strong demographic expansion and inferred higher population size estimates (N e = 20,978 and N anc = 6,842), but with 95% CIs that broadly overlapped with those of the ABC results. The timing of this event was dated more recently at ~ 14,274 years ago, but again, 95% CIs from both methods overlapped broadly, ranging from ~ 2,097-87,750 (ABC) and ~ 1,722-71,669 (IS) years. Considering the wide confidence ranges and the uncertainty concerning some parameters (e.g. mutation rate, generation time), it is not possible to clearly link this demographic expansion to any specific geological or climatic events, but it is likely a result of climate oscillations during the last glacial cycle of the late Pleistocene 55 . Cycles of global climate change have caused widespread changes in the distribution and abundance of flora and fauna, leaving distinct genetic signatures of glacial refugia and postcolonization events in contemporary populations [56][57][58] . Despite earlier assumptions that highly mobile marine species may have been less affected by fluctuating climate conditions, expansions of elasmobranch populations following the last glacial maximum (~ 26,500-19,000 years ago) are commonly reported [59][60][61][62] . Carcharias taurus is a coastal species, and migratory movements are restricted by limits to its thermal tolerance 63,64 . Hence, it is likely that severe changes in oceanographic dynamics during climate oscillations have led to demographic changes and re-distributions of C. taurus in the Southern Hemisphere. Although we cannot link the expansion observed here to a specific climatic phase, it is clear that this event is unrelated to human activities. Estimates of N e and the absence of evidence for a recent, anthropogenic decline support the notion that the South African population is relatively healthy, as previously indicated by mark-recapture and beach meshing catch data 34,36 . Estimates of N e found here closely approximate the mean estimate of the number of adult sharks obtained by mark-recapture data (16,700) 36 . This is in agreement with previous findings of high N e /N c ratios in the East Australian C. taurus population 65,66 and other elasmobranchs, such as the zebra shark, Stegostoma fasciatum 17 and the sandbar shark, Carcharhinus plumbeus 18 . The observed pattern can be explained by similar life history traits, especially age at maturity and adult lifespan, which were found to account for half of the variation found in N e /N c ratios across 63 animal and plant species 19 . Moreover, theoretical and empirical evidence show that high variability in reproductive success can strongly reduce N e /N c ratios [67][68][69] . Taking into account that the reproductive output in C. taurus is small and constant (i.e. a maximum of two pups per litter as a result of intrauterine cannibalism) 70 , the reproductive variance is expected to be very small and this should increase the N e /N c ratio compared to organisms with higher fecundity, such as teleost fish.
Exactly what size of N e is needed to retain evolutionary potential is still controversial. Guidelines suggest that an N e of at least 50 is required to prevent inbreeding depression 71 , and an N e greater than 1,000 is needed to avoid the accumulation of deleterious alleles 12 . In order to retain the long-term evolutionary potential of a species, it has been suggested that an N e greater than 500 or even 1,000 is necessary [72][73][74] . Although estimates of N e were derived with relatively large 95% CIs, it can be concluded that the South African C. taurus population is nowhere near the beginning of genetic erosion. This is in stark contrast to the East Australian population, where a recent study found N e to be as small as ~ 400 individuals 25 . The species gained a bad reputation in Australia due to its fierce appearance and was culled by spear and line fishers, resulting in severe population declines 75 . Despite recent conservation efforts, the Australian C. taurus populations are still under threat from incidental catch by recreational and commercial fishers, as well as beach meshing 75 . These threats, in combination with the species' extremely low reproductive output, seem to have prevented a recovery, particularly along the East Coast. Industrial fisheries started in South Africa around 1890, but the exploitation of marine resources commenced as long ago as the arrival of Dutch settlers in the mid-seventeenth century 76 . Carcharias taurus is susceptible to being caught as bycatch by small-scale coastal fisheries and the beach meshing program due its preference for shallow waters, but the species was never specifically targeted, and decommercialized through the Marine Living Resources Act of 1998 77 . This could explain the lack of evidence for a more recent population decline and indicates that anthropogenic pressure was not high enough to significantly reduce the size of this population. Given the critical status of this shark across much of its distribution, the South African C. taurus population may be the last remaining healthy population globally.

Data availability
The multilocus genotype data file generated and analysed in this study is available in the Supplementary Information online.