Genetic and historic evidence for climate-driven population fragmentation in a top cetacean predator: the harbour porpoises in European water

Recent climate change has triggered profound reorganization in northeast Atlantic ecosystems, with substantial impact on the distribution of marine assemblages from plankton to fishes. However, assessing the repercussions on apex marine predators remains a challenging issue, especially for pelagic species. In this study, we use Bayesian coalescent modelling of microsatellite variation to track the population demographic history of one of the smallest temperate cetaceans, the harbour porpoise (Phocoena phocoena) in European waters. Combining genetic inferences with palaeo-oceanographic and historical records provides strong evidence that populations of harbour porpoises have responded markedly to the recent climate-driven reorganization in the eastern North Atlantic food web. This response includes the isolation of porpoises in Iberian waters from those further north only approximately 300 years ago with a predominant northward migration, contemporaneous with the warming trend underway since the ‘Little Ice Age’ period and with the ongoing retreat of cold-water fishes from the Bay of Biscay. The extinction or exodus of harbour porpoises from the Mediterranean Sea (leaving an isolated relict population in the Black Sea) has lacked a coherent explanation. The present results suggest that the fragmentation of harbour distribution range in the Mediterranean Sea was triggered during the warm ‘Mid-Holocene Optimum’ period (approx. 5000 years ago), by the end of the post-glacial nutrient-rich ‘Sapropel’ conditions that prevailed before that time.


INTRODUCTION
Changes in the environment can affect the behaviour and fitness of individual organisms and the viability of populations. Populations can respond to changing habitats by adapting (through natural selection or phenotypic plasticity), moving (to avoid habitat of reduced suitability, or take advantage of emerging habitat), by adjusting population size or some combination of the above. Both natural selection and genetic drift can shape populations as they evolve in this context.
In the North Atlantic, studies have increasingly reported strong shifts in plankton and fish assemblages with the contemporaneous climate and ocean warming (Beaugrand et al. 2002(Beaugrand et al. , 2008Richardson & Schoeman 2004;Perry et al. 2005). These shifts in marine habitat and community structure are expected to drive major changes in the distribution, density and dispersal of apex predators such as marine mammals, with reduction in density in marginal parts of the habitat leading to fragmentation, local displacement or extinction if the habitat becomes unsuitable (Ryther 1969;Learmonth et al. 2006;O'Corry-Crowe 2008;de Bruyn et al. 2009;Marx & Uhen 2010). Such evidence is, however, hard to capture for highly mobile marine predators (Kintisch 2006) such as cetacean species (Learmonth et al. 2006;O'Corry-Crowe 2008). Key to answering how these apex predators will deal with current and future change in their environment lies in (i) understanding how cetaceans dealt with past climate changes and (ii) determining the impact of changing climate on cetaceans over contemporary time scales (O'Corry-Crowe 2008).
Past shifts in the Earth's climate have led to the alteration or loss of existing environments and the creation of new ones. These changes have been associated with speciation events, adaptive radiations and extinctions, population expansions and contractions, and changes in individual dispersal and breeding behaviour (Hewitt 2000(Hewitt , 2004Hofreiter & Stewart 2009). Recent studies supported such processes in cetacean species by providing evidence that past changes in the environment, and in particular global change in temperature and diatoms richness that are dominant marine primary producers at the basis of marine food web, were significantly linked to the species diversity of mysticete and odontocete whales (Pastene et al. 2007;O'Corry-Crowe 2008;Marx & Uhen 2010). The contemporaneous climate change is however unique with respect to its speed and to the fact that it embeds in global change where both natural and anthropogenic influences act in synergy on ecosystems (Brook et al. 2008;Pimm 2008Pimm , 2009). Understanding how these factors affect marine systems is therefore a crucial issue.
In this study, we investigated the demographic history of a small coastal cetacean widely distributed in the North Atlantic, the harbour porpoise Phocoena phocoena, with regards to recent variation in its habitat. At around 1.5 m in length and 50 kg in weight, it is the smallest cetacean in the North Atlantic (Read 1999), and often considered as a species living in the 'fast lane' (Read & Hohn 1995;Lockyer 2007). Reproductive costs of harbour porpoises are high (Lockyer 2007), as females are often gestating and lactating at the same time and parturition occurs shortly before mating (Lockyer 2003). Given their limited capacity to store energy, it is assumed that harbour porpoises must feed frequently without prolonged periods of fasting (Koopman et al. 1996(Koopman et al. , 2002. Some authors suggested that fasting periods exceeding as short as 3 days could affect body condition (Kastelein et al. 1997;Lockyer 2007). Thus, temporary shortages in prey availability can negatively impact on the animals and are likely to be responsible for changes in their distribution (Santos & Pierce 2003;Santos et al. 2004;Johnston et al. 2005;Herr et al. 2009). Relatively continuous accessibility to adequate prey is therefore critical. Any changes in prey availability may affect energy stores, and ultimately survival (MacLeod et al. 2007). Therefore, we expect the distribution of this species to be strongly tied to variation in the primary and secondary productivity that provides the basis for apex consumers (Lockyer 2007).
The harbour porpoise is currently distributed fairly continuously throughout cold to temperate coastal waters of the North Pacific, the North Atlantic and in the Black Sea, but it is absent from the Mediterranean (Read 1999). The southeastern part of the distribution range in the North Atlantic is of particular interest in regard to recent climate change, as it covers the biogeographic transition between arctic/boreal species and subtropical species (Southward et al. 1995). Changes in distribution patterns and abundance in this region are thus more evident than elsewhere in Europe for many marine species (Southward et al. 1995;Beaugrand et al. 2008). In this southeastern part of the North Atlantic, the harbour porpoise displays a disjunct distribution range. The absence of porpoises in the Mediterranean splits the Black Sea range from the Atlantic range. Black Sea porpoises are therefore recognized as a subspecies P. p. relicta distinct from the Atlantic one P. p. phocoena (Rosel et al. 1995), and from an evolutionary perspective, as an evolutionary significant unit (ESU; Moritz 1994Moritz , 2002 independent from the Atlantic one (Tolley & Rosel 2006;Viaud-Martinez et al. 2007). More recently, in a study covering the entire eastern Atlantic range, Fontaine et al. (2007) discovered a second disjunction, a population restricted to the cold-water upwelling zone along the Atlantic coasts of Iberia (IB) (Fiù za 1983) and distinct from porpoises in the northern Bay of Biscay (NBB), which are part of the main North Atlantic range, extending from French waters fairly continuously to the Arctic (figure 1). Fontaine et al. (2007) argued that oceanographic conditions, specifically those conditioning food availability, impose strong constraints on the dispersal, and thus the genetic structure of the species. If this hypothesis is correct, the existing genetic structure across the North Atlantic range will have a recent origin, as the oceanographic conditions in the eastern North Atlantic, in the Mediterranean Sea and in the Black Sea have changed dramatically over the course of the Holocene (Alley et al. 2003;Mayewski et al. 2004;Osborn & Briffa 2006;Springer et al. 2008;Rohling et al. 2009).
As a cold temperate species, porpoises are unlikely to favour the current conditions of the southern Bay of Biscay and the Mediterranean Sea, with deep, warm and oligotrophic (nutrient-poor) waters (Fontaine et al. 2007). However, in the recent past, these waters have at times been cold and nutrient-rich, and it seems likely that porpoise ranges have closely tracked this nutrient availability. This would suggest the Iberian porpoise range could have been continuous with that of North Atlantic waters during recent cold periods that affected European climate (Grove 2004). There have been a number of notable cold periods between the Last Glacial Maximum (LGM, 18 000 years before present (yr BP)) and the most recent cold period, known as the 'Little Ice Age' (LIA). The period of continuity between the Black Sea and Atlantic ranges is a much more long standing question (Frantzis et al. 2001). The straits of Bosphorus and Dardanelle have been periodically closed during the Quaternary, but have been open continuously since approximately 8400 years ago (Major et al. 2002;Rohling et al. 2009). Analyses of genetic polymorphism at the non-coding region of the mitochondrial genome (the control region, mtDNA CR) (Tolley & Rosel 2006;Viaud-Martinez et al. 2007) suggested, however, a much older splitting time between Atlantic and Black Sea mitochondrial clades (at least 175 000 years ago).
In order to investigate the splitting process among the three porpoise populations identified in the Northeast Atlantic, we extracted historical information from genetic polymorphism data using model-based Bayesian coalescent approaches, and correlated this genetic inference with palaeo-oceanographic and historical records. Specifically, we tested whether the genetic divergence observed among harbour porpoise populations in the Northeast Atlantic is compatible with the recent changes reported in the seascape and marine assemblages, or is more likely owing to ancient vicariance processes.

MATERIAL AND METHODS
(a) Microsatellite genotypes Individual genotypes at 10 microsatellite loci were from Fontaine et al. (2007). To avoid computation time becoming intractable for the coalescent approach described below, we restricted the analyses to 183 harbour porpoises at the borders of the genetic barriers evidenced in Fontaine et al. (2007): 74 animals from the northern Bay of Biscay, 31 collected along the Iberian coast and 78 from the Black Sea (figure 1). Details regarding descriptive statistics on the microsatellite dataset can be found in Fontaine et al. (2007).

(b) Data analysis
We analysed the datasets using the 'isolation with migration' (IM) model of population divergence (Hey & Nielsen 2004). The model uses coalescent simulations within a Bayesian inference framework to estimate marginal probability distributions for six demographic parameters scaled by the mutation rate m: the time since population-splitting (T ¼ t m), measures of neutral population genetic diversity of the two current and one ancestral population (u 1 , u 2 , u A ) proportional to their effective population sizes N e (u ¼ 4N e m) and bidirectional migration rates ( (Hey & Nielsen 2004).
We analysed population-splitting first comparing the IB with the NBB population. Then, we compared the two ESUs, i.e. harbour porpoises from the Black Sea (BS) with those in the Atlantic (AT, including IB and NBB). This procedure was adopted because previous works suggested that the Black Sea harbour porpoises separated from Atlantic populations long before IB porpoises separated from the NBB populations. This is supported by complete lineage sorting of mtDNA haplotypes between porpoises from the Black Sea and the North Atlantic (Tolley & Rosel 2006;Viaud-Martinez et al. 2007), in contrast with the lack of private mitochondrial and microsatellite alleles between the Iberian and northern Bay of Biscay (  Demo-genetic history of porpoises M. C. Fontaine et al. 2831 Atlantic populations may induce some bias in our inference regarding the Black Sea-Atlantic comparison, because the model assumes panmixia within each population. However, a recent simulation study showed that the model estimates were robust to such departure (Strasburg & Rieseberg 2010). On the other hand, comparing each pair of populations separately (IB -BS and NBB-BS), thus ignoring a population in the Atlantic (IB or NATL) that is more closely related to the other, profoundly violates the basic assumptions of the model (Hey & Nielsen 2004) and could have a more profound impact on the model estimates (Strasburg & Rieseberg 2010).
We conducted Markov Chain Monte Carlo (MCMC) simulations using the IM program (Hey & Nielsen 2004), assuming a stepwise mutation model of microsatellite evolution and uniform prior distributions over parameter ranges. Prior boundaries were empirically determined to ensure that the posterior distributions fell completely within the prior range. The modes of the posterior probability distributions were thus taken as the maximum-likelihood estimate (MLE) of the parameters. We estimated the credibility intervals as the 90 per cent highest probability density intervals (HPDIs) (i.e. the shortest span that includes 90% of the probability density of a parameter).
Convergence of such a kind of analysis and dataset is a tricky issue. To ensure reliable convergence towards the stationary distribution, we monitored multiple independent runs for the two datasets, each with 50-100 independent heated chains under Metropolis coupling to improve mixing (Hey & Nielsen 2004). Each run was initiated with a burn-in period of at least 300 000 updates, with each multi-chain simulation lasting several months (the longest runs were up to 3.1 Â 10 7 iterations, see the electronic supplementary material, table S1, for details regarding runs). Mixing properties of the MCMC were assessed by visual inspection of the parameter trend plots and by examining the level of autocorrelation between initial and final parameter values (using the effective sample size, ESS). Analyses were considered to have converged upon the stationary distribution if independent runs generated similar posterior distributions, with each having a lowest ESS of 50 for each estimated parameter as recommended in Hey & Nielsen (2004).
To convert posterior estimates of the model scaled by the mutation rate m (i.e. T and u) into demographic units (i.e. t and N e ), we used a value of m ¼ 5.0 Â 10 24 mutation per generation, which is considered as the average mutation rate over many species (Brinkmann et al. 1998;Estoup & Angers 1998;Ellegren 2000;Estoup et al. 2002;Sun et al. 2009). Assuming a generation time (G) for harbour porpoise ranging between 5 (G L , the average age of females at first parturition) (Lockyer 2003) and 7 years (G H , as suggested by Tolley & Rosel 2006), population-splitting time (T ) can be converted to calendar years (t) and estimates of population mutation rates (u 1 , u 2 and u A ) can be converted into estimates of effective population size parameters (N 1 , N 2 and N A , respectively, in a number of individuals) as described in Hey & Nielsen (2004). Migration parameters in the model (M 1 and M 2 ) can be used to obtain population migration rates (2N 1 m 1 and 2N 2 m 2 , i.e. the effective number of migrants per generation into population 1 and 2, respectively) using the MLE estimates of current u and of M (2N 1 m 1 ¼ u 1 M 1 /2 and 2N 2 m 2 ¼ u 2 M 2 /2) (Hey & Nielsen 2004).

RESULTS AND DISCUSSION
(a) Convergence of IM runs After several steps of trimming, the replicated IM simulations displayed good mixing properties of the coupled MCMCs as suggested by the parameter trend line plots (data not shown) with an ESS of at least 71 for the IB -NBB analyses and 53 for AT-BS analyses (details on these runs are provided in the electronic supplementary material, table S1). IM replicates provided consistent similar posterior probability densities for all parameters (electronic supplementary material, figures S1 and S2). All these indications suggested that we have sufficiently explored the parameter space. We selected runs that best converged (i.e. those displaying the highest ESS estimates) to convert them into demographic units (figure 2 and table 1) and to interpret the outputs.  (table 1).
The genetic inference thus indicates the barrier to gene flow located in the southern Bay of Biscay (Fontaine et al. 2007) rose during the past millennium and most probably during the few last centuries. It is now widely recognized that worldwide anomalies in climatic and oceanic conditions occurred during the past 2000 years (Osborn & Briffa 2006). The most significant according to Osborn & Briffa (2006) include the 'Medieval Warm Period' (1200 -800 yr BP), the 'LIA' (600 -150 yr BP) and the modern warm period. During the LIA that culminated 300 years ago, dramatic cooling in oceanographic and climatic conditions in the North Atlantic, and over Europe in particular, had a substantial impact on human activities, for which there are historical records, but also more generally on terrestrial and marine ecosystems (Grove 2004). Fishery records in Europe, which date back as far as the tenth century for some species (Alheit & Hagen 1997), report both significant increases in capture of cold-water fishes (e.g. herring Clupea harengus) as far south as in the Bay of Biscay and the southward retreat of warm-water species (e.g. sardine) during the LIA (Alheit & Hagen 1997). Given that cold-water fishes, like herring or sandeel (Ammoditidae sp.), are a major component of porpoise diet throughout its current range (Read 1999;Santos & Pierce 2003), it is highly probable that the harbour porpoise's past range was likewise extended southward into the Bay of Biscay. The IM demographic estimates in terms of timing, geneflow and the relative size of the subranges are consistent with such a scenario where harbour porpoises north and south of the Bay of Biscay formed a single range during the LIA cold conditions, and began to fragment with the onset of the current warm period 200 -300 yr BP (Osborn & Briffa 2006), together with substantial shifts in the marine community in the eastern North Atlantic and especially in the Bay of Biscay (Poulard & Blanchard 2005). The relatively cold nutrient-enriched conditions generated by coastal upwelling along the Iberian Atlantic coasts (Fiù za 1983) have allowed a small porpoise population and other cold temperate species to persist in these waters up to now (OSPAR Commission 2000). The predominant northward migration signal detected in the genetic data supports a northward exodus of Iberian porpoises during their isolation, as expected for a Northern Hemisphere temperate species in response to ocean warming.  (table 1 and figure 2f ). Previous studies suggested a much older separation time of Black Sea from Atlantic porpoises based on the analysis of genetic polymorphism at the non-coding region of the mtDNA CR (Tolley & Rosel 2006;Viaud-Martinez et al. 2007). We ran further IM simulations using published and new mtDNA sequences (see electronic supplementary material, methods, table S1, S2 and figure S3). Combined microsatellite and mtDNA data did not change inference compared with the results described above, but results from the mtDNA locus alone showed a strong discrepancy with the nuclear loci (electronic supplementary material, figure  S3; Viaud-Martinez et al. 2007). The Atlantic-Black Sea splitting time estimate based only on the mtDNA CR and converted into calendar years using the formerly admitted range of mutation rates (between 3.3 and 4.3 Â 10 28 substitutions per site and per year (s.s.yr 21 ) (Tolley & Rosel 2006;Viaud-Martinez et al. 2007) was up to three orders of magnitude greater than at nuclear loci. However, any estimate of divergence time is conditioned on the uncertainty regarding the mutation rate estimates, which is very difficult to assess and highly variable among species (Nabholz et al. 2008a). The IM algorithm accounts for variation in the mutation rate among loci by estimating the relative locus-specific mutation rate as a mutation rate scalar (Hey & Nielsen 2004). Provided we have a sound mutation rate estimate for at least one locus allows us to deduce the mutation rate of the others. Assuming that the average mutation rate for microsatellite loci is 1 Â 10 24 mutations per year (considering a generation time of 5 years), the mtDNA CR substitution rate can be estimated from our combined analysis at 1.3 Â 10 26 s.s.yr 21 (HPDI: 8.4 Â 10 27 -2.3 Â 10 26 ). This is at least two orders of magnitude faster than 'traditional' substitution rate estimates derived from interspecific phylogenetic datasets or the fossil record (Hoelzel et al. 1991;Henn et al. 2009) but, within the intervals of those found with the mean human mtDNA CR pedigree rate estimate derived from a meta-analysis (9.5 Â 10 27 s.s.yr 21 ) (Howell et al. 2003), with ancient DNA datasets in Adélie penguin Pygoscelis adeliae (9.6 Â 10 27 s.s.yr 21 ) (Lambert et al. 2002), in Southern elephant seal Mirounga leonina 9.8 Â 10 27 s.s.yr 21 HPDI: 1.7 Â 10 29 -2.1 Â 10 26 ) (de Bruyn et al. 2009) and in Steller sea lions Eumetopias jubatus 2.7 10 27 s.s.yr 21 (Phillips et al. 2009). Such a high intraspecific estimate is reminiscent to the time dependency of mtDNA mutation rates recently recognized as a general phenomenon applying to many species (e.g. primates, marine mammals, birds and fishes) (Lambert et al. 2002;Howell et al. 2003;Ho et al. 2005Ho et al. , 2007Bandelt 2008;de Bruyn et al. 2009;Endicott et al. 2009;Henn et al. 2009). Higher mutation rate estimates are obtained from events calibrated with more recent dates, such as pedigrees, and lower estimates result when events are calibrated with older dates.
The anomalous behaviour of the mitochondrion did not only concern the substitution rate, but also estimates of ancestral population genetic diversity (u) from mtDNA CR alone. mtDNA u estimates showed a diametrically opposed pattern to that of nuclear microsatellite loci, with estimates of ancestral population genetic diversity being much smaller than the current population diversity (electronic supplementary material, figure S3 and fig. 4 in Viaud-Martinez et al. (2007)). Under neutral evolution of the genetic markers used, historical variation in population size should affect genetic polymorphism at both nuclear and mitochondrial loci in a similar manner (albeit with larger drift effects on the mitochondrion). The opposing pattern we observed here could be the outcome of the large stochastic variance in the coalescence process affecting a single locus (Hey & Machado 2003), but it is worthwhile considering other potential causal factors. First, selective processes may have acted on the mtDNA genome (e.g. selective sweeps). Given the central role of mitochondria in heat and energy production, and the strong energetic constraints faced by porpoises which are small, warm-blooded marine carnivores with a demanding reproductive schedule, detecting a signal of selection on the mitochondrion would seem hardly surprising. Second, heterogeneous mutation rates along the mtDNA CR are poorly captured by the substitution model. Mutational hotspots resulting in high levels of homoplasy could provide a false signal of population expansion (Pakendorf & Stoneking 2005;Galtier et al. 2006;Strasburg & Rieseberg 2010). Finally, the discrepancy may be the result of complex and multiple determinants of genetic diversity on the mtDNA genome. A recent study examining mtDNA evolution in mammals showed that mitochondrial rates of evolution generally do not depart from neutral expectations, although they are highly stochastic among species, and essentially unpredictable from knowledge of species biology (Nabholz et al. 2008b;Galtier et al. 2009).
As harbour porpoises currently inhabit productive coastal waters of the North Atlantic and Black Sea (Read 1999), we expected that porpoises were present in the Mediterranean at a time when palaeoceanographic conditions were more productive than the present-day oligotrophic conditions. Geological records indicate that such conditions prevailed in the Mediterranean Sea throughout the deglaciation period following the LGM (approx. 18 000 yr BP) until the Mid-Holocene warm period (approx. 6000 yr BP) (Kerwin et al. 1999;Jimenez-Espejo et al. 2007;Rohling et al. 2009). The strongest contrast with current conditions was in the eastern Mediterranean basin, where organic-rich sediments, termed 'Sapropel', accumulated under strong hydrological changes generated by intense fresh water discharges and substantially enhanced primary productivity (Calvert et al. 1992;Rohling et al. 2009). At the end of this Sapropel episode, approximately 5500 years ago, primary productivity abruptly decreased toward present-day conditions. Demographic estimates from the IM model based on microsatellite data or combined microsatellite and mtDNA data are consistent with porpoise range fragmentation in the Mediterranean starting around 5500 years ago at the end of this Sapropel event, while productive conditions persisted in the northern Aegean Sea and the Black Sea. The posterior estimates of Atlantic -Black Sea time of separation based on microsatellite data firmly place this isolation event in the post-LGM period and point to the post-Sapropel period. The signal of migration from the Atlantic to the Black Sea captured in the IM analysis (figure 2f and table 1) probably reflects migration events during the early stage of the divergence process and is consistent with the proposal that the last remnants of Mediterranean porpoises (the 'admixed' link between Atlantic and Black Sea) ended up retreating to the Black Sea.
Post-glacial and very recent climate-driven environmental changes have triggered deep reorganization in Northeast Atlantic ecosystems, with changes detected in plankton and fish communities. Our study provides strong evidences of such a response for a small coastal cetacean that has to meet a particularly demanding energetic budget to survive. Harbour porpoise populations in the Northeast Atlantic waters have undergone recent dramatic changes in connectivity and distribution probably related to climatedriven variation in habitat suitability. Our inferences naturally depend upon a specific model of population evolution (IM model), which necessarily simplifies the population demographic history, and on a specific mutation rate. Departures from the model, such as departure from panmixia within population or change in the effective size of diverging populations, may induce some bias in our estimates. However, the model appeared quite robust to these departures (Strasburg & Rieseberg 2010). The biological congruence between our genetic inferences and historical and palaeoceanographic records suggests the model captures the essential aspects of this fragmentation process.
It seems inevitable that marine predators will need to adapt to a changing spatial distribution of primary and secondary production within pelagic marine ecosystems as they did in the past (Marx & Uhen 2010). This reaction to climate change could however be exacerbated on a much shorter time scale by the depletion of major commercial fish stocks induced by intense commercial fisheries (Worm et al. 2006(Worm et al. , 2009 and/or by detrimental interactions with human activities such as incidental catches, which probably constitute a much more serious short-term threat for small coastal cetaceans such as harbour porpoises (Stenson 2003;Hovelsrud et al. 2008). The combination of these effects raises serious concerns regarding the ability of top predators to adapt to such rapidly changing habitat conditions. Addressing these concerns will be a particular challenge for future research and management policy.
We are grateful to all fishermen, stranding networks and volunteers that contributed in the collection of samples used in the present study. We thank A. Jacomme, A. Estoup, C. Moritz, R. Cheddadi, E. J. Rohling, L. Crespin, R. Streiff, T. Malausa, S. Piry, S. J. Goodman, G. Roderick, J. Hey and N. Galtier for their valuable help and/or comments during the development of this research. This work was funded by the Belgian Science Policy (Contract SSTC EV/12/46A). M.C.F. was supported by research fellowships from the Belgian National Fund for Scientific Research (FRS-FNRS, mandate 'Aspirant'), and by the AGAPE Marie-Curie Fellowship programme, and K.A.T. by the Norwegian Research Council. This paper is an MARE publication no. 187.