First circumglobal assessment of Southern Hemisphere humpback whale mitochondrial genetic variation and implications for management

The description of genetic population structure over a species’ geographic range can provide insights into its evolutionary history and also support effective management efforts. Assessments for globally distributed species are rare, however, requiring significant international coordination and collaboration. The global distribution of demographically discrete populations for the humpback whale Megaptera novaeangliae is not fully known, hampering the definition of appropriate management units. Here, we present the first circumglobal assessment of mito chondrial genetic population structure across the species’ range in the Southern Hemisphere and Arabian Sea. We combine new and existing data from the mitochondrial (mt)DNA control region that resulted in a 311 bp consensus sequence of the mtDNA control region for 3009 individuals sampled across 14 breeding stocks and subpopulations currently recognized by the International Whaling Commission. We assess genetic diversity and test for genetic differentiation and also estimate the magnitude and directionality of historic matrilineal gene flow between putative populations. Our results indicate that maternally directed site fidelity drives significant genetic population structure between breeding stocks within ocean basins. However, patterns of connectivity differ across the circumpolar range, possibly as a result of differences in the extent of longitudinal movements on feeding areas. The number of population comparisons observed to be significantly differentiated were found to diminish at the subpopulation scale when nucleotide differences were examined, indicating that more complex processes underlie genetic structure at this scale. It is crucial that these complexities and uncertainties are afforded greater consideration in management and regulatory efforts.


INTRODUCTION
Multiple biotic and abiotic factors operating at various spatial scales affect connectivity and rates of gene flow, thereby influencing the genetic population structure of a species across its geographic range (Anderson et al. 2010). Sufficient sampling at the scale of a species' range is therefore necessary to test for the existence of, and explain, geographic patterns. Non-representative sampling may result in a disproportionate level of importance being afforded to local processes in terms of their influence on the spatial structuring of populations, whereas processes operating over broader spatial scales remain undetected (Petit 2008). For globally distributed species, hemisphere-wide genetic studies provide an unparalleled perspective on the distribution and connectivity of population units (e.g. Peacock et al. 2015). Assessments of population structure at the hemisphere scale are rare, however, as they require international coordination and collaboration between multiple governmental, managerial, and scientific entities.
Humpback whales Megaptera novaeangliae are globally distributed baleen whales that undertake annual long-distance migrations between warm lowlatitude breeding areas and cold, highly productive feeding areas near the poles (Mackintosh 1942). The International Whaling Commission (IWC) defines demographically discrete populations of humpback whales for management purposes as management units (or breeding stocks) that 'are often identified by [significant] differences in frequencies of [mitochondrial DNA] mtDNA haplotypes or nuclear alleles, regardless of the underlying phylogeny, i.e., evidence of limited gene flow (Moritz 1994), although they can also be drawn along non genetic lines, e.g. through geography or demography' (Jackson & Pampoulie 2012, p. 2). In the Southern Hemisphere, the species is managed by the IWC as 7 breeding stocks, termed A (BSA) through G (BSG), that migrate in the summer to feeding areas in the Southern Ocean (Fig. 1). An eighth non-migratory Arabian Sea humpback whale (ASHW) stock occurs in the northern Indian Ocean and is considered part of our analysis due to its complete geographical, biological, and evolutionary isolation from all other Northern Hemisphere breeding stocks (Pomilla & Amaral et al. 2014). Evidence from a number of regional genetic studies suggests more complex population structure than accounted for in the current stock designations, results that echo similar findings for other highly migratory cetacean species in both hemispheres (e.g. Mendez et al. 2011, Costa-Urrutia et al. 2013. BSA, located off Brazil, shows relatively high diversity and no genetic substructure within the stock in a data set of 9 microsatellites (Cypriano-Souza et al. 2010), despite direct connectivity (of 1 individual) recorded between BSA and the distant BSC in the western Indian Ocean and song similarity between BSA and BSB in the eastern South Atlantic (Darling & Sousa-Lima 2005, Stevick et al. 2011. In contrast, BSB shows evidence of complex population 552 Fig. 1. Sampling locations. Sampling size for each location is in parentheses. Breeding stock and subpopulation abbreviations defined in Table 1. Breeding stocks are denoted as BS in the main text (e.g. BSB1) substructure. Two substocks, BSB1, which breeds off Gabon, and BSB2, which seasonally feeds off and migrates past western South Africa (and which shows lower than expected direct connectivity with, and some genetic differentiation from, BSB1 based on the mtDNA control region and 10 microsatellites), are recognized (Rosenbaum et al. 2009, Carvalho et al. 2014, Kershaw et al. 2017. Fine-scale temporal genetic substructure has been observed within BSB, consistent with the temporal segregation of migrating whales on the basis of age, sex, and reproductive status (Carvalho et al. 2014). Such variation may influence the interpretation of the number of demographically distinct populations due to temporal sampling effects (Carvalho et al. 2014). BSC is currently divided by the IWC into 4 substocks (BSC1−BSC4), although evidence of considerable migrant exchange between 3 of the substocks has emerged based on mtDNA control region and microsatellite data as well as satellite telemetry (Rosenbaum et al. 2009, Ersts et al. 2011, Fossette et al. 2014, Kershaw et al. 2017, suggesting relatively fluid exchange between the Mascarene region (BSC4), the Madagascar ridge (BSC3), and the Mozambique channel (BSC2). Gene flow between BSB and BSC is also apparent (Best et al. 1998, Pomilla & Rosenbaum 2005, Rosenbaum et al. 2009, Kershaw et al. 2017, and further sampling is needed to explore the possibility of a westward directional bias (i.e. from BSC to BSB).
Discovery tag recoveries imply that humpback whales from western Australia (BSD) are likely to mix with individuals from eastern Australia (BSE1) on Antarctic feeding grounds (Chittleborough 1965, Dawbin 1966. Levels of genetic connectivity are little known, however, and a better understanding of relationships between these breeding stocks is needed; nonetheless, some low but significant levels of genetic differentiation have been shown based on mtDNA control region sequences and 10 microsatellite loci (Schmitt et al. 2014).
Populations in the western Pacific Islands of Oceania, comprising substocks in New Caledonia (BSE2), Tonga (BSE3), the Cook Islands (BSF[CI]), and French Polynesia (BSF[FP]), show limited demographic exchange and a high degree of fidelity to breeding areas based on mtDNA control region se quences (Olavarría et al. 2007) and extensive photo-identification data (Garrigue et al. 2002, 2011, Franklin et al. 2014, supporting the designation of distinct population units in this region. However, fully resolving the boundaries of these stocks and subpopulations will require the integration of acous tic, genetic, photo-identification, and satellite tele metry data (Garland et al. 2015).
BSG, located off northwestern South America, is thought to be relatively isolated but represents some mitochondrial haplotypes that are characteristic of North Pacific humpback whale populations, suggesting historic or ongoing trans-equatorial gene flow (Baker & Medrano-Gonzalez 2002). Finally, the ASHW population is considered to be small and extremely isolated and shows the greatest genetic differentiation from all the other Southern Hemisphere breeding stocks based on mtDNA control region sequences and microsatellite data (Rosenbaum et al. 2009, Pomilla & Amaral et al. 2014. Despite the information contributed by regional and local analyses in specific ocean basins, a comprehensive understanding of humpback whale population structure can only be obtained by considering all breeding stocks concurrently in the same hemisphere-wide analysis. Here, we combine new data not previously subject to peer review with published regional data sets, comprising more than 3000 individuals sampled across 14 locations, and undertake the first assessment of humpback whale mtDNA population structure across the entire species' range in the Southern Hemisphere and Arabian Sea. We explore genetic structure at the population (breeding stock) and subpopulation (substock) scale, with and without the a priori designation of sampling location, and interpret the results in the context of existing information from studies conducted at the regional and local scales. These analyses provide new insights into potential circumpolar drivers of humpback whale population structure that serve to assist the interpretation of previous findings.
Prior to their protection by the IWC in 1963, some 215 000 humpback whales were hunted in the Southern Hemisphere as a result of open-boat (i.e. 18 th and 19 th centuries) and more recent 20 th century modern whaling, including heavy illegal Soviet whaling (Rocha et al. 2014

Sample collection and DNA sequencing
We assembled a database of previously published mtDNA control region sequences collected from biopsy samples or sloughed skin from 3009 individuals across 14 sampling sites representing all the known breeding stocks in the Southern Hemisphere and Arabian Sea ( Fig. 1; no genetic samples were available from BSC4 at the time for inclusion, and so this substock is not included in this analysis). All research undertaken followed local regulations and guidelines.
Sample collection, preservation methods, total genomic DNA extraction, amplification, and sequencing of the mtDNA control region are described elsewhere (Olavarría et al. 2007, Rosenbaum et al. 2009, Anderson 2013, and sample information is detailed in Table 1. Rigorous error-checking procedures were carried out by each respective research group, enabling the integration of these data sets for the subsequent analysis (Morin et al. 2010). All 3009 sequences were aligned in Geneious Pro 4.6.5 (Biomatters) using MUSCLE 3.6 (Edgar 2004) with a maximum of 10 iterations. The consensus sequences for the contributing data sets (Olavarría et al. 2007: 470 bp;Rosenbaum et al. 2009: 486

Genetic diversity and differentiation
We computed haplotype and nucleotide genetic diversity for each breeding stock (Nei 1987) by grouping samples based on their sampling location using DnaSP v.5 (Librado & Rosaz 2009). Genetic diversity indices (number of haplotypes, haplotype diversity, nucleotide diversity with Jukes-Cantor correction, and average number of average pairwise nucleotide differences among sequences) were calculated in DnaSP. To account for the effect of different sample sizes on the number of haplotypes identi-fied, rarefaction curves were constructed using the function rarefaction (Jacobs 2011), which makes use of the function rarefy in the vegan package ) in R v.3.2.2 (R Core Team 2016). Rarefy calculates the expected number of species in a subsample of the specified size, as defined in (Oksanen 2016).
The spatial hierarchical structure of the data among breeding stocks and subpopulations was evaluated through an analysis of molecular variance (AMOVA) (Excoffier et al. 1992) in Arlequin v.3.5.1.2 (Excoffier et al. 2005). Pairwise genetic differentiation between stocks and substocks was calculated using fixation indices based on haplotype frequencies (F ST ; Weir & Cockerham 1984) and genetic divergence (Φ ST ; Excoffier et al. 1992). Statistical significance was evaluated using the null distribution generated from 10 000 non-parametric random permutations of the data at the 0.05 significance level. No correction for simultaneous tests was applied.
We tested for genetic isolation by geographical distance (isolation by distance, IBD; Wright 1943) by exploring correlations between F ST (and the Slatkin linear F ST ) statistics and the shortest geographical distance between the approximate centroids of all breeding stocks (Y1 matrix). We used Mantel tests with 10 000 permutations to evaluate significance of such spatial−genetic distance correlations, as implemented in Arlequin (Excoffier et al. 2005).
Sequences were collapsed to haplotypes, and the alignment was converted to Roehl data format (.RDF) using DnaSP v.5 (Librado & Rosaz 2009). Medianjoining haplotype networks, both with and without maximum parsimony post-processing, were calculated using NETWORK v.4.6.0.0 (Bandelt et al. 1999, www. fluxus-engineering.com) with ε = 0 (i.e. zero expected homoplasy showing all possible connections) and all variable sites weighted equally. Median-joining networks have been recommended over maximum parsimony approaches in intraspecific studies, as they capture a greater degree of ambiguity, thus enabling more realistic interpretations (Cassens et al. 2005).
Genetic structuring was further investigated using both a principal component analysis (PCA) and a spatial PCA (sPCA)  h: haplotype diversity; π: nucleotide diversity (i.e. average number of nucleotide differences per site between randomly chosen sequences); θ: sample-specific theta (genetic diversity accounting for population size and gene flow) estimated by the program MIGRATE. Standard error for H r25 , standard deviation for h and π, and 95% confidence intervals for θ are shown in parentheses. Because the consensus sequence of the combined data set is of a shorter length (311 bp) than each of the contributing data sets (470−474 bp), genetic diversity measures were lower than those reported in Olavarría et al. (2007), Rosenbaum et al. (2009), and Anderson (2013).
ASHW: Arabian Sea humpback whale lation genetic model and therefore offer a complementary approach to fixation indices. The 2 methods were implemented in the ade4 v.1.7-3 (Dray & Dufour 2007) and the adegenet v.2.0.0 (Jombart 2008) packages in R (R Core Team 2016), respectively. PCA is a general method for representing highly dimensional data in a smaller number of dimensions, thereby unveiling the main factors explaining the structure of genetic variation in large samples. sPCA explicitly incorporates spatial information in the investigation of genetic variability . We carried out the sPCA by constructing a connection network using a matrix of the inverse Euclidian distance between substock sampling locations. This network was used for the calculation of Moran's I. sPCA optimizes the product of the variance of individual scores, based on the allelic frequencies of substocks, and of Moran's I to summarize genetic variability in a spatial context. Tests for global and local spatial structure were also implemented with 9999 permutations .

Migration rates between breeding stocks
To estimate the degree and directionality of matrilineal gene flow between contiguous breeding stocks, we used the maximum likelihood procedures implemented in MIGRATE software (Beerli & Felsenstein 2001). MIGRATE provides estimates of M (m/µ) and θ (2N e µ), where M is the mutation-scaled immigration rate, θ is the mutation-scaled population size, m is the immigration rate, µ is the mutation rate, and N e is the effective population size. The product θ M results in the number of im migrants per generation 2N e m. We adopted a maximum likelihood analysis strategy and migration matrix model allowing for asymmetric migration rates and variable θ between breeding stocks and variable sub population sizes (without subsampling). We se lected an inheritance scaler selected for θ of 1.00. Our data were first tested with default starting values for the population size and M parameters in 20 in dependent replicates of the Markov chain scheme: each replicate employed a static heating scheme and comprised 15 short chains (dememorization: 10 000 genealogies; recorded genealogies: 100; sampling increment: 100) and 3 long chains (dememorization: 10 000 genealogies; recorded genealogies: 1000; sampling increment: 100). We monitored acceptance rates of genealogies (and other parameter updates), effective sample sizes, and autocorrelations between parameters throughout the run. Using as initial parameters the consis-tent resulting values from these 20 initial runs, and employing the same Markov chain scheme, we launched 100 new runs to estimate the degree and directionality of gene flow between our stocks of interest at the breeding stock level. Convergence of Markov chain Monte Carlo algorithms was assessed using the program Tracer (Rambaut & Drummond 2007), which plots the log probability of data given parameters across the run.

Genetic diversity and population structure
A total of 184 haplotypes were identified among the 3009 samples from the 14 sampling locations based on 74 polymorphic sites. Because the consensus sequence of the combined data set was of a shorter length (311 bp) than each of the contributing data sets (470−474 bp), genetic diversity measures were lower than those reported for region-specific analyses in Olavarría et al. (2007), Rosenbaum et al. (2009), and Anderson (2013). We found genetic diversity to vary broadly across breeding stocks at the regional level ( At the breeding stock scale, AMOVAs were statistically significant (F ST = 0.029, Φ ST = 0.166, p < 0.001 for both indices), and most fixation indices from pairwise breeding stock comparisons were statistically significant (p < 0.05; Table 2). All haplotype-based indices (F ST = 0.005−0.191) and most nucleotidebased indices (Φ ST = 0.003−0.016) were found to be statistically significant (p < 0.05), except for the comparisons between stocks BSA-BSB-BSC-BSD and BSC-BSD that resulted in non-significance (Table 2). All subpopulation comparisons based on haplotype information alone were statistically significant (F ST = 0.001−0.191, p < 0.05), whereas significance dropped markedly for various pairwise comparisons when nucleotide information was incorporated (Φ ST = 0.003− 0.016), particularly from BSB2 and BSC (Table 3). No evidence of IBD between breeding stocks was found [F ST = p(rY 1 rand > rY 1 obs = 0.39, Slatkin F ST = p(rY 1 rand > rY 1 obs = 0.37), where rand = random value, and obs = observed value].
Median-joining networks showed comparable re sults irrespective of whether or not maximum parsimony (MP) post-processing was included. As expec ted, the median-joining network without MP post-processing captured a larger number of inferred nodes and reticulations. However, as the fundamental relationships be tween haplotypes were not affec ted, only the more parsimonious network with MP post-processing is shown (Fig. 2). No geographic clustering of haplotypes was evident across the network, with all breeding stocks and subpopulations broadly distributed across the 184 haplotypes, in cluding for the highly isolated ASHW population (Fig. 2).   The first 2 principal components of the PCA explained 73.39% of the variance in allele frequencies among putative populations (Fig. 3). While there are some unexpected results (the lack of relative separation of ASHW, for example; Fig. 3), the PCA resolves some broad patterns reflected by the fixation indices. The first principal component shows a clear separation of the BSE subpopulations (Fig. 3). In addition, BSB1 and BSC3 are highly distinct, possibly due to a higher level of variance for these sub-populations as a result of relatively larger sample sizes (Table 1). The second principal component further shows broad separation of subpopulations by latitude, with the samples from Oceania and off Colombia more closely situated with the subpopulations of the South Atlantic and western and northern Indian Ocean. The subpopulations off eastern Australia (BSE1), New Caledonia, and Tonga (BSE3) again show clear separation from the other subpopulations (Fig. 3).
Haplotype network of mitochondrial control region sequences (population size, N = 3009, number of haplotypes = 186) created using a median-joining algorithm with maximum parsimony post-processing with ε = 0, where ε is a weighting parameter for genetic distance, and variable sites weighted equally. Nodes are shaded according to sampling location. Size of node corresponds to the frequency of that haplotype among sampled individuals. Internal white nodes represent reconstructed median haplotypes. Gray notches represent nucleotide differences between haplotypes. Haplotype frequencies are available on the Dryad Data Repository (doi:10.5061/dryad.8cs4f). Breeding stock and subpopulation abbreviations defined in Table 1 The sPCA investigated the sequence variability using allelic frequencies in a spatial context, which further resolved some of the inconsistencies de scribed for the PCA (Fig. 3). The first positive axis had the highest eigenvalue and was therefore retained. Both global and local tests for structure were not found to be statistically significant (simulated p-value: global = 0.958, local = 0.217). The first global structure that was retained in the sPCA shows low variance between the South Atlantic and western and northern Indian Ocean subpopulations, as supported by the non-significance observed for Φ ST (Table 3). BSC3, however, shows relatively strong differentiation (Fig. 4); this pattern might also be expected for BSB1 based on the results of the PCA (Fig. 3), but this was not reflected to the same degree in the sPCA (Fig. 4).
The ASHW population appears to be grouped more closely with BSC to the west than BSD to the east, a pattern not detectable from the fixation indices alone (Table 2, Fig. 4). A second grouping occurs for BSD and the sub- The neighboring network was based on a matrix of the inverse Euclidean distance between sampling locations. Each stock and subpopulation is mapped by color coding its sPCA lagged score as intensity of the given color channel (first axis: red). Breeding stock and subpopulation abbreviations defined in Table 1 Fig. 3. Subpopulation principal components on the first 2 axes. Each subpopulation is color coded and corresponds to the haplotype network key (Fig. 2). PC1: principal component 1; PC2: principal component 2. Breeding stock and subpopulation abbreviations defined in Table 1 populations of BSE and Oceania; however, there appears to be some localized variance between the subpopulations of BSE and BSF, reflecting the relatively high observed F ST values (Table 3, Fig. 4). BSG is highly distinctive and shows clear separation from breeding stocks to both the west and the east (Fig. 4).

Migration rates between breeding stocks
Matrilineal gene flow between contiguous breeding stocks was estimated to be asymmetric in all pairwise comparisons. Asymmetries ranged from a 1.03fold difference in eastward compared to westward migration in BSF-BSG to a 62-fold difference in BSC-ASHW (Fig. 5). BSA showed notable asymmetry, with very little gene flow into BSB and much greater gene flow from BSB to BSA. The largest directional gene flow magnitudes were from BSC to ASHW, followed by BSC to BSD and BSB to BSA. The ASHW population is the only stock in this study where no historical emigration is estimated to any stock, and only immigration from other stocks is exhibited. BSD, BSE1, and BSF exhibited intermediate degrees of gene flow (Fig. 5). BSE showed migration both to and from BSE1 to the west and BSF to the east (Fig. 5). Gene flow estimates between BSA and BSG were very low, with a very slight westward bias (Fig. 5).

DISCUSSION
Our analysis of the genetic structure of the humpback whale across its circumpolar range in the Southern Hemisphere, based on a 311 bp fragment of the mitochondrial control region, revealed low but significant patterns of population genetic structure. Population and subpopulation comparisons confirm previous observations of complex patterns of genetic differentiation and gene flow, supporting assertions regarding the hierarchical nature of humpback whale population genetic structure (Rosenbaum et al. 2009, Baker et al. 2013, Kershaw et al. 2017). In addition, the global context of the study provided an opportunity to consider potential drivers of circumpolar genetic structure, thus providing the large-  Table 1 scale geographic context required for the interpretation of existing regional and local analyses. All pairwise comparisons between breeding stocks (i.e. BSA−BSG) based on haplotype frequencies (F ST ) and the majority of pairwise nucleotide differences (Φ ST ) were statistically significant (p < 0.05) ( Table 2; haplotype frequencies are available on the Dryad Data Repository, doi: 10. 5061/ dryad. 8cs4f), indicating that population structure at this scale is driven by maternally directed site fidelity to breeding areas (Mackintosh 1942, Baker & Medrano-Gonzalez 2002, Baker et al. 2013. The dominant influence of maternally directed site fidelity at the regional scale has also been confirmed by comparative assessments of bi parentally inherited nuclear markers (Ruegg et al. 2013, Schmitt et al. 2014, Kershaw et al. 2017. Although significant genetic differentiation was observed at the regional scale, patterns of estimated gene flow between breeding stocks were found to be relatively high and asymmetric in all cases (Fig. 5).
Moreover, no clear phylogeogra phic patterns were found among the haplotypes analyzed, with all 184 haplotypes broadly distributed across the network, even for those breeding stocks known to have experienced significant levels of isolation (i.e. ASHW and BSG; Fig. 2). At local spatial scales within breeding stocks, more complex patterns of population structure have been observed using nuclear markers (Kershaw et al. 2017). A future comparative study employing microsatellite markers would therefore be useful to further elucidate circumpolar-scale population structure. The relatively high levels of historic gene flow for some breeding stocks (Fig. 5) and the absence of phylogeographic signals (Fig. 2) support the notion of a shared evolutionary trajectory through out the Southern Hemisphere. These findings support previous suggestions that the Southern Hemisphere has an independent evolutionary trajectory from the North Atlantic and North Pacific (Jackson et al. 2014).
The asymmetries in regional fidelity and mtDNA differentiation observed in this study are consistent with findings from North Pacific humpback whale populations (Baker et al. 2013), indicating that this may represent a general pattern for humpback whale population structure globally. The number of statistically significant comparisons between breeding areas is also comparable, although the magnitude of fixation indices is higher in many instances in the North Pacific (Table 2; Baker et al. 2013, their  Table 4). The lack of clear evidence of fidelity to feeding areas in the Southern Hemisphere (Schmitt et al. 2014, Amaral & Loo et al. 2016) and lower mag-nitudes of genetic differentiation in some instances (Table 2; Baker et al. 2013, their Table 4) may indicate that interchange of individuals between breeding areas may be greater in some parts of the Southern Hemisphere relative to the North Pacific, possibly facilitated by access to the continuous or near-continuous circumpolar feeding habitat of the Southern Ocean.

Connectivity between the South Atlantic and southern Indian Ocean
While significant genetic differentiation was maintained at the subpopulation level for haplotype frequencies, significance dropped markedly for nucleotide-based indices, particularly between BSA, BSB1 and BSB2, substocks BSC1 to BSC3, and BSD (Table 3). Statistical differences in haplotype frequencies but no statistical support for molecular distances suggest genetic divergence at the level of a demographic aggregation of breeding individuals with small evolutionary divergence (i.e. resulting in haplotypes that differ by only 1 or 2 mutations). The relatively high rates of migrant exchange observed between breeding stocks in the South Atlantic and southern Indian Ocean (Fig. 5), and the relatively lower levels of variance between BSA-BSB-BSC-BSD identified by the sPCA (Fig. 4), support a conclusion of shallow levels of genetic divergence compared to other regions (Tables 2 & 3 Regional-scale contact between, and interchange of, individuals from different breeding stocks (e.g. between BSB and BSC) is assumed to take place primarily on shared feeding areas in the Southern Ocean. Densities of krill, the primary prey of Southern Hemisphere humpback whales, are relatively lower in the Atlantic (3.9 million tons km −2 ) and the southwestern Indian Ocean (2.3 million tons km −2 ) sectors of the Southern Ocean compared to other areas (e.g. the Antarctic Peninsula = 131.0 million tons km −2 ; Nicol et al. 2000). The lower maximum krill density in the Atlantic and southwestern Indian Ocean sectors may provide an explanation for some of the individual long-range movements and broadscale connectivity observed across BSA to BSD, as individuals could be forced to undertake longdistance longitudinal movements to maximize feeding opportunities. Such movement may increase the likelihood of mixing with other stocks and/or switching between breeding grounds.
At the local scale (i.e. within stocks), a more nuanced interplay of processes appears to be influencing genetic subpopulation structure. Genetic evidence supports the existence of 2 demographically discrete substocks (i.e. BSB1 and BSB2) off West Africa; however, an alternative hypothesis proposes that each substock represents 2 temporal ends of a single population widely distributed in space and time (Carvalho et al. 2014, Rosenbaum et al. 2014. Recent analyses of 9 nuclear microsatellite loci indicate that connectivity between BSB2 and BSC1 may also be greater than previously supposed (Kershaw et al. 2017). Animals from BSB2 have been observed to leave the feeding area of Bouvet Island, directly south of the African continent, and disperse widely eastward and westward in the Southern Ocean into feeding areas generally associated with BSC (0° to 70°E; Seakamela et al. 2015). Divergence between the substocks of BSC is also relatively shallow, calling into question their current delineation (Ersts et al. 2011). For instance, there is increasing evidence that BSC2 is likely to represent a mixed stock of animals from BSC1 and BSC3 (Best et al. 1998, Kershaw et al. 2017). In addition, long-distance movements between northeastern Madagascar (BSC3) and the coasts of Kenya and Somalia in northern BSC1 appear more frequent than supposed (Fossette et al. 2014, Cerchio et al. 2016).
Significant structure between western Australia (BSD) and eastern Australia (E1; F ST = 0.030, Φ ST = 0.017, p ≤ 0.05; Table 2, Figs. 3 & 4), at levels a magnitude greater than those observed between the substocks of BSB and BSC (Olavarría et al. 2007, Rosenbaum et al. 2009), confirms the results of previous genetic analyses (Anderson 2013, Schmitt et al. 2014) and resighting data (Dawbin 1966). Gene flow between BSD and BSE1 has been considered to be moderate and bidirectional at the regional level (Anderson 2013); however, at the circumpolar scale, we estimate female-mediated gene flow to be low relative to most other pairwise stock comparisons (Fig. 5). Sporadic bidirectional interchange has been observed between these 2 stocks (Chittleborough 1965, Noad et al. 2000, potentially due to annual changes in prey distribution that result in the mixing of stocks on feeding grounds (Anderson 2013). How-ever, our results indicate that the westward bias predominates, at least in terms of long-term maternally directed gene flow (Fig. 5).
Eastern Australia (BSE1) and the breeding substocks of Oceania (i.e. New Caledonia, BSE2; Tonga, BSE3; Cook Islands, BSF[CI]; and French Polynesia, BSF [FP]) are significantly differentiated from one another at a magnitude greater than the substocks of BSB and BSC for F ST in most comparisons, despite being more geographically proximal in some instances (Tables 2 & 3, Fig. 4). These results are consistent with the relatively high site fidelity to each island and significant population genetic structuring found in previous studies (Dawbin 1966, Garrigue et al. 2006, Olavarría et al. 2007). The degree of differentiation between BSE2 (New Caledonia) and BSE3 (Tonga; F ST = 0.009, Φ ST = 0.006, p ≤ 0.05) is lower than that between BSE1 (eastern Australia) and BSE2 (F ST = 0.023, Φ ST = 0.009, p ≤ 0.05) or BSE1 and BSE3 (F ST = 0.030, Φ ST = 0.004, p ≤ 0.05; Table 2). Genetic and acoustic breaks have previously been identified that separate BSE1 and BSE2 from the rest of Oceania (Anderson 2013, Garland et al. 2015); however, photo-identification matches between BSE1 and BSE2 have been reported (Garrigue et al. 2002), indicating some level of exchange, and significant longitudinal displacement of whales from American Samoa, eastern Australia, and Tonga has also been observed (Robbins et al. 2011, Franklin et al. 2014. Gene flow between New Caledonia and Tonga has been found to be relatively high (Anderson 2013), suggesting a more recent link between these 2 stocks that is supported by photo-identification (Garrigue et al. 2006, Franklin et al. 2014, genotypic matches (Anderson 2013), and acoustical studies of song sharing (Garland et al. 2014). New Caledonia remains clearly differentiated from eastern Australia and Tonga in both the PCA (Fig. 3) and the sPCA (Fig. 4), however. Population-specific analyses suggest that New Caledonia may be an area of transience, with whales briefly passing through the study area and subsequently moving in a wide range of directions (Garrigue et al. 2015). A transient area would result in higher genetic variance relative to BSE1 and BSE3, as was detected by the PCA and sPCA (Figs. 3 & 4).
The only instance of non-significant pairwise genetic differentiation across Oceania was observed between Tonga (BSE3) and the Cook Islands (BSF[CI]) at the nucleotide level (Φ ST = 0.005, p > 0.05; Table 3). Resights in the Cook Islands based on photo-identification and genotypic matching are relatively low (Constantine et al. 2012), but satellite tele metry data have shown 6 individuals moving directly towards Tonga upon leaving the Cook Islands (Anderson 2013). While we did observe significant genetic differentiation for French Polynesia and the Cook Islands (Table 3), previous genetic studies have suggested that these 2 substocks represent a single panmictic population (Anderson 2013). The Cook Islands may also be a location where whales from both BSE and BSF mix, for example as a migratory corridor for whales wintering in Tonga (Constantine et al. 2012), or, alternatively, they may use this location in alternate years (Anderson 2013). The relatively high levels of gene flow we estimated between BSE1-BSE2-BSE3 and BSF(CI) and BSF(FP) support these observations of mixing across Oceania (Fig. 5). Resolution of feeding areas will assist in the elucidation of some of the processes driving connectivity across this region.

Genetic isolation of BSG and ASHW and fidelity to feeding areas
Genetic structure and gene flow estimates indicate that BSG in the eastern South Pacific and ASHW in the Arabian Sea are the most isolated (Tables 2 & 3 ,Figs. 4 & 5) and least genetically diverse populations (Table 1). BSG winters primarily off Colombia and Ecuador, extending northward to the coasts of Panama and Costa Rica. Given the northward extension of the breeding area, and the presence of private haplotypes found within BSG that may be of northern origin, levels of isolation for BSG observed here might decrease if Northern Hemisphere populations were also considered in the analysis (Anderson 2013). In the austral summer and fall, whales from BSG migrate south to feeding areas in the Fueguian Archipelago, west of the Antarctic Peninsula (Zerbini et al. 2006, Cypriano-Souza et al. 2010, and recently also to the Corcovado Gulf, in the northern Chilean Patagonian channels (Acevedo et al. 2013). Photoidentification indicates philopatry and an absence of movements among the 3 feeding areas (Acevedo et al. 2013), possibly due to high and stable densities of krill in the vicinity of the Antarctic Peninsula (Nicol et al. 2000). If the majority of individual interchange between stocks and substocks does indeed occur on feeding areas, the fidelity of BSG feeding may be sufficient to drive its isolation from BSF to the west and BSA to the east. One individual has been observed to move from BSG to BSA, indicating that long-distance movements are possible (Stevick et al. 2013). However, despite their geographic proximity, there is no evidence of overlap between the feeding areas of BSG and those of BSA, located in the Scotia Sea.
Gene flow estimates between the ASHW population and its contiguous breeding stocks are singular. The lack of evidence of gene flow from the ASHW population to any of its contiguous stocks (BSC to the southwest and BSD to the southeast) corresponds to the Northern Hemisphere timing its breeding season, resulting in it being offset by 6 mo from all other Southern Hemisphere stocks. In addition, the upwelling system of the Arabian Sea provides yearround foraging habitat, mitigating the need for feeding-related migration (Pomilla & Amaral et al. 2014). The very high estimates of gene flow into the ASHW population from BSC (Fig. 5) present an apparent contradiction with the population structure parameters and can be explained, at least partially, by the time scale at which these estimations are optimized. The migration rate estimation conducted by MIGRATE explicitly accounts for the phylogenetic signal in the data, likely pointing to the historical connectivity between these currently isolated stocks. This historic connectivity is also suggested by the results of the sPCA, which indicate that ASHW is more closely grouped with the western Indian Ocean substocks compared to BSD off western Australia (Fig. 4). Further testing of hypotheses regarding the origin of ASHW using samples from BSC and BSD is required (e.g. using approximate Bayesian computation analysis; Fontaine et al. 2014). Highly asymmetric gene flow from populations in East Africa to Arabia has also been observed for humpback dolphins (Mendez et al. 2011) and has been attributed in part to the dominant oceanographic current systems of these coasts, which facilitate northward connectivity but not southward movement. While recognizing that drivers of connectivity for a nearshore small cetacean and an offshore baleen whale may not be directly comparable, it may be pertinent to explore whether this represents a general pattern for cetaceans in the Indian Ocean.

CONCLUSIONS AND MANAGEMENT RECOMMENDATIONS
Collectively, these results indicate that patterns of population genetic structure and connectivity for Southern Hemisphere humpback whales across their circumpolar range cannot be generalized across spatial scales. While maternally driven fidelity to breeding areas is the primary influence on population genetic structure at regional scales, relatively high levels of connectivity in some regions indicate that long-distance movements by females should be afforded more consideration as a mechanism facilitating gene flow between both adjacent and non-adjacent stocks and substocks. At the circumpolar scale, levels of estimated gene flow appear to correlate in some instances with krill density on feeding areas, suggesting that an interplay of site fidelity to stable feeding areas for some stocks (e.g. BSG) versus the need for individuals from other stocks to undertake extensive longitudinal movements to fulfill their energetic requirements (e.g. BSB and BSC) may be one of the factors influencing humpback population genetic structure at this scale. More studies into prey density and distribution on feeding areas at multiple geographic scales, and particularly in under studied areas such as those south of Oceania, are needed to further assess this relationship.
Our study demonstrates that, at least based on the statistically significant differences at the haplotype and nucleotide level of a 311 bp sequence of the mtDNA control region, the breeding stocks (BSA− BSG) recognized by the IWC represent genetically differentiated population units that should be considered as separate units in humpback whale management and regulatory assessments. Our results present a strong rationale for the management prioritization of BSG and ASHW in light of their genetic distinctiveness at the circumpolar context. In addition, the breeding stocks and subpopulations of Oceania represent the least abundant breeding population in the Southern Hemisphere, with little detectable recovery from whaling (Constantine et al. 2012), and are in significant need of protection.
We recommend a precautionary approach to the use of current subpopulation delineations as a basis for management, however. Some subpopulations may represent multiple demographic populations (e.g. BSB2, BSC2, and BSE2), in some cases, for example, as a mixed-stock migratory corridor, and there are continuing efforts to appropriately define the boundaries of these populations on both breeding and feeding areas (Rosenbaum et al. 2009, Carvalho et al. 2014, Garland et al. 2015. Failure to incorporate underlying genetic structure into population status assessments can result in significant biases through over-or under-representation of discrete population units in management goals. Moreover, levels of genetic connectivity between many stocks and subpopulations show complex and asymmetric patterns, requiring further research. Some regulatory processes, such as the recent review conducted by the United States that led to the species' down-and delisting under the US Endangered Species Act (ESA; Federal Register 2016), recognize fewer genetically distinct population segments (DPSs) than the IWC, as the former body operates under a more stringent set of criteria. Specifically, the 14 DPSs defined under the ESA do not capture breeding substock structure in the Southern Hemisphere, even though studies indicate that some substocks are as genetically differentiated from one another as separate DPSs delineated in the Northern Hemisphere (Bettridge et al. 2015). It is worth noting that all of the DPSs in the Southern Hemisphere are listed as Not at Risk under the US ESA, possibly due to a lack of resolution at the substock level (Federal Register 2016). Our results demonstrate that the IWC stocks and subpopulations should represent the minimum number of discrete populations taken into consideration until further research has been conducted. Moreover, uncertainties in stock and subpopulation boundaries, and how these uncertainties affect estimates of recovery, need to be carefully and transparently accounted for in any regulatory review of this species.