Phylogeographic analysis and species distribution modelling of the wood frog Batrachyla leptopus (Batrachylidae) reveal interglacial diversification in south western Patagonia

Sample collection

Between 2009–2018 we collected 130 individuals and buccal swabs (most samples) from 19 localities throughout the distributional range of B. leptopus in south western Patagonia (Table 1; Fig. 1). Each sampling site was geo-referenced with a GPS Garmin GPSmap 76CSx. Eight individuals of B. taeniata were used as the outgroup. This study was carried out under supervision and approval of the Bioethics and Biosecurity Committee of the Universidad Austral de Chile (UACh, Resolutions No. 236/2015 and 61/15), and the Servicio Agrícola y Ganadero (SAG, Resolution No. 9244/2015).

DNA extraction, amplification, and sequence alignment

Whole genomic DNA was extracted either from liver tissues or buccal swabs according to Broquet et al. (2007), using the manufacturer’s recommended protocol for the Qiagen DNeasy tissue kit (Cat. No. 69506). We amplified three mitochondrial regions: a segment of the Control region (D-loop; Goebel, Donnelly & Atz, 1999) Cytochrome b (cyt b; Degnan & Moritz, 1992), and Cytochrome oxidase subunit I (coI; Folmer et al., 1994), and two nuclear regions: Propiomelanocortin (pomc; Gamble et al., 2008) and β Crystallin A1 (crybA1; Dolman & Phillips, 2004), via polymerase chain reaction (PCR). Reaction cocktails for PCR were according to Suárez-Villota et al. (2018). PCR products were sequenced in Macrogen Inc. (Seoul, Korea) and at the DNA Sequencing Center at Brigham Young University (Provo, USA). To transform sequence data to haplotypes we used PHASE v2.1.1 (Stephens & Donnelly, 2003) with the default model for recombination rate variation (Li & Stephens, 2003). We aligned sequences using the automatic assembly function in Sequencher v. 4.8 (Gene Codes Corp.), and inspected the aligned sequences by eye, and made corrections manually.

Molecular diversity and lineage structure

Haplotype and nucleotide diversity indexes (Nei, 1987) and their standard deviations, were estimated with DNASP v5.0 (Librado & Rozas, 2009) using all markers. The possibility of saturation in the rate of base substitutions was assessed by the method of Xia et al. (2003) using DAMBE v6.0 (Xia & Xie, 2001). Population genetic structure was examined using the package GENELAND v4 implemented in R v3.1.2 (Guillot, Mortier & Estoup, 2005), to infer the number of populations by giving a spatial model of cluster membership without prior designations. GENELAND was run with a model of uncorrelated allele frequencies for the mitochondrial locus (with all gene regions concatenated). We performed eight independent runs of 1.5 ×10⁷ iterations, with thinning set to 500 and a “burn in” of 20%. The number of possible clusters tested ranged from 1 to 19 (according to sampling locations). The level of population structure among the clusters obtained by GENELAND, was assessed by analysis of molecular variance (AMOVA; Holsinger & Weir, 2009) using ARLEQUIN v3.1 (Excoffier, Laval & Schneider, 2005) for mtDNA and nDNA separately. Also, using all loci we evaluated whether the sequences evolved under strict neutrality using Tajima’s D (Tajima, 1989), Fu & Li’s D (Fu & Li, 1993), and r² (Ramos-Onsins & Rozas, 2002) tests.

Phylogenetic trees reconstruction, split networks and divergence time estimates

Previous to phylogenetic analyses, evolutionary models and partitioning strategies were evaluated using Bayesian information criterion (BIC) scores (Schwarz, 1978) in PARTITIONFINDER v2.1.1 (Lanfear et al., 2017) (Table S1). The phylogenetic analyses were performed with the combined mitochondrial and nuclear matrix. Partitioned maximum likelihood analyses were conducted using GARLI 2.0 (Zwickl, 2006) with 200 replicates of nonparametric bootstrap for branch support. Bayesian analyses were performed using MRBAYES v3.2 (Ronquist et al., 2012) with four independent MCMC runs of 50 million generations, sampling every 2,000 generations. Posterior distributions for parameter estimates and likelihood scores to approximate convergence were visualized with the TRACER program v1.6.0 (Rambaut et al., 2014). The effective sample sizes (ESS) of each parameter (>200) allowed us to confirm that samples were adequate for all analyises. A maximum clade credibility tree was visualized with the program FIGTREE v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).

Posterior probability values >0.95 were taken as high statistical support for a clade being present on the tree (Huelsenbeck & Rannala, 2004). In order to obtain additional statistical support for the best tree obtained, topologies of different trees (ML and Bayesian) were compared with the use of the Shimodaira-Hasegawa (S-H) test (Shimodaira & Hasegawa, 1999) with resampling-estimated log likelihood (RELL), and bootstrapping of 1,000 replicates, using the program PAUP*.

We are aware that phylogenetic methods may not apply at the within-species level, due to multifurcating population genealogies in which descendant alleles coexist with ancestral ones, and recombination events may produce reticulate relationships (Posada & Crandall, 2001). To consider these caveats, we constructed unrooted phylogenetic networks using the method described by Huson & Bryant (2006), implemented in SPLITSTREE v4.14.4.

To determine when major clades and lineages diverged relative to Quaternary glaciation history, we estimated time since the most recent common ancestor (TMRCA), using the reconstructed species tree from concatenated mitochondrial and nuclear sequences. For this reconstruction, we used the multi-species coalescent module implemented in *BEAST of BEAST v1.8.4 (Drummond & Rambaut, 2007; Heled & Drummond, 2010), and the same models used for phylogenetic tree reconstruction found by PARTITIONFINDER. Because that it is not possible to date any of the nodes within B. leptopus, as there are no fossils or dated biogeographic events, we used prior Neobatrachian mutation rates of 0.291037% per million years for COI, 0.37917% per million years for each other mitochondrial markers (dloop and cytb), and a rate of 0.3741% per million years for pomc and crybA1 sequences, according to Irisarri et al. (2012). Bayes factor analysis (Li & Drummond, 2012) indicated that species trees reconstructed under a strict-clock model received decisive nodal support compared to uncorrelated exponential or uncorrelated lognormal relaxed-clock models. Markov chains in BEAST were initialized using the tree obtained by MRBAYES, to calculate posterior parameter distributions, including the tree topology and divergence times. We used BEAST to estimate divergence times with runs for 2 ×10⁷ generations, and sampling every 1000th generation. The first 10% of samples were discarded as “burn in”, and we estimated convergence to the stationary distribution and acceptable mixing (ESS >200) using TRACER v1.6.0.

Population-size dynamics through time

Hypotheses of historical demographic expansions and dynamics through geological time of the inferred lineages were estimated by “mismatch distributions” (Rogers & Harpending, 1992), and Extended Bayesian Skyline Plots (EBSP; Heled & Drummond, 2008) respectively. We used in a complementary way both approaches because small sample sizes apparently fail to provide enough power to Bayesian skyline plots to detect population expansion (Grant, 2015). The smooth, unimodal distributions typical of expanding populations can be readily distinguished from the ragged, multimodal distribution “signatures” of long-term stationary populations, by means of the ‘raggedness’ of these distributions (Rogers & Harpending, 1992). Confidence intervals for these estimates were obtained by simulations using the coalescence algorithm as implemented in DNASP v5.0. To estimate population-size dynamics through geological time, we reconstructed Extended Bayesian Skyline Plots (EBSP; Heled & Drummond, 2008) implemented in BEAST for the four lineages obtained with GENELAND. This coalescent-based, nonparametric Bayesian MCMC algorithm incorporates multi-locus data to reduce error estimates associated with single genes (e.g., traditional Bayesian Skyline Plots), and increases the power to resolve alternative demographic histories (Ho & Shapiro, 2011). For each EBSP, the appropriate model of nucleotide substitution was determined using PARTITIONFINDER. Genealogies and model parameters for each lineage were sampled every 1000th iteration for 2 ×10⁷ generations under a strict molecular clock with uniformly distributed priors, and a “burn in” of 2000. Demographic patterns for each analysis were plotted in EXCEL v14.7.7.

Test of phylogeographic hypotheses with ABC

A coalescent method was used to test phylogeographic hypotheses by constraining the genealogies to fit alternative evolutionary models, and assessing each model’s fit by comparing the observed genetic pattern with the range of simulated patterns. Competing phylogeographic hypotheses were compared using an approximate Bayesian computation method (ABC approach), as implemented in DIYABC v2.1 (Cornuet et al., 2014). We evaluated five demographic scenarios to test alternative divergence times and tree topologies of the four main lineages recovered in GENELAND and in the phylogenetic analyses. Furthermore, the divergence of the main four lineages occurred before the last glacial maximum so the divergence scenarios were proved in such range. The refugia hypotheses correspond to the points of coalescence, with four lineages strongly supported, the possibilities of coalescence are from 1 to 3, for which we tested the possibilities except scenarios with postglacial admixture, because there are no- paraphyly events between lineages. The prior coalescence points at time times t1 and t2 applied in the ABC correspond to those estimated by BEAST for the origin of Batrachyla leptopus (t2), and the divergence of the four lineages (t1, lower and higher range of the four lineages). Thus, all historically relevant scenarios differed only in the order of population divergence, and in the number and timing of demographic expansion events. These alternatives were: Scenario 1—the null model—all four lineages coalesced at t1 with equal divergence rates. Scenario 2—also a null model, but all four lineages coalesced at t2 with equal divergence rates. Scenario 3—the first coalescence of lineages A and B at t1, whose ancestor coalesced at t2 with lineage C and D. Scenario 4—the first coalescence of lineages C and D at t1, whose ancestor coalesced at t2 with lineage A and B. Scenario 5—one split event at t1 isolated the north (lineages A and B) from the south (lineages C and D) clades, and then a coalescence of both clades at t2 (see Results). We tested other scenarios whose divergence time were older and the probability was very low so they were not considered.

Prior values of Ne were set as 1,000–500,000 individuals with a uniform distribution, based on Ne calculated from MIGRATE-N v3.6 (Beerli, 2006). We performed maximum likelihood using 10 short chains of 1,000 steps, and two long chains of 10 000 steps, sampling each 100 steps, and a burn-in of 10 per cent. Ne was calculated using mitochondrial and nuclear rates reported by Irisarri et al. (2012).

Prior values for divergence of the ancestral populations were based on divergence times calculated here (see Results), and a generation time of 2–3 yr (Martin & Palumbi, 1993), using a uniform distribution. Divergence times were set at between 20,000–500,000 generations ago for t2, and 10,000–200,000 generations for t1.

Paleo-distribution and species distribution modelling

A total of 120 occurrence records were used for the species distribution models (SDMs) and paleo-distribution modelling. Records were obtained from peer-reviewed literature, our sampled sites, and online databases (GBIF: gbif.org, VetNet: vertnet.org, and iDigBio: idigbio.org). We modeled the SDMs using the standard 19 bioclimatic variables downloaded from Worldclim (Hijmans et al., 2005) for the current conditions (1960–1990), the Mid Holocene (Mid-Hol, ∼6,000 yrs BP), the LGM (∼22,000 yrs BP), and the Last Inter-glacial (LIG, ∼120,000–140,000 yrs BP). The variables were at 30 arc-sec for the current conditions and the LIG, and at 2.5 arc-min for the Mid-Hol and LGM. The Mid-Hol and the LGM variables were based on the Community Climate System Model (CCSM) and the Model for Interdisciplinary Research on Climate (MIROC), and while LIG conditions were based on Otto-Bliesner et al. (2016). We restricted the projection of the models by creating a buffer of 2° around the outermost occurrence records and the known distribution of B. leptopus. All SDMs were performed using MAXENT v3.4.0 (Phillips et al., 2017). To avoid model overfitting and account for the correlation between the variables and the presence of outliers encountered during data exploration, we reduced the number of variables to five. This was performed by retaining the variables with —rho— < 0.8 that contributed the most to ten cross-validated models, as shown by a Jackknife test. These variables were Bio2 = Mean Diurnal Range, Bio4 = Temperature Seasonality, Bio5 = Max Temperature of Warmest Month, Bio13 = Precipitation of Wettest Month, and Bio17 = Precipitation of Driest Quarter.

For all models, the equal training and sensitivity threshold rule was applied and the cloglog output was selected. Extrapolation was not used and clamping was applied when hind-casting the model of the current conditions to the past. The models for the Mid-Hol and for the LGM were overlaid within each time period to identify areas of agreement and disagreement between the models for each the time period. All models were transformed to binary using the selected threshold rule (Phillips et al., 2017).

Genetic structure using Bayesian clustering model

A total of 113 haplotypes were found when mitochondrial and nuclear markers were combined (Table 2). Saturation tests as a function of the genetic distance estimated under substitution model GTR showed non- significative saturation of DNA sequence alignments. The mtDNA Bayesian analysis with GENELAND yielded a modal number of four clusters (K = 4), recovered from all independents runs (Fig. 2A); this is based on the highest average posterior probability. The distribution of these four clusters, from north to south, was named as follows (Fig. 2B): lineage A (Los Queules), lineage B (Nahuelbuta), lineage C (Bahía Mansa, Cordillera Pelada, Máfil, and Pichirropulli), and lineage D (all remaining localities shown in Table 1). The highest estimate of haplotype diversity was found in lineage A, whereas the lowest values were found in lineage B (Table 2). Highest nucleotide diversity was detected in lineage D and the lowest one in lineage B (Table 2). Negative but non-significant values for Tajima’s D and Fu’s FS neutrality test indexes were found in all lineages except for B lineage. Here the non-significant positive values likely reflected the smallest sample size for this locality. Rozas’s r² were positive for all lineages but only significant for A and B lineages, suggesting recent expansion (Table 2).

The AMOVA results using the four lineages indicate a significant genetic structure (groups defined as Lineages A, B, C, and D): (1) variation among lineages = 35.64% and variation within lineages = 45.93% for mtDNA; (2) variation among lineages = 0% and variation within lineages = 51.12% for pomc; (3) variation among lineages = 4.49% and variation within lineages = 2.74% for crybA1 (Table 3).

Phylogenetic tree reconstruction, split networks, and lineage divergence time

The models selected for the ML and Bayesian analyses are described in Table S1. Because the Bayesian analyses recovered a maximum clade credibility tree similar to the best ML tree, and the Shimodaira-Hasegawa test showed that topological disagreements were restricted to “low-support” nodes, we show only the Bayesian tree (Fig. 3A). The same four lineages recovered by GENELAND were recovered in the phylogenetic reconstruction with moderate to high support values. Lineages A and B were recovered as sister groups (bootstrap = 96%; BPP = 0.99), as were the C and D lineages, but with moderate support (bootstrap 74%; BPP= 0.95). Mitochondrial phylogenetic analyses recovered similar results to concatenated datasets, but the nuclear phylogeny was highly polytomized (Fig. S1). The split networks (Fig. 3B) recovered the same four lineages obtained by the ML and Bayesian analyses. The fit index was 94.96, meaning that only 5.04% of the distances in the distance matrix are not represented by the network. Most of the internal splits have bootstrap support between 68 and 100%.

Divergence dating indicates that the A–B and C–D clades separated during the late Pleistocene, approximately 0.107 mya [95% confidence interval (CI) = 0.020–0.278 mya]. The divergence between A and B lineages would have occurred by the late Pleistocene (approximately 0.068 mya; 95% CI [0.036–0.147] mya) and divergence between C and D lineages was approximately 0.065 mya (95% CI [0.056–0.092] mya) (Fig. 3A).

Demographic patterns of the inferred clusters

Results of mismatch distribution analyses (Figs. 4A–4D) revealed a single primary peak for lineage A, but with a non-significant raggedness index (r = 0.0230, P > 0.1). Similarly, unimodal patterns were observed in lineages C (r = 0.003, P < 0.001) and D (r = 0.0008, P < 0.001). The small sample size (n = 7) and haplotype numbers (H = 2) precluded this analysis on lineage B. Reconstruction of the demographic histories by means of Extended Bayesian Skyline Plot (Figs. 4E–4H) suggested population expansions for all lineages except B (Fig. 4F). EBSP further resolved sequential demographic expansions from the oldest, lineage D (c. 18,000 years bp), then lineage A (c. 11,000 years bp), and most recently, lineage C (c. 5,000 years bp).

Hypothesis testing with ABC

Logistic regression analysis with DIYABC identified Scenario 4 as most strongly supported among the five tested (Fig. 5D), with a high posterior probability (0.93; Table 4); all other scenarios had much lower support (0.0–0.04). Moreover, the Type I and Type II error rates estimated for Scenario 4 were the lowest in both cases (0.14; 0.014–0.196; Table 4).

Scenario 4 placed the first divergence as the split between lineages A, B, and the ancestor of the southern clade (lineages C and D) at t2, and the second split between lineages C and D at t1 (Fig. 5D). The effective population size (Ne) and divergence time parameters, in terms of the number of generations (t), estimated for this divergence scenario (Table S2), corroborate the population expansions inferred by EBSP and mismatch distribution for lineages C and D.

SDMs and paleo-distribution models

The predicted distribution models of B. leptopus under four periods (last inter-glacial to current) are shown in Fig. 6. The model for the current conditions showed that the distribution of this species is mostly encompassed by the county-based current distribution mapped by the IUCN, with a high AUC value (0.961, SD = 0.013) (Fig. 6A). The two circulation models for the Mid-Hol showed different distributions: the MIROC-based model showed a continuous distribution similar to the current conditions model, while the CCSM-based model resolved disjunct distributions for the northern region of the current distribution model (Fig. 6B). Both LGM models showed clear disjunct distributions with concordance around Nahuelbuta mountain range (Fig. 6C). The MIROC based model indicated a distribution to the north of the current distribution model, and the CCSM resolved a highly fragmented distribution. The model for the LIG showed a disjunct distribution for the species at the northern portion of the current distribution of B. leptopus (Fig. 6D).

Biogeographic structure of the lineages

The two northernmost lineages (A and B) are currently separated by two major river systems (Itata and Bío Bío; Fig. 1); these boundaries also coincide with a region with sparsely distributed forest patches. These two lineages are genetically well-differentiated from each other (Fig. 2B); lineage A is restricted to the Los Queules Reserve, and lineage B is limited to a small area in Nahuelbuta range (Locality 2) near the Butamalal River. The evolution of distinct lineages or genetic clusters is often attributed to population isolation during glacial advances, which in combination with different selection pressures and/or genetic drift would drive population divergence (Hewitt, 2004). Similarly, low levels of current genetic diversity as observed in lineage B (Nahuelbuta; Table 2), suggest a recolonization history of founder effects, small population sizes, and genetic bottlenecks (Hewitt, 2004).

In contrast to these well-differentiated/low variability lineages, those in the southern distribution (lineages C and D) are geographically more heterogeneous. Lineage C is widely distributed from Máfil in the Los Ríos region (Locality 3 in Fig. 1) to Bahía Mansa, Los Lagos region (Locality 6 in Fig. 1), suggesting that certain landscape features, such as extensive forests, would have allowed dispersal among breeding groups over long timescales. Moreover, the encompassed areas of lineage C also include large mountains (e.g., Bahía Mansa and Cordillera Pelada; Fig. 1), suggesting that such orography could represent barriers to gene flow in B. leptopus, as reported in some co-distributed vertebrates and plants (Sérsic et al., 2011). The combined dataset suggests that lineage D is widespread throughout the rest of the species’ range, with few phylogeographical subdivisions.

Lineage divergence time

Divergence time estimates suggest that diversification of B. leptopus lineages may have occurred earlier than reported in other frogs such as the ground frog Eupsophus calcaratus (Nuñez et al., 2011), although co-distributed populations (e.g., Bahía Mansa) appear to have diverged later in time (0.025 mya for B. leptopus (Fig. 3A), and 0.065 mya for E. calcaratus (Nuñez et al., 2011; Fig. 2). Vidal et al. (2016) point out that some B. leptopus populations from Chiloé Island and the mainland (included in Lineage D in our study) diverged approximately 1.1 mya. Our results indicated that the initial split of B. leptopus into the North and South clades was during the Late Pleistocene (∼0.107 mya; Fig. 3). Moreover, the overall pattern suggests that B. leptopus has undergone several rounds of fragmentation, followed by successive radiations within each clade. Further, at least two more recent series of fragmentation events are inferred within each of these clades. Our calibrations place the split between lineages A and B at ∼0.068 mya, and between lineages C and D at ∼0.065 mya.

The discrepancy in divergence times between our results and those of the Vidal et al. (2016) may be due to the use of different mutations rates (0.8%) and a single marker (mitochondrial cyt b). It is well known from population genetics theory that the stochastic nature of the genealogical process implies a significant amount of variance associated with parameter estimation. In fact, Nabholz, Glémin & Galtier (2009) suggest that divergence inferences date should be based on statistical phylogenetic methods accounting for substitution rate variation across lineages.

Indeed, analysis of mtDNA sequence data can be enhanced in conjunction with nuclear sequences, which provide an independent estimate of phylogenetic relationships, mitigating the inherent stochasticity of genetic drift, and the variance associated with parameter estimation (Carstens & Dewey, 2010). Interpretations derived of the divergence time analyses also need to take into account the largely overlapping confidence intervals of the results for each lineage divergence. For example, if only the median values are considered, the results suggest coalescence of the lineages C and D at t1 (around 0.066 mya), and coalescence of the lineages A and B at t2 (around 0.107 mya). But if we consider the confidence intervals, the divergence of both clades, North and South could be closer in time than what is hypothesized in our reconstructions.

Hypothesized refugia and post-glacial expansions

Studies on past contraction-expansion climate cycles in Patagonian landscapes suggest that rapid population expansions should occur in the biota affected by these processes, as habitats became more available (Fraser et al., 2012). Despite the geomorphological differences in Patagonian landscapes, population genetics theory predicts that a population undergoing rapid expansion should be characterized by low genetic diversity, since each new founder population represents only a fraction of the ancestral population (Nichols & Hewitt, 1994; Hewitt, 2000; Hewitt, 2004; Waters, Fraser & Hewitt, 2013).

The last two Pleistocene glaciations in southwestern South America (180 kya and 20 kya) covered the Andes with large ice fields reaching the Pacific Ocean south to 39°S, where the ice sheet decreased in elevation to sea level, and extending further to the southern tip of South America (Rabassa, 2011). Consequently, Late-Pleistocene divergence time estimates (0.107 mya) for first diversification of B. leptopus populations separating the North and South clades (Fig. 3), are consistent with a hypothesis of Pleistocene isolation followed by interglacial dispersal. In fact, the divergence of lineages A and B suggests that this interglacial dispersal from the ancestral population occurred rapidly across the current range of the species.

The existence of two suitable areas for the species is supported by the SDM for the LIG (0.120–0.140 mya) in that Los Queules population (Locality 1, Table 1) showed the highest genetic diversity (Table 2). This is typical for refugial populations that have been stable over time (Fraser et al., 2012). This evidence and the agreement of the two circulation models for Los Queules area as a suitable habitat for the species during the LGM (Fig. 6C), suggests that it is highly probable that the present Los Queules location is a remnant of the northern refuge, derived from the last southern Patagonian glaciation (180 kya).

Demographic reconstruction in B. leptopus using an ABC framework also supports the hypothesis of two putative refugia at different times during the Pleistocene (Scenario 4, Fig. 5D, Table 4). This scenario suggests that the South clade populations (lineages C and D) are likely descended from a divergence event approximately 65 Kya. This scenario is concordant with the predicted patchy distribution of the species during the Mid-Holocene (∼6,000 years BP; Fig. 6B). On the other hand, the various demographic analyses show that the genetic structure of lineages C and D contains signatures of demographic expansion consistent with Pleistocene glacial retreat. In fact, strong support for recent population expansion is represented by significantly negative Fu’s Fs values, unimodal mismatch distributions with low raggedness indexes, and EBSP’s depicting rapid expansion following the retreat of the Patagonian ice sheet after 15 Kya (Table 2; Fig. 4). This hypothesis is also reinforced with SDMs produced by both circulation models, although when viewed separately, the MIROC circulation model is the only one predicting this scenario.

Past population dynamics and current conservation significance

During the Cenozoic climatic oscilations included multiple glaciations in southern South America (Rabassa & Clapperton, 1990; Coronato, Martínez & Rabassa, 2004). These geological events have been hypothesized as causes of the retreat and advance of temperate Nothofagus forests and conifers (Villagrán & Hinojosa, 1997; Premoli, Kitzberger & Veblen, 2000; Tremetsberger et al., 2009). Accordingly, several phylogeographic hypotheses suggest that Pleistocene glaciations had profound effects on the population genetic structure and variability of the Patagonian fauna. For example, glaciated populations of some fish species display molecular diversity (high haplotype diversity and low nucleotide diversity) significantly correlated with latitude (Ruzzante et al., 2008; Cosacov et al., 2010). The same genetic patterns have been observed in reptiles (Breitman et al., 2011; Fontanella et al., 2012), amphibians (Nuñez et al., 2011), and mammals (Himes, Gallardo & Kenagy, 2008; Lessa, D’Elía & Pardiñas, 2010).

In agreement with these previous studies, the location of glacial refugia and postglacial expansions identified here, indicate that the climatic niche of B. leptopus is likely to be related to an increase in the availability of suitable habitat in the southern part of its current distribution. Reconstruction of the potential distribution area of B. leptopus (Fig. 6) suggests that suitable habitats underwent expansions and contractions during glacial retreats and advances.

In addition to the inferences of past population dynamics, predictions about the ability of species to respond to future climate change play an important role in alerting potential risks to biodiversity. In fact, many studies have investigated the response of biodiversity to climate change, and most of them indicate that current and future rates of these changes may be too fast for ecological niches to evolve (Fraser et al., 2012; Rolland et al., 2018). This is particularly critical in species with low dispersal capacity such as amphibians, which makes them potentially less able to respond to changes induced by climate and, consequently, more vulnerable to extinction (Duan et al., 2016). In this context, the genetic impoverishment at northern area of the distributional range of B. leptopus is of great concern, given a climate change scenario based on increases in temperatures and aridity in central-southern Chile.