Insights Into Insular Isolation of the Bull Shark, Carcharhinus leucas (Müller and Henle, 1839), in Fijian Waters

The bull shark (Carcharhinus leucas) is a large, mobile, circumglobally distributed high trophic level predator that inhabits a variety of remote islands and continental coastal habitats, including freshwater environments. Here, we hypothesize that the barriers to dispersal created by large oceanic expanses and deep-water trenches result in a heterogeneous distribution of the neutral genetic diversity between island bull shark populations compared to populations sampled in continental locations connected through continuous coastlines of continental shelves. We analyzed 1,494 high-quality neutral single nucleotide polymorphism (SNP) markers in 215 individual bull sharks from widespread locations across the Indian and Pacific Oceans (South Africa, Indonesia, Western Australia, Papua New Guinea, eastern Australia, New Caledonia, and Fiji). Genomic analyses revealed partitioning between remote insular and continental populations, with the Fiji population being genetically different from all other locations sampled (FST = 0.034–0.044, P < 0.001), and New Caledonia showing marginal isolation (FST = 0.016–0.024, P < 0.001; albeit based on a small sample size) from most sampled sites. Discriminant analysis of principal components (DAPC) identified samples from Fiji as a distinct cluster with all other sites clustering together. Genetic structure analyses (Admixture, fastStructure and AssignPOP) further supported the genetic isolation of bull sharks from Fiji, with the analyses in agreement. The observed differentiation in bull sharks from Fiji makes this site of special interest, as it indicates a lack of migration through dispersal across deep-water trenches and large ocean expanses.


INTRODUCTION
Reduced genetic variation, population isolation and genetic drift may lead to genetic divergence and speciation among terrestrial organisms inhabiting islands (Grant, 1985). Assessing genetic structure among island demes can help in identifying small, genetically differentiated, and spatially isolated populations (Allendorf et al., 2008). While terrestrial organisms inhabiting oceanic islands generally show significant levels of genetic differentiation when compared to their continental conspecifics, marine taxa with large dispersal capabilities are expected to show little divergence as their habitats appear to be continuous (Roberts, 1997). However, growing evidence suggests that biogeographic barriers, produced for example by ocean currents (Hare et al., 2005), changes in sea surface temperature (Kelly and Eernisse, 2007), physical barriers such as oceanographic fronts (Galarza et al., 2009), and resource availability (Palumbi, 1994) can impede gene flow (Cowman and Bellwood, 2013) in some marine taxa (Cowen et al., 2000;Leis, 2002;Hauser and Carvalho, 2008). Spatial distribution of the genetic variation to establish connectivity patterns among marine taxa has been subject of extensive research in recent years (Cowen et al., 2000;Cowen and Sponaugle, 2009;Crandall et al., 2019), including sharks (e.g., Momigliano et al., 2017). However, data is limited when it comes to the quantification of population structure in coastal shark species inhabiting both island and continental coastal waters. Such studies can be vital for the identification of management units and metapopulations. For example, an analysis examining the population structure of the tiger shark on a global scale using a suite of mitochondrial and nuclear (microsatellites) genetic markers revealed population isolation for the Hawaiian Islands (Bernard et al., 2016).
Overexploitation has led to drastic declines in shark populations worldwide (Dulvy et al., 2014). In general, various studies indicate that many commercially exploited shark species have a low intrinsic rebound potential due to their low fecundity and slow growth rates, and have further been found to have low genetic diversity and smaller population sizes compared with other marine taxa (Martin et al., 1992;Daly-Engel et al., 2010;Karl et al., 2011;Vignaud et al., 2013). Additionally, coastal sharks are vulnerable to increased mortality rates due to their proximity to human populations, where fishing pressure is typically higher and habitat degradation is more pervasive. Therefore, identifying the genetic structuring of populations is fundamental for determining whether, and to what degree, population differentiation exists. Furthermore, this information is essential for understanding the state of subpopulations that need to be treated as separate management units (Palsbøll et al., 2007).
While global surveys examining broad scale patterns of genetic differentiation confirmed connectivity within and between ocean basins in some highly migratory pelagic sharks (Hoelzel et al., 2006;Castro et al., 2007;Sequeira et al., 2013), oceanic expanses can cause genetic population structure in coastal sharks (Schultz et al., 2008;Ahonen et al., 2009;Ovenden et al., 2009;Portnoy et al., 2010;Momigliano et al., 2017). For example, sicklefin lemon shark populations, Negaprion acutidens, in Australia and French Polynesia, separated by oceanic distances of at least 750 km across the Tonga trench with a maximum depth of over 10,000 m (van der Hilst, 1995), show moderate but statistically significant genetic differentiation (F ST = 0.070-0.087, p < 0.001), with South Pacific archipelagos probably serving as stepping stones for rare dispersal events (Schultz et al., 2008). In addition, to oceanic expanses, physical barriers such as fronts, eddies, gyres and other oceanographic factors can prevent dispersal. In some tropical and subtropical shark species, restricted gene flow is associated with cool thermal barriers such as the Benguela upwelling Keeney and Heist, 2006). For example, in the blacktip shark, Carcharhinus limbatus, the cold Benguela upwelling along the south-West coast of Africa restricts contemporary female-mediated gene flow between Atlantic and Indo-Pacific populations (Keeney et al., 2005;Keeney and Heist, 2006;Sodré et al., 2012). Contrastingly, the Benguela upwelling appears not to be a barrier to the temperate copper shark, C. brachyurus, as evidenced by a lack of genetic differentiation between Namibia and South Africa (Benavides et al., 2011). Additionally, coastal sharks may be limited by their habitat, exhibiting fidelity to discrete locations for feeding, mating, parturition, maturation and migratory routes (Castro, 1983;Simpfendorfer and Milward, 1993;Heupel et al., 2007;Speed et al., 2010;Brunnschweiler and Baensch, 2011). For example, increasing evidence suggests that many coastal sharks have strong philopatric behavior (Karl et al., 2011;Chapman et al., 2015). The degree of segregation of these sites, and site fidelity (philopatry) can directly affect the level of population subdivision and genetic divergence among geographic regions (Keeney et al., 2005). Reproductive philopatry has been described in scalloped hammerheads, Sphyrna lewini , blacktips (Keeney et al., 2005), and bull sharks, Carcharhinus leucas (Karl et al., 2011;Tillett et al., 2012). Furthermore, natal philopatry and long-term fidelity to parturition sites has been documented, for the first time in any shark species at Bimini, Bahamas (Feldheim et al., 2014). Restricted maternal gene flow between populations caused by female site fidelity may result in small local effective population size and thus increased vulnerability to exploitation (Duncan and Holland, 2006;Nance et al., 2011). Thus, overfishing may have detrimental effects on large coastal shark species due to limited exchange with other regional populations.
One of the most iconic coastal apex predators is the bull shark (C. leucas; Müller and Henle, 1839). This large shark is globally distributed in warm, temperate, tropical and subtropical waters (Compagno, 2001). Evidence based on highly restricted maternal gene flow suggests site-fidelity and reproductive philopatry (Tillett et al., 2012), while tracking studies documented largescale movements of up to 1,500 km in individual bull sharks (Carlson et al., 2010;Daly et al., 2014;Heupel et al., 2015;Lea et al., 2015). For example, a recent study documented that a substantial segment of a Great Barrier Reef bull shark population displayed long migrations of up to 1,400 km to other coral reefs and/or inshore coastal habitats along the Australian east coast (Espinoza et al., 2016). Microsatellite markers revealed, for the first time, genetic structuring in West Atlantic bull shark populations (Karl et al., 2011) and strong genetic differentiation between bull sharks from the West North Atlantic and Indo-Pacific, with significant genetic structuring in specimens from Fiji (Testerman, 2014). Similarly, a more recent study by Pirog et al. (2019) utilized both mitochondrial and microsatellites DNA markers in 11 sites across the Western Atlantic, Western Pacific and Western Indian Ocean and found clustering between the Western Atlantic and those from the Indian and Pacific Oceans. However, given that the bull shark is a coastal species (Brunnschweiler et al., 2010), it seems likely that barriers to gene flow exist within its range, which includes distant oceanic islands separated by hundreds of kilometers of open ocean.
New methodological advances such as co-dominant, genome-wide single nucleotide polymorphism (SNP) loci, have the potential to improve resolution in the estimation of population structuring, migration rates, dispersal, and population assignment (Morin et al., 2004;Benestan et al., 2015;Gagnaire et al., 2015;Selkoe et al., 2016). Here, we utilize SNP markers to test the hypothesis that barriers to dispersal created by large oceanic expanses and deep-water trenches result in a heterogeneous distribution of neutral genetic diversity between island bull shark populations compared to demes sampled in continental locations and connected through continuous coastline of continental shelves.
To provide a broad geographic framework for our study, and to place the findings in the context of existing knowledge on bull shark population structure, we obtained samples from South Africa, Indonesia, Western Australia, Papua New Guinea, eastern Australia, New Caledonia and Fiji. Specific aims of the study were: i) to compare inbreeding coefficients and levels of heterozygosity of bull sharks within the sampled range, and ii) to assess the population structure in bull sharks on a genomic scale using SNPs. This study contributes to the understanding of genetic patterns of insular -coastal shark populations that inhabit historically understudied regions.

Ethics Statement
All shipping procedures of shark tissue samples were conducted under the relevant import and export permits issued by The Australian Government, Department of Agriculture and Water Resources to Diversity Arrays Technology Pty Ltd., Canberra, Australia and the Ministry of Fisheries and Forests in Fiji. Bull shark tissues from Papua New Guinea came from the CSIRO/ACIAR/PNG-NFA FIS/2012/102 "Sustainable management of the shark resources of Papua New Guinea: socioeconomic and biological characteristics of the fishery." Handling of live shark specimens were approved under the Research Permit issued by the Fiji Ministry of Education, Heritage and Arts to Beqa Adventure Divers and performed in accordance with relevant guidelines and regulations.

Sampling Procedures
Bull shark samples (white muscle and fin clips, 1 cm 2 ) were obtained from fisheries independent surveys, observers on board commercial fleets, and recreational and artisanal fishers from 2004 to 2017. The samples of 215 individual bull sharks came from seven locations (defined as a priori populations, hereafter) and six countries across the Indian and Pacific Oceans (Figure 1, Table 1). All samples were stored in 95% ethanol until DNA extraction, library preparation and sequencing.

Extraction and Sequencing
Tissue samples were sent to Diversity Arrays Technology (DArT-Seq TM ) in Canberra Australia where genomic DNA was extracted from the tissues using standard robotic methods. DNA was processed for reduced representation library construction, sequenced and genotyped by DArT-Seq TM following previously developed and tested complexity reduction protocols for scalloped hammerhead sharks (Sphyrna lewini)   (Marie et al., 2019). Briefly, genome complexity reduction was achieved with a double restriction digest using a PstI and SphI methylation-sensitive restriction enzyme combination. Libraries were sequenced on an Illumina HiSeq 2500 platform, and raw reads obtained following sequencing were processed using Illumina CASAVA v.1.8.2 software for initial assessment of read quality and sequence representation. The DArT PL proprietary software pipeline, DArTtoolbox was used for further filtering and variant calling to generate the final genotypes set.

SNP Filtering
The dataset was filtered by excluding duplicate SNPs possessing identical Clone IDs and by removing loci with a call rate (proportion of individuals scored for a locus) <99%, read depth >7 and minor allele frequencies (MAF) <2%. Detection of loci under selection was implemented in BayeScan v.2.1 (Foll and Gaggiotti, 2008). The most conservative neutral model in BayeScan was used to minimize falsely detected loci under selection (Lotterhos and Whitlock, 2014). Runs consisted of 100,000 iterations with a burn-in length of 50,000 iterations (Foll and Gaggiotti, 2008;Foll, 2012). Once probabilities had been calculated for each locus, the BayeScan function plot_R was used in the R v.3.2.0 statistical package (Venables et al., 2003) to identify putative outlier loci. A range of false discovery rate (FDR) values from 0.01 to 0.20 were evaluated based on preliminary testing and recommendations by Gondro et al. (2013). Departure from Hardy-Weinberg Equilibrium was tested for each locus using the software Arlequin v.3.5.2.2 (Schneider et al., 2000) using an exact test with 10,000 steps in the Markov Chain method and 100,000 dememorizations. Tests for linkage disequilibrium was performed using the software PLINK v.1.07 (Chang et al., 2015). Related individuals were detected using TrioML with COANCESTRY v. 1.0.1.7 (Wang, 2011). To determine cut-off r values for categorizing relatedness coefficients produced by the various algorithms into full-sib and less-related brackets, COANCESTRY was again used to simulate 200 dyads in each group based on the same allele frequencies as the original population samples. These values became the cutoffs for delineating full siblings in the empirical dataset. One individual from each of the full sibling pairs was removed for the final neutral data set. This new data set was used in all the subsequent analyses.

Allelic Diversity and Population Structuring
Allelic diversity indices including average observed (H o ) and unbiased expected heterozygosities corrected for population sample size (H nb ) were computed in Genetix v.4.05.2 (Belkhir et al., 1996), together with the inbreeding coefficient (F IS ) values using GENEPOP v.4.6 (Rousset, 2008). Pairwise F ST estimates were calculated using Arlequin v.3.5.1.3 followed by correction of significance levels for multiple pairwise tests (Excoffier and Lischer, 2010). An analysis of molecular variance (AMOVA) was performed in GenoDive v.3.0 using 10,000 permutations to estimate F-statistics in order to detect population genetic partitioning between sampling locations (Mengoni and Bazzicalupo, 2002). P values were adjusted by using the Bonferroni correction (Rice, 1989). Multiple AMOVAs were performed: i) one overall analysis including all sampling locations as separate populations, ii) Fiji versus all other sites grouped as one population, iii) insular (New Caledonia, Fiji) versus continental locations, and iv) Pacific Ocean versus Indian Ocean. Visualization of broad-scale population structure and population assignment was carried out by performing a Discriminant Analysis of Principal Components (DAPC) in the R package adegenet 1.4.2 (Jombart, 2008;Jombart and Ahmed, 2011). DAPC was carried out for all neutral loci. An α-score optimization was used to determine the number of principal components to retain. Alpha-score optimization indicated that the first nine PCs (42.7% of total variation) should be retained for analysis. Given the wide temporal range for sampling in Fiji The determination of the optimal number of clusters (k) was tested using both fastStructure (Raj et al., 2014) and Admixture (Alexander et al., 2009). Both software packages were written for large autosomal SNP data sets and unlike other programs can test the likelihood of a k of one. Admixture utilizes a crossvalidation procedure when determining the optimal k. While, for fastStructure the full model was selected with a seed of 100. For both programs, k ranging from 1 to 10 was run. Additionally, fastStructure included the python script chooseK.py which was run after the ten k scenarios were finished and the optimal k was selected. Once the optimal k was determined the probability of individual population assignments was calculated and visualized using both Admixture and AssignPOP (Chen et al., 2018).

Filtering and Genotypic Diversity
Genotyping by sequencing using DArT-Seq TM resulted in 59,601 SNPs prior to quality control filtering (Carson et al., 2014), which in turn resulted in 1,494 loci using a 99% call rate, read depth >7, and 2% MAF. Among the 1,494 SNPs that passed all quality control filters for the 215 individuals, zero SNPs were identified as outlier loci putatively under positive selection (FDR 1%). No loci were detected to be out of Hardy-Weinberg Equilibrium and in addition no loci were determined to be linked.
A total of 28 full sibling pairs were identified in COANCESTRY (Table 1), with Fiji having the most pairs (12) followed by eastern Australia (9). Most of the other sites only had one or two pairs, while no full sibling pairs were identified within the West Australian samples. After removal of one individual from each of the pairs, 187 samples were retained.
Overall, population-level indices of genetic diversity, including H o , H nb , and F IS , were variable across a priori populations ( Table 2). The H o were lowest within Papua New Guinea (0.128 ± 0.162) and South Africa (0.157 ± 0.168) and highest in bull sharks from New Caledonia (0.182 ± 0.188) and Indonesia (0.214 ± 0.166) ( Table 2). Inbreeding estimates were highest in specimens collected in eastern Australia (F IS = 0.355) and lowest in bull sharks from New Caledonia (F IS = −0.064) and Western Australia (F IS = −0.049).

Population Structuring
Pairwise F ST values were highest (0.034-0.044) between Fiji and all other sites, while values were moderate (0.016-0.024) for New Caledonia and all other continental sites ( Table 3). The AMOVA results allocated the vast majority of the variation within individuals and when sampling sites were treated separately without groupings, roughly 2.7% of the variation was found among a priori populations ( Table 4). The analysis showed significant variation among groups, when Fiji and New Caledonia were grouped as one versus a second group consisting of all continental sites (Table 4). This reflects the pairwise F ST values, in that Fiji and New Caledonia were significantly different from all other sites.
The DAPC scatterplot showed a clear distinction between Fiji and the other locations. New Caledonia was not differentiated from the continental locations, which may be a result of the small sample size (Figure 2). The analysis of temporal samples within Fiji resulted in scatter plots that placed the individuals compared, from three different years, overlapping to   an extent that supports stability in the distribution of the neutral genetic diversity within this insular population. Furthermore, when sampled years were analyzed separately with the rest of populations, again, sharks strongly overlapped within the Indo-Pacific and within Fiji. This confirmed the patterns of differentiation observed irrespective of the sampling year in Fiji and ruled out any temporal substructure (Supplementary Materials; Supplementary Figure S1). Both fastStructure and Admixture reported an optimal k of two, and we arbitrarily labeled them population A and B. For Admixture the crossvalidation inflection point was at k = 2, while for fastStructure the python script chooseK.py selected a k of 2 based on marginal likelihood. The population assignment analyses for both Admixture and AssignPOP assigned the same individuals to the same populations, with all continental locations and New Caledonia belonging to population A, and Fiji belonging to population B (Figure 3). The best performing assignment model in AssignPOP was Support Vector Machine (SVM) with an assignment accuracy of 0.999 (s.d. = 0.001) for population A and 1.0 for population B.

DISCUSSION
This is the first population genetic structure study conducted in bull sharks using SNP markers. The major finding of this study is that the insular bull shark population from Fiji is genetically distinct from its continental counterparts. Therefore, the data presented here supports the alternative hypothesis that biogeographic barriers to dispersal result in a heterogeneous distribution of the neutral genetic diversity between island bull sharks compared to populations sampled from continental shelves across the Indo-Pacific region.

Genetic Diversity
The sampling regime and SNP marker set used herein possessed the ability to resolve relatively low levels of differentiation enough to define bull shark population structure (Table 3). Based on SNP data, levels of observed heterozygosity (H o = 0.128-0.214) in bull sharks are comparable to other tropically distributed sharks, such as the gray reef shark (Carcharhinus amblyrhynchos) (0.139-0.312) (Momigliano et al., 2017) and the Galapagos shark (C. galapagensis) (0.188-0.193) (Pazmiño et al., 2017).

Population Genetic Structure
All five analyses examining population structuring and connectivity (F ST , AMOVA, DAPC, Admixture and PopAssign) failed to detect any barriers to gene flow among the six continental a priori populations (from the Pacific and Indian Oceans) in this study. Initially this might be surprising, given the large distances between the sampled populations. However, a recent study by Pirog et al. (2019) incorporating 25 microsatellites and using DAPC and fastStructure analyses also detected large ocean wide clusters including all western Atlantic populations in one and western Indian Ocean and Pacific sites in the other. Additionally, to date, no ocean basin subdivision has been detected for bull sharks (Karl et al., 2011;Tillett et al., 2012;Testerman, 2014). This is not surprising since tracking studies indicate bull sharks are capable of long-distance movements of up to 1,506 km along coastlines (Carlson et al., 2010) and can traverse open-ocean expanses of 2,000 km (Lea et al., 2015). However, genetic connectivity does not necessarily equate to demographic connectivity, as the exchange of only a few individual migrants per generation is enough to homogenize the gene pool even of distantly located populations (Waples, 1998). Low levels of migration among locations in turn are likely to be insufficient to compensate localized declines in abundance, but rare excursions can have important consequences for population viability by prevention of inbreeding. However, the same five analyses all detected differentiation between all continental populations and Fiji. Significant pairwise F ST values ranged from 0.034 to 0.044 between all locations (including New Caledonia) and Fiji. While values for New Caledonia were significant, they were smaller ranging between 0.016 and 0.024 and no significant differentiation was detected between Indonesia (Banda Aceh) and New Caledonia (F ST = 0.001); locations which are separated by thousands of kilometers of ocean. These smaller values and conflicting patterns may be a result of the small sample size (N = 8) of the New Caledonia population and may not reflect the true genetic structure. The DAPC results clearly differentiated Fiji from the rest of the populations, with New Caledonia clustering well within the continental cluster. Further support comes from Admixture and fastStructure, which both identified an optimal k of two with all continental locations (including New Caledonia) in one cluster and Fiji samples in the other. Furthermore, given a k of two, PopAssign was able to assign accurately the same individuals to each of the two populations that Admixture allocated. Thus, the corroboration of these independent analyses clearly indicates genetic heterogeneity and divergence between the insular population of Fiji and the other continental demes.
Bull sharks have a continuous distribution across South Africa, Indonesia and Australia (Ebert et al., 2013). The presence of continuous continental shelf waters off those landmasses and the absence of physical barriers along the shoreline habitats of Indo-Australia may facilitate dispersal in this coastal shark species. This is also suggested by the differences between South Africa and Western Australia but not with eastern Australia (the F ST between SAF and WAS is the same as between SAF and EAS). Thus, this lack of migration to Fiji from the surrounding continental populations suggests that bull shark migration across wide expenses of deep oceanic waters is an extremely rare event. The lack of differentiation between SAF samples and the rest of the continental sites despite the deep waters of the Indian Ocean is likely due to the continuous coastal habitats extending from Africa to southern Asia to the Sunda Shelf in Indonesia, allowing for the avoidance of deep water trenches and a near continuous migration route.

Management
The genetic differentiation of the Fiji bull sharks is a robust finding and makes this population of special interest due to its genetic and geographic isolation. Indonesia and Western Australia appear to be sites of transition for continental associated bull sharks. Dispersal limitations due to vast ocean expanses has implications for the degree of genetic (and demographic) connectivity among populations, and therefore the spatial scale of their management (Waples and Gaggiotti, 2006). For example, genetically isolated populations such as presumably the Fiji bull sharks could benefit from local protection of large, fecund individuals (Gwinn et al., 2015). Also, uncovering potential mating and nursery areas of high conservation value (Heupel et al., 2007) and fine-scale genetic connectivity is key to further population genetic studies . Consequently, a comparison between information generated by nuclear DNA (i.e., from both parents, relevant for kinship studies) with results derived from mtDNA (maternal, relevant for studies on philopatry) would help to define populations relevant for management. Additionally, the previously reported female philopatry implies that subpopulations within oceanic basins are likely to be independent demographically over ecological timescales, and that further work to identify discrete matrilineal populations along continuous coastlines throughout the species range is imperative to formulating effective management policies.

DATA AVAILABILITY STATEMENT
The data presented in this manuscript are permanently archived at the Dryad Digital Repository: doi: 10.5061/dryad.0vt4b8gwf.

ETHICS STATEMENT
The animal study was reviewed and approved by University Research Ethics Committee of the University of the South Pacific. All shipping procedures of shark tissue samples were conducted under the relevant import and export permits issued by the Australian Government, Department of Agriculture and Water Resources to Diversity Arrays Technology Pty Ltd, Canberra, Australia and the Ministry of

ACKNOWLEDGMENTS
Great thanks are owed to all assisting research colleagues, fishermen and the staff from Beqa Adventure Divers in Fiji. We authors thank Dr. Will White (CSIRO) for the bull shark tissue samples from Papua New Guinea via collaboration with a PNG-NFA ACIAR study, and the observers from commercial fishing vessels and artisanal surveys who collected the samples. We thank Crispen Wilson and Les Noble for sharing samples collected in South Africa and Indonesia. KG thanks the Ausbildungs-Stiftung fuer den Kanton Schwyz und die Bezirke See und Gaster (Kanton St. Gallen), the Shark Foundation Switzerland and Pacific Scholarship for Excellence in Research & Innovation of the University of the South Pacific. We thank Mike Neumann and Ian Campbell for continued advice during the study.