Genomic Assessment of Global Population Structure in a Highly Migratory and Habitat Versatile Apex Predator, the Tiger Shark (Galeocerdo cuvier)

Abstract

Understanding the population dynamics of highly mobile, widely distributed, oceanic sharks, many of which are overexploited, is necessary to aid their conservation management. We investigated the global population genomics of tiger sharks (Galeocerdo cuvier), a circumglobally distributed, apex predator displaying remarkable behavioral versatility in its diet, habitat use (near coastal, coral reef, pelagic), and individual movement patterns (spatially resident to long-distance migrations). We genotyped 242 tiger sharks from 10 globally distributed locations at more than 2000 single nucleotide polymorphisms. Although this species often conducts massive distance migrations, the data show strong genetic differentiation at both neutral (F_ST = 0.125–0.144) and candidate outlier loci (F_ST = 0.570–0.761) between western Atlantic and Indo-Pacific sharks, suggesting the potential for adaptation to the environments specific to these oceanic regions. Within these regions, there was mixed support for population differentiation between northern and southern hemispheres in the western Atlantic, and none for structure within the Indian Ocean. Notably, the results demonstrate a low level of population differentiation of tiger sharks from the remote Hawaiian archipelago compared with sharks from the Indian Ocean (F_ST = 0.003–0.005, P < 0.01). Given concerns about biodiversity loss and marine ecosystem impacts caused by overfishing of oceanic sharks in the midst of rapid environmental change, our results suggest it imperative that international fishery management prioritize conservation of the evolutionary potential of the highly genetically differentiated Atlantic and Indo-Pacific populations of this unique apex predator. Furthermore, we suggest targeted management attention to tiger sharks in the Hawaiian archipelago based on a precautionary biodiversity conservation perspective.

genetic adaptation, genetic connectivity, elasmobranch, single nucleotide polymorphism

Elucidating the population dynamics of globally distributed shark species contributes to understanding the ecology and evolutionary biology of these upper trophic level predators, many of which, particularly oceanic species, are at risk of extinction (Dulvy et al. 2014; Pacoureau et al. 2021). The population genetic dynamics of sharks, including genetic diversity, connectivity, phylogeography, and mating systems, have been investigated for phylogenetically diverse species, revealing that many shark species demonstrate sex-biased dispersal and/or female-philopatry (reviewed in Chapman et al. 2015), and population structuring associated with historical or extant biogeographic barriers (e.g., Ovenden et al. 2009; Veríssimo et al. 2010; Mendonça et al. 2011; Daly-Engel et al. 2012; Portnoy et al. 2015; Bernard et al. 2017; Momigliano et al. 2017; Pazmiño et al. 2018; Dimens et al. 2019).

Although notable progress has been made in assessing shark population genetics over the past decade, assessments for the globally distributed, large-bodied, highly mobile, oceanic (pelagic habitat using) sharks are still limited. These species constitute a major component of exploited and incidentally captured sharks in pelagic fisheries and can traverse large distances across open ocean habitats (Block et al. 2011; Pacoureau et al. 2021), thus predicting the potential for long-distance gene-flow. With 2 exceptions (Junge et al. 2019, Carcharhinus obscurus; Kraft et al. 2020, Carcharhinus falciformis), knowledge of the population genetics of pelagic sharks has thus far been derived from use of traditional molecular markers such as sections of mitochondrial DNA sequences and relatively few nuclear microsatellites (e.g., Hoelzel et al. 2006; Vignaud et al. 2014; Oñate-González et al. 2015; Camargo et al. 2016; Bernard et al. 2016; Veríssimo et al. 2017; Bailleul et al. 2018; Corrigan et al. 2018; Pirog et al. 2019). Since these highly migratory, pelagic sharks travel thousands of kilometer distances (Block et al. 2011; Kohler and Turner 2019), even small amounts of gene flow may reduce the evolution of allele frequency differences between demographically independent populations (Slatkin 1987; Waples and Gaggiotti 2006). Consequently, the relatively small numbers of presumed neutral, traditional markers that have been used may not possess the necessary power or precision to detect subtle genetic differences occurring across smaller spatial scales.

Recent genomic approaches, including restriction site-associated DNA sequencing (RADseq; reviewed in Andrews et al. 2016), have allowed the development of thousands of neutral and non-neutral molecular markers across genomes (i.e., single nucleotide polymorphisms [SNPs]) from non-model organisms, allowing increased robustness for estimation of population genetic parameters and identification of population-level adaptive processes (Allendorf et al. 2010; Davey and Blaxter 2011; Crawford and Oleksiak 2016). Such large numbers of genomic-scale markers can provide the higher power and precision needed for identifying demographically independent genetic lineages that may have otherwise been undetected or overlooked using traditional markers (Kelley et al. 2016). Although the application of genome-wide SNP markers to assess the population genomics of widely distributed pelagic fishes is only recent, this approach is already revealing finer levels of population structuring than detected previously in the same species using traditional markers (e.g., Barth et al. 2017; Mullins et al. 2018; Pecoraro et al. 2018; Kraft et al. 2020; Mamoozadeh et al. 2020).

The tiger shark (Galeocerdo cuvier) is a large-bodied (commonly 3.5–4.5 m as adults), circumglobally distributed, highly mobile species inhabiting warm temperate and tropical waters, and an apex predator of conservation concern in parts of its range [“Near Threatened” on the IUCN Red List of Threatened Species (Ferreira and Simpfendorfer 2019)]. Across its global distribution, the tiger shark shows tremendous versatility in using very diverse habitats ranging from shallow coral reef to near coastal to distant pelagic environments. Telemetry studies have shown that tiger sharks exhibit complex habitat use and heterogeneous movement patterns including partial migration (Papastamatiou et al. 2013; Lea et al. 2015, 2018), site-fidelity (Meyer et al. 2010; Hammerschlag et al. 2012; Fitzpatrick et al. 2012), and very long-distance movements (sometimes repeated and seasonal migrations) of thousands of kilometers (Heithaus et al. 2007; Hammerschlag et al. 2012; Holmes et al. 2014; Ferreira et al. 2015; Lea et al. 2015, 2018; Ajemian et al. 2020). Thus, predicting a priori the levels of gene flow or the demographic connections linking tiger sharks across both large and small geographic scales remains difficult. Previous studies on tiger shark population genetics have identified various levels of genetic heterogeneity (matrilineal and/or biparental) across globally distributed populations, with the bulk of variation found between tiger sharks occupying different ocean basins (Bernard et al. 2016; Holmes et al. 2017; Carmo et al. 2019; Pirog et al. 2019). Within ocean basins, significant matrilineal population differentiation has been observed (Bernard et al. 2016; Carmo et al. 2019), while only more subtle, and in some cases conflicting, signals of nuclear microsatellite-based regional differentiation have been found (Bernard et al. 2016; Holmes et al. 2017; Pirog et al. 2019; Sort et al. 2021). These inconsistent findings among studies suggest that examination using additional and higher resolution, neutral and non-neutral nuclear markers is needed to obtain a complementary, and potentially further refined, view of the population genetics of the tiger shark. Here, we employ genome scale SNP markers to address the broad objective of obtaining a higher resolution view of global population genetic heterogeneity and diversity in this globally distributed apex predator that often migrates massive distances.

Methods

Sample Collection and DNA Extraction

Tiger shark fin clip or muscle tissues were collected from 285 individuals captured in fisheries or beach protective nets as well as during ecological studies between 1999 and 2013 (where collection date was available) from 14 globally distributed locations (Supplementary Table S1). All tissue samples were stored in 95–99% ethanol or a salt-saturated, 20% DMSO, 0.25 M EDTA solution, and genomic DNA (gDNA) was extracted using the Qiagen DNeasy Kit (QIAGEN, Inc., Valencia, CA) following manufacturer’s instructions.

Genotyping-by-Sequencing and SNP Data Filtering

Genotyping-by-sequencing (GBS) library preparation and sequencing were performed by the Genomic Diversity Facility at Cornell University (NY, USA) per protocols outlined in Elshire et al. (2011). Tiger shark genome complexity was reduced via digestion of gDNA using the restriction enzyme EcoT22I, and 100 bp of product was single-end (1X100SE) sequenced across 3 lanes of an Illumina Hiseq (Illumina, Inc.). Generated FASTQ files were analyzed using the Tassel 3.0.166 (Bradbury et al. 2007; Glaubitz et al. 2014) GBS non-reference Universal Network-Enabled Analysis Kit (UNEAK) pipeline (Lu et al. 2013). The resultant SNP dataset was quality filtered using VCFtools 0.1.13 (Danecek et al. 2011) with: 1) genotypes possessing less than 5× coverage coded as missing data; 2) loci filtered to retain only those tags genotyped in a minimum of 70% of sharks and possessing a minor allele frequency of at least 5%; 3) sharks filtered to retain only those samples possessing less than 30% missing data; 4) sites filtered to retain only those SNPs possessing a mean depth of 10× coverage across sharks; and 5) to reduce the possibility of paralogous sequences, sites were filtered to discard those loci whose mean read depth across individuals was more than 2× that of the global mean (mean of individual locus means) across sites. In addition, we tested for duplicate individuals using the R (R Core Team 2018) package strataG 2.0.2 (Archer 2017) and the function “dupGenoytpes”; if any pair of individuals possessed a shared genotype exceeding 85%, then the individual with the highest rate of missing data was discarded. Following quality filtering, we used the R package SNPRelate 1.16.0 (Zheng et al. 2012) to assess the degree of relatedness among genotyped individuals using the Manichaikul et al. (2010) framework of robust kinship-based inference for genome-wide association studies (KING-robust), which expressly allows for the existence of population structure when estimating relatedness among individuals. Pairs of individuals possessing kinship coefficients exceeding 0.125 were identified (at minimum) as potential second-degree relatives. As the biases associated with inclusion of related individuals within population genetic studies is a matter of debate (Peterman et al. 2016; Waples and Anderson 2017; Wang 2018), all downstream population genetic analyses were conducted on 2 tiger shark sample sets: 1) all tiger shark individuals (hereafter, “tiger-complete”), and 2) a pruned sample set obtained after removing a single shark individual from each related pair (hereafter, “tiger-kinship”); (see Results).

For each resultant sample set (i.e., tiger-complete and tiger-kinship), input files for downstream analyses were formatted as needed by the program PGDSpider 2.1.1.5 (Lischer and Excoffier 2012). Loci were tested for conformation to Hardy–Weinberg equilibrium (HWE) within each of the sampling locations (hereafter, a priori “subpopulations”; Figure 1 and see Results for description of final subpopulation groupings) using the “hw.test” function in the R package pegas 0.11 (Paradis 2010). Those loci deviating from HWE at P <0.01 in more than 2 subpopulations were discarded. After HWE filtering, PLINK 1.90b6.7 (Purcell et al. 2007) was used to prune the overall datasets to ensure a coefficient of disequilibrium (r²) ≤0.8 between locus pairs.

Figure 1.

Map of tiger shark (Galeocerdo cuvier) global sampling locations: BM, Bermuda; FLE, Florida East Coast; GOM, Gulf of Mexico; USVI, United States Virgin Islands; VZ, Venezuela-French Guiana; WSA, western South Atlantic; WSI, western South Indian; AS, Andaman Sea; ESI, eastern South Indian; CP, Central Pacific. Amazon Barrier marks the outflow of the Amazon River. The a priori subpopulations used in our analyses are shown as either a solitary black dot (i.e., ESI, AS, WSI, BM, USVI) which represents general sample collection site only, or regionally grouped black dots shown within the dashed rectangles (i.e., CP, GOM, FLE, VZ, WSA). Final shark sample sizes (tiger-complete sample set) used given in brackets.

Genetic Diversity and Population Structure Analyses

Molecular diversity estimates—observed (H_O) and expected heterozygosity (H_E), allelic richness (Ar), and the inbreeding coefficient (F_IS)—were estimated at ocean basin (western Atlantic Ocean, Indian Ocean, and Pacific Ocean) and subpopulation levels using the function “divBasic” as implemented in the R package diveRsity 1.9.90 (Keenan et al. 2013). Population genetic structure was estimated using population- and individual-based analytical approaches: A hierarchical analysis of molecular variance (AMOVA) was performed using the R package poppr 2.8.3 (Kamvar et al. 2014, 2015) with subpopulations grouped according to ocean basin. Significance was determined using 1000 permutations. Pairwise values of differentiation (F_ST; Weir and Cockerham 1984) were estimated among: 1) the individual subpopulations, and 2) pooled ocean-basin samples (i.e., western Atlantic vs. Indian vs. Pacific), using the R package strataG (“pairwiseTest” function) and the “pairwise.WCfst” function in the R package hierfstat 0.5–7 (Goudet and Jombart 2018). Significance of the strataG estimations were determined using 1000 permutations and adjusted using the false discovery rate (FDR; Benjamini and Hochberg 1995) and “p.adjust” function in the R package stats. Significance of hierfstat pairwise estimates of differentiation were inferred through the generation of 95% confidence intervals and bootstrapping 1000 iterations over loci using the function “boot.ppfst.”

Individual-based genetic cluster analyses were conducted using the program Admixture (Alexander et al. 2009) and the multivariate method “discriminant analysis of principal components” (DAPC; Jombart et al. 2010). First, Admixture was run implementing default settings and the cross-validation procedure (CV) (assuming K = 1–10) to determine K. Clustering outcomes at the K associated with the lowest CV were plotted and visualized in R. Second, DAPCs were performed using the R package adegenet 2.1.3 (Jombart 2008; Jombart and Ahmed 2011) and assumed each individual shark’s sampling location as its a priori defined subpopulation (or group membership). For the DAPCs, the optimal number of PCs to retain for analysis (lowest root mean squared error; RMSE) was identified using adegenet’s cross-validation function (“xvalDapc”), assuming a maximum number of 100 PCs, a 90% training dataset, and 1000 replicates. Cross-validation was performed at least twice per dataset to ensure consistency of results.

SNP Outlier Identification and Analysis

Given the potential for false positives with outlier detection methods (Narum and Hess 2011; Lotterhos and Whitlock 2015), 4 distinct methods were adopted to identify candidate outlier SNPs, that is, those loci deviating from selectively neutral expectations and potential candidates for selection: First, Excoffier et al.’s (2009) hierarchical approach implemented in the program Arlequin 3.5 (Excoffier and Lischer 2010) was used, as hierarchical population genetic structure is known to occur globally for tiger sharks (Bernard et al. 2016; Pirog et al. 2019; Carmo et al. 2019). Tiger shark subpopulations were grouped according to their respective ocean basins (western Atlantic, Indian, and Pacific), and FDIST (Beaumont and Nichols 1996) analyses within Arlequin were performed employing 50 000 coalescent simulations with 100 demes simulated per group, adopting a P value of 0.01. Second, BayeScan 2.1 (Foll and Gaggiotti 2008), which assumes an island migration model, was run using default settings and assuming a prior odds for the neutral model of 100, and an FDR of 5%. Third, the R package pcadapt 4.3.3 (Luu et al. 2017) was used to identify candidate loci for selection using a population PCA-based approach. The number of principal components (K) was selected for analysis using a scree plot and Cattell’s rule, assuming a maximum K of 20. Candidate outliers were identified where their respective q values were less than an alpha of 0.05. And fourth, the R package OutFLANK (Whitlock and Lotterhos 2015) was used, assuming 10 final subpopulations (see Results) and using default settings. Loci were deemed candidates for selection using OutFLANK if its resultant q value was less than 0.05. Herein, a tiger shark SNP locus was deemed a candidate outlier only if it was identified by at least 2 of the 4 methods (hereafter, “tiger-complete-outlier” and “tiger-kinship-outlier” sample sets).

To assess if genomic selection was influencing patterns of genetic structure, pairwise estimates of genetic differentiation (F_ST) and DAPC genetic analyses were repeated for both the outlier SNP loci dataset and the putatively neutral SNP loci dataset (i.e., where outlier SNP loci have been removed from the overall dataset; hereafter, “tiger-complete-neutral” and “tiger-kinship-neutral”).

Results

GBS library sequencing of 285 shark samples generated a total of ~818 million reads, with ~677 million passing initial quality standards. Each shark yielded at minimum 637 000 barcoded reads. Following the GBS UNEAK SNP discovery pipeline, 104 031 SNPs were identified across 285 genotyped sharks, and among these individuals, 256 were retained for downstream SNP quality filtering (Supplementary Table S1). Quality filtering yielded a dataset containing 2478 SNPs genotyped across 242 tiger sharks collected from 10 globally distributed subpopulations (Figure 1; statistics after each filtering step provided in Supplementary Table S2): that included 6 sampling locations within the western Atlantic Ocean (Bermuda, n = 19; Florida East Coast, n = 32; Gulf of Mexico, n = 27; United States Virgin Islands, n = 30; Venezuela-French Guiana, n = 11, western South Atlantic [south of Amazon River outflow], n = 8), 3 from the Indian Ocean (Andaman Sea, n = 24; eastern South Indian, n = 30; western South Indian, n = 29), and a single site from within the Pacific Ocean (Central Pacific [CP], n = 32). Within these 242 sharks, SNPRelate identified 11 putative related pairs of individuals with coefficients of relatedness ranging from 0.131 to 0.294 (Supplementary Table S3); each pair of related individuals was found to have been sampled in the same general location and captured within a 3-year period. The 2 final filtered datasets used for analyses consisted of: (1) tiger-complete (n = 242; 2478 SNPs), and (2) tiger-kinship (n = 231; 2478 SNPs). While all downstream population genetic analyses were performed in duplicate using both these datasets, in the main text we report only the results from the tiger-complete dataset except where noteworthy discrepancies or similarities were found between the 2 datasets. The tiger-kinship dataset analyses results are provided as Supplementary Materials S1 and associated tables and figures.

Following HWE testing, 59 loci were discarded from the tiger-complete dataset due to significant deviations at P <0.01 in more than 2 subpopulations and an additional 413 loci were pruned by PLINK through LD testing, generating a final tiger-complete dataset containing 2006 SNPs genotyped across 242 sharks.

Genetic Diversity and Population Structure

Diversity estimates were highly similar across the 10 subpopulations and among pooled ocean-level samples comparisons (Table 1), although the allelic richness was slightly higher overall in western Atlantic relative to Indian or Pacific Ocean tiger sharks.

Table 1.

Tiger shark global population genetic diversity statistics for 2006 SNPs genotyped across 242 sharks (tiger-complete sample set)

Location	n	Ar	HO	HE	FIS
Western Atlantic Ocean	127	1.95	0.25	0.26	0.00
BM	19	1.74	0.25	0.25	−0.03
FLE	32	1.76	0.25	0.25	−0.01
GOM	27	1.75	0.25	0.25	−0.01
USVI	30	1.76	0.26	0.25	−0.01
VZ	11	1.72	0.26	0.24	−0.06
WSA	8	1.69	0.26	0.24	−0.06
Indian Ocean	83	1.84	0.25	0.25	0.01
WSI	29	1.73	0.25	0.25	−0.01
AS	24	1.72	0.25	0.25	−0.01
ESI	30	1.73	0.25	0.25	0.00
Pacific Ocean	32	1.83	0.25	0.25	0.00
CP	32	1.73	0.25	0.25	0.00

n, sample size; Ar, allelic richness; H_O, observed heterozygosity; H_E, expected heterozygosity; F_IS, inbreeding coefficient. Location (subpopulation) abbreviations: BM, Bermuda; FLE, Florida East Coast; GOM, Gulf of Mexico; USVI, US Virgin Islands; VZ, Venezuela-French Guiana; WSA, western South Atlantic Ocean; WSI, western South Indian Ocean; AS, Andaman Sea; ESI, eastern South Indian Ocean; CP, Central Pacific Ocean.

All population structure analyses (AMOVA, F_ST, Admixture, and DAPCs) consistently identified strong genetic differentiation between western Atlantic and Indo-Pacific tiger sharks (Tables 2 and 3; Supplementary Table S4; Figures 2 and 3). The hierarchical AMOVA found significant differentiation between western Atlantic, Indian, and Pacific Ocean tiger sharks (16.29% variance; P = 0.003), as well as between subpopulations within ocean basin (P = 0.001) although the total variance explained by the latter was minimal (0.12%) (Table 2). Notably, upon hierarchical AMOVA analysis of the tiger-kinship data, the significant differentiation between subpopulations within ocean basin disappeared (0.03% variance, P = 0.176; Supplementary Table S5), indicating that any genetic differentiation within oceans is likely minimal. Akin to the AMOVA results, pairwise values of genetic differentiation (F_ST) between western Atlantic and Indian-Pacific subpopulations were large (F_ST = 0.180–0.193) and significant with both analysis methods (i.e., P values and CIs; Table 3), as were estimates when samples within oceans were pooled and then compared across ocean basins (Supplementary Table S4). In addition, very low, but statistically significant genetic differentiation between tiger sharks from the CP and those from all 3 Indian Ocean collection sites (i.e., WSI, AS, and ESI; F_ST = 0.003–0.005; Section II in Table 3) was identified via both permutation testing and bootstrapping. These findings were also consistent when samples from the Indian Ocean subpopulations were pooled and compared to those sharks from the lone Pacific Ocean collection site (Supplementary Table S4). Within ocean basins, few concordant indications of pairwise genetic population structure were found. While the within western Atlantic tiger-complete SNP data did yield several significant pairwise comparisons between subpopulations (Section I in Table 3), these intra-basin findings were largely inconsistent across statistical methods (i.e., P values and CIs), and also when compared to the pairwise results derived from the tiger-kinship data.

Table 2.

Tiger shark AMOVA assuming 3 hierarchical groupings [pops: (1) western Atlantic Ocean (WATL), (2) Indian Ocean (IND), (3) Pacific Ocean (PAC)] for 242 sharks genotyped at 2006 SNPs and 10 a priori defined subpopulations (subpops; tiger-complete sample set)

Sources of variation	DF	Sum Sq.	% variance	ϕ-statistic	P value
Between pops (WATL, IND, PAC)	2	14780.72	16.29	0.163	0.003
Between subpops within pops [WATL (BM, FLE, GOM, USVI, VZ, WSA), IND (WSI, AS, ESI), PAC (CP)]	7	1921.15	0.12	0.001	0.001
Between samples within subpops	232	59899.05	0.75	0.009	0.210
Within samples	242	61364.96	82.84	0.172	0.001
Total	483	137965.87	100.00	—	—

DF, degrees of freedom; Sum Sq., sum of squares; % variance, percent variance. Location (subpopulation) abbreviations: BM, Bermuda; FLE, Florida East Coast; GOM, Gulf of Mexico; USVI, US Virgin Islands; VZ, Venezuela-French Guiana; WSA, western South Atlantic Ocean; WSI, western South Indian Ocean; AS, Andaman Sea; ESI, eastern South Indian Ocean; CP, Central Pacific Ocean.

Table 3.

Global population-level pairwise values of genetic differentiation (F_ST) for tiger sharks genotyped at SNPs

Subpopulation comparison	FST 2006 SNPs, n = 242 tiger-complete (all)	FST 108 SNPs, n = 242 tiger-complete (outliers)	FST 1898 SNPs, n = 242 tiger-complete (neutral)
Section I: western Atlantic Ocean
BM vs. FLE	0.001	0.002	0.001
BM vs. GOM	0.000	−0.001	0.000
BM vs. USVI	0.001	−0.003	0.001
BM vs. VZ	*0.003	0.006	0.003
BM vs. WSA	0.003	−0.006	*0.004
FLE vs. GOM	0.000	0.000	−0.001
FLE vs. USVI	0.001	−0.002	0.001
FLE vs. VZ	0.002	*0.010	0.001
FLE vs. WSA	0.002	−0.009	0.003
GOM vs. USVI	0.000	−0.003	0.001
GOM vs. VZ	0.001	0.001	0.001
GOM vs. WSA	0.002	−0.011	*0.004
USVI vs. VZ	*0.003	0.006	*0.003
USVI vs. WSA	*0.004	−0.011	*0.006
VZ vs. WSA	*0.005	0.008	*0.005
Section II: Indian and Pacific Oceans
WSI vs. AS	0.001	−0.011	*0.002
WSI vs. ESI	0.001	−0.001	0.001
WSI vs. CP	*0.005	0.009	*0.005
AS vs. ESI	0.001	−0.001	0.001
AS vs. CP	*0.004	0.012	*0.004
ESI vs. CP	*0.003	−0.001	*0.003
Section III: western Atlantic Ocean vs. Indian and Pacific Oceans
BM vs. WSI	*0.187	*0.650	*0.136
BM vs. AS	*0.191	*0.622	*0.142
BM vs. ESI	*0.189	*0.647	*0.139
BM vs. CP	*0.188	*0.658	*0.137
FLE vs. WSI	*0.187	*0.593	*0.138
FLE vs. AS	*0.191	*0.570	*0.144
FLE vs. ESI	*0.189	*0.591	*0.141
FLE vs. CP	*0.187	*0.601	*0.138
GOM vs. WSI	*0.190	*0.623	*0.139
GOM vs. AS	*0.193	*0.598	*0.144
GOM vs. ESI	*0.192	*0.621	*0.141
GOM vs. CP	*0.190	*0.631	*0.139
USVI vs. WSI	*0.186	*0.604	*0.137
USVI vs. AS	*0.191	*0.579	*0.143
USVI vs. ESI	*0.188	*0.602	*0.140
USVI vs. CP	*0.187	*0.611	*0.137
VZ vs. WSI	*0.189	*0.740	*0.131
VZ vs. AS	*0.192	*0.712	*0.136
VZ vs. ESI	*0.191	*0.736	*0.135
VZ vs. CP	*0.189	*0.747	*0.132
WSA vs. WSI	*0.180	*0.755	*0.125
WSA vs. AS	*0.186	*0.727	*0.132
WSA vs. ESI	*0.184	*0.749	*0.129
WSA vs. CP	*0.182	*0.761	*0.127

Values in bold indicate statistical significance at P <0.05 (strataG) after false discovery rate (FDR) correction; values in italics indicate non-adjusted statistical significance at P <0.01 (strataG); asterisk (*) indicates values with 95% confidence intervals that do not overlap zero (hierfstat). Location (subpopulation) abbreviations: BM, Bermuda; FLE, Florida East Coast; GOM, Gulf of Mexico; USVI, US Virgin Islands; VZ, Venezuela-French Guiana; WSA, western South Atlantic Ocean; WSI, western South Indian Ocean; AS, Andaman Sea; ESI, eastern South Indian Ocean; CP, Central Pacific Ocean.

Figure 2.

Clustering outcome (biplot) of Admixture analysis reporting proportional ancestry (y axis) of 242 tiger sharks (x axis) sampled globally and genotyped at 2006 SNP loci (tiger-complete sample set). Location (subpopulation) abbreviations: BM, Bermuda; FLE, Florida East Coast; GOM, Gulf of Mexico; USVI, US Virgin Islands; VZ, Venezuela-French Guiana; WSA, western South Atlantic Ocean; WSI, western South Indian Ocean; AS, Andaman Sea; ESI, eastern South Indian Ocean; CP, Central Pacific Ocean.

Figure 3.

Biplots of the DAPCs performed for globally sampled tiger sharks (tiger-complete sample set; n = 242) genotyped at (a) 2006 SNP loci, (b) 108 candidate outlier SNPs, and (c) 1898 putative neutral SNPs. Analyses were performed by retaining the first (a) 20 PCs (30.2% of the conserved variance), (b) 30 PCs (83.3% of the conserved variance), and (c) 40 PCs (37.7% of the conserved variance). Inertia ellipses provide a graphical summary of the cloud of points. Location (subpopulation) abbreviations: IND, Indian Ocean; AS, Andaman Sea; ESI, eastern South Indian Ocean; WSI, western South Indian Ocean; PAC, Pacific Ocean; CP, Central Pacific Ocean; WATL, western Atlantic Ocean; BM, Bermuda; FLE, Florida East Coast; GOM, Gulf of Mexico; USVI, US Virgin Islands; VZ, Venezuela-French Guiana; WSA, western South Atlantic Ocean.

Admixture’s CV method identified the optimal value of K for the global tiger-complete shark SNP dataset as K = 2 (CV = 0.456; Supplementary Figure S1a). Visualization of clustering outcomes yielded 2 distinct groups consisting of: 1) western Atlantic, and 2) Indo-Pacific tiger sharks, with some minor degree of admixture between groups (Figure 2). Both hierarchical Admixture analyses (i.e., western Atlantic samples only or Indo-Pacific samples only) identified K = 1 as the most appropriate number of genetic clusters (Supplementary Figures S1b,c). For the DAPC, the cross-validation procedure for the global tiger-complete SNP dataset showed that the inclusion of 20 PCs in the analysis provided the lowest RMSE. Visualization of the DAPC results (after retaining 20 PCs; 30.2% of the conserved variance) showed strong differentiation of western Atlantic and Indo-Pacific tiger sharks, with some mild separation of CP sharks from all Indian Ocean tiger sharks when sample group membership (i.e., clustering) was identified a priori (Figure 3a).

SNP Outlier Identification and Analysis

Across all 4 outlier detection methods, up to 372 unique SNP loci were identified as potential candidates undergoing either diversifying or balancing selection (Arlequin: 8 loci; BayeScan: 372 loci; pcadapt: 108; OutFlank: 0 loci). Of those candidate markers, 108 were identified as undergoing divergent selection by both BayeScan and pcadapt. No other outlier detection method (i.e., Arlequin and OutFlank) identified candidate markers undergoing diversifying selection. Within ocean basins, although some significant pairwise genetic differentiation was found upon analysis of the tiger-complete-outlier dataset, these results were largely inconsistent (i.e., P values and CIs) across statistical methods (Table 3). DAPC outlier SNP results mirrored pairwise measures showing strong differentiation between western Atlantic versus Indian and Pacific Ocean sharks, and little population substructuring within these broad geographic regions (Figure 3b). Analysis of the strictly neutral SNP dataset (tiger-complete-neutral) yielded highly similar patterns of genetic clustering as the overall SNP dataset (Table 3; Figure 3c).

Discussion

Here, we provide the first genome-wide, SNP marker-based perspective on the global population genetic dynamics of the tiger shark, a highly migratory, immense-distance travelling, marine apex predator. The results reveal strong genetic differentiation between western Atlantic versus Indian and Pacific Ocean tiger sharks, a finding consistent with that of the nuclear microsatellite genetic assessment by Bernard et al. (2016). This direct comparison between studies is strengthened by the fact that our SNP study used approximately the same global distribution sample set used by Bernard et al. (2016). In this context, we note that the intra- and inter-ocean basin pairwise SNP F_ST values from the current study largely overlap with the corresponding microsatellite pairwise F_ST values in Bernard et al. (2016), providing further corroboration of the global population differentiation magnitude and dynamics in tiger sharks. In contrast, the intra-ocean basin, mitochondrial control region sequence F_ST values (from Bernard et al. 2016) are over 2 orders of magnitude higher than the corresponding SNP-based F_ST values. The consistently much higher mitochondrial than SNP F_ST values in a species known to migrate vast distances is concordant with female tiger shark philopatric behavior proposed by Bernard et al. (2016) and inferred from the high mitochondrial genetic differentiation in the western Atlantic reported by Carmo et al. (2019). Our SNP-based findings were also concordant with the high microsatellite differentiation between western Atlantic and Indo-Pacific tiger sharks reported by Holmes et al. (2017) and Pirog et al. (2019), although their studies were primarily Indo-Pacific focused but included 8 tiger shark samples from a location in the western South Atlantic.

Although the GBS methods used here sampled only a small fraction of the tiger shark’s nuclear genome, ~4.4% of the total SNPs surveyed were identified as candidate outliers (i.e., 108 SNPs common to 2 outlier detection methods). Subsequent population genetic analysis of these strictly candidate outlier SNP datasets indicated significant genetic differentiation between sharks from the western Atlantic and Indo-Pacific Oceans (Table 3 and Figure 3b), supporting the hypothesis of potentially adaptive genetic differences between tiger sharks inhabiting distinct ocean basins. Although little to no evidence of genetic differentiation was found among subpopulations within ocean basins according to these same candidate outlier SNP datasets, these data cannot be used to exclude the possibility that adaptive differences within basins exist. We note, however, that since outlier detection tests have the potential for false positives (Narum and Hess 2011; Lotterhos and Whitlock 2015); further study using whole genome sequence analysis is needed to test the hypothesis that adaptive genetic differences exist between these 2 oceanic populations of tiger sharks.

The handful of investigations into the population genomics of other shark species where outlier SNPs have been identified have also yielded evidence of potentially adaptive differentiation among subpopulations (e.g., Momigliano et al. 2017 [Carcharhinus amblyrhynchos]; Pazmiño et al. 2018 [Carcharhinus galapagensis]; Marie et al. 2019 [Sphyrna lewini]), including indications of latitude-associated selection in Sphyrna tiburo (Portnoy et al. 2015). These findings underscore the potential that even highly mobile shark species possessing expansive distributions across heterogenous environments may demonstrate local adaptation to specific habitats (Momigliano et al. 2017), including in some cases where little to no differentiation is seen at neutral markers [see Pazmiño et al. (2018) and Xu et al. (2017) for examples of sharks and marine teleosts, respectively]. As a complete and well annotated reference genome for the phylogenetically and physiologically enigmatic tiger shark (Castro et al. 2016; Swift et al. 2016) does not currently exist, we are unable to map the putative outlier SNP loci to a reference genome, and hence the underlying functions of any genes associated with these SNPs remains to be elucidated.

Population Connectivity—Western Atlantic

Western North Atlantic adult tiger shark movements include very long-distance seasonal migrations (Hammerschlag et al. 2012; Lea et al. 2015, 2018; Ajemian et al. 2020), with some population segments showing repeated migrations linking high latitude (~41°N) western North Atlantic temperate oceanic habitats to lower latitude (~10°N) warm water regions (Lea et al. 2015, 2018; unpublished data), indicating the potential for high genetic connectivity throughout at least the western North Atlantic. Furthermore, the broad temperature tolerance of tiger sharks (6.3–28 °C; Castro 2011; Vaudo et al. 2014) and the capability of this species to travel huge distances (even across the Atlantic; Afonso et al. 2017), makes it feasible that tiger sharks are capable of at least occasional demographic linkages between Atlantic northern and southern hemispheric subpopulations.

Despite the occurrence of such extensive migrations in this region, the 2 previous mitochondrial DNA surveys of western Atlantic tiger sharks have indicated the presence of significant matrilineal population structure, with the bulk of genetic differentiation occurring between sharks sampled from the southeastern USA Atlantic/Gulf of Mexico/The Bahamas/Panama waters and sharks occupying the more southern Brazilian coastline (Bernard et al. 2016; Carmo et al. 2019). In fact, Carmo et al. (2019) offered some evidence of multiple matrilineal subpopulations within Brazilian waters, and akin to Bernard et al. (2016), these authors identified extremely high levels of mitochondrial diversity within Brazilian tiger sharks. These findings complement the telemetry studies showing that some tiger sharks occupy relatively restricted activity spaces (Hammerschlag et al. 2012; Hazin et al. 2013; Vaudo et al. 2014), which coincides with the presence of some regionally distinct matrilineal subpopulations of these large predators.

The only other previous nuclear marker-based (survey of 10 microsatellite loci), population genetic study of tiger sharks in their western north to south Atlantic range (Bernard et al. 2016) found evidence of shallow (albeit inconsistent across analysis methods) population genetic differentiation between western north and south Atlantic tiger sharks. This signal of slight differentiation warranted further investigation (Bernard et al. 2016). Herein, a putatively higher resolution, complementary analysis of 2000+ SNP loci (all, neutral or candidate outlier loci) still provided mixed (i.e., inconsistent across datasets) evidence of genetic differentiation among western Atlantic tiger shark subpopulations, even between North and South of the Amazon Barrier and/or the Equator. Where significant pairwise differentiation was detected between subpopulations, the SNP results proffered very small magnitudes of differentiation (F_ST) similar to that seen by Bernard et al.’s (2016) microsatellite study. Per Bernard et al. (2016), we performed a post hoc analysis whereby tiger shark samples were pooled by hemisphere and tested for pairwise differentiation; while Bernard et al. (2016) found low, but statistically significant differentiation between North and South hemispheric collections (F_ST = 0.006, P = 0.050), our SNP hemispheric pooled results were equivocal (tiger-complete, F_ST = 0.003, P = 0.03; tiger-kinship, F_ST = 0.000, P = 0.57).

The wide-ranging distribution of these sharks within the western Atlantic (northern USA to Uruguay) and their highly individualistic and complex migratory behavior, placed in the context of 1) Carmo et al.’s (2019) findings of high matrilineal population structure along the Brazilian coastline, 2) Brazil sampled tiger sharks demonstrating an epicenter of high mitochondrial diversity (Bernard et al. 2016; Carmo et al. 2019), and 3) the equivocal SNP genotyping results obtained herein, suggest further genomic examination of the patterns of tiger shark gene flow within the western Atlantic, including larger sample sizes from WSA, may be beneficial. Given the heterogeneous oceanographic conditions spanning the extensive geographic range of tiger sharks in the western Atlantic and the discovery of a reasonably high number of candidate outlier SNP loci (108 outlier loci common to 2 methods, ~4.4% of total loci), the potential remains that further examination may reveal latitudinally adapted lineages in this range. We also note that no genetic study has yet included tiger sharks from the eastern Atlantic distribution of this species because of the difficulty of obtaining such samples; this broader geographic analysis is needed to obtain a through view of the population genetic dynamics of tiger sharks in the Atlantic.

Population Connectivity—Indo-Pacific

Genetic homogeneity across large portions of the tiger shark’s Indo-Pacific distribution has been consistently inferred across the 4 published nuclear microsatellite DNA studies (Bernard et al. 2016; Holmes et al. 2017; Pirog et al. 2019; Sort et al. 2021). Still, these studies showed some conflicting patterns of population structure, notably between tiger sharks from the easternmost Pacific location sampled in all 3 studies (i.e., the CP–Hawaii), and the other more westerly Indo-Pacific collection sites. The data of Bernard et al. (2016) suggested low but significant nuclear differentiation of Hawaii sampled tiger sharks (n = 65 animals) from the western Pacific and Indian Oceans, whereas the data of Holmes et al. (2017) and Pirog et al. (2019) [both studies used the same samples (n = 21) from Hawaii] did not support this differentiation. In the current study, we SNP genotyped a subset of Bernard et al.’s (2016) Hawaiian tiger sharks and found pairwise estimates of differentiation between Hawaii tiger sharks and those collected from 3 Indian Ocean sites (ESI, WSI, AS) were universally statistically significant, using both full and neutral SNP datasets and via both bootstrapping and permutation testing, with fixation indices estimated using robust sample sizes for SNP analyses (n = 29–32; Willing et al. 2012). The significant, pairwise values we obtained differentiating Hawaii sampled tiger sharks from the Indian Ocean subpopulations were very low (F_ST = 0.003–0.005) indicating weak differentiation between sites, and while subtle differentiation of Hawaiian tiger sharks was also found via DAPC (overall and neutral SNPs), the program Admixture did not reveal such differentiation. Notably, candidate outlier SNPs did not detect any evidence of adaptive differentiation between any of the Indo-Pacific tiger shark subpopulations (including Hawaii), indicating that genetic drift at neutral loci and/or a low mixing rate of sharks between CP and Indian Ocean waters may underly the shallow signals of differentiation seen here.

Within the Hawaiian archipelago, mature tiger shark females display parturition associated, long distance migrations from the Northwest Hawaiian Islands (NWHI) to the southern main Hawaiian Islands (distances of ~1200 km) (Lowe et al. 2006; Papastamatiou et al. 2013). Herein, our findings of 3 highly related tiger shark pairs (half-siblings or higher) from within the waters off Honolulu and Kaneohe Bay, are consistent with previous work indicating local reproduction within the Hawaiian archipelago (Whitney and Crow 2007; Papastamatiou et al. 2013; Meyer et al. 2018) and suggest some degree of site fidelity by these animals. These observations, along with the overall fairly high spatial residency displayed by these sharks around Hawaii (Meyer et al. 2018) suggest that tiger sharks inhabiting waters surrounding one of the earth’s most remote archipelagos have the potential for demographic independence, despite some individuals in the Pacific demonstrating clear capacity for long-distance, offshore excursions (Holmes et al. 2014; Werry et al. 2014; Ferreira et al. 2015; Meyer et al. 2018). We also note that even though several studies have investigated tiger shark movements in Hawaii (e.g., Holland et al. 1999; Papastamatiou et al. 2013; Meyer et al. 2009, 2010, 2014, 2018), to date, no tiger shark movements linking the Hawaiian archipelago to Australian or other coastal western Pacific waters have been observed.

Conservation Perspectives

The tiger shark is a phylogenetically and behaviorally distinctive species that likely plays a significant role in marine ecosystems as a globally distributed, large-bodied, apex predator with high versatility in habitat use and immense-distance migrations linking diverse ecosystems. Even though these particular natural history attributes would predict genetically homogeneous populations, this first genome-scale marker view of tiger shark population structure highlights the presence of at least 2 unambiguous, strongly differentiated populations (Atlantic and Indo-Pacific). With concerns about the health of marine ecosystems in an era of continued high exploitation of oceanic sharks in the midst of a rapidly changing environment, these findings make it imperative that tiger shark fishery management prioritize conservation of the evolutionary potential of both ocean basin genetic populations of this unique apex predator. Furthermore, the additional support provided by the genomics analysis for a putative, genetically differentiated population of tiger sharks in the waters surrounding Hawaii suggests that this population, which reproduces only triennially (Whitney and Crow 2007), should also receive targeted attention, at least from a precautionary biodiversity conservation perspective. A further genomic assessment that includes tiger sharks from its eastern and western Pacific Ocean range is warranted to clarify the population genetics of this species in this ocean basin.

Acknowledgments

Samples were generously provided by E. Brooks, D. Burkholder, D. Chapman, T. Daly-Engel, A. Danylchuk, P. Doukakis, K.M. Duncan, M. Grace, D. Grubbs, G. Harvey, J. Harvey, M. Heithaus, K. Holland, P. Mancini, C. Meyer, D. Pinhal, D. Reid, W. Robbins, B. Wetherbee, R. Wilborn, C. Wilson, S. Wintner, Panama City Observer Program (courtesy of S. Gulak and L. Hale), and the National Marine Fisheries Service. Samples obtained from ecological studies were conducted under protocols sanctioned by the Nova Southeastern University Institutional Animal Care & Use Committee (IACUC Authorization #DB1).

Funding

This work was supported by the Save Our Seas Foundation (SOSF157); Hai Stiftung-Shark Foundation (HF2019), and the Guy Harvey Ocean Foundation (GHOF19-1).

Data Availability

We have deposited the primary data underlying these analyses as follows in Dryad: doi:10.5061/dryad.xksn02vgn.