Biogeography in the deep: Hierarchical population genomic structure of two beaked whale species

display significant genetic structure between the Atlantic, Indo-Pacific and in Cuvier ’ s, the Mediterranean Sea. Within major ocean basins, clear differentiation is found between genetic clusters on the east and west sides of the North Atlantic, and some distinct patterns of structure in the Indo-Pacific and Southern Hemisphere. We infer that macroevolutionary processes shaping patterns of genetic diversity include biogeographical barriers, highlighting the importance of such barriers even to highly mobile, deep-diving taxa. The barriers likely differ between the species due to their thermal tolerances and evolutionary histories. On a microevolutionary scale, it seems likely that the balance between resident populations displaying site fidelity, and transient individuals facilitating gene flow, shapes patterns of connectivity and genetic drift in beaked whales. Based on these results, we propose management units to facilitate improved conservation measures for these elusive species.


Introduction
The deep sea is the largest habitat on Earth (Sutton and Milligan, 2019), yet less is understood about the species that live there, their distributions, population structure and connectivity than any other ecosystem.In the deep sea, underwater features, such as steep continental slopes, canyons, and oceanic islands, can create areas of upwelling, aggregating prey species on which predators forage (Worm et al., 2003).Biodiversity hotspots are often contained within well-defined biogeographical provinces (Costello et al., 2017), where barriers may be fully impassable (e.g., Isthmus of Panama), or permeable (e.g., shallow rises; Bianchi and Morri, 2000), providing the basis for hierarchical levels of population structure in the deep sea.
Beaked whales (Ziphiidae) comprise the most speciose group of deep sea marine mammals, yet their often elusive, offshore and deep-diving lifestyle means that most information on these taxa comes from occasional finds of dead, beach-stranded specimens and several species are listed as data deficient on the IUCN Red List.Blainville's and Cuvier's beaked whales (Mesoplodon densirostris (Blainville, 1817), and Ziphius cavirostris (G.Cuvier, 1823), respectively; 'Blainville's' and 'Cuvier's' hereafter) are among the most well-studied beaked whale species, based on the number of publications.Both have cosmopolitan distributions and occupy the deep seas of all major ocean basins, although only Cuvier's is regularly found in the Mediterranean and at higher latitudes (MacLeod et al., 2006).Both species show site fidelity to some nearshore areas with steep bathymetry (Baird, 2019), potentially leading to population genetic structuring, especially in the presence of biogeographical barriers.
While many natural populations are at risk of genetic erosion due to declining abundance from human impacts (Leroy et al., 2017), some marine mammals (and particularly cetaceans) may be amongst the most susceptible taxa (Cammen et al., 2016;Schipper et al., 2008).Beaked whales are known to strand in unusual mortality events (UMEs) driven by underwater noise from human activities (Bernaldo de Quirós et al., 2019;Cox et al., 2006;Simonis et al., 2020).The majority of UMEs involving beaked whales have been reported for Cuvier's and/or members of the Mesoplodon genus (D'Amico et al., 2009).However, the extent to which these mortalities might cause long-term, population-level impacts on these poorly understood species is currently unknown.Identifying the baseline patterns of population genetic diversity and structure of beaked whales, and understanding the underlying environmental drivers, would enable mortality to be linked with population units, thus facilitating future management efforts in the face of anthropogenic impacts.
To date, only mitochondrial DNA (mtDNA) has been analysed to investigate population structure and genetic diversity in Cuvier's and Blainville's (Dalebout et al., 2005;Morin et al., 2012).Cuvier's mtDNA control region and whole mitochondrial genome (mitogenome) data revealed significant differences in haplotype frequencies between ocean basins, but not reciprocal monophyly, suggesting limited gene flow or incomplete lineage sorting (Dalebout et al., 2005;Morin et al., 2012).In contrast, Blainville's mitogenomes formed two reciprocally monophyletic clades according to Atlantic and Pacific ocean basins, indicative of no recent (maternal) gene flow (Morin et al., 2012).Although these studies suggest ocean-basin-level genetic structure, the paucity of beaked whale data makes it difficult to confirm these hypotheses and extend analyses to a finer geographical scale.Indeed, other deep-sea cetaceans can exist in populations that are socially-driven (e.g., sperm whales, Physeter macrocephalus; Alexander et al., 2016), island-associated (e.g., rough-toothed dolphins, Steno bredanensis;Jefferson, 2018) or panmictic (e.g., Gray's beaked whale, Mesoplodon grayi; Westbury et al., 2021).
Here, we aim to shed light on the population genetic structure, phylogeography and diversity of Cuvier´s and Blainville´s beaked whale species across their cosmopolitan ranges to inform conservation efforts and understand the micro-and macroevolutionary processes driving genetic structure in deep-sea species.We improve the spatial resolution of previous genetic studies by Dalebout et al. (2005) and Morin et al. (2012) by identifying finer-scale structure using larger sample sizes, and analysing whole mitogenomes and reduced-representation nuclear genomic data (ddRAD; double digest restriction-site associated DNA).We propose macro-and microevolutionary drivers for the substantial hierarchical genetic structure identified, which also allows us to recommend evolutionarily significant units (ESUs; Waples, 1995) and demographically independent populations (DIPs; Martien et al., 2019;Palsbøll, Bérubé, and Allendorf, 2007) for conservation and management.

Sample collection and lab protocols
Tissue samples were selected from the newly established International Tissue Archive of Beaked Whales (ITABW) to cover the full geographical ranges of Cuvier's and Blainville's (Electronic Supplementary Material (ESM) 1, Supplementary Spreadsheet Table (SST) 1).All sampling was done under appropriate local permits and, if applicable, local animal ethics approvals.DNA was extracted from tissue samples (ESM1) and sequencing libraries were prepared following the method described by Peterson et al. (2012) and optimised for beaked whales (Carroll et al., 2016(Carroll et al., , 2021) (ESM2).To explore major phylogeographical patterns, and supplement the nuclear analyses, we used "shotgun" sequencing (Carøe et al., 2018); ESM2 for protocol) to generate and analyse complete mitogenomes from a subset of the samples used for ddRAD, as well as additional low-quality museum samples, and augmented with publicly available sequences (Morin et al., 2012;SST1).

Phylogeography
Global structure based on ddRAD SNPs was visualized by generating rooted neighbour-joining trees (BIONJ, Gascuel, 1997) for each species, using sequences from southern right whales (SRW, Eubalaena australis) as the outgroup (ESM5).Full mitogenomic sequences were analysed in a Bayesian phylogenetic framework implemented in BEAST v.2.5.2 (Bouckaert et al., 2019).For the mitogenome data, the 12srRNA, 16srRNA and 13 protein-coding genes were individually aligned, and the stop codons for the coding genes were removed in GENEIOUS PRIME v2020.1.2(https://www.geneious.com/).The final alignments from the two rRNA and 13 protein-coding genes were subsequently concatenated into a single alignment, and substitution models were inferred for each codon position for each gene with PARTITIONFINDER v.2.1.1 (Lanfear et al., 2017).Control region (CR) data were analysed separately (see below).The resultant partitions were incorporated into BEAUTI, the graphical user interface for generating BEAST XML files (SST2).
For the Cuvier's mitogenome phylogeny, a strict clock with a Yule tree model was applied, as closely related species are not expected to have deviating clocks.Calibration dates followed most recent common ancestor (MRCA) estimates inferred from McGowen et al. (2019).Baird's beaked whale (Berardius bairdii; GenBank accession: NC_005274.1)was used as an outgroup and the MRCA between all Cuvier's and the Baird's beaked whale outgroup was set to a mean age of 15 million years ago (mya; 95% probability range 12.3-17.6).The model was set to a 15,000,000 chain length, sampled every 1500 steps and replicated three times to assess model A.B. Onoufriou et al. convergence.Log and tree files were combined with LOGCOMBINER v.2.5.2 (Drummond and Rambaut, 2007) after a 10% burn-in for each run.All parameters in the combined log file had ESS values greater than 200.The Maximum Clade Credibility (MCC) tree was prepared in TREEANNOTATOR v.2.5.2 (Bouckaert et al., 2019) and visualised in FIGTREE v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).
For the Blainville's mitogenome phylogeny, a strict clock with a Yule tree model was also applied.Two outgroups were used as they are both closely related to Blainville's (Gray's beaked whale, GenBank accession: NC_023830.1;Stejneger's beaked whale, Mesoplodon stejnegeri, GenBank accession: NC_036997.1)and their nodes for the MRCA were also calibrated using estimates from McGowen et al. (2019).The calibration node for the Blainville's alignments and the two outgroups was set to the mean age of 5.19 mya (95% probability range 4.12-6.26mya).The calibration node for the Blainville's alignments with the Stejneger's beaked whale outgroup was set to the mean age of 4.59 mya (95% probability range 3.6-5.58mya).The Blainville's model needed fewer iterations to converge than the Cuvier's model and hence was set to a 10,000,000 chain length, sampled every 1000 steps and replicated three times to check  for model convergence.The tree and log files were combined, and the output visualized as described for the Cuvier's data.

Genetic population structure
Isolation by distance (IBD) for the ddRAD SNPs was calculated per species using a Mantel test in 'ade4' v1.7-16 in R (Dray and Dufour, 2007), to compare Euclidean genetic distance and geographical distances calculated via the least-cost (LC) path distance over seawater in 'marmap' v1.0.5 (Pante and Simon-Bouhet, 2013)

(ESM6).
We conducted a hierarchical clustering analysis incorporating ddRAD SNPs and geographical information in a spatially explicit, least-squares optimization approach to estimate ancestry in the R package 'tess3r ' (v1.1.0, Caye, Jay, Michel, and Francois, 2018).The most likely number of clusters was estimated from a range (K=2-10) with n = 20 repetitions and n = 200 maximum iterations of the optimisation algorithm.Cross-entropy scores were plotted against K values to infer the most likely number (ESM7).When hierarchical spatial structure is likely, genetic clustering algorithms tend to find structure at the highest level of divergence (Evanno et al., 2005).Previous genetic work on the species using mtDNA data (Dalebout et al., 2005;Morin et al., 2012) suggests ocean-basin-level   differentiation, so we carried out analyses both between and within ocean basins to allow for the detection of finer-scale structure.The genetic clusters determined using 'tess3r' were further explored using Discriminant Analysis of Principal Components (DAPC; ESM8; Jombart, Devillard, and Balloux, 2010).A simulation study suggests that an a priori approach to DAPC successfully resolves the composition of the genetic clusters and the genetic distance between clusters correlates well with the value of F ST used to generate the simulated data (Miller et al., 2020).On the other hand, de novo DAPC analyses perform poorly in situations that are common to marine species; low genetic differentiation and IBD (Jombart et al., 2010;Miller et al., 2020).Therefore, while we explored both de novo and a priori DAPC analyses, we only report on the a priori results.
In addition, for future comparisons with other studies and criteria recently suggested for species, subspecies and population delineations in cetaceans (Taylor et al., 2017), we estimated F ST (and associated p-value) and d A (net nucleotide divergence, Nei, 1987, equation 10.21) between ocean basins for full mitogenomes and just the CR with the software package DNAsp v6.12.03 (Trozas et al., 2017).

Phylogeography and genetic population structure
The Cuvier´s exhibited substantial phylogeographic and population genetic structure in both the nuclear ddRAD SNPs and the mitogenomes.The mitogenome phylogeny of 35 Cuvier's indicated a separation into three major clades with high levels of support (Fig. 3a); Clade 1 contained lineages sampled exclusively in the North Atlantic (including the Mediterranean), whereas Clades 2 and 3 contained individuals from the western North Atlantic, and divided Indo-Pacific lineages (and one from Argentina) broadly along the equator (Fig. 3a, see SST1 for individual metadata and clade assignments).Individuals from the North Atlantic (Clade 1) appear to have diverged from the mostly Indo-Pacific Clade 2 approximately 1.5 mya, while the individuals in the Mediterranean diverged from the rest of the North Atlantic approximately 0.5 mya (Fig. 3a).Unlike the paraphyletic mitogenome phylogeny, the BIONJ phylogeny based on the final ddRAD SNP data consisting of 33137 SNPs genotyped in 118 Cuvier's (ESM5) and 6 SRW outgroup showed clearer monophyly of Atlantic and Indo-Pacific Ocean clades (Fig. 3b).Like the mitogenome phylogeny, the Mediterranean forms a separate monophyletic clade in the ddRAD BIONJ phylogeny (Fig. 3b and ESM5 Figure S5.1).
In the 'tess3R' analysis based on ddRAD SNPs, the most likely number of genetic clusters for the 123 Cuvier's was K=3 (Fig. 1a and ESM7, Figure S7.1), clearly dividing samples from the Atlantic, Mediterranean and Indo-Pacific (Fig. 1a).This was supported by statistically significant levels of pairwise differentiation between samples from the Atlantic and the Indo-Pacific in both mtDNA (mitogenome F ST =0.178, CR F ST =0.241, p < 0.01, SST5) and nuclear markers (ddRAD F ST =0.018, p < 0.01, SST6).Substantial differentiation was also seen between samples from the Mediterranean and the Atlantic and Indo-Pacific in both the mtDNA (mitogenome F ST =0.623 and 0.375, CR F ST =0.531 and 0.347, respectively, SST5) and nuclear (ddRAD F ST =0.184 and F ST =0.197, p < 0.01, respectively; SST6) markers.Net nucleotide differences showed the same pattern as the CR (SST5).
In addition to ocean-level genetic structure, the Cuvier's ddRAD SNPs revealed finer scale structure within ocean basins where the most likely number of genetic clusters within the North Atlantic was K= 5 clusters, K= 3 in the Indo-Pacific and K= 2 in the Mediterranean (Fig. 1b-e and ESM7 Figure S7.1).In the Mediterranean (Fig. 1b-c), distinct structure was present between the west and east.Within the North Atlantic, there was little admixture and individuals grouped into five discrete clusters, reflecting their general

Table 1
Nuclear diversity statistics calculated from ddRAD SNPs (Cuvier's: n = 30479, Blainville's: n = 13988): Mean and 95% confidence intervals (in brackets) are provided for each ocean basin and genetic cluster of observed heterozygosity (H o ; Nei, 1987) gene diversity (H s ; Nei, 1987), inbreeding coefficient (F is ; Weir and Cockerham, 1984) and Tajima's D (100 kb sliding window).Mitogenomic diversity statistics are provided for whole mitogenomes and mtDNA control region (in brackets) for each ocean basin: no.segregating sites (S), no. of haplotypes (h), haplotype diversity (Hd) and nucleotide diversity (π).geographical sampling location (Fig. 1d).Intriguingly, four individuals from Puerto Rico and the Virgin Islands shared low levels of admixed ancestry (<30%), suggesting connectivity with the other North Atlantic clusters.One cluster consisted of only two individuals from Spain, and was not included in further analyses.Furthermore, all pairwise ddRAD SNP F ST values were significant (p < 0.05), except for between the France and Canary Island clusters in the Atlantic (SST6).Within the Indo-Pacific, Cuvier's subclusters also corresponded with the geographical origin of samples (Fig. 1e).Here, the broadest distribution of samples was found in the South Indo-Pacific (Indo_Sou).Three admixed individuals from the Philippines, Chile and Samoa shared most of their ancestry coefficient with the remaining South Indo-Pacific individuals (South Africa, Australia, and Aotearoa New Zealand), and were not included in the remaining analyses.
As DAPC and 'tess3r' analyses produced similar results, and because 'tess3r' is expected to perform better in the presence of admixture and IBD (reviewed in Francois and Waits 2016), we present the DAPC analyses primarily in the Supplementary Material (ESM8).The first axis in the DAPC separates the two Mediterranean clusters from the North Atlantic and Indo-Pacific ones (ESM Fig. S8.1).The second axis separates clusters from the Indo-Pacific from the North Atlantic, as well as those from the east and west Mediterranean.Finer-scale differentiation between clusters within ocean basins is observed when the 2nd vs 3rd (ESM Fig. S8.3) and 3rd vs 4th (ESM Fig. S8.4) axes are plotted.Samples were reassigned to the a-priori clusters with high probability, except for a few samples in the eastern North Atlantic and Indo-Pacific clusters (S8.1).
Isolation-by-distance was present and significant when all Cuvier's samples were combined and within ocean basins (ESM6, Table S6.1).

Genetic diversity and demography
Genetic diversity summary statistics (H o , H s , F is , and Tajima's D) and their 95% confidence intervals are presented for the 118 Cuvier's that were assigned a genetic cluster (i.e., not admixed, Table 1).Overall, the inbreeding coefficient F is ranged widely amongst Cuvier's clusters within ocean basins ( F is =0.077-0.141;Table 1).Mitogenome haplotype diversity was very high, with only one haplotype shared between two individuals from the Bahamas, and one shared between one Irish and two Canary Islands individuals (SST7).However, both the overall Mediterranean dataset and populations therein had positive values of Tajima's D (indicative of population contraction), the lowest levels of diversity and inbreeding (H o , H s , and F is ; Table 1), and the lowest mtDNA diversity (π and S, Table 1, SST5).Slightly greater overall diversity values were seen in the North Atlantic, compared to the Indo-Pacific, in both the ddRAD and mtDNA data (Table 1).All sample partitions outside of the Mediterranean had negative values of Tajima's D, indicative of population expansion (Table 1).

Genetic population structure and phylogeography
The mitogenome phylogeny of 23 unique Blainville's lineages showed clear division between the Indo-Pacific and North Atlantic, which occurred ~0.75 mya (Fig. 4a).One individual from Car Nicobar, Bay of Bengal, was nested within the North Atlantic clade, near the base of the group (Fig. 4a, SST1).No ddRAD data was available for this individual to compare locations in the phylogenies.Significant genetic clustering according to ocean basin (100% bootstrap support) was also found in the BIONJ phylogeny of 42 ddRAD Blainville's and the 6 SRW used as an outgroup (Fig. 4b and ESM5, figure S5.2).One individual sampled in South Africa was nested within the North Atlantic individuals.No mitogenome data were available to compare phylogenies for this individual.
In the 'tess3r' analyses based on ddRAD SNPs, the most likely number of genetic clusters of the n = 43 Blainville's was K= 2, corresponding with North Atlantic and Pacific Oceans (Fig. 2 and ESM7 Figure S7.2).Blainville's samples from South Africa also clustered with individuals from the Pacific, forming an Indo-Pacific cluster.The pairwise difference between the North Atlantic and Indo-Pacific clusters in Blainville's was substantial and statistically significant (ddRAD F ST =0.119, p < 0.001, SST6), and greater than the parallel differentiation between North Atlantic and Indo-Pacific Cuvier's clusters (ddRAD F ST =0.018, p < 0.01, SST6).Mitogenomic differentiation between the North Atlantic and Indo-Pacific clusters was also statistically significant (p < 0.01 for all comparisons) and extremely high (mitogenome F ST =0.711, CR F ST =0.328, SST5).Net nucleotide divergence between ocean basins was d A = 0.0097 (mitogenome) and d A = 0.0017 (CR) (SST5).
Like the Cuvier's, the ddRAD SNP data revealed considerable within-ocean-basin structure, with individuals grouping best into K=2 Atlantic and K=3 Indo-Pacific genetic clusters (Fig. 2a-c).Five individuals from across the North Atlantic (Atl_Oth, Fig. 2b) clustered together with equal amounts of admixture from the Bahamas and eastern North Atlantic clusters.Blainville's from South Africa (Indo_Afr) formed a cluster separate from the large South Pacific (Indo_Sou) cluster (Fig. 2c).One sample from New Zealand was equally admixed between the South African and Hawaiian clusters and was therefore removed from the remaining analyses.All pairwise F ST comparisons of clusters between ocean basins were highly differentiated and statistically significant (SST6).
In support of the 'tess3r' analyses, the first axis in the DAPC broadly separates the North Atlantic from the Indo-Pacific clusters, with the South African cluster spanning both sides (ESM8 Fig. S8.2).The second axis separates the east from the western North Atlantic clusters, while in the Indo-Pacific it divides the Hawaiian and southern Indo-Pacific individuals.Individuals were reassigned to the a priori clusters with high probability and few misassignments (ESM8 Fig. S8.1).
Finally, IBD was found for all Blainville's samples analysed together and within ocean basins (ESM6, Table S6.1), with two clearly separated patches of points in the Atlantic Blainville's suggesting genetic sub-structuring (ESM6, Fig. S6.2).

Genetic diversity and demography
Nuclear ddRAD diversity summary statistics and their 95% confidence intervals are presented for 42 Blainville's that assigned to a cluster, and mtDNA statistics are provided for whole mitogenomes and mtDNA CR data in Table 1.The overall Indo-Pacific dataset and clusters therein had higher levels of diversity (H o and H s ) and inbreeding (F is ) compared with the North Atlantic cluster, a pattern also reflected in the mitochondrial data (Table 1).One mitogenome haplotype was shared between two individuals sampled in Hawai'i and three pairs from the Bahamas shared unique haplotypes (SST7).Both oceanbasin clusters had negative Tajima's D values, though the North Atlantic value was very close to 0 (D=− 0.073) indicating long-term population size stability.Within ocean basins, all Indo-Pacific clusters and only the admixed North Atlantic (Atl-Other) cluster had negative Tajima's D values (Table 1), indicative of increasing population sizes.The remaining North Atlantic clusters showed positive Tajima's D values, signalling population declines (Table 1).

Discussion
Here we show that two deep-sea predators, the Cuvier's and Blainville's beaked whales, exhibit substantial hierarchical population structure on a global scale and substantial regional genetic differentiation on ocean basin scale.At a macroevolutionary scale, this suggests that patterns of spatial genetic structure can vary between deep-sea species depending on their evolutionary history as other similar predators demonstrate e.g., panmixia (Westbury et al., 2021;Winkelmann et al., 2013) or IBD (Amaral et al., 2012;Gonçalves da Silva et al., 2020).At a microevolutionary scale, we propose that patterns of gene flow and genetic drift are shaped by life history traits and site fidelity in Blainville's and Cuvier's.These microevolutionary processes may be similar to those thought to influence other deep-sea marine mammals, such as island preference (Albertson et al., 2017;Van Cise et al., 2017), social structure and prey specialisation (Foote et al., 2016;Martien, Taylor et al., 2019).
Recently, genetic structure has been identified in several other beaked whale species.At a global scale, recent work indicated that True's beaked whales (Mesoplodon mirus) in the Northern and Southern Hemispheres had been separated for 0.5-2 mya, resulting in the designation of a new species in the south (M.eueu, Carroll et al., 2021).At an ocean basinlevel, population structure or IBD was detected in Baird's beaked whale (Berardius bairdii, Morin et al., 2017) across the North Pacific, and the genetic divergence measured between the "black" and "gray" stocks using the mtDNA CR was found to represent separate species within the genus, B. minimimus and B. bairdii, respectively (Morin et al., 2017;Yamada et al., 2019).The beaked whale species with the most genetic studies is the northern bottlenose whale (Hyperoodon ampullatus), which is found exclusively in the North Atlantic (Dalebout et al., 2006;Dalebout et al., 2001;de Greef et al., 2022;Feyrer et al., 2019).Northern bottlenose whales show genetic structure and IBD across their range, likely through philopatry linked to bathymetric features like the "Gully" canyon (Dalebout et al., 2006), with resulting conservation management implications (COSEWIC, 2011).Thus, though notably absent from Gray's beaked whale (Thompson et al., 2016;Westbury et al., 2021), some degree of genetic structuring appears to be present in the beaked whale species studied at global and regional scales.We add to the sparse knowledge about beaked whale population structure, by discussing drivers of Cuvier's and Blainville's genetic structure at macro-and microevolutionary scales, and discuss the conservation management implications of this work.

Macroevolutionary drivers of genetic population structure in Cuvier's and Blainville's
Biogeographical barriers on both global and regional scales (Toonen et al., 2016) have been known to influence patterns of genetic structure in widely-distributed marine species (e.g., Clarke et al., 2015;Gaither et al., 2016;Leslie et al., 2018).Such barriers can be fully impassable (if only recently in geological history; such as the Isthmus of Panama), whereas others are likely physically passable features, such as shallow rises (e.g., the Siculo-Tunisian Strait; Bianchi and Morri, 2000;Mejri et al., 2009), and expanses of deep sea (e. g., East Pacific and Mid-Atlantic Ridge Barriers; Toonen et al., 2016).At a global scale, the genetic structuring found here in Cuvier's and Blainville's appears to be shaped by such biogeographical barriers previously highlighted in the marine realm, including the Isthmus of Panama, East-Pacific Barrier, Mid-Atlantic Ridge, Benguela Shelf and Sunda Shelf.These barriers vary in permeability between our focal species, primarily due to differences in the timing and geography of each species' origin and spread across the globe (e.g., before or after the formation of the Isthmus of Panama).
Blainville's demonstrate a significant genetic distinction between North Atlantic and Indo-Pacific clusters, supporting the previous mitogenome phylogeny showing reciprocal monophyly to ocean basin (Morin et al., 2012).In contrast, Cuvier's mtDNA CR previously showed no monophyly, but significant haplotype frequency differences between basins, indicating high and significant levels of differentiation (Dalebout et al., 2005;Morin et al., 2012).With additional mitogenomes and nuclear SNPs, we were able to extend these results for Cuvier's and show ocean-specific clades, greatly refining our understanding of genetic structuring in this species.This is similar to what was found in northern bottlenose whales, where additional population structure was identified when nuclear genomes were added to the analysis (Dalebout et al., 2006;de Greef et al., 2022).
The estimated divergence dates of the deepest splits that separate primarily North Atlantic and Indo-Pacific clades (Cuvier's: ~1.5-2 mya Fig. 3b) support a role of the formation of the Isthmus of Panama (O'Dea et al., 2016) in driving genetic structure in Cuvier's, as it has in other odontocete species (Amaral et al., 2012;da Silva et al., 2015;Van Cise et al., 2019).The later divergence date for Blainville's (~0.75 mya, Fig. 4b) suggests a link to the Mid-Pleistocene Transition (MPT), which is noted for its overall cooling of sea-surface temperatures, elongation and cooling of glaciation cycles, and northward shift of Antarctic polar fronts (Clark et al., 2006;Kemp et al., 2010).The resulting cooling of tropical waters is hypothesised to have driven the speciation between M. mirus and M. eueu as well (Carroll et al., 2021).Thermal constraints linked to the smaller body size of Blainville's compared to Cuvier's may determine the former's contemporary southerly limit, which coincides with the Southern Ocean Sub-Tropical Front (STF; Graham and De Boer, 2013;Pitman et al., 2020).A northward shift of the STF during the MPT would have further limited migration between the North Atlantic and Indo-Pacific via the southern tip of South America and potentially South Africa.Therefore, we posit that oceanographic conditions and thermal constraints may have reduced gene flow between Atlantic and Indo-Pacific Blainville's populations compared with those of Cuvier's.Given the well-described relationship in mammals between body size and lifespan (Blueweiss et al., 1978), we also hypothesise that genetic drift and subsequent differentiation could have occurred more quickly in Blainville's compared with Cuvier's, although we acknowledge that many factors can contribute to genetic drift (Leroy et al., 2017).
The East Pacific and Mid-Atlantic Ridge Barriers are large stretches of deep and open water that fall within the oligotrophic centres of the ocean and may act as "nutritional" barriers.Beaked whales prefer productive habitats with complex bathymetry, including underwater canyons and deep slopes (MacLeod and D'Amico, 2006), and no sightings records of Cuvier's nor Blainville's exist within the open areas of the East Pacific basin and the Mid-Atlantic (OBIS, 2021: although such surveys are difficult to conduct effectively; (Barlow et al., 2006).The Benguela and Sunda Shelf Barriers may also influence Cuvier's and Blainville's population structure.
However, small sample sizes from the Southern Ocean mean further sampling will be needed to clarify the relevance of these barriers as the Southern basins are connected via the Antarctic Circumpolar Current, and the resulting productive convergence zones could promote connectivity of species across wide geographical expanses, limiting genetic differentiation (e.g.Gray's beaked whales; Thompson et al., 2016;Westbury et al., 2021).
Genetic discontinuities are well documented across the Strait of Gibraltar that separates the North Atlantic and Mediterranean Sea in many marine taxa (Patarnello et al., 2007), including a number of odontocetes (e.g.Gaspari et al., 2007;Kraft et al., 2020).Here, we find genetic structuring not only between the North Atlantic and Mediterranean, but also within the Mediterranean, a pattern found in other odontocetes (Gaspari et al., 2007(Gaspari et al., , 2019;;Gkafas et al., 2017).A small amount of geneflow from individuals in the North Atlantic entering the Mediterranean may be driving higher genetic diversity in the western Mediterranean populations.Alternatively, sequential founder effects may have resulted in lower diversity in the eastern Mediterranean populations.Further investigation of the evolutionary history of Cuvier's in the Mediterranean requires samples from the Alborán Sea population (Cañadas and Vázquez, 2014).

Microevolutionary drivers of genetic population structure in Cuvier's and Blainville's
Physiology, life history, behaviour, and dispersal capabilities have also likely contributed to the population structure and genetic diversity patterns of Cuvier's and Blainville's observed in our work.Cuvier's are considerably larger than Blainville's (MacLeod, 2005), make deeper dives (Baird, 2019) and occupy a greater geographical range (including the circumpolar Southern Ocean; Macleod et al., 2005).Such factors increase the potential connectivity of Cuvier's populations in the Southern Hemisphere, though the substantial genetic structure between the North Atlantic and Indo-Pacific clusters indicates the low likelihood of individuals travelling great distances between ocean basins.
Blainville's and Cuvier's often live in small, resident populations around island groups (Baird, 2019;Hooker et al., 2019), where such philopatry could enhance inter-oceanic distinctiveness, regardless of the presence of biogeographical barriers.Such site fidelity is also seen in northern bottlenose whales from the Scotian Shelf, where their preference for the submarine canyons of the "Gully" has led to genetic distinctiveness from other populations (Dalebout et al., 2006;de Greef et al., 2022).Some genetic clusters from areas of known site fidelity observed via tagging and/or long-term photo ID studies such as Cuvier's from the Canary Islands (Reyes, 2018) and Blainville's from the Bahamas (Claridge, 2006(Claridge, , 2013)), have slightly lower diversity compared to the remaining North Atlantic clusters.In contrast, tagging studies indicate that Blainville's from Hawai'i belong to either a resident or oceanic population (Baird et al., 2011), and we found these to have greater genetic diversity and lower inbreeding coefficients than their island-associated Cuvier's counterparts.Improving abundance and residency estimates in areas without dedicated photo-identification studies will help to determine the role of site fidelity and ongoing gene flow in shaping genetic diversity in these species.

Proposed management units for conservation strategies
Genetic structure has now been found in most of the beaked whale species studied to date (see discussion above), suggesting that genetically structured populations may be a common feature of the ziphiids.Many of these species surveyed to date are known to occur at low abundances in small, localized populations, which given their substantial site-fidelity and population structure, may render them particularly vulnerable to human disturbances (Hooker et al., 2019).
Maintenance of genetic diversity within populations is a priority for conservation and management measures, and small populations, such as those observed in many beaked whales (Hooker et al., 2019), are more susceptible to genetic erosion (Leroy et al., 2017).Consequences of genetic erosion include inbreeding depression, accumulation of deleterious alleles, maladaptation, and reduced adaptive potential to climate change, and hence monitoring of genetic diversity is imperative for the management and conservation of wild populations (Leroy et al., 2017).In particular, multiple populations of Blainville's and Cuvier's have been subject to repeated UMEs, including in the Ligurian Sea where more than 35 Cuvier's mortalities have occurred as a result of UMEs since the 1960s, corresponding to a loss of approximately 1/3 of the current estimated population size (Podestà et al., 2006(Podestà et al., , 2016)).Similar incidences have impacted Cuvier's in Puerto Rico and the Virgin Islands, coinciding with US Navy activity between the 1960s and 2000 (Mignucci-Giannoni, 1996;Mignucci-Giannoni and Rosario-Delestre, 1998).In the eastern North Atlantic, UMEs between 2008 and 2018 involved nearly 200 deep-diving whales, around half of which were Cuvier's (Brownlow et al., 2014(Brownlow et al., , 2018;;Dolman et al., 2010).
Given that only a small proportion of dead animals reach land during UMEs, and there are estimated to be less than 2300 Cuvier's in the eastern North Atlantic (Rogan et al., 2017), these events may have had substantial impacts on the species.The classification of taxonomic and management units below species is therefore necessary to determine how these human activities are impacting populations with respect to conservation objectives, and investigate the potential evolutionary consequences such activities may have (Taylor et al., 2017).
In Table 2, we propose both ESUs; important evolutionary components of a species which are necessary for maintaining genetic variability and successfully responding to future environmental changes (Taylor et al., 2010;Waples, 1995), and DIPs; population units that depend on internal demographic dynamics rather than immigration (Martien et al., 2019;Palsbøll et al., 2007;Taylor et al., 2010;Waples and Gaggiotti, 2006).
Nearly reciprocal monophyly to ocean basin in both the mitogenomic and ddRAD phylogenies and minimal amounts of admixture between ocean basins in the cluster analysis provide clear and strong evidence that populations are on separate evolutionary trajectories, consistent with distinct ESUs.To quantify and further demonstrate the distinctiveness of ocean-basin-level populations, we used mean F ST and d A as evidence of population isolation (Taylor et al., 2017;Waples et al., 2018).
For both species, nuclear SNP based F ST (Cuvier's F ST =0.018-0.197,Blainville's F ST =0.119) values fall within the range of F ST values calculated from nuclear SNPs in recognized odontocete subspecies (spinner dolphins (Stenella longirostris): F ST =0.0045-0.012;pantropical spotted dolphins (S. attenuata): F ST =0.055; harbour porpoise (Phocoena phocoena): F ST = 0.187 − 0.207; and short-finned pilot whale (Globicephala macrorhynchus): F ST = 0.1 − 0.4; Lah et al., 2016;Leslie and Morin, 2018;Van Cise et al., 2019).The F ST of ocean-basin level populations using full mitogenomes (Cuvier's F ST =0.178-0.623,Blainville's F ST =0.711) are much higher than mitogenomic F ST values calculated in recognized odontocete subspecies (spinner dolphins: F ST =0.013, pantropical spotted dolphins: F ST =0.013, fin whales (Balaenoptera physalus): F ST = 0.005 − 0.018; Archer et al., 2019;Leslie et al., 2018).The F ST of the ocean-basin populations calculated using the mtDNA CR (Cuvier's F ST = 0.241-0.531,Blainville's F ST =0.628) is also greater than the range of F ST values in recognized cetacean subspecies (F ST =0.013-0.209;Rosel et al., 2017).In contrast, our estimates of d A did not reach the previously described minimum subspecies threshold of 0.004 (Taylor et al., 2017).The exception to this was the Indo-Pacific and Mediterranean Cuvier's comparison (Cuvier's d A =0.0045, F ST = 0.531).Beaked whales show a contrasting pattern to other cetaceans, where the CR is less diverse compared to other regions of the mitochondrial genome (Dalebout et al., 2004;Morin et al., 2012).It is possible that d A measurements using only the CR are not as appropriate for beaked whales as for other cetaceans, and additional evaluation of mitogenomic divergence across known taxonomic groups would be useful to further resolve divergence levels and taxonomy within and among beaked whale species.
The proposed ESUs show different patterns of expansion and contraction and may therefore respond differently to UMEs and other human activities (Taylor et al., 2017).For example, both proposed Mediterranean ESUs have low genetic diversity and historical demography indicative of population contraction, in contrast to Cuvier's in other regions.Furthermore, in both species, individuals from the Southern Hemisphere showed genetic divergence from the Northern Hemisphere and while sample sizes in the former are small, Southern Hemisphere Cuvier's and Blainville's should be managed as distinct ESUs until confirmation of genetic distinctiveness with further sampling.
The designation of DIPs stems from multiple lines of evidence (Martien et al., 2017).Based both on the statistically significant genetic differentiation shown here and the results of long-term studies on resident populations (recently reviewed by Hooker et al., 2019), we suggest the following DIPs for future conservation and management initiatives: Bahamas, Canary Islands and Hawai'i in both Cuvier's and Blainville's, and separate Mediterranean Cuvier's populations in the Ligurian and Ionian seas (Table 2).We acknowledge the limitations of our ability to identify DIPs given the available dataset, as increasing the number and distribution of samples will likely increase the number of identifiable DIPs.However, Martien et al. (2019) advise that for rare or elusive species such as beaked whales, the sample size used to infer management units should not only be assessed with respect to statistical power and representativeness, but the length of time it could take to increase the sample size.Given that the samples in the current study span 140 years, we recommend application of the precautionary principle by incorporating the management units suggested here into conservation and management activities in light of the anthropogenic pressures that these species face (Moore and Barlow, 2013;Parsons, 2016).

Conclusions
We show that deep-sea predators can have substantial genetic structuring across geographical scales, likely driven by a combination of macroevolutionary (e.g., biogeographical boundaries, body size and thermal preference) and microevolutionary (e.g., philopatry and dispersal) processes.This contributes to the growing body of knowledge of drivers of genetic patterns in deep sea species.The dataset compiled here is the largest for any published analysis of Cuvier's or Blainville's, and is unlikely to be matched soon due to the inherent difficulties in sampling beaked whales (Baird et al., 2020).Indeed, this study would not have been possible without this global collaboration providing samples from a variety of sources.We hope to encourage future cooperation between researchers, Indigenous groups, conservation organisations, stranding networks, and governments to continue collecting and archiving tissue from as many elusive, deep sea species as possible to facilitate future studies.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Sampling locations and genetic clustering of Cuvier's beaked whales (Ziphius cavirostris, n = 123 individuals, n = 30479 SNPs).Main map shows sample locations plotted on the Spilhaus projection, to highlight the connectivity and continuity of the world's oceans for a globally distributed species, with names and locations of relevant major biogeographical barriers.Bar plots of co-ancestry coefficients generated in 'tess3r' are presented for (a) all samples, (c) Mediterranean (location shown more detail in (b) map inset), (d) Atlantic and (e) Indo-Pacific.The three Indo-Pacific samples enclosed in a red box represent the admixed individuals and are plotted on the map as red points.

Fig. 2 .
Fig. 2. Distribution and genetic clustering of Blainville's beaked whales (Mesoplodon densirostris, n = 43 individuals, n = 13988 SNPs).Main map shows sample locations plotted on the Spilhaus projection, to highlight the connectivity and continuity of the world's oceans for a globally distributed species, with names and locations of relevant major biogeographical barriers.Bar plots of co-ancestry coefficients generated in 'tess3r' are presented for (a) all samples, (b) Atlantic and (c) Indo-Pacific.

Fig. 3 .
Fig. 3. Phylogenetic relationship among globally distributed and sampled Cuvier's beaked whales (Ziphius cavirostris).(a) Bayesian mitogenome phylogeny generated from 32 unique Cuvier's mitogenome haplotypes with posterior estimates (as a proportion) with divergence shown in million years ago and the purple bars represent the 95% probability range of divergence dates (HPD).(b) BIO neighbour-joining (BIONJ) phylogeny of 118 Cuvier's beaked whales with bootstrapped support for each node > 50% (n = 1000 bootstraps) generated using 33137 nuclear SNPs.The six southern right whales (Eubalaena australis) used to root the tree were removed prior to plotting (see ESM5 Fig. S5.1 for tree with outgroup).

Fig. 4 .
Fig. 4. Phylogenetic relationship among globally distributed and sampled Blainville's beaked whales (Mesoplodon densirostris).(a) Bayesian mitogenome phylogeny generated from 23 unique Blainville's mitogenome haplotypes with posterior estimates (as a proportion) with divergence shown in million years ago (ma) and the purple bars represent the 95% probability range of divergence dates (HPD).(b) BIO neighbour-joining (BIONJ) phylogeny of 42 Blainville's beaked whales with bootstrapped support for each node > 50% (n = 1000 bootstraps) generated using 29904 SNPs.The six southern right whales (Eubalaena australis) used to root the tree were removed prior to plotting (for tree with outgroup see ESM5, Fig. S5.2).

Table 2
Sampling locations, sample size and 'tess3r' population ID across and within ocean basins for Cuvier's (Ziphius cavirostris, n = 123 individuals, n = 30479 SNPs) and Blainville's beaked whales (Mesoplodon densirostris, n = 43 individuals, n = 13988 SNPs).Proposed ESUs formed distinct genetic clusters in 'tess3r' that corresponded to distinct geographic regions and significant ddRAD F ST > 0.01.Proposed DIPs were based on resident populations investigated by long-term photo ID and telemetry studies and differentiated by significant F ST .Abbreviations include Med for Mediterranean and S. Hem for Southern Hemisphere, and cardinal compass points abbreviated to their first letter (e.g., west is W).