Population genetics of four heavily exploited shark species around the Arabian Peninsula.

The northwestern Indian Ocean harbors a number of larger marine vertebrate taxa that warrant the investigation of genetic population structure given remarkable spatial heterogeneity in biological characteristics such as distribution, behavior, and morphology. Here, we investigate the genetic population structure of four commercially exploited shark species with different biological characteristics (Carcharhinus limbatus, Carcharhinus sorrah, Rhizoprionodon acutus, and Sphyrna lewini) between the Red Sea and all other water bodies surrounding the Arabian Peninsula. To assess intraspecific patterns of connectivity, we constructed statistical parsimony networks among haplotypes and estimated (1) population structure; and (2) time of most recent population expansion, based on mitochondrial control region DNA and a total of 20 microsatellites. Our analysis indicates that, even in smaller, less vagile shark species, there are no contemporary barriers to gene flow across the study region, while historical events, for example, Pleistocene glacial cycles, may have affected connectivity in C. sorrah and R. acutus. A parsimony network analysis provided evidence that Arabian S. lewini may represent a population segment that is distinct from other known stocks in the Indian Ocean, raising a new layer of conservation concern. Our results call for urgent regional cooperation to ensure the sustainable exploitation of sharks in the Arabian region.

Introduction

Understanding the spatio-temporal patterns of gene flow among geographically separated populations has long been a major focus in ecology. Limited genetic differentiation over broad spatial scales is often associated with the high dispersal capacities of marine organisms, resulting from either a highly dispersive larval phase affected by ocean currents or the active movements of juvenile and adult specimens in animals lacking a planktonic larval stage. Yet, there are numerous well-known examples of barriers to gene flow within and among populations that result in higher than expected genetic structure, even in species with presumed high levels of vagility (e.g., dolphins: Andrews et al. 2010; Möller et al. 2011; killer whales: Foote et al. 2011; sharks: Blower et al. 2012; tuna: Dammannagoda et al. 2008; Kunal et al. 2013).

Patterns of genetic population structure in sharks are not uniform across species, but range from localized genetic subdivision (e.g., leopard shark: Lewallen et al. 2007; nurse shark: Karl et al. 2012; zebra shark: Dudgeon et al. 2009) and population structuring on relatively small geographic scales (e.g., blacktip reef shark: Vignaud et al. 2014a; bull shark: Karl et al. 2011; dusky shark: Benavides et al. 2011; grey nurse shark: Ahonen et al. 2009; lemon shark: Schultz et al. 2008; sandbar shark: Portnoy et al. 2010), to population differentiation detectable only across ocean basins (e.g., shortfin mako shark: Schrey and Heist 2003; whale shark: Castro et al. 2007; Schmidt et al. 2009; Vignaud et al. 2014b) and nearly global panmixia (basking shark: Hoelzel et al. 2006). Genetic subdivision in sharks is commonly facilitated by geographic dispersal barriers, such as large oceanic expanses (lemon shark: Schultz et al. 2008; spot-tail shark: Giles et al. 2014) or environmental gradients along continuous landmasses extending across different geographic regions (blacktip shark: Keeney and Heist 2006). In addition, the degree of species- and/or location-specific genetic differentiation is typically reflected by a combination of individual vagility, foraging habits, habitat preferences, reproductive mode, and sensitivity toward natural and anthropogenic influences (Dudgeon et al. 2012). The wide range of life histories and movement patterns exhibited by even closely related shark species hence hampers the a priori inference of spatial population structure.

There is compelling evidence to investigate the genetic population structure of sharks in the water bodies surrounding the Arabian Peninsula, that is, the Arabian Sea, the Gulf of Oman and two semi-enclosed bodies of water, the Red Sea, and the Arabian/Persian Gulf (hereafter “the Gulf”) (Fig.1). First, a number of resident marine vertebrate taxa display remarkable heterogeneity in biological aspects, such as distribution, behavior, morphology, and population genetics. The Arabian Sea off the Oman coast, for instance, harbors the world's most isolated and most distinct population of nonmigratory humpback whales, Megaptera novaeangliae (Pomilla et al. 2014). Hawksbill turtles in the Gulf are significantly smaller than those in Omani waters (Pilcher et al. 2014), and sea snakes, which are abundant and diverse in the Gulf and present in the Arabian Sea, are entirely absent from the Red Sea (Sheppard et al. 1992). In addition, barriers to gene flow have been indicated between the Red Sea and the western Indian Ocean for several invertebrates (crabs: Fratini and Vannini 2002; sponges: Giles et al. in press) and some reef fishes (DiBattista et al. 2013), but not for others (Kochzius and Blohm 2005; DiBattista et al. 2013). In the Gulf, the large and highly mobile sailfish, Istiophorus platypterus, was described as phylogeographically isolated (Hoolihan et al. 2004), while another epipelagic predator, the Spanish mackerel, as well as the fiddler crab, does not appear to exhibit genetic subdivision between the Gulf and the Arabian Sea (Hoolihan et al. 2006; Shih et al. 2015).

Figure 1

Map of the Arabian Sea region, displaying collection locations (circles) of Carcharhinus limbatus, C. sorrah, Rhizoprionodon acutus, and Sphyrna lewini. Numbers indicate fish markets or landing sites in Saudi Arabia, Oman, the United Arab Emirates, and Bahrain from where samples were obtained. (1) Jeddah, (2) Salalah, (3) Mirbat, (4) Masirah, (5) Sur, (6) Muscat, (7) Seeb, (8) Barka, (9) Sohar, (10) Shinas, (11) Dibba, (12) Khasab, (13) Ras Al Khaimah, (14) Sharjah, (15) Dubai, (16) Abu Dhabi, (17) Bahrain. See Table S1 for number of tissue samples obtained from each landing site or fish market. Triangles display other main landing sites in Saudi Arabia from which sharks are transported to the main fish market in Jeddah. Geographical color codes refer to haplotypes in Fig.2.

Figure 2

Mitochondrial control region haplotype networks for Carcharhinus limbatus (A), C. sorrah (B), Rhizoprionodon acutus (C), and Sphyrna lewini (D) constructed by statistical parsimony in TCS 1.21 (Clement et al. 2000). Circles are sized in proportion to the number of individuals with that haplotype. Each connecting line represents a single mutation. Black dots represent inferred mutational steps. Ocean basins are indicated by colors: The study region is color coded by geographical regions displayed in Fig.1, dark blue (Red Sea), green (OAB). Haplotypes sampled in previous studies are indicated by red (Atlantic), yellow (Indian), light blue (Pacific), yellow fading to blue (shared Indian Pacific), gray (South-East Asia), purple (Australia), salmon (New Caledonia) and are numbered to match their designations in those studies. (A) CL5–CL7 represent novel haplotypes discovered in this study. Haplotypes sampled in previous studies are indicated by ovals (Keeney et al. 2003, 2005; Keeney and Heist 2006) and rectangles (Sodré et al. 2012). CL1–CL4 are identical to Indian Ocean and Indo-Pacific haplotypes discovered by Keeney and Heist (2006). CL1 = H33; CL2 = H24, H26, H27, and H35; CL3 =  H31; and CL4 =  H32. (B) CS4–CS10 and CS12 represent novel haplotypes. Haplotypes sampled by Giles et al. (2014) are represented by ovals. Haplotype CS1 is identical to H5, CS2 to H36, CS3 to H11, CS11 to H12, CS13 to H6, CS14 to H26, and CS15 to H38 in Giles et al. (2014). (D) SL1–SL5 represent novel haplotypes. Haplotypes sampled in previous studies are indicated by ovals (Duncan et al. 2006), rectangles (Chapman et al. 2009), and triangles (Nance et al. 2011).

Second, existing studies suggest variation in distributional and morphological patterns within Arabian elasmobranch species. Several elasmobranch species in the Arabian region have highly localized known distributions (e.g., C. leiodon: Moore et al. 2011) with a number of species endemic to the Red Sea (e.g., H. bentuviai: Baranes and Randall 1989) and the Gulf (e.g., H. randalli: Last et al. 2012). In addition, a large number of common elasmobranch species, which are reliably reported from the Gulf of Oman and the Gulf of Aden, have not been reported in the Gulf and the Red Sea, respectively (Moore 2011; Spaet et al. 2012). Furthermore, significant morphological differences between Gulf elasmobranchs and “typical forms” were suggested (Moore 2011), and a number of Gulf and Red Sea taxa still remain undescribed (unpublished data).

Recent global genetic studies of elasmobranchs have identified the Arabian region as one of four regions harboring a substantial proportion of taxa that are genetically distinct from their closest relatives in neighboring regions (Naylor et al. 2012). Moreover, global and range-wide studies on several species that included samples from ocean basins in the Arabian region demonstrated substantial genetic differentiation between this region and widely separated Indo-Pacific locations, as well as a strong separation between Indo-Pacific and Atlantic clades for blacktip reef (Vignaud et al. 2014a), silky (Clarke et al. 2015), spot-tail (Giles et al. 2014), and whale sharks (Schmidt et al. 2009; Vignaud et al. 2014b). Yet, in spite of the evident ecological distinctiveness of this region, no study to date has specifically focussed on the genetic population structure of elasmobranchs or indeed any other large vertebrate species around the Arabian Peninsula.

Despite its ecological relevance, the Arabian region features an alarming fisheries situation. Traditional and industrial shark fisheries exist throughout most of the region and for several countries have reached unsustainable exploitation levels (Bonfil 2003; Moore 2011; Jabado et al. 2014a; Spaet and Berumen 2015). Nonetheless, management strategies for shark resources are found in only a fraction of these countries, and proper enforcement of fisheries laws is essentially nonexistent (Bonfil 2003; Moore 2011; Spaet and Berumen 2015). In addition to an apparent general lack of concern toward the conservation of sharks in this region (Bonfil 2003; Spaet and Berumen 2015), the proper assessment and management of elasmobranch stocks has so far been hampered by insufficient information on the biology, ecology, and fisheries of exploited species (Moore 2011; Spaet et al. 2012). Only recently, efforts have been made to bridge this gap, contributing to our knowledge on country-specific fisheries and species-specific biological characteristics (Bonfil 2003; Henderson et al. 2006, 2007, 2009; Moore 2011; Spaet et al. 2011; Moore et al. 2012; Moore and Peirce 2013; Jabado et al. 2014a; Spaet and Berumen 2015).

Patterns of dispersal and population structure can vary significantly even among closely related species in shared habitats (Toonen et al. 2011; DiBattista et al. 2012). Therefore, in this study, we aimed to assess the genetic population structure of four shark species within the Arabian region with different biological, ecological, and life-history characteristics: (1) the blacktip shark, Carcharhinus limbatus (Müller & Henle, 1839); (2) the spot-tail shark, Carcharhinus sorrah (Müller & Henle, 1839); (3) the milk shark, Rhizoprionodon acutus (Rüppell, 1837); and (4) the scalloped hammerhead shark, Sphyrna lewini (Griffith & Smith, 1834).

Carcharhinus limbatus and S. lewini are found in coastal and semi-oceanic waters worldwide, although several studies suggest that undescribed diversity exists within both species (e.g., Zemlak et al. 2009; Naylor et al. 2012). Carcharhinus sorrah is found on continental and insular shelves, in the tropical and subtropical Indo-West Pacific, and R. acutus occurs along the continental shelf across the eastern Atlantic and Indo-West Pacific (Compagno 2001). Rhizoprionodon acutus is the smallest of the four species and reaches maximum total lengths (TL) of 98 cm in the study region while C. sorrah, C. limbatus, and S. lewini can reach 196 cm, 287 cm, and 303 cm TL, respectively (R. W. Jabado unpubl. data). Carcharhinids and Sphyrnids are placental livebearers with typically low intrinsic rates of increase. Although S. lewini exhibits the highest fecundity of all four study species, (12–41: White et al. 2008 cf. 1–11: Carrier et al. 2012 (range of the other three species)), resilience to exploitation is low due to the species' late age at maturity (10–30 years: Baum et al. 2007 cf. 2–7 years Compagno 1984). Based on International Union for Conservation of Nature (IUCN) Red List criteria, R. acutus is globally categorized as Least Concern, C. limbatus and C. sorrah are classified as Near Threatened, and S. lewini is listed as Endangered. Except for R. acutus, dispersal capacities for all species are considered very high. Tagging studies of C. sorrah demonstrated an individual maximum travel distance of 1116 km, although almost half of the tagged specimens were recaptured within 50 km of the tagging location (Stevens et al. 2000). Movements of up to 2148 km were observed for C. limbatus (Kohler et al. 1998), and an individual S. lewini specimen has reportedly traversed 1600 km of deep ocean habitat (Kohler and Turner 2001). Although no movement studies are available for R. acutus, the smaller body size of this species implies lower vagility compared to the three larger species, potentially indicating greater genetic subdivision.

We use a combination of mitochondrial (control region (CR)) and nuclear (microsatellites) markers. Congruence between both types of markers has been shown to yield a high degree of intraspecific resolution, providing a useful tool for the delineation of marine lineages and populations (e.g., Nance et al. 2011; Ovenden et al. 2011). Moreover, contrasting nuclear and mitochondrial data have been used successfully to identify sex-biased dispersal patterns in different elasmobranch species (e.g., Pardini et al. 2001; Portnoy et al. 2010; Daly-Engel et al. 2012). By combining two kinds of genetic markers over four species with variable biology, life-history characteristics, and vagility, we intend to resolve intraspecific spatial genetic patterns representative of a range of elasmobranchs in this region. We discuss the implications of our findings in light of fisheries management and conservation in the Arabian Peninsula.

Materials and Methods

Sample collection and DNA extraction

Tissue samples of C. sorrah and R. acutus were collected between 2010 and 2013 from whole sharks at fish markets and landing sites in Saudi Arabia (Red Sea coast), Oman, the United Arab Emirates (UAE), and Bahrain; C. limbatus and S. lewini were collected from all locations except Bahrain (site 17, Fig.1), where these species were uncommon or absent in a previous landings survey (Moore and Peirce 2013). Details of species-specific sample numbers per landing site are given in Table S1.

Animals were initially identified based on morphological characteristics. Saudi Arabian samples were obtained from one fish market only (Jeddah), but landings at this site originated from fishing grounds spanning the country's entire Red Sea coast (Spaet and Berumen 2015) (Fig.1). Samples from the UAE were collected from landing and market sites in Abu Dhabi, Dubai, Sharjah, and Ras Al Khaimah as described in Jabado et al. (2014a, 2015). Samples from Oman were collected directly from landing sites along the Omani coast; samples from Bahrain were obtained at the wholesale market of the capital, Manama (Fig.1; Table S1).

At all collection sites, special care was taken to avoid inclusion of specimens for which catch location data were unavailable. This was achieved by interviewing fishermen and traders onsite and verifying the obtained information by a thorough assessment of license plates and origin information of transport trucks used. Based on interpreted assisted fishermen interviews, in the Red Sea, 90% of all four species originated from the five main landing sites displayed in Fig.1, from where they were transported to the Jeddah market by trucks. The remaining 10% originated from smaller landing sites along the Saudi Arabian Red Sea coast. Based on the limited operating range of fishing vessels in Saudi Arabia, all fishing grounds were assumed to lie within a 1–30 km radius of the landing sites. While hence no exact catch location data were available, all samples from Jeddah could definitely be assigned to the Red Sea Basin. The operational range of vessels landing into sites in Oman tends to be small, generally limited to within a few kilometers of the landing site (Henderson et al. 2007). Fishermen in the UAE remain in Gulf waters, yet they are known to travel up to 130–185 km from their landing sites to find productive fishing grounds (Jabado et al. 2014b). The majority of Bahrain specimens were caught in local Bahraini waters (Moore and Peirce 2013) although some may have come from nearby Saudi Arabian or Qatari waters. Despite extensive efforts to determine exact catch locations for more detailed seascape genetic analyses, it was not always possible to assign the origin of samples to their respective landing site regions with 100% certainty. As a precautionary approach, all genetic analyses were hence run with pooled data for the two main geographic groups, combining all samples obtained from the Red Sea into one group (Red Sea) and all samples obtained from outside the Red Sea into a second group representing other Arabian basins (OAB), that is, the Arabian Sea, the Gulf of Oman, and the Gulf (Fig.1).

At all market locations, small fin clips or gill tissue were collected from each specimen and preserved in 99% ethanol. Total genomic DNA was extracted from 10 to 20 mg of preserved tissue using the Macherey-Nagel Genomic DNA from tissue extraction kit (Bethlehem, PA) following the manufacturer's instructions and subsequently stored at −80°C until further analysis.

Microsatellites – laboratory methods and data analysis

Shark samples were genotyped at 8 to 12 microsatellite loci (C. limbatus, 12 loci; C. sorrah, 9 loci; R. acutus, 8 loci; S. lewini, 12 loci). Microsatellite loci were adopted from Feldheim et al. (2001), Keeney and Heist (2003), Ovenden et al. (2006), and Nance et al. (2009) and were directly applied to target species or cross-amplified in nontarget species. Between two and three multiplex PCRs were performed per individual for all species. PCRs were performed in 11 μL total volume containing 2 μL genomic DNA, 5 μL Qiagen Multiplex PCR Master Mix, 3.5 μL H20, and 0.5 μL of primer mix (each primer at 2 μmol/L). Thermal profiles consisted of a denaturation step at 95°C for 15 min, followed by 30 cycles of 30 sec at 94°C, annealing for 90 sec at loci-specific temperatures between 55°C and 60°C (Table S2), and an extension of 60 sec at 72°C, with a final extension of 30 min at 60°C. Fragment analysis was conducted in an Applied Biosystems 3730 XL genetic analyzer, and microsatellite alleles were scored using GENEMAPPER software (v4.0 Applied Biosystems, Foster City, CA). The null hypothesis of Hardy–Weinberg equilibrium (HWE) was tested using GENEPOP on the Web (v4.2 Rousset 2008). MICRO-CHECKER (v2.2.3 van Oosterhout et al. 2004) was used to determine likely causes for deviations from HWE. GENEPOP was also used to characterize genetic diversity (expected (HE), observed (HO) and unbiased (UHE) heterozygosity, allelic richness, and mean number of alleles.

STRUCTURE (v2.3.4 Pritchard et al. 2000) was used to infer the number of putative discrete populations in all samples. We set K = 1–10 for each run, assuming prior population information and an admixture model allowing for mixed ancestry of individuals. Each run was repeated three times with independent allele frequencies, 100,000 steps, and a burn-in of 10,000 steps. We used STRUCTURE Harvester (Earl 2012) to determine which K best describes the data according to the highest averaged maximum-likelihood score and Evanno's delta K (Evanno et al. 2005). We then re-ran STRUCTURE with pooled data for the two main geographic groups, combining all samples obtained from the Red Sea into one group (Red Sea) and all samples obtained from the OABs into a second group (Fig.1). A hierarchical analysis of molecular variance (AMOVA) implemented in ARLEQUIN (v3.5 Excoffier and Lischer 2010) and FST (Weir and Cockerham 1984) values were calculated using ARLEQUIN. All microsatellite FST values were corrected ( in Hedrick (2005)) using SMOGD (v1.2.5 Crawford 2010) to compensate for the downward bias in FST associated with highly variable microsatellites.

Mitochondrial DNA – laboratory methods and data analysis

For each species, we examined genetic subdivision based on sequence variation in the mtDNA CR. Approximately 1120 base pairs (bp) of the 5′end of the mtDNA CR was amplified for C. limbatus, C. sorrah, and R. acutus using the forward primer ProL2 and the reverse primer PheCacaH2 (Pardini et al. 2001). A different primer set was used for S. lewini to identify potential specimens of the recently described cryptic species S. gilberti (Quattro et al. 2013). The forward primer CRF6 and the reverse primer CRR10 (Pinhal et al. 2012) were shown to clearly distinguish between the two species and were hence used in our study to amplify approximately 700 bp of the initial portion of the mtDNA CR for all S. lewini specimens. Amplification protocols were the same for both primers and followed those described in Spaet and Berumen (2015). For S. lewini and R. acutus, 700 bp and 1021 bp of the CR were sequenced in the forward and reverse direction, respectively. For C. limbatus and C. sorrah, approximately 600 bp of the CR, respectively, was sequenced in the forward direction only, but to ensure accuracy of nucleotide designations, rare and questionable haplotypes were sequenced in both directions. The program Codon Code Aligner (v4.7.2 CodonCode Corporation, Dedham, MA) was used to assemble, check, manually edit, and subsequently align sequences using the MUSCLE algorithm. Aligned sequences were exported to FaBox (Villesen 2007) and collapsed into haplotypes. Initial species identifications based on morphological characters during market sampling were confirmed by comparison with reference CR sequences in the GenBank database through BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi). In the case of R. acutus, no reference sequences were available prior to this study. Therefore, to validate initial species identification, all R. acutus samples were amplified for the COI gene using the primer combination Fish F1 and Fish R1 (Ward et al. 2005). The PCR protocol used was identical to the one used for the CR locus. PCR products were purified and sequenced following Spaet and Berumen (2015). Resultant COI sequences were compared to reference sequences in GenBank (http://www.ncbi.nlm.nih.gov) for species recognition. If sequence data did not match the original identification, respective specimens (C. limbatus: three, C. sorrah: four, S. lewini: eight, R. acutus: seven) were excluded from the data set.

For C. limbatus, C. sorrah, and S. lewini, haplotype networks were constructed to explore the relationships among intraspecific haplotypes. Published haplotypes were sourced from Keeney et al. (2003, 2005), Keeney and Heist (2006), and Sodré et al. (2012) for C. limbatus; Giles et al. (2014) for C. sorrah; and Duncan et al. (2006), Chapman et al. (2009), and Nance et al. (2011) for S. lewini, aligned with novel haplotypes for each species, trimmed to one length (C. limbatus: 554 bp, C. sorrah: 455 bp, S. lewini: 534 bp), and subsequently assessed using a statistical parsimony network constructed in TCS (v1.21 Clement et al. 2000). For R. acutus, a parsimony network was constructed based on the haplotypes recorded in this study.

An AMOVA under the Tamura–Nei (TN) model of sequence evolution, which was individually selected as the most appropriate model for all four species in jModelTest (v2.1.4. Darriba et al. 2012), was used to assess population genetic structure in ARLEQUIN. ARLEQUIN was also used to describe the genetic variation between the Red Sea and OAB sampling regions by haplotype and nucleotide diversity (h and π, respectively).

Ramos-Onsins and Rozas (2002) demonstrated that Fu's Fs neutrality test (Fu 1997) has the greatest power to detect population expansion for non-recombining regions, such as mtDNA, under a variety of different circumstances, when population sample sizes are large (>50). We hence calculated Fu's Fs to assess deviations from selective sequence neutrality that could be attributed to selection and/or population size changes. Significance was tested with 100,000 permutations. Recent population expansion is indicated by negative (and significant) Fs values. The time since the most recent population expansion was estimated by fitting the population parameter τ (Rogers and Harpending 1992) for both sampling regions and each species. Mutation rate estimates were available from previous studies for S. lewini: 0.8% divergence between lineages per million years or 0.4 × 10−8 mutations per site per year (Duncan et al. 2006) and for C. limbatus: 0.43% or 0.215 × 10−8 (Keeney and Heist 2006); no species-specific mutation rates were available for C. sorrah and R. acutus. For those species, we hence used the averaged mutation rate (0.62%) reported for other shark species (Galván-Tirado et al. 2013). Generation time estimates were available from previous studies for all four species, C. limbatus: 10 years, C. sorrah: 4.3 years (Cortés 2002), and R. acutus: 2.5 years (Simpfendorfer 2003). Generation time estimates for S. lewini are controversial and vary among ocean basins (e.g., Branstetter 1987; Liu and Chen 1999). As no estimates were available for the Indian Ocean, we used the generation time estimated for the closest ocean region for which an estimate was available, the west Pacific: 16.7 years (Cortés 2002). We estimated population expansion times assuming a constant molecular clock and rates using the Mismatch Calculator tool developed by Schenekar and Weiss (2011).

Results

Genetic diversity and summary statistics

Microsatellites

Microsatellite indices of genetic diversity, that is expected (HE), observed (HO), and unbiased (UHE) heterozygosities, allelic richness, and mean number for each locus and species within each sample region are provided in Table S2. No signs of linkage disequilibrium were detected among any pairs of loci after correction for multiple comparisons.

In all species, several microsatellite loci showed deviations from HWE in one or both of the putative populations and signs of null alleles (Table S2). To test whether significant differences between expected vs. observed heterozygosities at some loci could confound population level analyses, we removed all those loci and re-ran AMOVA analyses. A comparison of FST values calculated from the subset of loci in HWE and from the full data set was not significant for any of the species (paired t-tests calculated in JMP P > 0.6 in all species). To ensure that the pattern of microsatellite structure (or lack thereof) was not being driven by a single locus, we conducted locus-by-locus AMOVA analyses (data not shown), which gave consistent results across all except one locus (Cli118, C. sorrah). This locus was solely responsible for the observed pattern of significant population structure and was subsequently removed from the analysis.

Mitochondrial DNA

Low haplotype (h) and nucleotide (π) diversities were found for C. limbatus and S. lewini, while C. sorrah and R. acutus showed slightly higher h and π values (Table1).

Table 1

Mitochondrial DNA control region sample size and genetic diversity indices for Carcharhinus limbatus, C. sorrah, Rhizoprionodon acutus, and Sphyrna lewini across both sampling regions. Haplotype (h) and nucleotide (π) diversities, neutrality statistics (Fu's FS), and estimates of times since last population expansion are shown. Population expansion ranges (below expansion times) are given for 95% confidence intervals of tau

Species	n (Red Sea)	n (OAB)	Time since expansion Yrs (Red Sea)	Time since expansion Yrs (OAB)	h ± SD Red Sea	h ± SD OAB	π ± SD Red Sea	π ± SD OAB	Fu's FS Red Sea	Fu's FS OAB
C. limbatus	172	115	182,772 (96,129–273,696)	269,498 (4.2–1,273,612)	0.3490 ± 0.0387	0.3054 ± 0.0525	0.000724 ± 0.00748	0.000755 ± 0.000769	−2.90484	−1.89320
C. sorrah	159	216	214,134 (11–656,541)	178,331 (8551–478,324)	0.3270 ± 0.0479	0.4606 ± 0.0412	0.001030 ± 0.000893	0.001314 ± 0.001408	−7.23390	−7.43393
R. acutus	77	217	190,831 (89,096–358,124)	177,561 (48,023–378,345)	0.7365 ± 0.0345	0.6599 ± 0.0306	0.001397 ± 0.000956	0.001265 ± 0.000881	−3.57008	−12.17461
S. lewini	82	151	151,245 (71,174–249,110)	139,679 (64,724–194,172)	0.4998 ± 0.0318	0.4661 ± 0.0317	0.000086 ± 0.000232	0.000116 ± 0.00027	0.73386	−0.13352

Numbers in bold are significant, P < 0.02 for Fu's FS estimates, Fu (1997).

Fu's Fs statistics were negative for all four species and both sampling regions, yet significant only for C. sorrah for both regions (Fu's Fs Red Sea = −7.23; P = 0.002; Fu's Fs OAB = −7.43; P = 0.006) and for R. acutus for the OAB region only (Fu's Fs = −12.17; P = 0.01) (Table1). The range of τ values yielded estimates of time since last population expansion with very similar expansion time estimates in all four species and both regions (139.679–269.498 years, Table1).

Genetic structure

The results obtained from all STRUCTURE runs yielded K = 1, indicating no differentiation among tentative populations.

FST values were small and nonsignificant for mtDNA analyses in all four species. Very low, yet significant genetic population subdivision was found using microsatellite allele frequencies for C. limbatus (0.012; P = 0.00), R. acutus (0.002; P = 0.04), and S. lewini (0.006; P = 0.001) (Table2).

Table 2

FST results and associated P-values for both regions, characterizing spatial structure with both mtDNA and microsatellites. Standardized FST values (, Hedrick 2005) are shown in brackets

	mtDNA	Microsatellites
Carcharhinus limbatus	0.0025; P = 0.236	0.012; P = 0.00 (0.0128)
C. sorrah	0.0057; P = 0.099	0.000; P = 0.58 (0.000)
Rhizoprionodon acutus	0.0608; P = 0.583	0.002; P = 0.04 (0.000864)
Sphyrna lewini	0.0130, P = 0.050	0.006; P = 0.001 (0.009604)

Carcharhinus limbatus

A 554-bp sequence was obtained for 287 C. limbatus individuals. A total of seven haplotypes (GenBank Accession Numbers: KR232982-KR232988) were defined, characterized by five polymorphic sites composed of five transitions (Table S3A). Except for three singletons, all haplotypes were found in both putative populations and matched known Indian Ocean and Indo-Pacific mtCR haplotypes from the global data set of Keeney and Heist (2006). One haplotype (CL1) clearly dominated the sample set and was found in both populations in almost identical numbers (Red Sea: n = 100; OAB: n = 131). Two singletons were unique to the Red Sea and one was unique to the OAB. Novel haplotypes were very closely related to Indian Ocean and Indo-Pacific haplotypes reported in Keeney and Heist (2006) and at least nine mutational steps away from any Atlantic haplotypes (Fig.2A).

Carcharhinus sorrah

A 455-bp sequence was resolved for 375 individuals and resulted in 15 mtDNA haplotypes (GenBank Accession Numbers: KR232989-KR233003), characterized by 12 polymorphic sites composed of 10 transitions, one transversion, and one deletion (Table S3B). All common haplotypes were observed in both putative populations. One haplotype clearly dominated the sample set (CS1). Seven haplotypes matched haplotypes from the range-wide data set of Giles et al. (2014). All novel haplotypes were closely related to Indian Ocean and South-East Asian haplotypes reported in Giles et al. (2014) and formed a lineage distinct from Australian and New Caledonian haplotypes (Fig.2B).

Rhizoprionodon acutus

Variation in a 1021-bp fragment of 294 R. acutus specimens defined 25 haplotypes (GenBank Accession Numbers: KR232957-KR232981) characterized by 22 polymorphic sites composed of 18 transitions, five transversions, and two deletions (Table S3B). All common haplotypes were separated by two mutational steps at most. Three singletons and one haplotype, recorded from two individuals only, were separated from the cluster of common haplotypes by up to 10 mutational steps. Except for one (RA7) that was unique to the OAB, all common haplotypes were shared in both sampling regions. Haplotype RA17 dominated the sample set and was found in more than half of all OAB samples (Fig.2C).

Sphyrna lewini

A 562-bp sequence revealed low levels of diversity for 233 S. lewini specimens: five haplotypes (GenBank Accession Numbers: KR232952-KR232956), characterized by four polymorphic sites composed of two transitions and two transversions (Table S3) that differed by no more than one mutational step from each other. Two haplotypes clearly dominated the sample set (Fig2D). All five haplotypes were novel, that is, not present in the global data set of Duncan et al. (2006) or in any of the regional data sets by Chapman et al. (2009), Nance et al. (2011), and Castillo-Olguín et al. (2012). The parsimony network provided evidence that the haplotypes discovered in this study form a Distinct Population Segment (DPS).

Discussion

This study is the first to assess the population structure of large mobile marine vertebrates between the Red Sea and all other Arabian Ocean Basins. Our analyses were based on a comparatively large number of samples (total n = 1189) of four different shark species, from collection locations spanning across over 5000 km of coastline genotyped at two types of genetic markers (mtDNA and nuclear DNA). Contrary to previous findings of significant population genetic structure across the region in different taxa, our results indicate that dispersal of sharks around the Arabian Peninsula is not limited by any obvious barriers to gene flow. Furthermore, ecological, morphological, and life-history differences among the investigated species do not appear to significantly influence their patterns of population structure. Divergent haplotypes in one of our study species (S. lewini), however, are suggestive of an Arabian population that is genetically distinct from others in the Indian Ocean.

Several previous studies have shown the existence of historical, oceanographical, and ecological barriers to gene flow resulting in genetic subdivision in a range of marine organisms among ocean basins surrounding the Arabian Peninsula. Our analyses did not provide compelling evidence for more than one Arabian Sea genetic stock for C. limbatus, C. sorrah, R. acutus, or S. lewini. There was slight evidence of genetic structure between the Red Sea and the OAB for C. limbatus, R. acutus, and S. lewini based on microsatellite allele frequencies; however, FST values were low (0.002–0.012) and not consistent among different statistical tests. These inconsistencies might stem from the high number of null alleles in our data set, which might have been caused by (1) cross-species rather than species-specific loci used in this study due to the limited availability of microsatellite loci for all investigated species and/or (2) species-specific loci, which were developed for specimens sampled in other ocean regions. For future studies on elasmobranch species from regions that have not previously been included in samples used for the design of microsatellite markers, we hence recommend designing species-specific markers based on samples originating from the targeted study region.

The homogenous population structure observed here was not unexpected, given the contiguous shelf habitat around the Arabian Peninsula and the high potential mobility of our study organisms. While previous regional and range-wide studies on C. limbatus, C. sorrah, and S. lewini demonstrated restricted dispersal across deep ocean habitats, genetic structure along continental margins was shown to be relatively minor (Duncan et al. 2006; Keeney and Heist 2006; Nance et al. 2011; Daly-Engel et al. 2012; Giles et al. 2014). Studies on all four species across spatial scales similar to this study in Australia and Indonesia demonstrated heterogeneous population structure in C. sorrah and R. acutus, but not for S. lewini and C. limbatus between central Indonesia and northern Australia based on nuclear and mtDNA markers (Ovenden et al. 2009, 2010, 2011). The observed subdivision in the two smaller, less vagile species was suggested to arise from the Timor Trough acting as a deep water dispersal barrier (Ovenden et al. 2009). While large expanses of deep water dividing shallow habitats are absent in our study region, potential oceanographic barriers to gene flow may still exist, for example, regional upwelling systems or local turbid water regions (Schott 1983). Present-day oceanic currents and habitat heterogeneity in the study area have recently been suggested to inhibit gene flow in teleost larval dispersal (DiBattista et al. 2013; Nanninga et al. 2014). Sharks, however, are lacking the dispersive larval phase of most teleost fish, and based on our results, their swimming capacities as juveniles and especially as adults are likely too strong to be influenced by ocean currents characteristic of the Arabian region. Intermittent historical barriers like the ones created by Pleistocene glacial cycles have also reportedly impacted gene flow in teleost species between the Red Sea and the Indian Ocean (Klausewitz 1989; DiBattista et al. 2013). A potential significant reduction in population size during this period was demonstrated by negative and significant indices of neutral evolution (Fu's Fs test) for C. sorrah and R. acutus, indicating recent population expansion events between approximately 178,000 and 214,000 years ago (Table1). Those events likely followed substantial bottleneck events that were caused by re-occurring limitations of inflow and exchange of surface water between the Red Sea and the Indian Ocean (Siddall et al. 2003). The decrease in sea water level during those periods likely caused increased evaporation, raising temperatures, and salinity levels beyond the tolerance limits of most marine fauna (Biton et al. 2008). Another reason for the observed excess of low-frequency haplotypes might be caused by positive selection. However, to unambiguously discern between the effects of natural selection and demographic population expansion would necessitate an analysis of several unlinked loci in the genome, because selection only acts on specific loci (Akey et al. 2004).

Additionally to the apparent homogenous population structure, we also found no indication of differences between male and female dispersal in any of the study species. This finding stands in contrast to previous studies describing marked philopatric behavior (Feldheim et al. 2014) in C. limbatus and S. lewini based on contrasting mitochondrial and nuclear data (Keeney et al. 2003, 2005; Daly-Engel et al. 2012). We suggest that long-shore movements of both males and females along the continuous coastline stretching from the Red Sea all the way into the Gulf cause panmixia over large spatial scales across the region.

Genetic diversity for C. limbatus and S. lewini was relatively low, with only seven and five haplotypes, respectively. Yet, this pattern appears to be typical for both species throughout their global range (Duncan et al. 2006; Keeney and Heist 2006) and hence may not necessarily be a function of overexploitation. While all our samples of C. limbatus and C. sorrah matched previously published Indian Ocean and Indo-Pacific mtDNA haplotypes, all haplotypes discovered for S. lewini were novel.

There are two possible explanations for the observed genetic separation between S. lewini specimens sampled around the Arabian Peninsula and specimens sampled in other, nearby Indian Ocean regions (e.g., the Seychelles and Madagascar, Duncan et al. 2006). First, Arabian Seas S. lewini may have evolved to breed differently from conspecifics outside this area. Estuaries have repeatedly been reported as an important nursery habitat for S. lewini elsewhere in the world (e.g., Clarke 1971; Snelson 1981; Simpfendorfer and Milward 1993; Duncan and Holland 2006). Due to the desertification of the Arabian region in the past few thousand years, permanent estuaries are now entirely absent for several thousand kilometers of continental coastline from Iraq to Somalia, an area that encompasses our study region. Suggested nursery areas and breeding grounds for S. lewini, however, exist near Djibouti City (Bonfil 2003) and in habitats in the central Saudi Arabian Red Sea (J. L. Y. Spaet, unpubl. data), suggesting that the species may not depend on estuarine habitat in these areas. Arabian S. lewini may thus have evolved to no longer require estuaries as breeding/nursery grounds, eliminating the need for reproductive migrations and thereby reducing gene flow with other populations. Such scenarios may also explain why C. limbatus and C. sorrah, which are not reported as being strongly dependent on estuary nurseries, are genetically well connected to other Indian Ocean populations. Second, regional oceanography and upwelling zones may form temporary barriers between Arabian and other Indian Ocean populations. The Somali Current, for instance, which only operates between June and September (Schott 1983), may coincide with key migration/breeding periods of S. lewini, but not with those of C. limbatus and C. sorrah. In this case, mixing with south Indian Ocean populations might be inhibited for S. lewini, but not for the other two species. Additional data on migration routes, migration times and breeding cycles, however, are needed to confirm either of these hypotheses.

The fact that Arabian S. lewini, which comprise a large amount of the commercial harvest in the Arabian region (Jabado et al. 2015; Spaet and Berumen 2015), might represent a DPS raises a new layer of conservation concern and may warrant species-specific conservation actions under the Convention on International Trade in Endangered Species (CITES), Convention on the Conservation of Migratory Species of Wild Animals (CMS), and a re-evaluation of its IUCN Red List conservation status. Future research should focus on the identification of broader scale genetic breaks by sampling all four species further to the west and east of our sampling locations. In addition, research on dispersal mechanisms based on nongenetic techniques, for example, tagging or parasite studies coupled with molecular methods would provide interesting insights into the actual dispersal mechanisms underlying the observed homogenous population structure.

Conclusions

Molecular studies on a diverse range of elasmobranch species have done much to illuminate issues that complicate fisheries management and conservation (see Dudgeon et al. 2012 for a review). Here, we provide the first multispecies analysis of population structure between Red Sea and Arabian Sea, Gulf of Oman, and Gulf elasmobranchs indicating that dispersal of four different shark species is not limited by any obvious barriers to gene flow in the waters surrounding the Arabian Peninsula. Three broad conclusions are apparent:

Existing contemporary barriers such as regional upwelling systems and ocean currents are likely not influencing long-shore or stepping-stone connectivity even in smaller, less vagile shark species like R. acutus.

A comparison of novel S. lewini haplotypes with published western Indian Ocean haplotypes revealed the possibility of a S. lewini population in the Arabian region that is distinct from other Indian Ocean populations.

Similar dispersal patterns in sharks with contrasting ecological, morphological, life-history, and distributional patterns indicate that populations of other shark species are likely to also function as common stocks across all ocean basins surrounding the Arabian Peninsula.

Overall, our results call for urgent regional cooperation on the management of shark stocks in all countries surrounding the Arabian Peninsula to ensure a sustainable future for this vital component of the marine biodiversity in the western Indian Ocean. Regulations on the exploitation of only one part of the stock will not suffice and management arrangements need to be implemented, enforced, and coordinated among all responsible authorities. Given current harvesting levels and the apparent connectedness of stocks, unregulated exploitation in one or several countries is likely to cause uniform depletion across the entire stock.