Tracing glacial refugia of Triturus newts based on mitochondrial DNA phylogeography and species distribution modeling
Frontiers in Zoology volume 10, Article number: 13 (2013)
The major climatic oscillations during the Quaternary Ice Age heavily influenced the distribution of species and left their mark on intraspecific genetic diversity. Past range shifts can be reconstructed with the aid of species distribution modeling and phylogeographical analyses. We test the responses of the different members of the genus Triturus (i.e. the marbled and crested newts) as the climate shifted from the previous glacial period (the Last Glacial Maximum, ~21 Ka) to the current interglacial.
We present the results of a dense mitochondrial DNA phylogeography (visualizing genetic diversity within and divergence among populations) and species distribution modeling (using two different climate simulations) for the nine Triturus species on composite maps.
The combined use of species distribution modeling and mitochondrial phylogeography provides insight in the glacial contraction and postglacial expansion of Triturus. The combined use of the two independent techniques yields a more complete understanding of the historical biogeography of Triturus than both approaches would on their own. Triturus newts generally conform to the ‘southern richness and northern purity’ paradigm, but we also find more intricate patterns, such as the absence of genetic variation and suitable area at the Last Glacial Maximum (T. dobrogicus), an ‘extra-Mediterranean’ refugium in the Carpathian Basin (T. cristatus), and areas where species displaced one another postglacially (e.g. T. macedonicus and western T. karelinii). We provide a biogeographical scenario for Triturus, showing the positions of glacial refugia, the regions that were postglacially colonized and the areas where species displaced one another as they shifted their ranges.
The Quaternary Ice Age (~2.59 Ma-present) is associated with large climatic oscillations [1–4]. Long, cold and dry glacial cycles are alternated by relatively short, warm and wet interglacials. The transition between the two takes place in a geological blink of an eye [5, 6]. The climatic oscillations have a major impact on species distribution patterns . Generally, impacts of glacial cycles are more extreme further away from the equator (and at higher elevations): areas at a higher latitude (and elevation) become inhospitable, whereas climate at a lower latitude (and elevation) remains habitable [1, 7]. Populations at higher latitude may cope with climate change in situ, through adaptation or phenotypic plasticity, or alternatively they may track suitable habitat [6, 8]. However, a more likely outcome is that such populations go extinct [8–10]. At lower latitudes, populations can endure glacial periods relatively unimpaired, in so called glacial refugia [2, 10, 11]. During subsequent interglacials, species can reclaim their former distribution, by rapidly recolonizing the large tracts of still uninhabited, but freshly habitable land from glacial refugia .
The fluctuating climate during the Quaternary has left its mark on patterns of genetic diversity [1–3, 5, 12]. Populations persisting in glacial refugia have a relatively long and stable demographical history compared to those in areas claimed postglacially. As a result, populations in glacial refugia are characterized by high levels of genetic diversity, whereas populations established after the most recent glacial cycle typically show little genetic variation. This is the concept – devised from a northern hemisphere perspective – of ‘southern richness and northern purity’ . Furthermore, species displacement after postglacial secondary contact regularly coincides with introgression of genetic material (especially of cytoplasmic DNA) [12, 13]. By uncovering the spatial structuring of genetic lineages between glacial refugia and along recolonization routes, and by detecting mismatches between genetic markers and species boundaries, phylogeographical surveys provide insights in past distribution rearrangement .
Independent from genetics, species distribution modeling can be applied to answer historical biogeographical questions . Species distribution modeling involves the approximation of the ecological requirements of a species, based on the range of environmental conditions experienced at known localities . The constructed model can then be extrapolated on current climate layers, to determine the species’ potential distribution. Similarly, the model can be projected on climatological reconstructions of the past. Niches evolve over time, questioning the validity of predicting past distributions based on present day models. However, considering the relatively short time period spanning the shift from the last glacial cycle to the present day, niche conservatism is a realistic assumption [16, 17]. Comparison of present and past potential distribution provides information on range shifts.
The distribution of amphibians is tightly linked to environmental conditions and thus has to follow suit in the face of rapid climate change [18–21]. The marbled and crested newts (genus Triturus) are a group of nine closely related species [22, 23] (Table 1). Triturus newts are distributed across most of Europe and adjacent Asia (Figure 1). They are found in regions generally regarded as important glacial refugia, such as the Iberian, Italian and Balkan Peninsulas, Anatolia, Caucasia and the southern Caspian basis [3, 9, 10]. On the other hand, Triturus newts occupy large tracts of land which would have been uninhabitable during glacial periods, particularly temperate Eurasia . Thus, within this single model system, we expect to observe varying responses to glaciation.
We first conduct a dense phylogeographical survey: employing mitochondrial DNA, we determine geographical genetic structuring (i.e. diversity within and divergence among populations) for each of the different Triturus species. We also delineate areas where mitochondrial DNA has been asymmetrically introgressed. Subsequently, we approximate the distributions of the different Triturus species at the Last Glacial Maximum (~21Ka) using species distribution modeling. The outcome of the two independent approaches is visualized on composite maps, which show the nine different species side by side. Our hypothesis is that those areas predicted suitable at the Last Glacial Maximum based on species distribution models are also the ones showing the highest genetic diversity. Vice versa, we expect areas suggested to have been unsuitable at the time are genetically depleted. We summarize the signatures of past distribution dynamics as inferred from mitochondrial DNA phylogeography and species distribution modeling on a map.
Generally the identification of Triturus newts is straightforward based on geography (Figure 1). For individuals occurring near the contact zone more care should be taken. The Triturus individuals included were in the first place identified based on morphology [26, 27] and for a subset nuclear genetic data confirming their identity  (Arntzen et al., submitted; Wielstra et al., unpublished data).
We included genetic data (658 bp of subunit 4 of the NADH dehydrogenase gene complex; ND4) for 2470 Triturus newts, representing 493 populations (Figure 2 and Additional file 1). A large proportion of these individuals (n = 1795) was taken from previous studies (published [25, 28–32] or submitted [Arntzen et al.]; see Additional file 1 for details). The remainder (n = 675) was newly sequenced for the current paper, using the primers in Table 2 and following the protocol outlined in . Sequences were manually aligned and identical ones merged into haplotypes using MacClade 4.08 . Our purpose was 1) to test for mitochondrial DNA introgression between species, and 2) to infer the spatial genetic variation within each species.
To detect which individuals retain their ‘original’ mitochondrial DNA and which possess introgressed mitochondrial DNA (derived from one of the other species), we first assigned mitochondrial DNA haplotypes to species by conducting a Neighbor Joining analysis in MEGA 5.05 . Calotriton asper was used as outgroup (sequence taken from ; GenBank accession number: GU982378).
To infer interspecific geographical structuring, we first excluded introgressed mitochondrial DNA (as it does not properly reflect the evolutionary history of either ‘host’ or ‘donor’). Subsequently, we determined – for each species – 1) the genetic diversity within populations and 2) the genetic distance among populations. The mean number of pairwise differences among haplotypes (π), as determined with Arlequin 3.5 , was used as a measure of genetic diversity within populations. To determine the genetic distance among populations we used Alleles in Space 1.0 . Alleles in Space connects the populations in a network, based on Delaunay triangulation. Subsequently, the program produces values of average genetic distance among the populations that are connected by the network (the proportion of mismatched nucleotide sites; Z i ). These values are positioned at the midpoints of the connections in the network. For the measure of genetic diversity within populations we only included populations for which more than one sequence was available (π will always be zero for populations with only one sequence).
We interpolated the values for π and Z i across geographical space using inverse distance weighting in the spatial analyst extension of ArcGIS (http://www.esri.com) . The output for each species was cropped according to its distribution range (Figure 1). For both π and Z i , we compiled a single composite map from the nine crops (i.e. one for each Triturus species). We used a color scheme running from blue to red to reflect low to high values of π and Z i . Using a single scale for all Triturus species facilitates comparison among them, but runs the risk that variation in genetically poor species is overshadowed by that of genetically rich species. Hence, we additionally apply a species specific scale to better express intraspecific structuring. This way we managed to visualize 1) regions of relatively high and low intraspecific genetic diversity and 2) the relative genetic divergence among populations within species.
Species distribution modeling approach
We composed a database of 4532 Triturus localities (incorporating the 493 populations included in the genetic analysis; see Figure 2 and Additional file 2). This database is based on  and updated in the light of the taxonomic developments  (see Table 1). For most species the majority of localities concern accurate locations, but in particularly for the wide ranging T. cristatus we mostly only had the name of the nearest village of a locality or digitized atlas data to our disposal. Species distribution models are still expected to perform relatively accurately if spatial error in localities is introduced and the method we use (Maxent, see below) is robust against such errors . We used the bioclimatic variables available at a 2.5 arcminute resolution (c. 5 × 5 km) from the WorldClim database version 1.4  to calibrate our species distribution models. In order to prevent model overfitting, which would negatively influence transferability [16, 42], we minimized multicolinearity among data layers by selecting a subset that showed a Pearson’s correlation of r < 0.7. Furthermore, focusing on climate layers that are deemed biologically meaningful based on life history knowledge of the model system yields the most appropriate species distribution models . Variation in the availability of standing water during the breeding season appears to have driven the ecological radiation among Triturus newts . Taking these two points into consideration, we included a set of layers that encompasses seasonal variation in evaporation and precipitation, in casu 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. Bioclimatic variables are also available from the WorldClim database for the Last Glacial Maximum (~21 Ka). These data are derived from the Paleoclimate Modelling Intercomparison Project phase 2 ; [http://pmip2.lsce.ipsl.fr/] and based on two climate simulations: the Model for Interdisciplinary Research on Climate version 3.2 (MIROC) and the Community Climate System Model version 3 (CCSM) .
Species distribution models were created with Maxent 3.3.3 k . We restricted the feature type to hinge features as this produces smoother model fits, so forcing models to be focused on key trends rather than potential idiosyncrasy in the data. This approach facilitates extrapolation to a different time (or place) . The environmental range covered by the pseudo-absence data, used to discriminate presence data from background, should neither be too narrow nor too broad [48–50]. Too narrow a range results in complex models that do not generalize well, whereas too broad a range results in too simple models that focus on coarse-scale and neglect fine-scale variation. A practical solution is to focus on area that is potentially accessible to the species of interest if it were not for abiotic factors, i.e. an area where spread would not be hampered by major physical barriers. However, competition could lead to further exclusion of the target species, even though the area would be suited in the absence of such biotic interactions . Taking these considerations into account, we restricted the area from which pseudo-absence was drawn to the distribution of the entire genus Triturus. This area was broadly defined as a 200 km buffer zone  around known Triturus localities (Additional file 2).
The species distribution models were tested for statistical significance against a null model derived from random localities . For each tested species distribution model, we created a null distribution of 99 AUC values. These AUC values were derived from species distribution models, 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 100th value and deemed statistically significant if it ranked higher than the 95th value (i.e. above the 95% confidence interval). Random point data were created with ENMTools 1.3 ; [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 .
The species distribution models were projected on the current and Last Glacial Maximum climate layers. Composite maps were created for each set of climate layers (in the same way as explained for the genetic approach). Maxent provides predicted probability values between zero and one and we used a color scheme running from blue to red to reflect these values.
The 2470 Triturus sequenced newts comprise 315 haplotypes (see Additional file 1 for details and Additional file 3 for GenBank accession numbers). The Neighbor Joining phylogeny shows that haplotypes group in nine reciprocally monophyletic lineages, corresponding to species (Figure 3). A large number of individuals (n = 408; c. 16.5%) that belong to one species, possess mitochondrial DNA characteristic of another (Figure 2 and Additional file 1). This mitochondrial DNA introgression is mostly restricted to near the contact zones; only T. pygmaeus, T. macedonicus, T. cristatus and central T. karelinii have extended ranges in which foreign mitochondrial DNA is present (derived from T. marmoratus for the first one and from western T. karelinii for the other three).
The number of populations included to determine genetic diversity within (π) and genetic divergence among (Z i ) populations was 44 and 48 out of 52 populations for T. carnifex, 60 and 88 out of 104 for T. cristatus, 29 and 40 out of 44 for T. dobrogicus, 55 and 79 out of 79 for western T. karelinii, 13 and 16 out of 24 for central T. karelinii, 27 and 31 out of 31 for eastern T. karelinii, 32 and 37 out of 72 for T. macedonicus, 34 and 36 out of 36 for T. marmoratus and 49 and 50 out of 51 for T. pygmaeus (populations containing asymmetrically introgressed mitochondrial DNA were excluded completely and populations represented by only a single individual were additionally excluded in the calculation of π). Measures of genetic diversity within and genetic divergence among populations can be found in Additional file 1 and 4. Composite maps visualizing genetic structuring in each species are depicted in Figures 4 and 5. Because Alleles in Space places the values for genetic divergence among populations at the midpoint of the network connecting the populations, information falling outside the current range is lost in Figure 5 (e.g. in the case of two allopatric populations). Uncropped maps for each species, showing a more comprehensive picture per species but at the cost of conciseness, can be found in Additional file 5. Triturus species generally show considerable spatial variation in their genetic composition (reflected by ‘warm’ and ‘cold’ areas in Figures 4 and 5).
Species distribution models perform significantly better than random (Additional file 6). Composite maps depicting the predicted suitability of each species’ distribution range at the Last Glacial Maximum, based on the two different climate simulations (MIROC and CCSM), are provided in Figure 6. As species were not necessarily bound to their current ranges through time, projections for each species on a wider area (roughly the distribution of the entire genus Triturus) are provided in Additional file 7. The ranges of all Triturus species are predicted to have been restricted at the Last Glacial Maximum. A composite map showing species distribution models projected on present day climate layers is provided in Figure 6. Predicted suitability of the distribution range of each species under current conditions shows a considerable overlap with the sketched outlines of their ranges (Figure 1).
Before we discuss our results, we first consider where care should be taken when interpreting genetic or species distribution modeling results in isolation and explain how interpretation of the results of the two independent techniques together at least partially negates these problems. Then we identify the probable glacial refugia and postglacially colonized area for each Triturus species. Finally, we summarize our scenario regarding the biogeographical response of Triturus newts to the climate shift from the previous glacial period to the current interglacial on a map.
Considerations and combination of the genetic and species distribution modeling results
A high value for genetic diversity within populations (Figure 4) can result from two distinct processes: postglacial secondary contact from distinct source populations or long term presence in situ throughout glacial-interglacial cycles (potentially enforced by a low connectivity between populations due to e.g. xeric environment) . To separate these two processes, phylogenetic relationships among haplotypes provide some insight. For example, the obvious ‘hotspot’ found in T. cristatus (SE Poland) likely reflects the comingling of haplotypes belonging to the two distinct clades present in this species (Figure 3; Additional file 1). On the other hand, the hotspot found in eastern T. karelinii (SE Azerbaijan) likely represents standing genetic variation, given the lack of phylogeographic structure among the haplotypes involved. Species distribution modeling can further help to distinguish both processes by establishing whether an area became suitable only recently or whether it has continuously been so on the long term.
Maps showing genetic divergence among populations (Figure 5) reflect a different aspect of geographical genetic structuring. On the one hand, areas with a high genetic overturn among populations (e.g. the southern part of the range of T. macedonicus) show many hotspots. Such areas would be expected to have been habitable on the long term. On the other hand, genetically homogenous areas, due to high levels of gene flow and/or demographic expansion stand out as cold areas (e.g. both of these processes likely determine the absence of genetic structuring among T. dobrogicus populations). These two processes are more difficult to disentangle: whereas demographic expansion would particularly be expected in areas that became suitable postglacially, gene flow could reduce genetic structuring both inside and outside of glacial refugia. In some cases geographical barriers to gene flow stand out, such as the Greater Caucasus mountain range separating eastern T. karelinii on its northern and southern side. Such barriers reflect climatically unsuitable regions and can be identified as such in the species distribution models.
The reliability of genetic diversity values for populations depends on sample size, and some populations are relatively poorly sampled in our study. The genetic turnover between populations provides additional insight that counters the adverse effects of low sample sizes. If populations within an area are genetically diverse, but sample sizes per population are low, genetic diversity per population might be underestimated due to chance effects in sampling. However, these effects are random per population. Therefore, the more populations that are being compared within that area, the less likely it is that actual genetic diversity remains hidden by chance. For this same reason determining the genetic turnover between populations is less susceptible to sample size to begin with: low sample size per population adds up to a big sample sizes within a region. So, in short, if multiple populations within a region are found to be genetically similar, even though samples sizes per population are low, together they provide a good indication that genetic diversity in the region really is low. We thus recommend exploring genetic diversity within populations and genetic divergence among populations together to improve the estimate of spatial genetic variation, especially when sample sizes are low.
We present the results of a single (mitochondrial) genetic marker. Although we consider our results indicative (and use the species distribution modeling results as an independent check), it should be taken into account that a single marker may not accurately reflect species history . In fact, we explicitly exploit the mismatch between species history and mitochondrial DNA in the case of asymmetrically introgressed mitochondrial DNA. Asymmetrically introgressed mitochondrial DNA may suggest that species displaced one another (e.g. [29, 55]). For Triturus, the similarity of the introgressed mitochondrial DNA to that present in the ‘donor’ species suggests a recent (postglacial) transference across the species boundary. Given the climatic oscillations during the Quaternary Ice Age, the Triturus newt species are presumed to have been in periodic contact through time. However, we did not identify ancient mitochondrial DNA introgression events (as in e.g. [56–58]). A potential explanation is that introgressed mitochondrial DNA was limited to areas of postglacial expansion and never penetrated those areas that repeatedly acted as glacial refugia. In effect, any mitochondrial DNA that became introgressed during an interglacial would be erased by the next glacial period. Species distribution modeling allows us to test this hypothesis.
A disadvantage of species distribution models is that they cannot distinguish between potential and realized distribution: areas that had suitable conditions at the Last Glacial Maximum may not have actually been inhabited. Genetic information can narrow down the position of glacial refugia (e.g. both France and the Carpathian Basin are suggested to have had suitable conditions during the Last Glacial Maximum for T. cristatus but only the latter shows the high genetic variation expected to be present in a former glacial refugium). Another issue is that species distribution models are based on the climatic conditions currently experienced by a species, but the actual conditions under which they could occur might be broader (e.g. no suitable conditions were predicted to have been present at the Last Glacial Maximum for T. dobrogicus, but yet it did survive to the present day).
Glacial refugia and area postglacially colonized for each Triturus species
The two marbled newt species illustrate the different effect of glaciations on species with a relatively southern and relatively northern distribution. The northern T. marmoratus is predicted to have withdrawn its range at the Last Glacial Maximum more extensively than the southern T. pygmaeus (especially according to the CCSM climate simulation; Figure 6). In T. marmoratus genetic structuring is highest in the south of its range, but for T. pygmaeus it is not (Figures 4 and 5). The presence of T. marmoratus mitochondrial DNA in the northern part of the range of T. pygmaeus suggests that in this area the former was replaced by the latter [31, 59, 60]. This is confirmed by the distribution models projected on Last Glacial Maximum climate data outside of the current species ranges (Additional file 7): at the time, area predicted suitable for T. marmoratus (but not for T. pygmaeus) was present in the northern part of the area currently occupied by T. pygmaeus. Striking is that both T. marmoratus and T. pygmaeus mitochondrial DNA occurs throughout this area (Figure 2). Although a species will often simply show foreign mitochondrial DNA where it replaced another  see below for examples in Triturus, exceptions similar to the T. marmoratus-T. pygmaeus case are known (e.g. ).
The species distribution models based on the two climate simulations agree that most of the range of T. cristatus was unsuitable at the Last Glacial Maximum (Figure 6). Suitable area was suggested to be present at the Last Glacial Maximum in both the Carpathian Basin and in France. The genetic data suggest that only the Carpathian Basin acted as a glacial refugium at the Last Glacial Maximum as genetic diversity here is high (Figures 3, 4, 5). This patterns corresponds to the geographical distribution of genetic variation based on nuclear DNA markers found in previous studies (MHC genes, allozymes and microsatellites [62, 63]). This also suggests that the remainder of T. cristatus’ current range (including France) was colonized post-glacially, as the species is genetically depleted outside of the Carpathian Basin. Although T. cristatus shows asymmetrically introgressed western T. karelinii mitochondrial DNA south of the Danube river, relatively distinct T. cristatus haplotypes also occur here (Figures 2 and 3, Additional file 1), suggesting a complex pattern of both range expansion at the cost of western T. karelinii but also long term presence south of the Danube.
For T. carnifex there is a discrepancy between the two climate simulations, with MIROC showing a more extensive predicted distribution at the Last Glacial Maximum than CCSM. Particularly CCSM suggests a more considerable range reduction and no suitable area in the current northern Balkan range (Figure 6). The genetic data better correspond to the modeled distribution based on MIROC, as it shows distinct genetic clusters distributed on both sides of the Adriatic Sea and high genetic diversity in Italy (Figures 3, 4, 5), supporting long term presence in the Balkans and a stable population in Italy (see  for more detail). Both climate simulations agree that T. carnifex colonized the part of its range east and north of the Alps postglacially. The T. carnifex mitochondrial DNA introgressed into T. cristatus where the two species meet suggests that Balkan T. carnifex were involved in this range expansion.
Compared to the other Triturus species, intraspecific genetic diversity (i.e. phylogenetic depth) is low in T. dobrogicus (Figure 3). Triturus dobrogicus shows a ‘starburst’ phylogeny – a large number of similar haplotypes – suggesting a demographic expansion after a bottleneck. At the Last Glacial Maximum, the range of T. dobrogicus is predicted to have been unsuitable (Figure 6). This is in line with a population bottleneck. However, the species distribution models do not make it possible to pinpoint a potential refugium: the entire currently occupied range was predicted unsuitable at the Last Glacial Maximum. Suitable area was also not available outside the current range (Additional file 7). Evidently, T. dobrogicus must have had a refugium somewhere and given its specialization on river floodplains (e.g. ) this was likely positioned within its current range. Furthermore, compared to the other Triturus species T. dobrogicus does not show reduced genetic diversity of nuclear DNA (allozymes ). This suggests that perhaps the actual ecological niche of T. dobrogicus is broader than is currently realized (and used to construct the species distribution models). Another option would be that the niche of T. dobrogicus has evolved in a brief time interval to adapt to the environmental shift . We consider this less likely as T. dobrogicus must have survived multiple glacial-interglacial cycles in situ.
The western part of the current range of T. macedonicus is predicted to have been suitable during the Last Glacial Maximum (Figure 6). This does not fully correspond to mitochondrial DNA data, which shows structuring in the southern part of the range, including the southeast, but not in the northwest (Figures 4 and 5). Genetic variation in T. macedonicus is the highest of all Triturus newts (Figure 3). In a large part of its range, T. macedonicus possesses western T. karelinii mitochondrial DNA (Figure 2). This can be explained by T. macedonicus displacing western T. karelinii there after the Last Glacial Maximum as the area only became unsuitable for T. macedonicus after the Last Glacial Maximum (Figure 6; see  for a detailed analysis).
In western T. karelinii, a basal split separates a clade occurring in European Turkey and extreme southeastern Bulgaria (Figure 3). The species distribution modeling approach suggests suitable area was present here during the Last Glacial Maximum (Figure 6). The genetic distinction of this clade is expressed by the ‘barrier’ in Figure 5 surrounding the populations in which it is represented. The other clade shows extensive structuring in Asiatic Turkey, suggesting the region acted as a refugium. Suitability at the Last Glacial Maximum is confirmed by the species distribution modeling approach. A single sub-clade derived from (i.e. nested within) Asiatic Turkey has recently colonized the Balkan part of the range; it shows a starburst phylogeny in line with demographic growth during a range expansion . The species distribution modeling approach agrees that the area only became hospitable after the Last Glacial Maximum.
Central T. karelinii shows a fragmented distribution at the Last Glacial Maximum (Figures 6 and 7). This is reflected by an east west divide (the ‘barrier’ in Figure 5) between two distinct clades separated by a basal split (Figure 3), which are currently in secondary contact (the ‘warm’ spot in Figure 4). Central T. karelinii contains western T. karelinii mitochondrial DNA in northwestern Asiatic Turkey (Figure 2), suggesting postglacial species displacement (see  for a detailed analysis).
For eastern T. karelinii, the southern Caspian basin and western Georgia (Colchis) are suggested to have harbored suitable area at the Last Glacial Maximum according to species distribution models (Figure 6). Genetic diversity is particularly high along the southern Caspian Sea shore (Figures 3, 4, 5). Relatively recently a single nested clade colonized the Caucasus and Crimea from the Caspian part of the range . However, this clade is still distinct (reflected by a relatively deep coalescence; see Figure 3), suggesting that colonization happened during a previous interglacial and thus that the clade persisted here in a glacial refugia during one or several glacial cycles.
Historical biogeographical scenario
We provide a historical biogeographical scenario for the genus Triturus in Figure 7, summarizing glacial refugia, areas of postglacial expansion and regions where species displaced one another. The genetic and species distribution modeling results provide insight into the distribution dynamics of Triturus in response to climate change during the Quaternary Ice Age. Triturus conforms to the general pattern of the Iberian, Italian and Balkan Peninsulas functioning as glacial refugia [1, 67]. Furthermore, several regions that are increasingly appreciated as glacial refugia are also identified as such for Triturus: Anatolia, the southern Caspian Basin and the Caucasus region (e.g. [10, 68]). The novel idea of extra-Mediterranean refugia, positioned more to the north , also applies to Triturus: the glacial refugium of T. cristatus was situated in the Carpathian Basin.
Similarly, areas typically identified as having been postglacially colonized from glacial refugia are predicted as such for Triturus: T. marmoratus, T. carnifex and (to a lesser extent) ‘eastern T. karelinii’ extended their ranges into temperate Europe from their southerly positioned refugia. However, postglacial expansion is best exemplified by T. cristatus. The contemporary range of T. cristatus outside its glacial refugium in the Carpathian Basin encompasses a huge (c. 4.75 million km2) postglacially acquired area, stretching all the way from western Europe to Scandinavia and central Russia. This gives an indication of the speed with which postglacial colonization can be accomplished.
The combination of phylogeography and species distribution modeling aids locating of glacial refugia . The two approaches to visualize intraspecific geographical genetic structuring together illustrate which areas are genetically rich and thus suspected to reflect long term inhabited area and which areas are genetically poor and thus probably only recently became habitable (Figures 4 and 5). Additionally, asymmetrically introgressed mitochondrial DNA suggests where species displaced one another upon post-glacial secondary contact (Figure 2). Independently, the two reconstructions of potential distribution at the Last Glacial Maximum (based on the MIROC and CCSM climate simulations) provide an indication of which areas were habitable at the time and which were unsuitable for Triturus (Figure 6). We provide a qualitative comparison of the results of the two independent techniques (summarized in Figure 7). We anticipate that future developments in the synthesis of phylogeography and species distribution modeling will make it possible to integrate the two approaches into a single analysis even further.
Hewitt G: The genetic legacy of the Quaternary ice ages. Nature. 2000, 405: 907-913. 10.1038/35016000.
Keppel G, Van Niel KP, Wardell-Johnson GW, Yates CJ, Byrne M, Mucina L, Schut AGT, Hopper SD, Franklin SE: Refugia: identifying and understanding safe havens for biodiversity under climate change. Global Ecol Biogeogr. 2012, 21: 393-404. 10.1111/j.1466-8238.2011.00686.x.
Schmitt T: Molecular biogeography of Europe: Pleistocene cycles and postglacial trends. Front Zool. 2007, 4: 11-10.1186/1742-9994-4-11.
Cane MA, Braconnot P, Clement A, Gildor H, Joussaume S, Kageyama M, Khodri M, Paillard D, Tett S, Zorita E: Progress in paleoclimate modeling. J Climate. 2006, 19: 5031-5057. 10.1175/JCLI3899.1.
Hofreiter M, Stewart J: Ecological change, range fluctuations and population dynamics during the Pleistocene. Curr Biol. 2009, 19: R584-R594. 10.1016/j.cub.2009.06.030.
Hof C, Levinsky I, Araújo MB, Rahbek C: Rethinking species’ ability to cope with rapid climate change. Glob Change Biol. 2011, 17: 2987-2990. 10.1111/j.1365-2486.2011.02418.x.
Ohlemüller R, Huntley B, Normand S, Svenning J-C: Potential source and sink locations for climate-driven species range shifts in Europe since the Last Glacial Maximum. Global Ecol Biogeogr. 2011, 21: 152-163.
Hoffmann AA, Sgro CM: Climate change and evolutionary adaptation. Nature. 2011, 470: 479-485. 10.1038/nature09670.
Provan J, Bennett KD: Phylogeographic insights into cryptic glacial refugia. Trends Ecol Evol. 2008, 23: 564-571. 10.1016/j.tree.2008.06.010.
Stewart JR, Lister AM, Barnes I, Dalen L: Refugia revisited: Individualistic responses of species in space and time. P R Soc Lond B-Biol Sci. 2010, 277: 661-671. 10.1098/rspb.2009.1272.
Bennett KD, Provan J: What do we mean by ’refugia’?. Quaternary Sci Rev. 2008, 27: 2449-2455. 10.1016/j.quascirev.2008.08.019.
Excoffier L, Foll M, Petit RJ: Genetic consequences of range expansions. Annu Rev Ecol Evol Syst. 2009, 40: 481-501. 10.1146/annurev.ecolsys.39.110707.173414.
Currat M, Ruedi M, Petit RJ, Excoffier L: The hidden side of invasions: Massive introgression by local genes. Evolution. 2008, 62: 1908-1920.
Svenning J-C, Fløjgaard C, Marske KA, Nógues-Bravo D, Normand S: Applications of species distribution modeling to paleobiology. Quaternary Sci Rev. 2011, 30: 2930-2947. 10.1016/j.quascirev.2011.06.012.
Kozak KH, Graham CH, Wiens JJ: Integrating GIS-based environmental data into evolutionary biology. Trends Ecol Evol. 2008, 23: 141-148. 10.1016/j.tree.2008.02.001.
Peterson AT: Ecological niche conservatism: a time-structured review of evidence. J Biogeogr. 2011, 38: 817-827. 10.1111/j.1365-2699.2010.02456.x.
Wiens JJ, Graham CH: Niche conservatism: Integrating evolution, ecology, and conservation biology. Annu Rev Ecol Evol Syst. 2005, 36: 519-539. 10.1146/annurev.ecolsys.36.102803.095431.
Buckley LB, Jetz W: Linking global turnover of species and environments. Proc Natl Acad Sci USA. 2008, 105: 17836-17841. 10.1073/pnas.0803524105.
Vieites DR, Nieto-Román S, Wake DB: Reconstruction of the climate envelopes of salamanders and their evolution through time. Proc Natl Acad Sci USA. 2009, 106: 19715-19722. 10.1073/pnas.0902956106.
Araújo MB, Nogués-Bravo D, Diniz-Filho JAF, Haywood AM, Valdes PJ, Rahbek C: Quaternary climate changes explain diversity among reptiles and amphibians. Ecography. 2008, 31: 8-15. 10.1111/j.2007.0906-7590.05318.x.
Araújo MB, Thuiller W, Pearson RG: Climate warming and the decline of amphibians and reptiles in Europe. J Biogeogr. 2006, 33: 1712-1728. 10.1111/j.1365-2699.2006.01482.x.
Arntzen JW, Wielstra B: Where to draw the line? A nuclear genetic perspective on proposed range boundaries of the crested newts Triturus karelinii and T. arntzeni. Amphibia-Reptilia. 2010, 31: 311-322. 10.1163/156853810791769509.
Wielstra B, Baird AB, Arntzen JW: A multimarker phylogeography of crested newts (Triturus cristatus superspecies) reveals cryptic species. Mol Phylogenet Evol. 2013, 67: 167-175. 10.1016/j.ympev.2013.01.009.
Hewitt GM: Genetic consequences of climatic oscillations in the Quaternary. Philos Trans R Soc Lond B Biol Sci. 2004, 359: 183-195. 10.1098/rstb.2003.1388.
Wielstra B, Arntzen JW: Unraveling the rapid radiation of crested newts (Triturus cristatus superspecies) using complete mitogenomic sequences. BMC Evol Biol. 2011, 11: 162-10.1186/1471-2148-11-162.
Arntzen JW: Triturus cristatus Superspecies - Kammolch-Artenkreis (Triturus cristatus (Laurenti, 1768) - Nördlicher Kammolch, Triturus carnifex (Laurenti, 1768) - Italienischer Kammolch, Triturus dobrogicus (Kiritzescu, 1903) - Donau-Kammolch, Triturus karelinii (Strauch, 1870) - Südlicher Kammolch). Handbuch der Reptilien und Amphibien Europas Schwanzlurche IIA. Edited by: Grossenbacher K, Thiesmeier B. 2003, Wiebelsheim: Aula-Verlag, 421-514.
Arntzen JW, Wallis GP: Geographic variation and taxonomy of crested newts (Triturus cristatus superspecies): Morphological and mitochondrial data. Contrib Zool. 1999, 68: 181-203.
Wielstra B, Espregueira Themudo G, Olgun K, Poyarkov NA, Arntzen JW, Güclü Ö: Cryptic crested newt diversity at the Eurasian transition: The mitochondrial DNA phylogeography of Near Eastern Triturus newts. Mol Phylogenet Evol. 2010, 56: 888-896. 10.1016/j.ympev.2010.04.030.
Wielstra B, Arntzen JW: Postglacial species displacement in Triturus newts deduced from asymmetrically introgressed mitochondrial DNA and ecological niche models. BMC Evol Biol. 2012, 12: 161-10.1186/1471-2148-12-161.
Ivanović A, Wielstra B, Olgun K, Litvinchuk SN, Kalezić ML, Arntzen JW, Üzüm N: Is mitochondrial DNA divergence of Near Eastern crested newts (Triturus karelinii group) reflected by differentiation of skull shape?. Zool Anz. 2013, 252: 269-277. 10.1016/j.jcz.2012.08.005.
Espregueira Themudo G, Nieman A, Arntzen JW: Is dispersal guided by the environment? A comparison of interspecific gene flow estimates among differentiated regions of a newt hybrid zone. Mol Ecol. 2012, 21: 5324-5335. 10.1111/mec.12026.
Canestrelli D, Salvi D, Maura M, Bologna MA, Nascetti G: One Species, three Pleistocene evolutionary histories: Phylogeography of the Italian crested newt. Triturus carnifex. PLoS One. 2012, 7: e41754-10.1371/journal.pone.0041754.
Arntzen JW, Espregueira Themudo G, Wielstra B: The phylogeny of crested newts (Triturus cristatus superspecies): Nuclear and mitochondrial genetic characters suggest a hard polytomy, in line with the paleogeography of the centre of origin. Contrib Zool. 2007, 76: 261-278.
Maddison DR, Maddison WP: MacClade 4: Analysis of phylogeny and character evolution, version 4.08. 2005, Sunderland (Massachusetts): Sinauer Associates
Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S: MEGA5: Molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011, 28: 2731-2739. 10.1093/molbev/msr121.
Excoffier L, Lischer HEL: Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010, 10: 564-567. 10.1111/j.1755-0998.2010.02847.x.
Miller MP: Alleles In Space (AIS): Computer software for the joint analysis of interindividual spatial and genetic information. J Hered. 2005, 96: 722-724. 10.1093/jhered/esi119.
Fijarczyk A, Nadachowska K, Hofman S, Litvinchuk SN, Babik W, Stuglik M, Gollmann G, Choleva L, Cogălniceanu DAN, Vukov T: Nuclear and mitochondrial phylogeography of the European fire-bellied toads Bombina bombina and Bombina variegata supports their independent histories. Mol Ecol. 2011, 20: 3381-3398. 10.1111/j.1365-294X.2011.05175.x.
Wielstra B, Beukema W, Arntzen JW, Skidmore AK, Toxopeus AG, Raes N: Corresponding mitochondrial DNA and niche divergence for crested newt candidate species. PLoS One. 2012, 7: e46671-10.1371/journal.pone.0046671.
Graham CH, Elith J, Hijmans RJ, Guisan A, Townsend Peterson A, Loiselle BA, The Nceas Predicting Species Distributions Working G: The influence of spatial errors in species occurrence data used in distribution models. J Appl Ecol. 2008, 45: 239-247.
Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A: Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005, 25: 1965-1978. 10.1002/joc.1276.
Arntzen J: From descriptive to predictive distribution models: a working example with Iberian amphibians and reptiles. Front Zool. 2006, 3: 8-10.1186/1742-9994-3-8.
Rödder D, Schmidtlein S, Veith M, Lötters S: Alien invasive slider turtle in unpredicted habitat: A matter of niche shift or of predictors studied?. PLoS One. 2009, 4: e7843-10.1371/journal.pone.0007843.
Braconnot P, Otto-Bliesner B, Harrison S, Joussaume S, Peterchmitt JY, Abe-Ouchi A, Crucifix M, Driesschaert E, Fichefet T, Hewitt CD: Results of PMIP2 coupled simulations of the Mid-Holocene and Last Glacial Maximum - Part 1: experiments and large-scale features. Clim Past. 2007, 3: 261-277. 10.5194/cp-3-261-2007.
Collins WD, Bitz CM, Blackmon ML, Bonan GB, Bretherton CS, Carton JA, Chang P, Doney SC, Hack JJ, Henderson TB: The Community Climate System Model Version 3 (CCSM3). J Climate. 2006, 19: 2122-2143. 10.1175/JCLI3761.1.
Phillips SJ, Anderson RP, Schapire RE: Maximum entropy modeling of species geographic distributions. Ecol Model. 2006, 190: 231-259. 10.1016/j.ecolmodel.2005.03.026.
Elith J, Kearney M, Phillips S: The art of modelling range-shifting species. Methods Ecol Evol. 2010, 1: 330-342. 10.1111/j.2041-210X.2010.00036.x.
VanDerWal J, Shoo LP, Graham C, Williams SE: Selecting pseudo-absence data for presence-only distribution modeling: How far should you stray from what you know?. Ecol Model. 2009, 220: 589-594. 10.1016/j.ecolmodel.2008.11.010.
Stokland JN, Halvorsen R, Støa B: Species distribution modelling - Effect of design and sample size of pseudo-absence observations. Ecol Model. 2011, 222: 1800-1809. 10.1016/j.ecolmodel.2011.02.025.
Elith J, Phillips SJ, Hastie T, Dudík M, Chee YE, Yates CJ: A statistical explanation of MaxEnt for ecologists. Divers Distrib. 2011, 17: 43-57. 10.1111/j.1472-4642.2010.00725.x.
Barve N, Barve V, Jiménez-Valverde A, Lira-Noriega A, Maher SP, Peterson AT, Soberón J, Villalobos F: The crucial role of the accessible area in ecological niche modeling and species distribution modeling. Ecol Model. 2011, 222: 1810-1819. 10.1016/j.ecolmodel.2011.02.011.
Raes N, ter Steege H: A null-model for significance testing of presence-only species distribution models. Ecography. 2007, 30: 727-736. 10.1111/j.2007.0906-7590.05041.x.
Warren DL, Glor RE, Turelli M: ENMTools: A toolbox for comparative studies of environmental niche models. Ecography. 2010, 33: 607-611.
Ballard JWO, Whitlock MC: The incomplete natural history of mitochondria. Mol Ecol. 2004, 13: 729-744. 10.1046/j.1365-294X.2003.02063.x.
Zieliński P, Nadachowska-Brzyska K, Wielstra B, Szkotak R, Covaciu-Marcov S, Cogălniceanu D, Babik W: No evidence for nuclear introgression despite complete mtDNA replacement in the Carpathian newt (Lissotriton montandoni). Mol Ecol. 2013, 22: 1884-1903. 10.1111/mec.12225.
McGuire JA, Linkem CW, Koo MS, Hutchison DW, Lappin AK, Orange DI, Lemos-Espinal J, Riddle BR, Jaeger JR: Mitochondrial introgression and incomplete lineage sorting through space and time: Phylogenetics of crotaphytid lizards. Evolution. 2007, 61: 2879-2897. 10.1111/j.1558-5646.2007.00239.x.
Bryson JRW, De Oca AN-M, Jaeger JR, Riddle BR: Elucidation of cryptic diversity in a widespread Nearctic treefrog reveals episodes of mitochondrial gene capture as frogs diversified across a dynamic landscape. Evolution. 2010, 64: 2315-2330.
Liu K, Wang F, Chen W, Tu L, Min M-S, Bi K, Fu J: Rampant historical mitochondrial genome introgression between two species of green pond frogs. Pelophylax nigromaculatus and P. plancyi. BMC Evol Biol. 2010, 10: 201-10.1186/1471-2148-10-201.
Espregueira Themudo G, Arntzen JW: Newts under siege: range expansion of Triturus pygmaeus isolates populations of its sister species. Divers Distrib. 2007, 13: 580-586. 10.1111/j.1472-4642.2007.00373.x.
Arntzen JW, Espregueira Themudo G: Environmental parameters that determine species geographical range limits as a matter of time and space. J Biogeogr. 2008, 35: 1177-1186. 10.1111/j.1365-2699.2007.01875.x.
Krosby M, Rohwer S: A 2000 km genetic wake yields evidence for northern glacial refugia and hybrid zone movement in a pair of songbirds. P R Soc Lond B-Biol Sci. 2009, 276: 615-621. 10.1098/rspb.2008.1310.
Babik W, Pabijan M, Arntzen JW, Cogalniceanu D, Durka W, Radwan J: Long-term survival of a urodele amphibian despite depleted major histocompatibility complex variation. Mol Ecol. 2009, 18: 769-781. 10.1111/j.1365-294X.2008.04057.x.
Wallis GP, Arntzen JW: Mitochondrial-DNA variation in the crested newt superspecies: Limited cytoplasmic gene flow among species. Evolution. 1989, 43: 88-104. 10.2307/2409166.
Vukov TD, Sotiropoulos K, Wielstra B, Džukić G, Kalezić ML: The evolution of the adult body form of the crested newt (Triturus cristatus superspecies, Caudata, Salamandridae). J Zoological Syst Evolutionary Res. 2011, 49: 324-334. 10.1111/j.1439-0469.2011.00633.x.
Vörös J, Arntzen JW: Weak population structuring in the Danube crested newt, Triturus dobrogicus, inferred from allozymes. Amphibia-Reptilia. 2010, 31: 339-346. 10.1163/156853810791769518.
Losos JB: Phylogenetic niche conservatism, phylogenetic signal and the relationship between phylogenetic relatedness and ecological similarity among species. Ecol Lett. 2008, 11: 995-1003. 10.1111/j.1461-0248.2008.01229.x.
Hewitt G: The structure of biodiversity - insights from molecular phylogeography. Front Zool. 2004, 1: 4-10.1186/1742-9994-1-4.
Tarkhnishvili D, Gavashelishvili A, Mumladze L: Palaeoclimatic models help to understand current distribution of Caucasian forest species. Biol J Linn Soc. 2012, 105: 231-248. 10.1111/j.1095-8312.2011.01788.x.
Schmitt T, Varga Z: Extra-Mediterranean refugia: The rule and not the exception?. Front Zool. 2012, 9: 22-10.1186/1742-9994-9-22.
Waltari E, Hijmans RJ, Peterson AT, Nyári ÁS, Perkins SL, Guralnick RP: Locating Pleistocene refugia: Comparing phylogeographic and ecological niche model predictions. PLoS One. 2007, 2: e563-10.1371/journal.pone.0000563.
D. Canestrelli generously shared at the time unpublished sequence data. M.P. Miller provided advice on the use of Alleles in Space. J. van Alphen, W. Babik, W. Beukema, S. Bogaerts, S. Branković, Henrik Bringsøe, D. Canestrelli, D. Cogălniceanu, A. Crottini, G. Džukić, G. Espregueira Themudo, E. Froufe, A. Geraldes, Ö. Güclü, I. Haxhiu, D. Jovic, M. Kalezić, H.G. Kami, A. Loureiro, A. Maletzky, Đ. Milošević, N. Moin, O. Mozaffari, P. Mikulíček, B. Naumov, K. Olgun, M. Pabijan, N.A. Poyarkov, B. Safaei, J.F. Schmidtler, M. Slijepcevic, A. Urošević, N. Üzüm, J. Vörös, B. Vrenozi and A. Zuiderwijk aided with the collection of samples and locality data. The PMIP 2 Data Archive is supported by CEA, CNRS and the Programme National d’Etude de la 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. JCI was supported by a DAPTF Seed Grant 2003 and Grant 173025 of the Ministry of Education and Science of the Republic of Serbia. TV was supported by Grant 173043 of the Ministry of Education and Science of the Republic of Serbia.
The authors declare that they have no competing interests.
BW, AKS, AGT and JWA designed the study. BW, JCI, SNL, BR, KS, NT, TV and JWA collected data. BW analyzed the data. BW drafted the manuscript. All authors read and approved the final manuscript.
Electronic supplementary material
Additional file 2:Locality data. Details on the Triturus populations used for species distribution modeling. (XLS 522 KB)
Additional file 4:Genetic distances. Values of average genetic distance among populations (Zi) within each Triturus species. (XLS 110 KB)
Additional file 5:Full genetic similarity maps The genetic (dis)similarity among populations within the different Triturus species, not cut according to the current species ranges. (PDF 1 MB)
Additional file 6:Model performance. The AUC values for the species distribution models of each Triturus species, tested against null models. (XLS 42 KB)
Additional file 7:Full species distribution models. The species distribution model of each Triturus species projected for Last Glacial Maximum and current climate conditions, not cut according to the current species ranges. (PDF 2 MB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
About this article
Cite this article
Wielstra, B., Crnobrnja-Isailović, J., Litvinchuk, S.N. et al. Tracing glacial refugia of Triturus newts based on mitochondrial DNA phylogeography and species distribution modeling. Front Zool 10, 13 (2013). https://doi.org/10.1186/1742-9994-10-13
- Contact zone
- Historical biogeography
- Ice Age