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.

Introduction

The size of a population is a determining factor in how evolutionary and demographic processes affect its long-term survival. In conservation biology, small populations are of concern as they are typically more susceptible to demographic, environmental and stochastic genetic events[1-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 (Ne) 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 Ne is often reduced in comparison to actual, or census population size (Nc)[7-9].

Several studies have attempted to establish the ratio of Ne and Nc, 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 Ne. Reported average Ne/Nc ratios across a range of species from different taxa range from approximately 0.1 to 0.5[11-13]. Ne/Nc ratios in teleost fish with high fecundity, high juvenile mortality and sweepstake recruitment (type III survivorship), can be as low as 10–5[14,15] (but see Waples et al. 2018[16]). In contrast, in elasmobranchs (sharks, rays and skates), Ne closely approximates Nc[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]. Nevertheless, considerable variation exists among species, as shown by Chevolot et al. 2006[20], where Ne/Nc 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 Ne 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 Ne and demographic changes through time can provide valuable information to guide effective management and conservation of vulnerable elasmobranch species and populations[24-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 FST and Dest 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 Ne.

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 Ne. 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 × 106 data sets) to determine posterior distributions of Ne 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 Ne 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). 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 × 106 iterations was simulated (1 × 106 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 et al. 2008[48]. A detailed description of the evaluation of model fit and reliability of model choice is provided in the Supplementary Information (‘DIYABC: model evaluation’).

Figure 1

Demographic scenarios for Carcharias taurus tested using approximate Bayesian computation. N1 is the current effective population size and t is the time (in generations) at which the size of the population changed. N2 and N3 represent smaller and larger historical effective population sizes relative to N1, respectively. Prior distributions of each parameter are listed in Supplementary Table S3.

Likelihood-based importance sampling

Likelihood-based importance sampling (IS)[43,49,50] in MIGRAINE version 0.5.4.1 (https://kimura.univ-montp2.fr/~rousset/Migraine.htm) was used as a complementary method to explore the demographic history of C. taurus. A combination of strict likelihood and PAC-likelihood (products of approximate conditionals[51]) computations (PACanc option) was used to infer parameters of interest under a model of a single population undergoing an exponential change in size, starting at some point in the past and continuing to the present. The generalized stepwise mutation model (GSM) was used with the parameter P describing the geometric distribution of stepwise mutations. Point estimates and corresponding 95% CIs were inferred for the following scaled parameters: current population size θ = 2Nµ, ancestral population size θanc = 2Nancµ, where µ is the mutation rate per locus per generation and N is the number of gene copies per deme, and D = Dg/2N, an estimator of the time when population size started to change. The extra composite parameter Nratio = θ/θ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 Nratio 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 (Ne and Nanc) 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 FST and Dest 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.

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 last ~ 400 years did not provide conclusive posterior distributions as much of the genetic signal present in the data was excluded, estimates of demographic parameters were derived from scenario 2 simulated in analysis 1. Numerical values of posterior distributions with corresponding 95% highest posterior density intervals (95% HPDI) for all parameters estimated from scenario 2 are shown in Supplementary Table S4. A population size change occurred approximately 4,120 generations ago (95% HPDI: 669–9,420), corresponding to ~ 32,960 years, where Ne increased from 1,100 (95% HPDI: 134–5,810) to the current effective size of 15,900 (95% HPDI: 9,810–21,600). In addition, simulation-based evaluation demonstrated high power to distinguish between competing demographic models investigated here and validated this model as the most likely demographic scenario. A full description of these results is provided in the Supplementary Information (‘DIY model evaluation’, Table S4, Fig. S3, Fig. S4).

Figure 2

Estimation of relative posterior probabilities of three demographic scenarios based on a logistic regression-based estimate with (a) prior for t set to 1–5 × 103 (analysis 1) and (b) prior for t restricted to 1–50 (analysis 2). Shaded areas indicate 95% confidence intervals. Y-axis is posterior probability, x-axis is number of simulated data sets closest to the observed data.

Inferred pairwise likelihood-ratio profiles showed clear peaks for all parameters, indicating that the inference algorithm performed well (Fig. 3). Current and ancestral population size were estimated at θ = 42.24 (95% CI: 25.76–88.38) and θanc = 12.55 (95% CI: 5.25–21.41). Using a mutation rate of 4.23 × 10–4, these estimates correspond to Ne = 24,964 (95% CI: 15,224–52,234) and Nanc = 7,417 (95% CI: 3,103–12,653), respectively. Nratio was estimated at 3.37 (95% CI: 1.55–7.10), indicating a significant population expansion. This demographic event appears quite recent as the timing since the beginning of this expansion was dated at Dg/2N = 0.0159 (95% CI: 0.0012–0.0830), which corresponds to T = 14,289 years ago (95% CI: 1,078–74,592), assuming a generation time of 9 years.

Figure 3

Pairwise likelihood-ratio profiles for South African Carcharias taurus inferred by MIGRAINE for θanc (2Nancµ): ancestral population size, θ (2Nµ): current population size, D (Dg/2N): timing of demographic event. The profiles show a significant demographic expansion in the population. Axes are on a logarithmic scale. See Results for point estimates and associated 95% confidence intervals.

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 Ne 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 (Ne = 20,978 and Nanc = 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 post-colonization events in contemporary populations[56-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-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 Ne 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 Ne 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 Ne/Nc 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 Ne/Nc ratios across 63 animal and plant species[19]. Moreover, theoretical and empirical evidence show that high variability in reproductive success can strongly reduce Ne/Nc ratios[67-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 Ne/Nc ratio compared to organisms with higher fecundity, such as teleost fish.

Exactly what size of Ne is needed to retain evolutionary potential is still controversial. Guidelines suggest that an Ne of at least 50 is required to prevent inbreeding depression[71], and an Ne 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 Ne greater than 500 or even 1,000 is necessary[72-74]. Although estimates of Ne 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 Ne 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.

Supplementary Information 1.

Supplementary Information 2.