Postglacial species displacement in Triturus newts deduced from asymmetrically introgressed mitochondrial DNA and ecological niche models

Background If the geographical displacement of one species by another is accompanied by hybridization, mitochondrial DNA can introgress asymmetrically, from the outcompeted species into the invading species, over a large area. We explore this phenomenon using the two parapatric crested newt species, Triturus macedonicus and T. karelinii, distributed on the Balkan Peninsula in south-eastern Europe, as a model. Results We first delimit a ca. 54,000 km2 area in which T. macedonicus contains T. karelinii mitochondrial DNA. This introgression zone bisects the range of T. karelinii, cutting off a T. karelinii enclave. The high similarity of introgressed mitochondrial DNA haplotypes with those found in T. karelinii suggests a recent transfer across the species boundary. We then use ecological niche modeling to explore habitat suitability of the location of the present day introgression zone under current, mid-Holocene and Last Glacial Maximum conditions. This area was inhospitable during the Last Glacial Maximum for both species, but would have been habitable at the mid-Holocene. Since the mid-Holocene, habitat suitability generally increased for T. macedonicus, whereas it decreased for T. karelinii. Conclusion The presence of a T. karelinii enclave suggests that T. karelinii was the first to colonize the area where the present day introgression zone is positioned after the Last Glacial Maximum. Subsequently, we propose T. karelinii was outcompeted by T. macedonicus, which captured T. karelinii mitochondrial DNA via introgressive hybridization in the process. Ecological niche modeling suggests that this replacement was likely facilitated by a shift in climate since the mid-Holocene. We suggest that the northwestern part of the current introgression zone was probably never inhabited by T. karelinii itself, and that T. karelinii mitochondrial DNA spread there through T. macedonicus exclusively. Considering the spatial distribution of the introgressed mitochondrial DNA and the signal derived from ecological niche modeling, we do not favor the hypothesis that foreign mitochondrial DNA was pulled into the T. macedonicus range by natural selection.


Background
Speciation involves the evolution of reproductive isolating mechanisms, preventing diverging gene pools from merging [1]. If, however, interspecific hybridization still occurs, the stage is set for the transfer of segments of DNA across the species boundary [2][3][4][5]. In the majority of eukaryotes, mitochondrial DNA is passed clonally through the female line only and therefore either gets transmitted to the next generation in its entirety or not at all. Mitochondrial DNA does, hence, not 'dilute' via recombination over the generations, as is the case with nuclear DNA. [6,7]. Through an initial hybridization event and subsequent backcrossing of female progeny with the paternal species, mitochondrial DNA can spread on the foreign species background (Figure 1a).
Parapatric species can show geographically asymmetric introgression, with one species showing genes typical of the other over a significant part of its range [reviewed in ref. 4]. This could be caused by natural selection, pulling favorable foreign genes into a species' range [8]. However, empirical and simulation studies show that if displacement of one species by another coincides with hybridization, introgression of neutral genes from the native to the invading species is the expected outcome [4,5]. Furthermore, markers showing reduced intraspecific gene flow (such as mitochondrial DNA) are particularly susceptible to introgression, as their dispersal limitation reduces the chances of being genetically swamped [4,5].
We explain the principle of introgression due to species displacement using mitochondrial DNA as an example (Figure 1b). The rationale is as follows: mitochondrial DNA of a common, native 'donor' species is captured by an initially rare invader, through introgressive hybridization. As the invader expands its range at the expense of the native species, its population increases and the introgressed mitochondrial DNA is co-amplified. Consequently, the introgressed mitochondrial DNA surfs the wave of advance and is left in its wake. Initial introgression does not have to occur at a high frequency to result in a mismatch between species identity and mitochondrial DNA type spanning an extensive area.
The climate oscillations related to the glacial/interglacial cycles of the Quaternary ice age repeatedly forced species to shift their distributions [9]. Distribution rearrangement of closely related species after postglacial secondary contact is considered a likely driver of introgression [3,4]. To explore the role of postglacial displacement in introgression, ecological niche modeling -the approximation of the ecological requirements of a species based on the range of environmental conditions experienced at known localities [10] -can be used. By extrapolating ecological niche models on both current and past climate layers, distribution shifts can be identified [11]. By determining which part of the current distribution was uninhabitable at the Last Glacial Maximum (~21 Ka), the area colonized post-glacially can be approximated. Introgression of neutral genes after species displacement would likely be restricted to this particular region, whereas such a geographical limitation would not apply to introgression by positive selection.

Model system: Triturus newts
We use the two parapatrically distributed crested newt species Triturus macedonicus and T. karelinii as a model system. The contact zone of these newts is situated on the Balkan Peninsula in southern Europe (Figure 2a). Note that T. karelinii comprises three distinct mitochondrial DNA lineages, which might better be regarded as different species [12,13]. When we refer to T. karelinii, in this paper, we mean the lineage that occurs in the Balkans and western Turkey, in parapatry with T. macedonicus. Their parapatric distribution suggests the two crested newt species, at least locally, are able to competitively exclude one another. This is likely to be driven by different abiotic preferences: T. macedonicus is presumed to be relatively more aquatic then T. karelinii [13,14].
Triturus macedonicus and T. karelinii are not sister groups: their most recent common ancestor is the ancestor of all crested newt species (the T. cristatus superspecies), estimated to have lived c. 10-11 million years ago [13,15,16]. Analysis of a suite of allozyme markers has shown the nuclear gene pools of T. macedonicus and T. karelinii to have remained isolated [17,18, Arntzen et al.  Figure 1a depicts the transfer of mitochondrial DNA across the species boundary via introgressive hybridization. Large circles reflect the nuclear DNA composition of individuals, small ones their mitochondrial DNA type. There is an initial hybridization event between the members of two species, a red female and a green male. The F1 offspring contain a mix of red and green nuclear DNA, but only red mitochondrial DNA (due to mitochondrial DNA's matrilineal transmission). Subsequent backcrossing of admixed females with green males over the generations dilutes the red nuclear DNA out, in effect reconstituting the green species' nuclear genome, but with red mitochondrial DNA (ensured by mitochondrial DNA's clonal transmission). Figure 1b shows how the outcompeting of one species by another can result in asymmetric mitochondrial DNA introgression. Small circles now reflect the spatial distribution of mitochondrial DNA type and the background the geographical nuclear DNA composition of the population. At the top, the ranges of a red and a green species are in allopatry, but green expands its range towards red. In the middle, green and red have become parapatric. The species hybridize at the contact zone, where red mitochondrial DNA introgresses into the green species (as in Figure 1a). At the bottom, green shifts its distribution further to the right, at the expense of red. As the members of the green species leading the expansion contain red instead of green mitochondrial DNA, only red mitochondrial DNA spreads in the region where the green species displaces the red species. This figure is based on [3,5]. submitted]. The mitochondrial DNA of the two species is highly distinct [13]. Furthermore, T. macedonicus and T. karelinii can be distinguished by morphology, including meristic, morphometric and qualitative characters [14,19].   Figure 2a shows the distribution of T. macedonicus (shaded green) and T. karelinii (blue). Note that the northwestern part of the range of T. karelinii is disconnected from the main range (i.e. it constitutes an enclave). Ranges are limited by sea (white), uninhabited land (grey) and other Triturus taxa (purple). Green dots represent T. macedonicus localities containing original T. macedonicus mitochondrial DNA, blue dots T. karelinii localities containing T. karelinii mitochondrial DNA, red dots T. macedonicus localities containing T. karelinii mitochondrial DNA and two orange stars T. macedonicus localities containing both original and introgressed mitochondrial DNA. Large dots represent localities for which mitochondrial DNA was sequenced and small dots additional localities included in ecological niche modeling. Localities where genetic admixture of the nuclear genome [based on allozyme data presented in Arntzen et al. submitted] is absent are denoted with a thick black border and those where at least one individual is present showing signs of nuclear genetic admixture are denoted with a thick grey border. Localities for which the composition of the nuclear genome could not be consulted were identified based on diagnostic morphological characters. Figure 2b shows the geographical distribution of mitochondrial DNA based on Thiessen polygons, where each polygon covers the area that is closer to its corresponding locality than to another one. Color codes are the same as in Figure 2a. The red and orange polygons were combined to delimit the zone of asymmetric mitochondrial DNA introgression. enclave of T. karelinii, cut off from the main distribution range [19]. These factors would facilitate asymmetric introgression. Indeed, mitochondrial DNA typical of T. karelinii has been found in T. macedonicus, at localities removed from the current contact zone [20,21]. Here we further explore the introgression of T. karelinii mitochondrial DNA into T. macedonicus. We first delimit the zone of T. karelinii mitochondrial DNA asymmetrically introgressed into T. macedonicus through a dense phylogeographic survey. We then build ecological niche models for both species and project these on current, mid-Holocene (~6 Ka) and Last Glacial Maximum (~21 Ka) climate layers. Inferred distribution shifts across the current zone of asymmetric mitochondrial DNA introgression suggest post-glacial displacement of T. karelinii by T. macedonicus.

Results
As T. macedonicus and T. karelinii mitochondrial DNA differs considerably (D XY = 0.091), sequenced individuals could be unambiguously assigned to mitochondrial DNA type (details in Additional file 1). Out of 71 of the T. macedonicus localities for which sequence data are available, 35 have only original T. macedonicus mitochondrial DNA, 34 have only T. karelinii mitochondrial DNA and two localities show syntopy of both mitochondrial DNA types ( Figure 2a). There is no T. macedonicus mitochondrial DNA found in any of the 86 T. karelinii localities for which sequence data is available. Based on the known geographical distribution of mitochondrial DNA we interpret mitochondrial DNA type for the remaining localities for which sequence data are not available ( Figure 2a). By employing Thiessen polygons, we delimit an extensive (ca. 54,000 km 2 ) area of asymmetric mitochondrial DNA introgression, in which T. macedonicus contains T. karelinii mitochondrial DNA ( Figure 2b).
The T. karelinii mitochondrial DNA shows extensive structuring and we discuss three (not necessarily monophyletic) spatio-temporal groups ( Figure 3). A 'Balkan ancestral' clade, restricted to the extreme south-east Balkan Peninsula, is distinct from the remaining T. karelinii haplotypes (D XY = 0.035). A genetically diverse group of haplotypes is distributed in western ' Asiatic Turkey' (Figure 3). A 'Balkan derived' clade is nested within, and closely related to, the ' Asiatic Turkey' haplotypes (separated by a single substitution). The 'Balkan derived' haplotypes show a starburst pattern: they are genetically similar, comprising a few common and a large number of rare haplotypes ( Figure 3). For the T. karelinii mitochondrial DNA, only the 'Balkan derived' clade shows signs of demographic expansion ( Figure 4). The original T. macedonicus mitochondrial DNA is genetically diverse (Additional file 2) and does not indicate demographic growth ( Figure 4). All introgressed T. karelinii mitochondrial DNA belongs to the 'Balkan derived' clade ( Figure 3). Three introgressed haplotypes are also found in T. karelinii (two of which represent the most frequent haplotypes found in both species) and 12 have only been identified in T. macedonicus.
The ecological niche models of both crested newt species perform statistically significantly better than expected by random chance (Additional file 4). The ecological niche models are most affected by the mean temperature of the coldest quarter in both species and precipitation in the wettest quarter plays an important additional role in T. macedonicus (Additional file 5). When projected on Last Glacial Maximum (~21 Ka) data layers, the models suggest that the ranges of both crested newt species were restricted at the time; the present day introgression zone was uninhabitable ( Figure 5). At the time of the mid-Holocene (~6 Ka) the two species would have been able to occupy much of their current range, but the predicted suitability maps differ from the current situation ( Figure 5). The suitability of the area where the introgression zone is situated generally increased for T. macedonicus since the mid-Holocene, whereas it decreased for T. karelinii. The northwestern part of the current introgression zone is predicted to have been continuously unsuitable for T. karelinii ( Figure 5).

A scenario of species displacement
We confirm that the range of T. karelinii is bisected by T. macedonicus (Figure 2a), corroborating a previously hypothesized T. karelinii enclave [19]. Considering that crested newts are surface-bound and have limited dispersal capabilities, the presence of the T. karelinii enclave is best explained by it having been cut off from the main T. karelinii range by expansion of T. macedonicus. However, the mitochondrial DNA data suggest that a considerably larger area is involved in the geographical displacement of T. karelinii by T. macedonicus than that required to explain the enclave per se (Figure 2b). How was this introgression zone established? The structuring of the mitochondrial DNA of the two crested newt species and predicted distributions based on ecological niche modeling provide useful insight.
The ranges of both crested newt species were reduced at the Last Glacial Maximum ( Figure 5). For T. karelinii, the southwestern margin of its current Balkan range (extending into land currently occupied by T. macedonicus) is predicted to have been hospitable. However, the mitochondrial DNA data suggest that the 'Balkan ancestral' clade was likely restricted to the south-east, whereas the 'Balkan derived' clade originates from a recent colonization from Asiatic Turkey ( Figure 3). Indeed, Turkey's west coast was also predicted suitable at the Last Glacial Maximum and mitochondrial DNA indicates a stable demographical history here ( Figure 5). Furthermore, the 'Balkan derived' clade shows the signature of demographic growth, in line with a range expansion ( Figure 4). The west side of the current range of T. macedonicus was also predicted to have been suitable at the Last Glacial Maximum, corresponding to long-lasting presence in situ as suggested by the mitochondrial DNA data (Figures 4 and 5).
The area where the introgression zone is situated was uninhabitable for both crested newt species at the Last Glacial Maximum and secondary contact between the two was established after glacial conditions had where T. macedonicus contains T. karelinii mitochondrial DNA is outlined with a white stippled line. For each group the nucleotide diversity (π) is determined. In the phylogenetic tree, red arrows signify haplotypes that are (also) found as introgressed mitochondrial DNA in T. macedonicus. Statistically significantly supported clades (with a Bayesian posterior probability of 0.95 or higher) are denoted with an asterisk (*). In the haplotype network, pie diameter expresses haplotype frequency (Additional file 3) and bars the number of mutations along a branch if more than one. Blue pies reflect T. karelinii mitochondrial DNA found in T. karelinii and red pies T. karelinii mitochondrial DNA present in T. macedonicus. If T. karelinii haplotypes are shared between T. karelinii and T. macedonicus, frequencies are reflected by the size of the slices. The numbers in the phylogeny and haplotype network refer to T. karelinii mitochondrial DNA haplotypes as coded in Additional file 3. alleviated ( Figure 5). Furthermore, the T. karelinii haplotypes introgressed into T. macedonicus are similar or identical to those found in T. karelinii itself. Taken together, this suggests that introgression happened after the Last Glacial Maximum. We did not detect more ancient mitochondrial DNA introgression (as found in other taxa by e.g. [22][23][24][25]). Given the climatic oscillations during the Quaternary Ice Age, the two newt species are presumed to have been in periodic contact through time during previous interglacials. However, as subsequent glacial periods would cause their ranges to recede again, any mitochondrial DNA that became introgressed in the area outside of their glacial refugia would be erased.
The climate changed considerably between the Last Glacial Maximum and the present, but we have an intervening snapshot at the mid-Holocene. Compared to the mid-Holocene, the suitability of the location of the present day introgression zone generally increased for T. macedonicus, whereas it decreased for T. karelinii ( Figure 5). Life history differences between the two species provide some suggestion of how T. macedonicus was subsequently able to outcompete T. karelinii, under the changing abiotic environment. The two crested newt species appear to differ in the time annually allocated to an aquatic and a terrestrial life style, with adult T. macedonicus presumed to spend one more month in the water than T. karelinii [13,14]. This suggests that T. macedonicus is adapted to relatively wetter conditions. Bioclimatic variable contribution to the ecological niche models is in line with this hypothesis: the occurrence of T. macedonicus correlates with a higher precipitation than that of T. karelinii (Additional file 5). The climate layers used in the environmental niche modeling provide further insight: since the mid-Holocene, precipitation seasonality in the area where the introgression zone is positioned decreased; whereas precipitation during the wettest season decreased, it increased during the driest season (Additional file 6). Further research is required to understand the adaptive advantage of T. macedonicus over T. karelinii in the introgression zone.
For T. karelinii, the northwestern part of the current introgression zone is predicted to be unsuitable currently, but also in the two past time periods ( Figure 5). Note that T. karelinii probably never colonized this part of the introgression zone: T. karelinii mitochondrial DNA probably reached this region via its T. macedonicus host. The reasoning is as follows: T. macedonicus newts that colonized the northwestern part of the introgression zone derived from the stock in the south of the introgression zone. As this T. macedonicus stock already possessed T. karelinii mitochondrial DNA, no matter if there were T. karelinii present in the northwestern part of the introgression zone, T. macedonicus would spread T. karelinii mitochondrial DNA here, not T. macedonicus mitochondrial DNA. We summarize the hypothesized distribution dynamics of the two crested newt species in Figure 6.

Reinforcement of asymmetric mitochondrial DNA transmission
Geographical species displacement involves a demographical inequality: at the leading edge, the invader is rare whereas the species being displaced is common. This demographic inequality is in itself sufficient to explain asymmetrical mitochondrial DNA introgression [4]. However, the pattern can be reinforced by male biased dispersal and through the operation of Haldane's rule [26]. Males dispersing into territory occupied by a species with which they can hybridize would only encounter females of the indigenous species to mate with. Most of their genetic material can potentially be transmitted to future generations, but their offspring would always contain indigenous mitochondrial DNA, as this is transmitted matrilinearily only. Whether there is male-biased dispersal in crested newts is not known. Haldane's rule [27] states that if one sex is absent, rare, or sterile in F1 offspring, that sex is the heterogametic sex. In the case of crested newts, this would be the male [28]. In effect, markers transmitted via the female line (i.e. mitochondrial DNA) would be particularly prone to introgress [26]. Male biased dispersal and Haldane's rule could act in combination while assortative mating would further amplify asymmetric mitochondrial DNA introgression [4]: if members of the native species prefer to mate among themselves, the invading males and the excess of hybrid females would be more likely to mate with each other.
Even when the parental species occurs in the same frequency, an asymmetry in mitochondrial DNA transmission could be strengthened by factors favoring the relative frequency (prezygotic factors) or viability (postzygotic factors) of one hybrid class over the other [2][3][4]. Examples of pre-and postzygotic selection are known in nature. In the tree frogs Hyla cinerea and H. gratiosa, differences in behavior (courting at the edge versus inside the breeding pond and presence versus absence of satellite mating behavior) cause H. cinerea males to more often intercept H. gratiosa females than vice versa, resulting in biased mitochondrial DNA transmission [29]. For the crested newt T. cristatus and the marbled newt T. marmoratus, F1 hybrids make up about four percent of the adult population where the ranges of the two species overlap [30]. Uneven viability of hybrids causes the frequency of T. marmoratus-mothered hybrids to drop from fifty percent in embryos to ten percent in adults. The T. cristatus-mothered hybrids show an excess of females (i.e. a Haldane's effect). However, the T. marmoratus-mothered hybrids consist of males only, suggesting an incompatibility between the T. cristatus X chromosome and T. marmoratus cytoplasm [30].
In the present case, a higher frequency or reproductive success of T. karelinii-mothered hybrids over T. macedonicus-mothered hybrids would have acted like a filter, hampering the spread of T. macedonicus mitochondrial DNA into the introgression zone. Research into the presence of such asymmetries in T. macedonicus -T. karelinii hybridization is a prospect for future study.

What about positive selection?
Asymmetric mitochondrial DNA introgression due to the positive selection of foreign mitochondrial DNA [6,8,31] is characterized by a distinctly different pattern of geographical spread than if caused by species displacement: instead of the contact zone between two species moving across the initial mitochondrial DNA boundary, which is the case with geographical displacement (Figure 1b), positive selection results in mitochondrial DNA of one species being 'pulled' into the range of the other. A selective sweep, where foreign mitochondrial DNA replaces the original mitochondrial DNA, has been demonstrated experimentally in Drosophila flies [32]. Positive selection has also been suggested to cause asymmetric mitochondrial DNA introgression in nature [33,34]. However, when encountering asymmetric mitochondrial DNA introgression, it is difficult to choose between a scenario of species displacement or natural selection. For example, in Lissotriton newts the original mitochondrial DNA of L. montandoni has been fully replaced by that of L. vulgaris [35]. However, this exchange of mitochondrial DNA did not have a single origin but involved multiple, independent introgression events of well-differentiated L. vulgaris mitochondrial DNA. The observed pattern could either reflect repeated geographical replacement during multiple interglacials, when the more mountainous L. montandoni shifted its range to a lower elevation, or several uptakes of foreign mitochondrial DNA that subsequently spread in L. montandoni by natural selection.
Why do we in the present case favor the geographical displacement scenario over the positive selection scenario? The bisecting of the T. karelinii range by T. macedonicus provides positive evidence in favor of T. macedonicus expanding its range at the expense of T. karelinii. This displacement is further supported by T. karelinii mitochondrial DNA having spread into a region where we infer T. karelinii was never presentthis would require a T. macedonicus host that was expanding its range after having captured T. karelinii mitochondrial DNA. Furthermore, foreign T. karelinii mitochondrial DNA is only present in T. macedonicus in postglacially colonized areas. If positive selection had caused the asymmetric mitochondrial DNA introgression, foreign mitochondrial could also have penetrated into the region that acted as a glacial refugium for T. macedonicus. Only if an adaptive advantage of foreign mitochondrial DNA was somehow limited to the area where the introgression zone is currently positioned, it would be indistinguishable from a scenario involving species displacement. However, we know that T. macedonicus has at least partially displaced T. karelinii, based on the well-supported T. karelinii enclave. Considering that species displacement occurred at least in part of the introgression zone, the most parsimonious explanation is that species displacement alone caused the pattern of asymmetric mitochondrial DNA introgression observed, instead of a combination of species displacement and localized positive selection of mitochondrial DNA. To test this, future research should employ a large battery of nuclear DNA sequences. Selectively neutral nuclear genes would be expected to show a comparable geographical pattern of introgression compared to mitochondrial DNA when asymmetrical mitochondrial DNA introgression was caused by species displacement [4]. This would, however, not be the expectation in the case of positive selection of mitochondrial DNA.

Conclusions
The T. macedonicus -T. karelinii case provides a unique insight into asymmetric mitochondrial DNA introgression. The fact that T. macedonicus bisects the range of T. karelinii provides a clear indication of directionality: T. macedonicus likely invaded T. karelinii territory instead of the other way around. Based on a combination of phylogeography and ecological niche modeling, we suggest that the asymmetric introgression of mitochondrial DNA from T. karelinii into T. macedonicus is due to postglacial species displacement (see the summarized scenario in Figure 6). As these crested newt species expanded their ranges in response to climate change after the conclusion of the Last Glacial Maximum, they came into spatial contact. Subsequently, the newly established hybrid zone between the two moved across the landscape, as T. macedonicus outcompeted T. karelinii. However, the former distribution of T. karelinii, before T. macedonicus started to invade its range, could be inferred based on the asymmetrically introgressed T. karelinii mitochondrial DNA. The asymmetrical introgression of mitochondrial DNA has been regularly documented in a wide range of taxa and can provide key insights into historical biogeography [4,36].

Distribution and genetic data
We included 139 T. macedonicus and 135 T. karelinii localities ( Figure 2 and Additional file 1). Species identity was primarily based on morphological characters (throat and belly pattern and morphological measurements); for a subset (conveniently arranged along transects across the contact zone) diagnostic allozyme markers have been published [17,18,Arntzen et al. submitted]. Note that at five localities the species are found in syntopy. We obtained a 658-bp segment of ND4 for 600 newts (70 previously published in, and 530 newly sequenced following the protocol of [12]) from 71 of the T. macedonicus and 86 of the T. karelinii localities (Figure 2 and Additional file 1). Among them, the sequenced individuals contained 83 haplotypes (Additional file 3). We delimited the zone of asymmetrically introgressed mitochondrial DNA by Thiessen polygons, using ArcGIS (www.esri.com). Each polygon represents a locality and contains the area that is closer to that particular locality than to another one [as in 37].

Genetic analyses
Phylogenetic trees were constructed with MrBayes 3.1.2 [38], employing two, four-chain, twenty million generation runs, with a sampling frequency of 0.001 and a heating parameter of 0.1. MrModeltest [39] identified GTR + I, HKY + I + G and HKY + G as the most suitable models of sequence evolution for codon positions 1, 2 and 3. Data partitions were unlinked. Tracer 1.5 [40] was used to check for stabilization of overall likelihood within and convergence between runs. The first quarter of the sampled trees was discarded as burn-in. As an outgroup we used T. marmoratus (GenBank accession number GU982379). Minimum spanning haplotype networks were created with HapStar 0.5 [41], based on distance matrices produced with Arlequin 3.5 [42]. DnaSP5 [43] was used to determine the nucleotide diversity (π) within and the average number of nucleotide substitutions (D XY ) between groups of haplotypes. To test for signals of demographic expansion, we compared the observed mismatch distribution with that expected under population expansion [44]. We interpreted the differences based on the Ramos-Onsins and Rozas' R 2 and Fu's Fs statistics [45], for which statistical significance was obtained with a 1,000 coalescent simulations in DNAsp, and based on the sum of square deviations, for which statistical was evaluated with a 1,000 bootstrap replicates, in Arlequin 3.5 [42].

Environmental data
For climate layers we used bioclimatic variables at 2.5 arcminute resolution, available from WorldClim 1.4 [46]. To obtain realistic and transferable ecological niche models, it is recommended to mirror the physiological limitations of the study species and minimize the effects of multicolinearity among data layers [47][48][49][50]. Crested newt species differ in the length of their annual aquatic period [13,14]. Therefore, we include a set of layers that likely reflects the availability of water bodies during the breeding season, i.e. seasonal variation in evaporation and precipitation: bio10 = mean temperature of warmest quarter, bio11 = mean temperature of coldest quarter, bio15 = precipitation seasonality, bio16 = precipitation of wettest quarter, and bio17 = precipitation of driest quarter. These layers show a Pearson correlation < 0.7. Bioclimatic variables are also available for the mid-Holocene (~6 Ka) and the Last Glacial Maximum (~21 Ka), being derived from the Paleoclimate Modelling Intercomparison Project phase 2 [51; http://pmip2.lsce.ipsl.fr/]. We used datalayers for the former based on simulations using the ECHAM3 model (Model Max Planck Inst. für Meteorologie) and the latter based on simulations using the Community Climate System Model version 3 (CCSM) [52].

Ecological niche modeling
Ecological niche models were created using Maxent 3.3.3e [53]. We restricted the feature type to hinge features as this produces smoother model fits, so forcing models to be more focused on key trends rather than potential idiosyncrasy in the data, which would hamper extrapolation to a different time or place [54]. We restricted the area from which pseudo-absence was drawn to the distribution of the entire crested newt T. cristatus superspecies. This area was broadly defined as a 200 km buffer zone [55] around known crested newt localities [Wielstra et al., submitted]. The ecological niche models were tested for statistical significance against a null model derived from random localities, to test whether they perform better than expected by random chance [56]. We created a null distribution of 99 AUC values, based on as many random localities as used for the tested species distribution model. The AUC value of the tested species distribution model was treated as a 100 th value and deemed statistically significant if it ranked higher than the 95 th value (i.e. above the 95% confidence interval). Random point data were created with ENMTools 1.3 [57; http://enmtools.blogspot.com/]. The null model approach prevents interpreting model quality based on an arbitrary AUC threshold and precludes the requirement to set aside part of the localities for model testing [56]. Models were projected on the current and Last Glacial Maximum climate layers. We explored variable importance and response curves and conducted a jackknife test to determine how ecological niche models of the two species differed.
Dynamique du Climat (PNEDC). We acknowledge the international modeling groups for providing their data for analysis and the Laboratoire des Sciences du Climat et de l'Environnement (LSCE) for collecting and archiving the model data.