Thirteen dorsal-fin mounted satellite tags were deployed in continental shelf and slope waters of southern Australia between 2008 and 2013. These included position-only Sirtrack KiwiSat 202 and K2F161A tags, Wildlife Computers™ Smart Position or Temperature (SPOT) tags, and data-collecting Argos SPLASH and Mk10A tags. Sirtrack 202 and WC SPOT tags were programmed to transmit daily signals, whereas the SPLASH and Sirtrack K2F161A tags were duty-cycled to transmit at a 2-day frequency to optimize battery life. Capture and satellite tagging techniques are described in Rogers et al. (2015b). Shark total length, sex, and tag deployment locations were recorded (Table S1).

Satellite tags transmitted signals to the low polar orbiting environmental satellite network receiver stations, which were forwarded to Argos centers in France and the USA. Argos position estimates were accessed using Telnet and Tera Term Pro software, downloaded in seven location classes ranging from the highest to the lowest quality between 3, 2, 1, 0, A, B, and Z with predicted accuracies of 3 = < 250 m, 2 = 250–500 m, 1 = 500–1,500 m, classes 0 – B = >1,500 m, and Z = no position (http://www.argos-system.org). Raw data were mapped in MapInfo v. 11.5. to remove positions on land. Argos data were filtered by estimating locations using a Kalman filter under a continuous-time state-space framework using the (C)orrelated (RA)andom (W)alk (L)ibrary “CRAWL” package in R v. 2.15.2 (Johnson et al., 2008; R Core Team, 2013). Locations were interpolated along each filtered track to reduce sampling bias due to irregular transmission of Argos location data. To establish a set of spatial scale-based movement parameters, we estimated mean rate of movement (ROM) per day, minimum cumulative distance traveled and distal displacement distances (linear distance between tagging location and most distant location) for each individual (Table S1).

Tissue samples were obtained from 389 shortfin mako collected opportunistically from recreational and commercial fisheries catches or through collaboration with international research organizations. Samples were collected from six regions throughout the Southern Hemisphere (N = 275: Indo-Pacific, eastern Australia, southern Australia, Western Australia, New Zealand, and South Africa; Figure 1). Two regions from the Northern Hemisphere (N = 114: Northern Atlantic and northern Indian Oceans) were also sampled to assess trans-equatorial connectivity. Individuals were sampled at several locations within these broad regions to ensure that fine-scale geographic structure would be detected, if present (Figure 1). Tissue was preserved in either 95% ethanol or salt-saturated 20% dimethyl sulfoxide, and genomic DNA was extracted using a modified salting out protocol (Sunnucks and Hales, 1996).

The mitochondrial control region and portions of the flanking tRNAs (~1142 bp) were amplified by the Polymerase Chain Reaction (PCR) using primers Shark tPheR 5′-TYTCATCTTAGCATCTTCAGTGC-3′ and Shark tProF 5′-AGCCAAGATTCTGCCTAAACTG-3′. Reactions were conducted in 25 μL volumes comprised of 15–30 ng template DNA, 2.5 mM MgCl₂, 1 × MangoTaq™ reaction buffer (Bioline, Taunton USA), 0.25 mM dNTPs, 30 pmol forward and reverse primers and 1.25 U MangoTaq™ DNA polymerase (Bioline, Taunton USA). PCR cycling consisted of initial denaturation at 94°C, followed by “touchdown” cycling of 30 s denaturation at 94°C, 45 s annealing, and 1 min extension at 72°C. Annealing temperatures began at 59°C and decreased by 2°C at each touchdown, stabilizing at 51°C for 30 cycles. Amplified products were purified using ExoSAP-IT (Affymetrix USB® Products, Affymetrix, Inc., Cleveland USA), according to the manufacturer's protocol. Sanger sequencing was performed bi-directionally using internal primers mako405F 5′-GCCCGCTAGTTCCCTTTAATG-3′ and mako572R 5′-CCTTTCAGTTATGGTCAACTTGACAATC-3′, and BigDye® Terminator chemistry on an ABI 3730xl genetic analyzer (Applied Biosystems®, Life Technologies, Grand Island USA). DNA sequences were edited and aligned using Geneious® Pro v. 6.1.7 (Biomatters Ltd, Auckland New Zealand http://www.geneious.com). Sequences were cropped to 791bp for downstream analyses due to variability in sequence quality on either end.

Ten microsatellite loci were amplified using primers Iox-12 and Iox-30 (Schrey and Heist, 2002) Iox-M01, Iox-M110, Iox-M115, Iox-M192, Iox-M36, Iox-M59, Iox-B3, and Iox-D123 (via GenBank accession numbers KJ454433, KJ454434, KJ454435, KJ454436, KJ454437, KJ454438, KJ454439, KJ454440, respectively). Forward primers were tailed with a fluorescently labeled M13 tag (Schuelke, 2000). Reactions were conducted in 5 μL volumes comprising 15–30 ng template DNA, 3 mM MgCl₂, 1 × MangoTaq™ reaction buffer (Bioline, Taunton USA), 0.1 mM each dNTP, 0.1 pmol M13 tailed forward primer, 0.3 pmol reverse primer, 0.1 pmol fluorescently labeled M13 primer, 0.5 μg bovine serum albumin, and 0.25 U MangoTaq™ DNA polymerase (Bioline, Taunton USA). PCR cycling consisted of initial denaturation at 94°C, followed by “touchdown” cycling of 30 s denaturation at 94°C, 45 s annealing, and 1 min extension at 72°C. Annealing temperature began at 65°C and decreased by 2°C at each touchdown, stabilizing at 57°C for 30 cycles. Products were separated on an ABI 3730xl genetic analyzer (Applied Biosystems®, Life Technologies, Grand Island USA). Reference samples for each locus were included in all PCR programs and during capillary separation of fragments to ensure consistent genotype calling. Microsatellite alleles were visually inspected, binned, and sized according to the GeneScan™ 500 LIZ™ size standard (Applied Biosystems®, Life Technologies, Grand Island USA) using the Third Order Least Squares algorithm in the microsatellite plugin for Geneious® Pro v6.1.7 (Biomatters Ltd, Auckland New Zealand. http://www.geneious.com). Genotypes were checked for signatures of possible scoring errors due to null alleles, short allele dominance, and stutter peaks using Microchecker v. 2.2.3 (Van Oosterhout et al., 2004).

To avoid biases associated with limited sampling, samples from Western Australia and southern Australia, the Indo-Pacific and eastern Australia, were pooled for all frequency-based analyses of both mitochondrial and microsatellite data. Analysis of molecular variance (AMOVA) did not indicate any significant difference among these sampling locations, confirming the validity of this pooling scheme.

Mitochondrial DNA sequence variation and the extent of population differentiation were explored in Arlequin v. 3.5.1.2 (Excoffier and Lischer, 2010). Number of haplotypes, haplotypic, and nucleotide diversities were calculated assuming the Jukes and Cantor model of DNA substitution (Jukes and Cantor, 1969). An exact test of population differentiation was performed and population pairwise estimates of the parameters F_ST and Φ_ST were calculated. Significance was assessed via permutation (100,000 permutations) and interpreted following non-parametric Bonferroni correction (Rice, 1989). Hierarchical AMOVA was conducted partitioning total variance into within population, among population, and among regional covariance components (Cockerham, 1973) and testing for significance via permutation (10,100 permutations). A median-joining network (Bandelt et al., 1999) was constructed in Network v. 5.0 (Fluxus Technology Ltd) and illustrated in Network Publisher v. 2.0.0.1 (Fluxus Technology Ltd). Epsilon was set to 0 and hyper-variable sites were down weighted.

Microsatellite diversity was characterized using GenAlEx v. 6.5 (Peakall and Smouse, 2012) by calculating allele frequencies, number of alleles, effective number of alleles, and average observed, expected and unbiased expected heterozygosities per sampling location. Allelic richness was calculated in FSTAT v. 2.9.3.2 (Goudet, 2001). Genepop v. 4.2 (Raymond and Rousset, 1995) was used to assess whether the data conformed to expectations under Hardy-Weinberg and linkage equilibrium models. Bonferroni corrections for multiple comparisons were applied prior to interpretation.

Powsim v. 4.1 (Ryman and Palm, 2006) was used to determine the alpha error and statistical power with which significant genetic differentiation could be determined using our microsatellite dataset. Datasets were simulated with the same sample size, number of loci, and average allele frequencies as our observed and populations allowed to drift for a user-specified number of generations in order to attain a pre-defined level of differentiation (F_ST = 0.0005 to 0.05, 500 replicates per value). Statistical power was determined as the proportion of simulations for which Fisher's exact and Chi-square tests showed significant differentiation. Statistical α (type I) error was assessed by calculating the probability of rejecting H₀ when it is true for simulations omitting the drift step (i.e., F_ST = 0).

Population differentiation was investigated in GenAlEx v. 6.5. Pairwise fixation indices, Nei's G_ST, and Hedrick's G_ST^′′, were calculated following Meirmans and Hedrick (2011). Allelic differentiation, D_EST, was calculated following Jost (2008). Arlequin v. 3.5.1.2 was used to conduct an AMOVA of microsatellite data, partitioning total variance into within population, among population, and among regional covariance components. Significance was assessed via permutation (10,100 permutations).

Model-based clustering of genotypic data was conducted using Structure v. 2.3.4 (Pritchard et al., 2000). Since mobility is high in shortfin mako, allele frequencies were assumed to be similar across populations (Falush et al., 2003) and individuals were assigned using the admixture model of ancestry. Prior information regarding sampling location was allowed to inform ancestry in order to assist clustering (Hubisz et al., 2009). Inference was conducted over one million iterations (100,000 burn-in). Five independent runs were performed for each value of K, which varied from one to the number of sampled localities. Priors for the average and standard deviation of F (drift within populations) were set to 0.01 and 0.05, respectively, following Falush et al. (2003). A uniform prior (0, 10) on α (the parameter shaping the distribution of admixture proportion) was assumed. Following Evanno et al. (2005), ΔK (the second order rate of change of the log probability of the data given K (Ln P(X|K)) was calculated using Structure Harvester v. 0.6.93 (Earl and Vonholdt, 2012) and used to guide inference regarding the number of populations. Replicate clustering analyses were aligned using CLUMPP v. 1.1.2 (Jakobsson and Rosenberg, 2007) and visualized using distruct v. 1.1 (Rosenberg, 2004).

Analyses of sex-biased dispersal were conducted on a reduced dataset consisting only of individuals for which sex data were available (152 females (F) and 151 males (M) total; northern Indian 41 F: 40 M, South Africa 34 F: 57 M, eastern Australia 28 F: 20 M, southern Australia 21 F: 22 M, and New Zealand 28 F: 12 M).

The likelihood that an individual originates from its sampled location was calculated following Paetkau et al. (1995) using GeneClass2 v.2.0 (Piry et al., 2004). Log transformed likelihood values were corrected for population effects following Favre et al. (1997) resulting in corrected Assignment Indices (AI_c) that averaged zero per population and whereby negative values indicate lower than average probability of being born locally (migrants). The distributions of AI_c were calculated and compared for males and females, with the expectation that the more dispersive sex would show a more negative frequency distribution (Favre et al., 1997; Mossman and Waser, 1999). Following Goudet et al. (2002), the parameters F_ST, F_IS, and the mean (mAI_c) and variance (vAI_c) of AI_c were calculated and compared among sexes by taking the difference between the more dispersive and philopatric sex for F_IS (F_ISd – F_ISp), the difference between the more philopatric and dispersive sex for mAI_c and F_ST (mAI_cp – mAI_cd, F_STp – F_STd); or the ratio of the more dispersive to philopatric sex for vAI_c (vAI_cd / vAI_cp). Significant bias was detected using a randomization approach in FSTAT v. 2.9.3.2.

Following Banks and Peakall (2012), multivariate spatial autocorrelation analyses (Smouse and Peakall, 1999; Peakall et al., 2003) were compared across sexes to look for any sex-bias in fine-scale spatial patterns of genetic structure. Pairwise genetic distances were calculated following Peakall et al. (1995) and Smouse and Peakall (1999). Autocorrelation coefficients (Smouse and Peakall, 1999) were calculated across a range of distance classes that varied so as to incorporate comparisons within sampling localities, among adjacent localities, and more distant comparisons. Confidence intervals (95% CIs) about r were calculated by bootstrapping (Peakall et al., 2003) and the null hypothesis of no sex-bias was accepted if there was overlap in the CIs between sexes. The alternative hypothesis predicts that r values are significantly greater in the more philopatric sex. Heterogeneous autocorrelation across sexes was also assessed using single- (t²) and multi-distance (ω) class criteria as implemented in the non-parametric heterogeneity tests described by Smouse et al. (2008). These analyses were conducted in GenAlEx v. 6.5 and assessed for significance using 10,000 permutations and 10,000 bootstrap replicates.

Shortfin mako exhibited fidelity to the neritic waters of the Great Australian Bight, Bass Strait, southern Western Australia and the broad oceanic area to the west of Tasmania along the Sub-Tropical Front (Figure 2). Some sharks exhibited oceanic transit phases, leaving continental shelf waters to travel northward into the tropical waters of the northeastern Indian Ocean. During these migrations, three sharks traveled via the Perth and Carnarvon Canyons to the Bartlett and Karma Seamounts, located to the south of Indonesia. Other long-distance movements included four sharks that traveled southward to the Sub-Tropical Front. One shark traveled to the Coral Sea and another individual crossed the Tasman Sea to New Zealand shelf waters, followed by a northward migration to tropical waters near New Caledonia. A single individual moved as far west as the Crozet Plateau in the Indian Ocean.

Thirteen individuals were tracked for 249–672 days (mean 418 ± 37 per individual). Six individuals were tracked for more than one year. Figure 2 shows the spatial scale of movements by all individuals according to CRAWL model fits to the Argos data. Total cumulative distances traveled by shortfin mako ranged from 8,776 km in 262 days to 24,213 km in 551 days (Table S1). Distal displacement distances ranged from 1,500 to 7,520 km (mean 3356 ± 509 km per individual), with 69% of individuals exhibiting distal displacements greater than 2,000 km and 38% of individuals moving more than 4,000 km from their tagging locations. Many of these movements however, represented return migration events (Figure 2). There were no apparent sex-biases in scale of movement. The two longest (M5, M8) and shortest (M9, M12) movements were undertaken by both a male and a female. Males and females traveled an average of 40.8 and 37.5 km, respectively, per day.

Mitochondrial DNA data suggest matrilineal substructure across hemispheres, while nuclear DNA data indicate shortfin mako may constitute a globally panmictic population. There was generally high genetic connectivity within Australian waters.

The mitochondrial control region was sequenced for 365 shortfin mako resulting in 48 unique haplotypes, defined by 31 polymorphic sites, sampled across eight broad geographic regions (Table 1, Figure 1). Overall, haplotypic diversity was high (0.894 ± 0.013) while nucleotide diversity was very low (0.004 ± 0.003). Diversity metrics including sample size, number of haplotypes, haplotypic, and nucleotide diversity are shown in Table 1.

The haplotype network (Figure 3A) is characterized by a single, abundant haplotype that was sampled from all locations and in approximately 30% (108/365) of individuals. Other haplotypes are closely related, mostly separated by only a single substitution. Three, or fewer, substitutions were required to link any two haplotypes using parsimony. The network does not indicate any apparent geographic partitioning of haplotypes. Frequencies differ across sampling sites but most haplotypes are found at several, often geographically disparate, locations. One third of haplotypes (16/48) were unique to a single location, 13 of which were singletons.

Pairwise values of F_ST and Φ_ST based on mitochondrial data were low to moderate. There was significant differentiation among both Northern Hemisphere locations (Northern Atlantic and Northern Indian Ocean) and all Southern Hemisphere localities. There was significant differentiation between South Africa and all other locations based on exact tests of population differentiation. Weak but significant differentiation was also detected between South Africa and southern Australia based on F_ST, but this result was not corroborated by Φ_ST estimates, which showed no significant differentiation between South Africa and any of the Australasian locations (southern Australia, eastern Australia, or New Zealand). Within Australasia, significant differentiation was detected between southern Australia and New Zealand based on F_ST, Φ_ST and exact tests of population differentiation (Table 2A).

The results from AMOVA based on F_ST and Φ_ST were qualitatively similar. Interpretations presented herein are thus based on Φ_ST only. The global Φ_ST estimate was low, but significant (Φ_ST = 0.080; P = 0.000). Total variation in the dataset could be separated into five major regions: the northern Atlantic, the northern Indian, South Africa, western Australasia (western and southern Australia), and eastern Australasia (Indo-Pacific, eastern Australia, and New Zealand). While most of the total variation in the dataset was found within populations (91.5%, Φ_ST = 0.085, P = 0.000), among region variance accounted for a significant 8.2% (Φ_CT = 0.082, P = 0.009).

Ten microsatellite loci were genotyped for 355 shortfin mako sampled across the eight broad geographic locations (Figure 1, Table 1). There was no evidence of scoring errors, although Iox-12 and Iox-D123 showed evidence of null alleles in samples from a single location. The frequency of null alleles was low overall (< 10%) and all loci and populations conformed to Hardy-Weinberg expectations following Bonferroni correction. Linkage disequilibrium was detected between Iox-M110 and Iox-B3, Iox-12 and Iox-30, and Iox-M192, and Iox-D123, also in samples from a single location. All loci were therefore included in final analyses.

Genetic diversity at microsatellite loci was moderate to high. The number of alleles per locus ranged between 9 and 30, with means per population ranging from 7.1 to 15.8. The effective number of alleles per locus ranged between 5.3 and 9.4 across populations. Allelic richness was relatively consistent across populations, ranging between 6.7 and 7.2. Observed heterozygosity ranged from 0.75 to 0.86 (unbiased expected heterozygosity = 0.80–0.87; Table 1).

The microsatellite dataset has good statistical power with a 100% probability of detecting a true F_ST as low as 0.0025, and a high probability (65–70%) of detecting an F_ST as low as 0.001. The alpha error was ≤ 5%. The majority of microsatellite variation was within populations (99.8%) and the global multilocus F_ST estimate was very low (F_ST = 0.002), but significant (P = 0.020). This result was driven by significant F_ST values at just two of the 10 loci (Iox-M192, F_ST = 0.005, P = 0.004; Iox-M36 F_ST = 0.009, P = 0.001). Population pairwise estimates of G_ST, G_ST^′′, and D_EST were low and only a single pairwise comparison, South Africa vs. southern Australia, indicated significant differentiation (Table 2B, Table S2). The among-region variance component of AMOVA was not significant (F_CT = 0.004, P = 0.060), accounting for < 1% of total variation in the dataset.

High connectivity among sampling locations was also supported by model-based clustering analyses (Figure 3B). The mean estimated log probability of the data was highest for K = 1, while the modal value of the distribution of ΔK (Evanno et al., 2005) suggested that two clusters could be identified in the data. The ΔK metric cannot be estimated for K = 1 and so panmixia could not be assessed as a possible scenario using this approach. Further, this metric does not take into account the scale of ΔK. Observed values were two orders of magnitude smaller than is typical of cases of real structure and bar plots of the estimated cluster membership coefficients for each individual did not support K = 2 (Figure 3B). There was considerable variance in parameter estimates across runs for each individual K, suggesting non-convergence of the analysis despite running for a sufficient length of time. Together these observations are consistent with there being no signal of population structure in the data.

Patterns of molecular variation across sexes trended toward a signal of male-biased dispersal, however, this was not statistically supported.

Pairwise fixation indices (F_ST) based on microsatellite markers were low overall, but higher in females (F_ST = 0.003) than males (F_ST = 0.000). This difference bordered on significance (P = ~0.050), however the observed values of the test statistics for these parameters were within the range of the null distribution that dispersal is not biased by sex (Figure S1). F_IS was higher in males (F_IS = 0.009) than females (F_IS = 0.001), but this difference was not statistically significant (P = 0.203) and the observed value of the test statistic was also within the null distribution. Corrected assignment (AI_c) values ranged between −8.0 and 7.9 for males and −6.2 and 10.8 for females. The frequency distributions of AI_c values for males and females were largely overlapping and both sexes showed a similar proportion of negative values (54% for females and 52% for males). The mean and variance of AI_c were higher for females (mAI_c = 0.19, vAI_c = 11.87) than for males (mAI_c = −0.19, vAI_c = 9.11), however, these differences were not statistically significant (P = 0.158 and P = 0.838, respectively). The observed value of the test statistics for both mAI_c and vAI_c fell within the range of the null distribution representing the probability that dispersal is not biased by sex (Figure S1).

Spatial patterns of genetic structure were similar across sexes. The male and female 95% bootstrap confidence intervals about r overlapped in all distance classes (Figure 4). The single distance class t²-tests were all non-significant, as was the multi-class ω test of overall correlogram heterogeneity (ω = 6.2, P = 0.411). Low but significant positive autocorrelation among genotypes was detected for both males and females at small (100 km) distance classes (Table S3; Figure 4).

Satellite telemetry data indicated that shortfin mako in Australian waters exhibit multiple movement phases. Periods of residency in neritic habitats are probably indicative of time spent foraging. These contrasted with highly directional transitory movements, within neritic waters, and across vast oceanic expanses including among seamounts, ridges, and adjoining basins (Figure 2). For example, we hypothesize that features such as the eastward flowing South Indian Current, Sub-Tropical Front, and east–west running bathymetric features such as the Naturaliste Plateau, Diamantina Fracture Zone, and Broken Ridge may have facilitated the long-distance movement into the Indian Ocean that was recorded for one individual (Figure 2).

Overall, the geographical scale of the telemetry dataset spanned over 10,700 km from east to west. Observed movement patterns aligned closely with those documented in a previous study of juveniles (Rogers et al., 2015b) and both males and females exhibited similar scales of movement. While movement at these spatial scales indicate that shortfin mako in the Australian region are among the most HMS of pelagic sharks (Benavides et al., 2011; Rogers et al., 2013a,b; Holmes et al., 2014), it is important to also note that many of these movements were return events and some individuals exhibited fidelity to particular areas for extended periods.

Notably, no satellite tagged individuals traversed the lower tropical latitudes, nor passed through equatorial thermal frontal systems. Observed northernmost turning points of directional migrations aligned with surface water temperatures of 28–29°C. Southernmost latitudinal turning points coincided with the Southern subtropical frontal zone and were generally characterized by 9–11°C surface water temperatures. Although a single individual tagged with a standard tag as part of the New South Wales Department of Primary Industries Game Fish Tagging Program was recaptured after having apparently traversed the equator from the east coast of Australia to the Philippines (Rogers et al., 2015a), such transequatorial migration events appear to be uncommon in this species. Long-term telemetry studies of shortfin mako in the northeast Pacific Ocean also reported tropical thermal fronts aligned with turning points during similarly vast return migrations to shelf waters of the California Current ecosystem (Block et al., 2011). This apparent thermal preference has been documented by other studies (Holts and Bedford, 1993; Abascal et al., 2011; Musyl et al., 2011; Rogers et al., 2015b), providing further evidence that warm water may act as a potential barrier to dispersal among hemispheres.

We found considerable mitochondrial DNA diversity in shortfin mako. Haplotypic diversity was close to, or higher than, 0.9 at most sampling sites (Figure 3, Table 1), which is similar to that observed in previous studies of this species (Heist et al., 1996; Michaud et al., 2011; Taguchi et al., 2011, 2015) and toward the higher end of the range typically observed for elasmobranchs (Hoelzel et al., 2006; Dudgeon et al., 2008; Schultz et al., 2008; Chabot and Allen, 2009; Benavides et al., 2011; Blower et al., 2012; Corrigan et al., 2016). Also consistent with previous work, our mitochondrial DNA data showed evidence of trans-equatorial matrilineal substructure. Although haplotype sharing was observed among all locations, both Northern Hemisphere sampling locations (North Atlantic and northern Indian) showed significant differentiation from all other sampling sites (Table 2A). This indicates reduced matrilineal gene flow between hemispheres, consistent with the observation that trans-equatorial migration events appear to be infrequent according to tracking data.

Based on observed haplotype sharing between the Atlantic ocean basin and Australia/New Zealand, Michaud et al. (2011) hypothesized that gene flow between Pacific and Atlantic populations of shortfin mako occurs primarily via the Indian Ocean. Taguchi et al. (2011) were similarly unable to distinguish western Indian Ocean sampling sites from those in the eastern Indian or Pacific Oceans. Our analyses provide some support for this hypothesis. Pairwise Φ_ST estimates between South Africa and the Australasian sampling sites were low relative to comparisons between Northern and Southern Hemisphere sampling sites and not statistically significant, indicating gene flow across the Indian Ocean basin (Table 2A). A distal displacement distance of 7,520 km was recorded for one tracked individual tagged off southern Australia, representing a return westward movement to a location about 2,000 km east of South Africa. This suggests that suitable oceanic migratory pathways exist that could facilitate trans-Indian Ocean linkages between Australasian and South African populations (Figure 2, Table S1). However, exact tests of population differentiation indicated significant differentiation between Australasian and South African sampling sites, and a single pairwise comparison between South Africa and southern Australia was also statistically significant based on F_ST (Table 2A).

Taguchi et al. (2011) reported that the eastern Indian Ocean sample was highly differentiated from most other sampling sites, although this was based on limited sampling from the region. We are unable to comment on the validity of their finding because, despite extensive efforts, we too obtained only few samples from the eastern Indian Ocean off Western Australia. Taguchi et al. (2011) also indicated possible population structure between the eastern and western coasts of Australia. The Bassian Isthmus in southern Australia is a well-characterized biogeographic barrier that is thought to have promoted bicoastal population subdivision in several marine taxa (Teske et al., 2017). For example, Blower et al. (2012) reported matrilineal subdivision between eastern and southwestern coastal regions of Australia in the white shark, Carcharodon carcharias, a close relative of shortfin mako (Naylor et al., 2012). In contrast, we did not find any evidence of matrilineal population structure in shortfin mako sampled from around the Australian continent (Table 2A). Interestingly, however, the comparison between southern Australia and New Zealand indicated significant differentiation. The Φ_ST estimate between southern Australia and New Zealand is lower than those observed between Northern and Southern Hemisphere sampling sites, indicating that gene flow between these locations is less constrained than across the equator, but nevertheless restricted enough to represent significant differentiation (Table 2A). It is possible that matrilineal gene flow occurs in a “stepping stone” fashion throughout the region whereby southern Australia and New Zealand are connected via the east coast of Australia. Alternatively, the statistical significance of this single comparison may be artefactual (discussed below).

Further investigations into the extent of connectivity between Australian and South African waters would benefit from tracking information from additional adult individuals, as this may uncover links between neritic habitats on either coastline and reveal how these animals may use bathymetric features during transoceanic movements. Genetic data from an improved sampling of individuals from the southeastern Indian Ocean region would also help clarify the extent of gene flow across the Indian Ocean Basin and between the east and west coasts of Australia.

Schrey and Heist (2003) report very weak evidence of population structure between the North Atlantic and North Pacific Oceans according to their analysis of microsatellite DNA. Based on an analysis of a larger number of microsatellite markers, Taguchi et al. (2015) report that shortfin mako lack differentiation across their Pacific Ocean range. Sampling from Australasia and the Indian Ocean were limited in both of their studies, allowing little prior inference regarding nuclear genetic structure across the region.

Similar to Schrey and Heist (2003) and Taguchi et al. (2015), we inferred little evidence of population structure based on our microsatellite data. Only a single pairwise comparison, South Africa vs. southern Australia, indicated significant differentiation (Table 2B). The model-based clustering analysis suggested only subtle differences in allele frequencies across regions (Figure 3) and no apparent population structure.

The biological relevance of significant pairwise comparisons of fixation indices should generally be interpreted with caution given their observed small magnitude. F_ST and analogs measure the effects of gene flow on population differentiation and are thus influenced by both population size and migration rate. The magnitude of F_ST and analogs decreases non-linearly with increasing migration rate, such that a similar signal of weak to no genetic differentiation can be produced under a scenario of panmixia as well as when populations are large but sufficiently independent to warrant separate management (Waples and Gaggiotti, 2006; Waples et al., 2008; Gagnaire et al., 2015). This makes it difficult to precisely estimate these parameters, and interpret the significance of their magnitude, when population sizes are large and dispersal potential is high. This is likely the case of most marine species, including shortfin mako. Moreover, restricted sampling from a highly diverse set of genotypes can mean that some low estimates of pairwise differentiation are spuriously rendered statistically significant due to minor allele frequency differences (Waples, 1998; Waples et al., 2008).

Given that a high percentage of our tracked individuals showed long distance movements and that there was no evidence of genetic structure based on clustering analyses, it seems plausible that the statistically significant differentiation we detected based on fixation indices between Southern Australia and New Zealand in the mitochondrial data, and South Africa and Southern Australia in the microsatellite data, are artefactual (Waples, 1998; Waples et al., 2008; Gagnaire et al., 2015). However, it is worth noting that many long-distance movements by shortfin mako are return events and this could potentially promote genetic differentiation at smaller geographic scales than their mobility predicts. Additionally, the South Africa vs. Southern Australia comparison for the microsatellite data was also statistically significant based on Jost's D_EST, a complementary measure of population structure that quantifies allelic differentiation rather than nearness to population fixation and is less influenced by demographic variables (Jost et al., 2017). This may indicate some potentially biologically relevant partitioning of allelic diversity between these areas, although the magnitude of this parameter was still low.

Methods for collecting genomic-scale data are becoming readily available and the possibility to analyze data from thousands of markers may allow better estimates of low values of genetic differentiation in the future. However, this will likely remain challenging for many marine species with large population sizes (Waples et al., 2008; Gagnaire et al., 2015; Waples, 2015). Quantifying adaptive divergence between locations offers a solution for assessing differentiation even when population sizes are large and gene flow is high. Unlike genetic drift, selection counteracts the homogenizing effect of migration more efficiently in large populations. Genome scans of large marker datasets make it easier than ever to identify loci under selection for this purpose. Applying such methods to studying shortfin mako in the future may offer the possibility to delineate locally adapted units that require independent management even though they may be highly connected (Allendorf et al., 2015; Gagnaire et al., 2015; Jost et al., 2017). Contrasting patterns of neutral vs. adaptive variation may also be particularly relevant to make future predictions regarding how populations may respond to changing environmental conditions or fishing pressure.

Understanding biases in dispersal patterns can reveal ecologically important areas such as feeding or breeding grounds. This information can also guide fisheries management to avoid selective overharvest of a more philopatric sex (Hueter et al., 2005; Chapman et al., 2015). Male-biased dispersal has been demonstrated in a number of elasmobranch species (Schultz et al., 2008; Daly-Engel et al., 2012; Portnoy et al., 2015). This includes other Lamniformes such as white sharks (Pardini et al., 2001; Blower et al., 2012) in which both sexes are known to undertake oceanic scale movements (Bonfil et al., 2005; Bruce et al., 2006). Skewed sex ratios in shortfin mako catches indicate regional and seasonal sexual segregation (Mucientes et al., 2009; Francis, 2013) and Schrey and Heist (2003) propose male-biased dispersal as a possible mechanism to explain differing patterns of matrilineal vs. nuclear genetic structure in this species.

Under a scenario of sex-biased dispersal, allele frequencies should be more similar across sampling sites among individuals of the dispersing sex than those of the more philopatric sex. Because of this, the dispersing sex will show greater variance in assignment index and lower probability of local assignment, while F_ST will be higher among the more philopatric sex (Goudet et al., 2002). Observed values for these parameters were consistent with these expectations, suggesting that dispersal may be male-biased in shortfin mako. However, the difference in these parameters between sexes was not statistically significant (Figure S1).

There are several caveats to the interpretation of these results. These tests lack power when dispersal occurs at intermediate rates and sex-bias is subtle (< 80:20; Goudet et al., 2002). Tracking data, low pairwise fixation indices and the clustering analysis based on genetic data suggest that both male and female shortfin mako are highly mobile. Given their mobility and pelagic habit, it seems more likely that any female philopatry will be weak, perhaps at the scale of hemispheres given that we detected a signature of trans-equatorial matrilineal substructure and warm water at the equator appears to represent a physical barrier to dispersal. These tests are also only applicable if sampling occurs during the dispersed phase (Goudet et al., 2002). This assumption is likely violated here given that we sampled multiple cohorts of mostly juveniles and sub-adults. It is possible that violation of these assumptions is masking any signal of sex-bias in these particular analyses, although the trend indicates a male bias.

We also did not detect any statistically supported differences in spatial genetic structure across sexes based on spatial autocorrelation analysis. Detecting a sex bias using this method similarly requires large sample sizes and the development of strong spatial genetic structure in the more philopatric sex (Banks and Peakall, 2012). Banks and Peakall (2012) stress the importance of sampling and concentrating pairwise data points at the scale at which dispersal is restricted in the more philopatric sex. This analysis and our inferences regarding sex-biased dispersal in general, would thus benefit greatly from more information regarding the movements and mating behavior of adult individuals of both sexes, particularly identifying regions that are used for mating and parturition. Satellite tracking together with genetic analysis of a large sample of mature sharks collected during the breeding season from both hemispheres will be required.

Inferences of high connectivity based on analyses of our long-term telemetry and molecular datasets spanning six key regions support defining shortfin mako as a pelagic HMS in Australia and neighboring regions of the Southern Hemisphere. Although highly mobile, molecular data presented herein and elsewhere (Michaud et al., 2011; Taguchi et al., 2011, 2015) indicate separation between the Northern and Southern Hemispheres and weaker evidence of separation within the Southern Hemisphere. There appears to be distinct populations in the southeastern and southwestern Pacific (Michaud et al., 2011; Taguchi et al., 2015), and potentially southern Australia and the western Indian Ocean, though connectivity across the Indian Ocean is somewhat complicated to interpret. From a management perspective, it is most important to determine whether inferred differentiation is biologically meaningful such that units warrant management as independent stocks. Only a small number of migrants are required to homogenize allele frequencies across regions (Spieth, 1974; Mills and Allendorf, 1996). Significant spatial partitioning may occur despite high genetic connectivity and the number of migrants per generation required to allow stock rebuilding may be much higher than is required to produce genetic homogeneity (Waples, 1998; Waples et al., 2008).

Tagging data support the separation of Northern and Southern Hemisphere populations of shortfin mako, with only one tagged shark known to have crossed the Equator (Sippel et al., 2011; Rogers et al., 2015a,b; Holdsworth and Saul, 2017). Most movements of tagged mako have occurred within half-ocean basins (e.g., the southwest Pacific). Although Australasian mako frequently make long-distance movements, they also often return to near their tagging location, and importantly, show fidelity to specific areas of continental shelf and slope over several to many months (Rogers et al., 2015a,b; M. Francis unpublished data; data presented herein). Hundreds of mako tagged by gamefishers with standard tags in Australian and New Zealand waters all remained within the southwest Pacific Ocean (Sippel et al., 2011; Rogers et al., 2015a; Holdsworth and Saul, 2017). Fourteen mako tagged with electronic tags in New Zealand and tracked for 34–588 days (mean 251 days, all but one longer than 120 days) spent 42–100% (median 77%) of their time in the New Zealand Exclusive Economic Zone (EEZ; M. Francis, unpublished data). Those results, in combination with our own telemetry data from southern Australia that also showed fidelity to Australian waters, indicate that mako do not wander randomly across the globe, but instead may be resident in comparatively small areas for extended periods. Mako do cross international boundaries and the high seas, however, such that management at the scale of Regional Fisheries Management Organizations is important. But the propensity for mako to spend extended periods within national EEZs means that the homogenizing effect of large-scale movements likely occurs at a rate that is too slow to combat differing levels of fishing mortality across the entire genetic stock. This means that effective fisheries management of shortfin mako must occur at national as well as international levels.

JGP was employed by Pepperell Research and Consulting Pty Ltd. LM was employed by Stick Figure Fish Illustration. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.