Postglacial Colonization of Northern Coastal Habitat by Bottlenose Dolphins: A Marine Leading-Edge Expansion?

Oscillations in the Earth’s temperature and the subsequent retreating and advancing of ice-sheets around the polar regions are thought to have played an important role in shaping the distribution and genetic structuring of contemporary high-latitude populations. After the Last Glacial Maximum (LGM), retreating of the ice-sheets would have enabled early colonizers to rapidly occupy suitable niches to the exclusion of other conspecifics, thereby reducing genetic diversity at the leading-edge. Bottlenose dolphins (genus Tursiops ) form distinct coastal and pelagic ecotypes, with finer-scale genetic structuring observed within each ecotype. We reconstruct the postglacial colonization of the Northeast Atlantic (NEA) by bottlenose dolphins using habitat modeling and

phylogenetics. The AquaMaps model hindcasted suitable habitat for the LGM in the Atlantic lower latitude waters and parts of the Mediterranean Sea. The time-calibrated phylogeny, constructed with 86 complete mitochondrial genomes including 30 generated for this study and created using a multispecies coalescent model, suggests that the expansion to the available coastal habitat in the NEA happened via founder events starting ~15 000 years ago (95% highest posterior density interval: 4 900-26 400). The founders of the 2 distinct coastal NEA populations comprised as few as 2 maternal lineages that originated from the pelagic population. The low effective population size and genetic diversity estimated for the shared ancestral coastal population subsequent to divergence from the pelagic source population are consistent with leading-edge expansion. These findings highlight the legacy of the Late Pleistocene glacial cycles on the genetic structuring and diversity of contemporary populations.
Subject areas: Molecular systematics and phylogenetics, Population structure and phylogeography Key words: genetic diversity, habitat modeling, Last Glacial Maximum (LGM), multispecies coalescent, phylogenetics, timedependency During the Late Quaternary period (1 Ma to present), the Earth's climate was governed by a series of glacial and interglacial events and temperature fluctuations that occurred at approximately 100 000-year intervals (Shackleton 2000). These glacial cycles are thought to have played an important role in shaping the current distribution and genetic structuring of species and populations. In the Northern Hemisphere, this is more pronounced at high latitudes where the presence of ice sheets and Arctic temperatures during cold stadial periods restricted the available habitat to warmer refugia for many temperate-adapted species (Hewitt 2000). The distribution of temperate-adapted species at high latitudes is thought to have been characterized by cyclical range contractions and expansions throughout the Pleistocene (Hewitt 2000), and present-day populations represent lineages descended from relict populations that survived in glacial refugia (Hofreiter and Barnes 2010). However, the capacity to adapt to environmental change varies among species and populations through their dispersal ability, genetic diversity, and generation time (Stewart et al. 2010;Montgomery et al. 2014;Younger et al. 2016) and has likely affected their present-day distribution. Hewitt (1999Hewitt ( , 2000 proposed that colonization of terrestrial species during the warm interglacial periods typically occurred via long-range dispersal events, or "leading-edge range expansions". During these events, a few individuals from the source population, carrying a subset of the genetic diversity, first colonized and then expanded their range to fill emerging geographic and ecological niches leading to the exclusion of secondary waves of colonizers and thereby reducing genetic variability at this leading-edge. Unlike most terrestrial mammal species, marine mammals have a relatively low cost of movement (Tucker 1975;Williams et al. 1992;Williams 1999) and few geographic barriers to dispersal, so the leading-edge colonization model may, or may not, be realistic for such highly mobile organisms that are able to move thousands of kilometers within a few weeks (Gabriele et al. 1996;Mate et al. 1997). High mobility certainly makes them capable of undertaking long-distance dispersal events, but at the same time, also increases the potential to retain ongoing gene flow and migration between the leading-edge and source populations. Compared to terrestrial species, little is known about the post-glacial colonization history of Europe's marine fauna, but during recent years a growing number of studies have explored the dynamics of postglacial range expansion in marine taxa (e.g., de Bruyn et al. 2009;Fontaine et al. 2010Fontaine et al. , 2014Catchen et al. 2013;Foote et al. 2013). Genetic analyses of carbon-dated subfossils indicate that temperate climate-adapted species such as gray whales Eschrichtius robustus (Lilljeborg, 1861) and North Atlantic right whales Eubalaena glacialis (Müller, 1776) replaced cold-adapted species such as bowhead whales Balaena mysticetus (Linneaus, 1758) in mid-latitude European waters during the Late Pleistocene-Early Holocene transition (Foote et al. 2013;Alter et al. 2015). However, the refugial distribution during the Last Glacial Maximum (LGM) of the source populations for these recolonizing temperate marine species is largely unknown.
Here, a multispecies coalescent model was used together with information on species and population structure to reconstruct the postglacial colonization in a highly mobile marine genus Tursiops. Specifically, the aim was to determine the timing of the colonization into the northeastern Atlantic coastal habitats and to investigate whether the colonization occurred via leading-edge expansion. The genus Tursiops currently consists of two separate species recognized by the International Union for the Conservation of Nature, T. truncatus (Montagu, 1821; common bottlenose dolphin) and T. aduncus (Ehrenberg, 1833; Indo-Pacific bottlenose dolphin), although T. australis (Burrunan dolphin) inhabiting the coastal waters in southern Australia has been proposed as a separate species (Charlton-Robb et al. 2011). For simplicity, we refer to the latter as T. australis throughout this paper but acknowledge the unofficial status of this putative species. Species in the Tursiops genus are found in coastal inshore waters, continental shelf regions and open ocean environments (Wells and Scott 2009) in tropical and temperate waters (Leatherwood and Reeves 1989). Throughout much of their range, they exhibit a hierarchical population structure, with the greatest genetic divergence found between pelagic and coastal populations (Hersh and Duffield 1990;Curry and Smith 1998;Hoelzel et al. 1998;Natoli et al. 2005;Tezanos-Pinto et al. 2009;Louis et al. 2014b;a;Lowther-Thieleking et al. 2015;Nykänen et al. 2018). The broad-scale population structuring between coastal and pelagic ecotypes in the Northeast Atlantic (NEA) (Louis et al. 2014b) likely reflects colonization of emerging or available inshore habitats after the LGM from a pelagic source population followed by divergence in allopatry (Moura et al. 2013;Louis et al. 2014a). It has been suggested that the colonization of coastal habitats in the NEA was possibly achieved via a founder event by a small number of individuals (Louis et al. 2014a). However, the possible role of leading-edge colonization has not been fully investigated in this genus by quantifying historic population sizes. Therefore, to better understand the climatic, temporal, and spatial context of the evolutionary processes that gave rise to the present-day diversity and habitat use of bottlenose dolphins in the NEA, mitogenomic data were analyzed in combination with habitat mapping, using newly available tools that can robustly infer ancestral effective population sizes and therefore give insights into colonization processes.

Models of Suitable Habitat
Following the methods in Foote et al. (2013), Alter et al. (2015), and Morin et al. (2015), the AquaMaps modeling approach (Ready et al. 2010;Kaschner et al. 2011) was used to predict the distribution of suitable habitat for common bottlenose dolphins during the present, for the year 2100 (under the Institut Pierre Simon Laplace climate scenario SRES A2), and for the LGM (~20 000 ya). AquaMaps is a bioclimatic model that combines large sets of occurrence data (e.g., from visual observations) with available expert knowledge on species preference and tolerance to different environmental parameters and ultimately generates large-scale predictions of the probability of occurrence for different marine species. This way the preferred habitat of a species can be estimated based on a predefined set of environmental parameters including water depth, sea surface temperature, salinity, primary production, sea ice concentration, and proximity to land, and projected into geographic space as relative probability of occurrence. For the purpose of this study, a slightly modified version of the original AquaMaps model (which is available on www.aquamaps.org) was used; specifically, primary production was excluded from the model, as there are no available data for this parameter for the Pliocene or Pleistocene. Current environmental conditions were assumed to be representative of the entire Holocene as the variability in conditions during this time period has been small compared to the differences between glacial cycles (Folland et al. 2002). AquaMaps has been previously used to hindcast suitable habitat predictions for bowhead whales, gray whales, and killer whales, Orcinus orca (Linnaeus, 1758), during the LGM (Foote et al. 2013;Alter et al. 2015;Morin et al. 2015) by using mean annual environmental conditions during the LGM based on the GLAMAP project data set (Schäfer-Neth and Paul 2003).

Sample Collection
Tissue samples from 12 individuals from T. truncatus were obtained by biopsy sampling free-ranging animals (see Krützen et al. 2002) in coastal waters of Ireland and France in 2005-2012. Additionally, skin samples were collected from 21 individuals that stranded along the coast of Ireland, France, and the United Kingdom in 1991-2010 ( Figure 1). Samples were specifically selected for sequencing based on previous genetic assignment to distinct ecotypes/populations   Mirimin et al. 2011;Louis et al. 2014b) and the lack of evidence for recent migration between the populations (Louis et al. 2014b). For consistency, we refer to individuals not assigned to a coastal population as belonging to a "pelagic" ecotype. However, we acknowledge that little is known about the ecology of these individuals, and they may inhabit both neritic and oceanic waters. In addition to the 30 samples successfully sequenced for the study, 57 mitochondrial haplotypes from the genus Tursiops from other regions of the world, along with 17 delphinid sequences, were downloaded from GenBank (Appendix S1, Supplementary Material).

DNA Extraction, Amplification, and Sequencing
DNA was extracted from 33 tissue samples using the Qiagen DNeasy (Qiagen DNeasy, Valencia, CA) kit following the manufacturer's guidelines. DNA yield was quantified using a Qubit (Life Technologies) and ranged between 10 and 300 ng/μL for all samples. The DNA samples were then sheared to fragments of ~150-200 bp using a Diagenode Bioruptor NGS run with 20 cycles of 30 s on and 30 s off.
To generate mitochondrial genome sequences, we employed a simple shotgun sequencing method (see Tilak et al. 2015). Illumina sequencing libraries were built on the sheared DNA extracts using NEBNext (Ipswich, MA) DNA Sample Prep Master Mix Set 1 following Meyer and Kircher (2010). Libraries were subsequently index amplified for 15 cycles using Phusion High-Fidelity Master Mix (Finnzymes) in 50-μL reactions following the manufacturer guidelines. The libraries were then purified using a MinElute PCR purification kit (Qiagen, Hilden, Germany). Concentrations of amplified libraries were initially checked using a Qubit (Life Technologies) and fragment size distribution was visualized on agarose gel, before a 1/10 diluted aliquot was run on Agilent Bioanalyser 2100 (Palo Alto, CA) to determine molarity and concentration and facilitate equimolar pooling of index amplified libraries. The index amplified libraries were then sequenced in subpartitions of single channels on an Illumina HiSeq 2500 ultra-high-throughput sequencing platform using single read (SR) 100-bp chemistry.
Conversion of Illumina's *.bcl files to fastq, and demultiplexing were performed using Illumina's CASAVA (version 1.8.2) software allowing for no mismatch in the 6-nucleotide indices used for barcoding. Sequencing reads within the generated fastq files were processed with ADAPTER-REMOVAL (Lindgreen 2012) to trim residual adapter sequence contamination and to remove adapter dimer sequences as well as low-quality stretches at 3′ ends (i.e., consecutive stretches of N's and of bases with a quality score of 2 or lower). Reads that were ≤30 bp following trimming were discarded. Filtered reads with average phred-scores ≥30 were then mapped to a reference T. truncatus mitogenome sequence (KF570351.1) using BWA version 0.5.9 (Li and Durbin 2009), requiring a mapping quality of Q ≥ 30 for all bases. The reference mitochondrial sequence was modified as per Morin et al. (2015) to improve assembly coverage at the "ends" of the linearized mitogenome by adding 40 bp from each end to the opposite end (so that reads could map across the artificial break point of the linearized sequence). The influence of potential NUMTs was excluded given the absence of stop codons or frame shifts in the aligned protein encoding genes, and the lack of regions of excessive coverage due to the mapping of both mt and nuDNA. Clonal reads were collapsed using the rmdup program of the SAMTOOLS (version 0.1.18) suite ). Consensus mitogenome sequences were then reconstructed using bam files, which were assembled and visualized in GENEIOUS (Biomatters Ltd.), allowing indels and unique variants to be visually verified in the BAM files, and only including nucleotide positions with a read depth of ≥3× coverage; ambiguous bases were assigned as N's.

Time-Calibrated Phylogeny of Delphinids and the Genus Tursiops
There were insufficient data for the reconstruction of complete mitochondrial genomes for 3 sequenced samples, which were therefore excluded from further analyses. The remaining 30 whole mitochondrial T. truncatus genomes were used in the subsequent analyses. The inference of time-calibrated phylogenies followed a 2-step methodology following Morin et al. (2015) with the purpose of the first step (delphinid phylogeny) to estimate a calibration point (divergence of the 2 most evolutionarily distant Tursiops spp. haplotypes) for the second step (Tursiops spp. phylogeny).
The 30 whole mitochondrial T. truncatus genomes sequenced for this study were first aligned with 56 sequences from the genus Tursiops and a rough-toothed dolphin (Steno bredanensis) sequence downloaded from GenBank (see Appendix S1, Supplementary Material for sequence accession numbers) that was used as an outgroup using the MUSCLE algorithm in MEGA 6 (Tamura et al. 2007). A topology tree in MRBAYES (Huelsenbeck and Ronquist 2001;Ronquist and Huelsenbeck 2003) was then built with these sequences after initial model selection for the best substitution scheme in JMODELTEST2 (Posada 2008;Darriba et al. 2012), and the resulting consensus tree was inspected to find the two most divergent Tursiops haplotypes (SABD1, representing a haplotype of T. australis obtained from South Australia, and SA99 representing a haplotype of T. aduncus collected from South Africa (see Moura et al. 2013). The 13 protein coding genes of the mitochondrial genome from these 2 most divergent Tursiops sequences were then aligned with the additional 18 delphinid sequences downloaded from GenBank (Appendix S1, Supplementary Material). We used the package ClockstaR (Duchêne et al. 2014) in R 3.3.3 (R Core Team 2016) to investigate whether the different mitochondrial genes had different mutation rates or whether a single rate could be applied in further time-calibrated models.
The "greedy" search in PARTITIONFINDER (v1.1.0) (Lanfear et al. 2012) was used to find the best partitioning scheme for each gene, and a time-calibrated phylogenetic tree was built in BEAST2 (Drummond et al. 2012) with three different data partitions for nucleotide substitution models (Appendix S2, Supplementary Material). We used a single clock model based on the ClockstaR results and estimated a single tree due to the nonrecombining properties of mtDNA to find the time to most recent common ancestor (TMRCA) between the 2 most divergent Tursiops haplotypes. The TMRCA of Delphinidae (mean = 10.08 Myr, SD = 1.413) (McGowen et al. 2009) was used to calibrate the root of the tree and a Calibrated Yule prior was used for the branching rate. Both uncorrelated log-normal relaxed clock (Drummond et al. 2006) and strict clock models were tested with 40 000 000 Markov Chain Monte Carlo (MCMC) steps, 25% pre-burn-in, and a sampling frequency of 4 000. Each model was run twice and the convergence of chains and the Effective Sample Size (ESS) values relating to the model parameters were checked in TRACER (Rambaut et al. 2014). After verifying convergence, LOGCOMBINER and TREEANNOTATOR (Drummond et al. 2012) were used to combine and summarize the trees, respectively. Model selection between the two clock models was done by inspecting the ucldStdev parameter in the relaxed clock model; a standard deviation close to zero (<0.1) in this parameter indicates negligible variation in the substitution rates across branches and a better fit of the strict clock model (Drummond and Bouckaert 2015).
We tested for potential saturation (multiple substitutions) of the third codon positions by plotting the uncorrected pairwise genetic distances and model-corrected (Tamura-Nei substitution model, Tamura and Nei 1993) genetic distances in R. To investigate the potential effect of incomplete purifying selection, we compared the node TMRCA estimates and clock rates derived using all codon positions (as described above) to a time-calibrated phylogeny created using only the third codon positions, as most mutations that occur in third codon positions are silent and therefore code for the same amino acid (Lagerkvist 1978). This makes them less likely to be directly targeted by purifying selection (but does not account for purifying selection upon a linked site), thereby maintaining a more constant substitution rate over evolutionary time. The approach of using only the third codon positions has been recently used in a study resolving killer whale (Orcinus orca) mitochondrial lineages (Morin et al. 2015).
The best nucleotide partitioning scheme for the coalescent analysis of the genus Tursiops was determined for the data set, consisting of 13 protein coding gene regions of the 30 mitochondrial genomes sequenced for this study supplemented by 56 Tursiops haplotypes downloaded from GenBank (Appendix S1, Supplementary Material), using PARTITIONFINDER (v1.1.0) (Lanfear et al. 2012) (Appendix S2, Supplementary Material). These samples/haplotypes from T. truncatus, T. aduncus, and T. australis had been determined to have originated from coastal or pelagic regions in different parts of the world in previous studies and had been delineated into separate populations or ecotypes (Figure 1; Appendix S1, Supplementary Material). This information on population structure was used in this study together with mitochondrial gene regions to create a timecalibrated phylogeny for the genus Tursiops and to estimate the divergence times for the different populations of T. truncatus in the North Atlantic.
As the multispecies coalescent (MSC) model has been found not only to delimit species but also to identify genetic structuring between putative populations (Sukumaran and Knowles 2017), this approach was chosen to construct the time-calibrated phylogeny using StarBEAST2 (Ogilvie et al. 2017), an extension of BEAST2 (Drummond et al. 2012). This method is thought to provide accurate inferences of species tree topologies, divergence times, and substitution rates even in the presence of incomplete lineage sorting, that is, when individual gene trees do not correspond well with the species tree (Ogilvie et al. 2017). For this reason, information on population structure obtained from previous studies on nuclear microsatellite markers were included in the MSC model. To estimate historical effective population sizes (N e ) for the populations sequenced for this study, the MSC model was run with the population model "Linear with Constant Root Populations" assuming that population sizes change over time but keeping the size of the root (the ancestor of all species/populations) constant. The scaled population sizes for pelagic Northeast Atlantic, Coastal North, and Coastal South produced by the StarBEAST2 models were adjusted by the average generation time of 21 years estimated for T. truncatus and T. aduncus (Taylor et al. 2007) to convert the estimates to population sizes. As with the delphinid phylogeny, to investigate the possible effect of incomplete purifying selection, StarBEAST2 was run with two different models: a model including all codon positions and a second model including only the third codon positions. Both models were run with an uncorrelated log-normal relaxed clock (Drummond et al. 2006) and a strict clock and with 200 000 000 MCMC steps, 10% burn-in, and a sampling frequency of 20 000. Based on the results of ClockstaR analysis (Duchêne et al. 2014) on delphinids, a single clock model was used for all the data partitions. We applied a constraint on the root age obtained in the delphinid phylogeny (split of the 2 most divergent Tursiops haplotypes; 3.662 Myr [95% highest posterior density interval,  with the all codon model, 3.692 Myr [95% HPDI: 2.567-4.153] with the third codon model) to calibrate the trees. The convergence of chains and model performance was inspected in TRACER (Rambaut et al. 2014) after running each model twice, and LOGCOMBINER and TREEANNOTATOR (Drummond et al. 2012) were used to combine the log-and tree-files and summarize the trees, respectively. The resulting summary trees were plotted using phytools (Revell 2012) in R.
Using DnaSP 5.0 (Librado and Rozas 2009), we calculated the number of segregating sites (k) and 4 different indices of historical population demographic changes (bottlenecks or expansions); Tajima's D (Tajima 1989), Fu and Li's D* and F* (Fu and Li 1993), and Ramos-Onsins & Rozas' R2 (Ramos-Onsins and Rozas 2002), for the pelagic, the Coastal North and the Coastal South populations. The examination of deviation from the neutral theory model at equilibrium between mutation and genetic drift by these indices was based on 10 000 coalescent simulations. Expectations of these statistics are ~zero in a constant-size population, whereas significant negative or positive values can indicate a sudden expansion or a recent population bottleneck, respectively. In addition, differences in nucleotide and haplotype diversity (π nuc and π hap , respectively) were compared between pelagic and the Coastal North and Coastal South populations with the "genetic_diversity_diffs" function (Alexander et al. 2016) implemented in R and using 10 000 resampling permutations, in order to test whether the observed differences in diversity between the populations were greater than expected by chance. A haplotype network of the whole mtDNA sequences (16 390 bp) and including all the North Atlantic T. Tursiops samples was created in popART (Leigh and Bryant 2015) using the median-joining network (Bandelt et al. 1999).

Models of Suitable Habitat
According to the AquaMaps model, the core suitable habitat for common bottlenose dolphins during the LGM ranged from ~42°N (approximately the latitude of the state of Connecticut) in the Northwest Atlantic southwards along the coast to lower latitudes ( Figure 2; Appendix S3, Supplementary Material). In the Northeast Atlantic, their coastal range was restricted to approximately 40°N, corresponding to the south of the Iberian Peninsula, and to an area around the Alboran Sea in the Mediterranean. In addition, suitable habitat was hindcasted to cover an area in the pelagic North Atlantic around the Mid-Atlantic ridge and around the islands on both sides of the Atlantic. Note that AquaMaps hindcasted suitable habitat also in the Black Sea, however, this remained closed and thus inaccessible for the dolphins until after 10 kya (Kerey et al. 2004).

DNA Sequencing
A total of 292 × 10 6 sequencing reads were generated from the samples; >1 million reads per individual for 28 of the 33 individuals included on the sequencing lane. Following QC filtering, removal of duplicate reads, mapping, and only including nucleotide positions with a read depth of ≥3× coverage, complete mitochondrial genomes were assembled for 30 individuals with a mean sequence coverage ranging from 10× to >100×. Control region sequences were compared with those generated for the same individuals by Sanger sequencing in previously published studies (Islas-Villanueva 2010; Mirimin et al. 2011;Louis et al. 2014b), in some cases the same individual had been previously sequenced independently by all 3 studies. In all cases, our high-throughput sequencing generated data agreed 100% with the Sanger sequence data. Twenty new haplotypes from the samples sequenced for this study were identified (Appendix S1, Supplementary Material), which are deposited in GenBank (KT601188-KT601207).

Time-Calibrated Delphinid Phylogeny
The approximately linear relationship between the uncorrected and model-corrected genetic distances indicated no evidence of saturation of the third codon position of the protein coding genes (Appendix S4, Supplementary Material). Low standard deviation (<0.1) of the clock rate parameter (ucldStdev) in both delphinid models, the all codon model and the third codon model, run with an uncorrelated log-normal relaxed clock indicated low variation in substitution rates between lineages; therefore, the models run with a strict clock were considered to be a better fit to the data. Average substitution rate for the strict clock all codon model was estimated as 6.663 × 10 −3 substitutions/site/Myr (95% HPDI: 4.872 × 10 −3 -8.865 × 10 −3 ), and for the third codon model as 0.017 substitutions/site/Myr (95% HPDI: 0.012-0.022). Both models produced similar topology and node age estimates for the clade including the Tursiops and Stenella species, Sousa chinensis, Steno bredanensis, and Delphinus capensis (Appendix S5, Supplementary Material). However, Lagenorhynchus albirostris and Cephalorhynchus heavisidii were placed in different clades in the trees depending on the model used and the models also estimated the divergence of killer whales (Orcinus orca) differently. The placement of these species has been uncertain in previous studies (LeDuc et al. 1999;Caballero et al. 2008;McGowen et al. 2009;Steeman et al. 2009;Vilstrup et al. 2011) and is thought to reflect rapid radiation in the Delphinidae family (Vilstrup et al. 2011). We compared the time to most recent common ancestor (TMRCA) estimates for the common nodes derived using the 2 models and found them to be significantly correlated (Pearson's r 2 = 0.99, P < 0.001; Appendix S6, Supplementary Material) with a slope of 1.01, implying that incomplete purifying selection was not an issue in the delphinid phylogeny. The TMRCA of the two most divergent Tursiops samples was estimated as 3.662 Myr (95% HPDI: 2.567-4.725) with the all codon model and as 3.692 Myr (95% HPDI: 2.567-4.772) with the third codon only model. The Effective Sample Size (ESS) values for the different parameters were all >500 in the third codon only model (most of them >2000), whereas in the all codon model they were >200 (most of them >2000), except for one parameter for the gamma shape and one parameter for the proportion of invariant sites, which were both >150 (Appendix S7, Supplementary Material).

Time-Calibrated Tursiops Phylogeny
Very low standard deviation (<0.01 in both all codon and third codon only models) in the substitution rates estimated with the uncorrelated log-normal relaxed clock indicated negligible variation in substitution rates between lineages; therefore, in StarBEAST2, a strict clock was used to estimate the time-calibrated phylogeny of the genus Tursiops. The ESS values for the parameters in both models were all >300 (most of them >6000), indicating no sign of autocorrelation between samples and a good convergence of chains (Appendix S8, Supplementary Material). Both models resulted in identical species/population topologies and similar node age estimates and posterior probabilities (Appendix S9, Supplementary Material). Similar to the Delphinid phylogeny, the node ages were significantly correlated with a slope of 1.01 (Pearson's r 2 > 0.99, P < 0.001; Appendix S10, Supplementary Material), indicating no evidence of incomplete purifying selection affecting the Tursiops tree.
The summary consensus species tree, inferred using the third codon only model and showing only the NEA and Mediterranean and Black Sea T. truncatus populations (Figure 3), indicates that the branch including haplotypes from the Mediterranean and Black Seas split from the NEA branch ~ 35 900 years ago (95% HPDI: 11 200-60 500). The coastal NEA populations, on the other hand, diverged from the pelagic population around 14 800 years ago (95% HPDI: 4 900-26 400) with a further more recent divergence of the northernmost coastal population (Coastal North) occurring ~7 200 years ago (95% HPDI: 0-14 200). However, the node posterior probabilities of the pelagic and coastal NEA populations are quite low (Figure 3), likely due to incomplete lineage sorting. Indeed, incomplete lineage sorting, especially of the coastal North Atlantic T. truncatus sequences, is evident in the consensus gene tree consisting of all samples from the genus Tursiops (green tips in the tree; Appendix S12, Supplementary Material), and is accompanied by older node coalescence times compared to the time-calibrated species tree (Figure 3). The consensus gene tree also indicates that the coastal bottlenose dolphin population currently occupying the northern parts of the NEA (i.e., Coastal North) originates from at least 2 separate mitochondrial lineages.
The species tree depicting the effective population sizes (N e ) through Late Pleistocene-Holocene (Figure 3; see also Appendix S11, Supplementary Material, for population size estimates) indicates that, unlike the pelagic population, which is estimated to have retained a relatively stable effective size throughout the Late Pleistocene and Holocene, all of the present-day coastal populations, and especially the Black Sea and the Coastal North populations, are inferred to have undergone a genetic bottleneck around 7 000-8 000 years ago followed by expansion towards the present day ( Figure 3; Appendix S11, Supplementary Material). The significantly lower nucleotide (π nuc ) and haplotype diversities (π hap ) in the Coastal North population compared to the pelagic population, and the significantly lower π nuc in the Coastal South population compared to the Pelagic East Atlantic (all P < 0.05, based on 10 000 permutations) are further indications of ancestral population bottlenecks in the coastal populations, as are the reduced number of segregating sites (k) (see Table 1). Conversely, the nonsignificant values of Tajima's D, Fu and Li's D* and F*, and Ramos-Onsins & Rozas' R2, albeit the first 3 indicators being negative for the Coastal South population (Table 1), fail to provide support for our inference of postglacial demographic expansion from a small founder population and indicate that there is no significant excess of low frequency variants in the Coastal North and Coastal South populations. This could be an artifact of lack of power to detect significance, or it could be due to a violation of the assumption of panmixia inherent in these statistics. For example, the Coastal North population consists of demes that are interconnected by limited and varying gene flow (Louis et al. 2014b;Nykänen et al. 2018). This could increase the number of rare variants, but it could also inflate intermediate variants; for example, variants that are at high frequency in 1 or 2 demes would be identified as intermediate frequency variants in our sample set. The possible unaccounted population structuring, in turn, would inflate the effective population sizes, which could help explain the unrealistically high recent time N e values estimated for the coastal populations. Furthermore, the fact that the divergence of the coastal populations from the pelagic population was estimated to have occurred already ~15 000 years ago could mean that the rare variants expected to arise in a population after a bottleneck or founder event have had sufficient time to drift to extinction or intermediate frequency.
The haplotype network shows the 2 source lineages of the haplotypes shared among different NEA coastal populations (Figure 4). One of the haplotypes is shared between 6 individuals from 2 different populations, 5 from the Coastal North and 1 from the Coastal South, and another haplotype is private to the Coastal North population and is shared by 5 individuals. Thus, 71% of the Coastal North individuals sampled in this study shared just these 2 haplotypes; the remainder have haplotypes recently derived from these ancestral lineages. This is consistent with the demographic expansion of the Coastal North population from just 2 founding lineages and a probable ancestral origin in the pelagic population. Similarly, the Coastal South population may have originated from only 1 or 2 pelagic founding lineages; however, sequencing more samples from this population will be required to determine this with certainty. Nevertheless, it is likely that the pelagic population gave rise to coastal lineages with current populations located in the NEA as well as in the Mediterranean and Black Seas.

Discussion
Our analyses indicate that bottlenose dolphins have colonized the northernmost range of their distribution through a leading-edge expansion, in which the early colonists have expanded to fill the available coastal niches, to the exclusion of others. Due to its strict maternal inheritance, the mitochondrial genome tracks only the matrilineal population history, a very stochastic coalescent process, and generally has less power to detect ancestral demographic changes than estimates based on analyses of multiple nuclear markers (Ho and Shapiro 2011). For this reason, a MSC model that utilized information on population structure from previous studies on biparentally inherited nuclear DNA markers (microsatellites) was used to infer divergence times of distinct species and populations of bottlenose dolphins. Incomplete lineage sorting and rapid radiation  Only the sequences generated for this study (N = 30) have been included in the analyses. k = number of segregating sites, i.e., polymorphisms. π nuc /π hap = nucleotide/haplotype diversity. a An outgroup (Steno bredanensis) was used in the calculations of the indices. b Significant difference at α = 0.05 level was detected between the pelagic and both of the coastal populations. c Significant difference at α = 0.05 level was detected between the pelagic and the Coastal North population.
among delphinids (e.g., McGowen et al. 2009;Amaral et al. 2012;Moura et al. 2013;Louis et al. 2014b;Morin et al. 2015;Foote and Morin 2016) can complicate the estimation of the timing of speciation events and population divergence in the genus Tursiops over phylogenetic time scales (Moura et al. 2013). One advantage of using MSC model, such as StarBEAST2 (Ogilvie et al. 2017), compared to traditional methods which assume that evolutionary histories of gene trees match the evolutionary history of the species tree, is that the branch lengths and divergence times are not overestimated (Ogilvie et al. 2017). Another advantage of StarBEAST2 is the possibility to infer historical effective population sizes while accounting for population structuring, as this structuring combined with uneven sampling can lead to bias in N e estimates inferred with traditionally used Bayesian Skyline methods (Heller et al. 2013) that assume panmixia. Additionally, time-dependency of increasing mutation rate towards the present caused by incomplete purifying selection-the time lag between the appearance of a mutation and its removal due to purifying selection-has the potential to further bias divergence time estimates (e.g., Ho et al. 2005;. The use of sites under strong selective constraint extends terminal branches potentially leading to overestimation of older node ages He et al. 2019). In this study, based on the comparison of time-calibrated phylogenies inferred using all codon positions and only the third codons, we found no evidence of incomplete purifying selection affecting the divergence time estimates in the delphinid family or in the phylogeny consisting of Tursiops spp. sequences and haplotypes. However, this may have been an artifact of the small number of common nodes (N = 8) included in the comparison or the relatively young node ages of the branches consisting of T. truncatus samples, and the issue of time-dependency in Tursiops phylogenies should be reexamined with more samples collected from various populations. Nevertheless, third codon positions are thought to be less constrained by purifying selection compared to first and second codon positions, and they are therefore expected to evolve in a more clock-like fashion, accumulating mutations at a more constant rate and thus minimizing the effect of time-dependence ). The average substitution rate for the third codon positions across all branches for the samples from the genus Tursiops was estimated as 0.017 substitutions/site/Myr. This estimate corresponds well with previously estimated rates of 0.024 for third codon positions for cetaceans (Ho and Lanfear 2010) and 0.021 substitutions/site/Myr for bottlenose dolphin mitogenome-wide rates (Moura et al. 2013). Similarly, the estimated divergence date for the second youngest node for the genus Tursiops in this study, represented by the divergence between Mediterranean and Black Sea bottlenose dolphins (around 7 900 years ago), is concordant with that derived previously by Moura et al. (2013), and is consistent with the geological timeframe of the opening of the Bosphorus Strait as well as the climatic context linking to the LGM. In contrast, the average root age for the genus Tursiops in this study was estimated as 3.7 Myr, which is considerably older than 1.1 Myr estimated by Moura et al. (2013). Possible reasons for this incongruence may be a result of the different modeling strategies applied to different parts of the mitogenome (mtDNA protein coding regions were used in this study whilst Moura et al. (2013) used the whole mtDNA sequence). However, a more likely reason is the fact that these authors used a relatively narrow prior for their biogeographical calibration [uniform prior of 0.003-0.01, supplementary table S5 in Moura et al. (2013)] that overpowered the fossil calibrations with much wider priors, and this ultimately resulted in poor convergence of the older node age posterior estimates based on low (<100) ESS values [supplementary table S6 in Moura et al. (2013)].
Pleistocene climatic oscillations are thought to have played a major role in shaping species distribution and divergence and in promoting speciation (Avise et al. 1998). The divergence of the Northwest Atlantic coastal population from other Atlantic T. truncatus populations around 600 000 years ago coincides with the end of a cooler period during the mid-Pleistocene when the SST started to gradually rise again up to ~10 °C in the Northern Hemisphere (Clark et al. 2006). Consistent with previous studies (Natoli et al. 2005;Moura et al. 2013;Louis et al. 2014a), the phylogenetic tree estimated in this study for the genus Tursiops (Figure 3) suggests that the colonization of coastal habitats in the NEA occurred from a pelagic source population. According to our results, this happened via an ancestor that occupied lower latitude Atlantic waters. The AquaMaps model also indicated that bottlenose dolphins may have inhabited parts of the Mediterranean Sea during the LGM, as some areas were highlighted as suitable habitat ( Figure 2B). Even if the northern parts of the Mediterranean Sea (e.g., Tyrrhenian Sea) were too cold or otherwise unsuitable for the bottlenose dolphins during the LGM, a refugial population may have existed in the areas around the warmer Alboran Sea, with models of SST showing consistently warmer annual temperatures above 10 °C for this area throughout the last ice age (Cacho et al. 2001). Similar refugia at lower latitudes have been documented, for example, with Atlantic salmon, Salmo salar (Consuegra et al. 2002;Finnegan et al. 2013) and European eel, Anguilla anguilla (Kettle et al. 2011). Suitable habitat for bottlenose dolphins during the LGM was also hindcasted with AquaMaps in the Black Sea. However, the land bridge between the Black Sea and the Mediterranean collapsed only ~10 000 years ago (Kerey et al. 2004). Moreover, an accumulation of species found in the Black Sea sediment that are indicative of marine, rather than brackish, conditions was dated to 7 000-6 400 ya (the Kalamitian transgression [Yanko-Hombach et al. 2007]). This makes the Black Sea an unlikely habitat for bottlenose dolphins until ~7 000 years ago at the earliest, and the divergence time between the eastern Mediterranean and Black Sea populations estimated as ca. 7 900 ya in this study fits quite well within this timeframe.
Temperature changes after the LGM have likely also played a role in the colonization of the northern parts of the NEA. The post-LGM divergence of the coastal populations in the NEA from the pelagic ~15 000 years ago coincides with the rising temperatures in the lower latitudes of the North Atlantic Ocean. The subsequent recent split of the Coastal North population from its southern counterpart currently occupying lower latitude European waters was estimated to have occurred around 8 000 years ago, and the timing of the split is supported by 2 T. truncatus subfossils found in the Dutch Southern Bight and radio-carbon aged to 7 390 and 8 135 years, respectively (Post 2005). It thus appears that these cladogenetic events are correlated with periods of temperature changes, with warmer temperatures leading to an increase in sea-level in coastal areas and the subsequent release of available habitat. Furthermore, it is also likely that suitable coastal habitat became available gradually in the NEA coastal areas after the LGM and the onset of the warmer Holocene leading to sequential colonization and divergence.
According to Hewitt (1999Hewitt ( 2000, the colonization of suitable niches in high latitudes would occur rapidly via long-range dispersal events. Being highly mobile, cetaceans are capable of undertaking long distance movements on relatively short time scales (Forcada 2009;O'Brien et al. 2009;Robinson et al. 2012). Coupled with this potential for long-range dispersal, the fact that bottlenose dolphins form small local resident populations in coastal European waters (e.g., Ingram and Rogan 2002;Cheney et al. 2013;Louis et al. 2015) makes postglacial leading-edge colonization of the northern NEA a likely dispersal scenario. The demographic contraction and subsequent expansion observed in the coastal populations of this study, as well as the significantly reduced nucleotide and haplotype diversities, suggest that bottlenose dolphins indeed expanded their range northwards to coastal areas via founder events during leading-edge expansion. A similar founder event has recently been suggested to give rise to the present-day resident-type killer whale communities in the Russian Pacific (Filatova et al. 2018), and rapid demographic expansions after the LGM have been documented across other marine taxa (e.g., Klimova et al. 2014;Jenkins et al. 2018).
Our results indicate that the colonization and expansion to fill the available coastal habitat around the British Isles happened via serial founder events after the LGM by a small leading-edge founder group, that comprised as few as 2 maternal lineages, from the source pelagic population. This finding is supported by previous publications of larger datasets of 136 mitochondrial control region sequences from coastal populations around the British Isles, which comprise the same 2 mitochondrial haplogroups (Mirimin et al. 2011;Louis et al. 2014b). The inferred process of postglacial expansion into emerging habitat reflects the model inferred by Hoelzel et al., (2007) for colonization of the coastal regions of the North Pacific and the formation of structure in killer whales. This process may thus reflect a commonality in the demographic history of high latitude social odontocetes. More broadly, serial founder events by a small subset of the parent population, leading to sequential loss of genetic diversity, have been linked to the phylogeographic history in a number of species, including modern humans (Homo sapiens) during the migration out of Africa (DeGiorgio et al. 2009;Deshpande et al. 2009;Henn et al. 2012). Still today the greatest nuclear and haplotype diversity among modern humans is found in Africa (reviewed in Campbell and Tishkoff 2008), and it is generally accepted that Africa was the source of all current modern human populations. Similarly, we identify a loss of genetic diversity in bottlenose dolphins at the northern edge of their coastal range relative to that found in the large pelagic population (e.g., Quérouil et al. 2007;Mirimin et al. 2011;Louis et al. 2014b).
In summary, our results suggest that the present day geographic and genetic structuring of bottlenose dolphins in the North Atlantic were shaped by the climatic cycles of the Late Quaternary, and the most northerly coastal population represents a leading-edge expansion by 2 founding mitochondrial lineages. We estimate these lineages originated from a diverse pelagic source population, and then rapidly spread throughout and retained the available coastal territory and associated resources around the British Isles. However, the responses to climate change are likely to vary even among mobile species depending on their tolerance to environmental conditions and the derived habitat preference (Foote et al. 2013;Sydeman et al. 2015). Although it seems likely that the amount of available suitable habitat for bottlenose dolphins will increase and expand northwards in the next hundred years with the ongoing warming of the oceans (Appendix S13, Supplementary Material), the complexity of ecosystems that bottlenose dolphins live in, combined with a range of social behaviors that they exhibit, makes predicting the consequences of the climate change to populations an extremely difficult task.

Supplementary Material
Supplementary data are available at Journal of Heredity online.

Funding
This work was supported by Thomas Crawford Hayes Scholarship (School of Biological, Earth and Environmental Sciences, University College Cork) to M.N.; National Parks and Wildlife Service, The Department of Culture, Heritage and the Gaeltacht, Ireland, to E.R. and S.I.; Fondation Fyssen to M.L.; Fondation Total; Agence de l'Eau Seine-Normandie; Fonds de Dotation pour la Biodiversité; Agence des Aires Marines Protégées; Direction Régionale de l'Environnement, de l'Aménagement et du Logement; Ministère de l'Ecologie, du Développement Durable et de l'Energie; Conseil Général de la Manche; and Communauté d'Agglomération de la Ville de La Rochelle.