# Multiple glacial refugia and contemporary dispersal shape the genetic structure of an endemic amphibian from the Pyrenees

Federica Lucati, Manon Poignet, Alexandre Miró, Audrey Trochet, Fabien Aubret, Laurent Barthe, Romain R. Bertrand, Teresa Buchaca, Olivier Calvez, Jenny Caner, et al.

► To cite this version:

Federica Lucati, Manon Poignet, Alexandre Miró, Audrey Trochet, Fabien Aubret, et al. Multiple glacial refugia and contemporary dispersal shape the genetic structure of an endemic amphibian from the Pyrenees. Molecular Ecology, 2020, 29 (15), pp.2904-2921. 10.1111/mec.15521. hal-02991153

HAL Id: hal-02991153

https://ut3-toulouseinp.hal.science/hal-02991153v1

Submitted on 7 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
MS. FEDERICA LUCATI (Orcid ID: 0000-0001-9925-7137)

DR. MATHIEU DENOEL (Orcid ID: 0000-0002-3586-8323)

Article type : Original Article

# Multiple glacial refugia and contemporary dispersal shape the genetic structure of an endemic amphibian from the Pyrenees

Running title: Linking phylogeography and recent dispersal

Federica Lucati*<fnref n="1" />,<fnref n="2" /> , Manon Poignet*<fnref n="3" />, Alexandre Miró<fnref n="2" />, Audrey Trochet<fnref n="3" />,<fnref n="4" />, Fabien Aubret<fnref n="3" />, Laurent Barthe<fnref n="5" />, Romain Bertrand<fnref n="3" />, Teresa Buchaca<fnref n="2" />, Olivier Calvez<fnref n="3" />, Jenny Caner<fnref n="2" />, Elodie Darnet<fnref n="3" />, Mathieu Denoël⁶, Olivier Guillaume<fnref n="3" />, Hugo Le Chevalier<fnref n="3" />, Albert Martínez-Silvestre⁷, Marc Mossoll-Torres⁸,⁹, David O’Brien¹⁰, Víctor Osorio<fnref n="2" />, Gilles Pottier<fnref n="5" />, Murielle Richard<fnref n="3" />, Ibor Sabás<fnref n="2" />, Jérémie Souchet<fnref n="3" />, Jan Tomàs<fnref n="2" /> and Marc Ventura<fnref n="2" />

*Both authors contributed equally to this work

This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process, which may lead to differences between this version and the Version of Record. Please cite this article as doi: 10.1111/MEC.15521

This article is protected by copyright. All rights reserved

<footnote n="1">¹ Centre for Ecology, Evolution and Environmental Changes (cE3c), Faculty of Sciences, University of Lisbon, Lisbon, Portugal</footnote>

<footnote n="2">²Center for Advanced Studies of Blanes (CEAB-CSIC), Blanes, Spain</footnote>

<footnote n="3">³ CNRS, Station d'Ecologie Théorique et Expérimentale (SETE), Université Paul Sabatier, Moulis, France</footnote>

<footnote n="4">⁴ Société Herpétologique de France, Muséum National d'Histoire Naturelle, Paris, France</footnote>

<footnote n="5">⁵ Association Nature En Occitanie, Maison de l'Environnement de Midi-Pyrénées, Toulouse, France</footnote>
**Abstract**

Historical factors (colonization scenarios, demographic oscillations) and contemporary processes (population connectivity, current population size) largely contribute to shaping species' present-day genetic diversity and structure. In this study, we use a combination of mitochondrial and nuclear DNA markers to understand the role of Quaternary climatic oscillations and present-day gene flow dynamics in determining the genetic diversity and structure of the newt *Calotriton asper* (Al. Dugès, 1852), endemic to the Pyrenees. Mitochondrial DNA did not show a clear phylogeographic pattern and presented low levels of variation. In contrast, microsatellites revealed five major genetic lineages with admixture patterns at their boundaries. Approximate Bayesian computation analyses and linear models indicated that the five lineages likely underwent separate evolutionary histories and can be tracked back to distinct glacial refugia. Lineage differentiation started around the Last Glacial Maximum at three focal areas (western, central and eastern Pyrenees) and extended through the end of the Last Glacial Period in the central Pyrenees, where it led to the formation of two more lineages. Our data revealed no evidence of recent dispersal between lineages, whereas borders likely represent zones of secondary contact following expansion from multiple refugia. Finally, we did not find genetic evidence of sex-biased dispersal. This work highlights the importance of integrating past evolutionary processes and present-day gene flow and dispersal dynamics, together with multilocus approaches, to gain insights into what shaped the current genetic attributes of amphibians living in montane habitats.

**Keywords:** *Calotriton*, genetic structure, phylogeographic history, Pyrenean brook newt, recent dispersal, Pyrenees

<footnote n="6">⁶ Laboratory of Ecology and Conservation of Amphibians (LECA), Freshwater and Oceanic science Unit of reSearch (FOCUS), University of Liege, Liege, Belgium</footnote>

<footnote n="7">⁷ Catalonia Reptile and Amphibian Rescue Center (CRARC), Masquefa, Spain</footnote>

<footnote n="8">⁸ Pirenalia, Encamp, Andorra</footnote>

<footnote n="9">⁹ Bomosa, Les Escaldes, Andorra</footnote>

<footnote n="10">¹⁰ Scottish Natural Heritage, Scotland, UK</footnote>

<footnote>□ Corresponding author. Email: federalucati@hotmail.com</footnote>
## Introduction

Unveiling the mechanisms driving species genetic diversity and structure is of crucial interest in phylogeography (Avise, 2000). The extent of genetic structure of a species is regarded to result primarily from the interplay of historical factors (e.g. colonization scenarios, demographic oscillations) and current population connectivity, namely gene flow (Hewitt & Butlin, 1997; Nichols & Beaumont, 1996). Unravelling the phylogeographic history of species and populations is important to understand their present-day and future distribution, genetic structure and adaptations (Hewitt, 2004). Historical processes are largely dependent on past climatic conditions and geological events. Such climatic and geological changes have significantly contributed to laying the genetic foundations of contemporary populations, which can be used to make inferences on their past dynamics (Cabrera & Palsbøll, 2017; Hewitt & Butlin, 1997). In addition, dispersal, which can include gene flow, is a significant component of metapopulation structure and dynamics and can counteract both neutral and selective processes (Johnson & Gaines, 1990; Ronce, 2007; Tallmon, Luikart, & Waples, 2004). A reduction in connectivity will ultimately result in a lack of dispersal among populations, increasing the risk of genetic variability loss (Ronce, 2007; but see Orsini, Vanoverbeke, Swillen, Mergeay, & De Meester, 2013). For this reason, dispersal is deemed crucial for the long-term survival of populations under changing conditions (Saccheri et al., 1998). In some circumstances, other processes might explain the genetic variability of populations, such as isolation by environment (reduction in gene flow among ecologically divergent habitats as a result of local adaptation) and by colonisation (reduction in gene flow among all populations in the landscape caused by local genetic adaptation following colonisation; Orsini et al., 2013). The literature on historic vs. contemporary mechanisms shaping the genetic attributes of species is mostly focused on either landscape genetics or dispersal processes alone, or tackle temporal dynamics dealing with the relatively recent past (Chiucchi & Gibbs, 2010; Epps & Keyghobadi, 2015; Noguerales, Cordero, & Ortego, 2017; Zellmer & Knowles, 2009). An integrative approach that combines the study of past evolutionary and phylogeographic processes and present-day gene flow and dispersal dynamics is required to shed light on the mechanisms underlying spatial patterns of contemporary genetic diversity and population structure, which can ultimately help to predict their responses to ongoing or future environmental changes.

In Europe, Quaternary climatic oscillations played a major role in shaping the geographic distribution and genetic constitution of species (Hewitt, 2000, 2004). Glacial and interglacial periods caused repeated changes in species' distributions, leading to events of contraction and
expansion and, consequently, to periodic waves of colonization or recolonization. Mountain ranges across Europe are regarded as biodiversity cradles, where diversification is promoted during periods when species' ranges are restricted to geographically isolated glacial refugia (Hewitt, 2000; Schmitt, 2009). As glaciers repeatedly advance and retreat, species are displaced outside or to the margin of mountain systems into lowland and peripheral areas, respectively, or survive in nunataks, namely areas above glaciers not covered with ice (Holderegger & Thiel-Egenter, 2009). Mountain ecosystems are home to many endemisms that still carry genetic imprints of these past dynamics, and thus represent excellent models with which to study the influence of climatic fluctuations on the diversification and postglacial colonization of species (Schmitt, 2009).

As one of the major European mountain ranges and separating the Iberian Peninsula from the rest of continental Europe, the Pyrenees played a considerable role in limiting postglacial dispersal routes of numerous temperate species (Taberlet, Fumagalli, Wust-Saucy, & Cosson, 1998). During glacial periods, the Pyrenees were largely covered with ice (Calvet, 2004; González-Sampériz et al., 2006). Nevertheless, it is suggested that some species could have survived glaciations in ice-free areas along the chain, such as nunataks and peripheral lower areas that served as glacial refugia (Bidegaray-Batista et al., 2016; Charrier, Dupont, Pornon, & Escaravage, 2014; Liberal, Burrus, Suchet, Thebaud, & Vargas, 2014; Mouret et al., 2011). Following the end of glacial periods, deglaciation allowed recolonization along routes spreading from these refugia and this ultimately sculptured a complex genetic structure in the Pyrenees (Hewitt, 1999; Taberlet et al., 1998). However, there has been little attempt to identify the geographic location of putative refugia where Pyrenean endemics survived glaciations, and to trace back their postglacial recolonization routes.

Dispersal capability is a crucial trait affecting the genetic composition of species and populations (Clobert, Le Galliard, Cote, Meylan, & Massot, 2009; Ronce, 2007; Tallmon et al., 2004), implying that variation in vagility generally leads to clear differences in genetic patterns. Good dispersers are likely to present less structured metapopulations than low vagility organisms (Allentoft, Siegismund, Briggs, & Andersen, 2009; Burns, Eldridge, & Houlden, 2004; Kraaijeveld-Smit, Beebee, Griffiths, Moore, & Schley, 2005; Vos, Antonisse-De Jong, Goedhart, & Smulders, 2001). Amphibians are generally regarded as low vagility and philopatric species (Gill, 1978) but this is being confuted in a number of studies (Denoël, Dalleur, Langrand, Besnard, & Cayuela, 2018; Smith & Green, 2005, 2006). Selective pressures favouring or restraining
dispersal may act differently on males and females and result in sex-specific dispersal strategies (Li & Kokko, 2019). Accordingly, sex-biased dispersal has been identified in a number of species, including newts (Denoël et al., 2018; Trochet et al., 2016). Furthermore, orographic features such as ridges and valleys can act as either barriers or bridges to dispersal and thus drive genetic structuring (Caplat et al., 2016; Noguerales, Cordero, & Ortego, 2016). Although it is deemed important to better understand the processes underlying genetic differentiation in natural populations, the combined influence of sex differences and orographic features on dispersal has rarely been studied (Roffler et al., 2014; Tucker, Allendorf, Truex, & Schwartz, 2017). Indeed, males and females may have different dispersal abilities and therefore orographic features may differently affect them, resulting in contrasting patterns of gene flow between sexes in mountain regions (see Cayuela et al., 2020 for a review).

The genus *Calotriton* (Gray, 1858) includes two species restricted to north-eastern Iberian Peninsula (Carranza & Amat, 2005). Speciation within the genus has been dated to the beginning of the Pleistocene (Carranza & Amat, 2005) but how these species endured Quaternary glaciations is still uncertain. The Pyrenean brook newt (*C. asper* Al. Dugès, 1852) is a small-bodied amphibian endemic to the Pyrenees (Bosch et al., 2009). It is a largely aquatic montane species that inhabits brooks, alpine lakes and caves between 250 and 2,500 m a.s.l. (Clergue-Gazeau & Martínez-Rica, 1978; Martínez-Rica & Clergue-Gazeau, 1977). As expected for many amphibian species, *C. asper* is believed to have low dispersal ability (Milá, Carranza, Guillaume, & Clobert, 2010; Montori, Llorente, & Richter-Boix, 2008), although little attention has been paid to this aspect. Following metamorphosis, a juvenile dispersal phase of at least 2 years is described before reaching the adult stage (Montori & Llorente, 2014), but it remains unclear how far individuals can disperse.

So far, few studies have analysed the genetic differentiation of the Pyrenean brook newt in a geographic context. Analysis of allozymes (Montori, Llorente, & García-París, 2008) and mitochondrial DNA (mtDNA; Milá et al., 2010; Valbuena-Ureña, Amat, & Carranza, 2013) revealed low levels of genetic variation. Higher levels of genetic differentiation and population structuring were detected using genome-wide amplified fragment length polymorphism (AFLP; Milá et al., 2010) and microsatellite markers (Valbuena-Ureña et al., 2018). However, these studies were either based on small numbers of populations and markers with low variability (Montori, Llorente, & García-París, 2008; Valbuena-Ureña et al., 2013), did not characterize the entire range and habitat types of the species (Milá et al., 2010), or addressed specific questions
targeting the role of geographic gradients and habitat type in shaping the current genetic attributes
of the species (Valbuena-Ureña et al., 2018). Furthermore, the timing of lineage divergence and
the relative importance of phylogeographic processes versus contemporary dispersal have not been
studied in *C. asper*.

Here, we employ a multilocus approach aimed to disentangle major historical and
contemporary processes that contributed to shaping the present genetic constitution of *C. asper*
over most of its distribution range. We combine comprehensive sample collection across all
habitat types with coalescent model frameworks and dispersal analyses to shed light on the
evolutionary history of the species and determine the degree of connectivity of present-day
populations and habitats. Specifically, we explore the effect of Quaternary climatic oscillations on
the evolutionary diversification of lineages and the formation of postglacial colonization routes.
Furthermore, we describe contemporary patterns of dispersal and investigate whether sex-specific
dispersal strategies, orography or geography played a role in determining the species' current
genetic structure.

**Materials and Methods**

*Sampling and DNA extraction*

Sampling was conducted in the period 2004-2017 across the whole Pyrenees (Figure 1; Table S1), encompassing most of the species range. DNA was sampled via buccal swab or toe clipping of metamorphosed individuals. Samples were preserved in EDTA or absolute ethanol and stored at -20°C until DNA extraction. The collection of samples was approved by the corresponding authorities: as for the French sampling, by the Conseil Scientifique Régional du Patrimoine Naturel (CSRPN, DREAL) of the Region of Occitanie; as for the Andorran sampling, by the Principality of Andorra; as for the Spanish sampling, by the Departament d'Agricultura, Ramaderia, Pesca, Alimentació i Medi Natural of the Catalan Government and the Instituto Aragonés de Gestión Ambiental of the Aragonese Government. Procedures followed guidelines established by the Association for the Study of Animal Behaviour and complied with current French, Andorran and Spanish regulations.

Genomic DNA was extracted using QIAGEN DNeasy Blood and Tissue Kit (Qiagen™, Hilden, Germany) according to the manufacturer's protocol, or following the HotSHOT method (Montero-Pau, Gómez, & Muñoz, 2008), in a total volume of 100 µl.
## Mitochondrial DNA sequencing and microsatellite screening

A fragment of the cytochrome *b* (cyt-*b*) gene was sequenced from 258 individuals from 59 sampling sites (Table S1). We amplified a fragment of 374 bp using primers Cytb1EuprF and Cytb2EuprR (Carranza & Amat, 2005). Amplification conditions were those described in Carranza, Arnold, Mateo, and López-Jurado (2000). Sequences were aligned using the ClustalW algorithm in MEGA 7 (Kumar, Stecher, & Tamura, 2016).

A total of 1,299 individuals from 96 sampling sites were genotyped for a set of 17 microsatellite loci combined in three multiplexes (Table S1; Drechsler et al., 2013). Fragments were sized with LIZ-500 size standard and binned using either GeneMapper v4.0 (Applied Biosystems) or Geneious 11.0.5 (Kearse et al., 2012). Only individuals that could be scored in a reliable manner for at least 15 loci were included in the analyses.

## Mitochondrial DNA analysis

Gene genealogy networks were generated using Haploviewer (Salzburger, Ewing, & Von Haeseler, 2011). jModelTest 2.1.3 (Darriba, Taboada, Doallo, & Posada, 2012) was run to determine the appropriate nucleotide-substitution model, under the Akaike Information Criterion (AIC). Phylogenetic reconstructions among haplotypes were estimated using a maximum likelihood approach as implemented in RAxML 7.7.1 (Stamatakis, 2006), and the best generated tree was used to estimate the haplotype network. The program was run with a GTRCAT model of rate heterogeneity and no invariant sites, applying 1,000 bootstrap replicates. Haplotype network reconstruction was implemented in Haploviewer, based on all sequences available from GenBank and this study. Overall number of haplotypes (H) and polymorphic sites (S), as well as haplotype (Hd) and nucleotide (Π) diversity indices were calculated in DNASP 6.11.01 (Rozas et al., 2017).

## Microsatellite analysis

The presence of potential scoring errors, stuttering, large allele dropout and null alleles was tested using MICRO-CHECKER 2.2.3 (Van Oosterhout, Hutchinson, Wills, & Shipley, 2004). The frequency of null alleles for each locus and population was further investigated using the expectation maximization algorithm implemented in FreeNA (Chapuis & Estoup, 2006). The same program was used to calculate global $F_{ST}$ values corrected for null alleles following the Excluding Null Alleles (ENA) correction method. Bootstrap 95% confidence intervals (CI) were calculated using 1,000 replicates over loci. We tested for linkage disequilibrium between loci and for
deviations from Hardy-Weinberg equilibrium (HWE) in each population and for each locus in
GENEPOP 4.2 (Rousset, 2008). Significance levels for multiple comparisons were adjusted using
the Bonferroni correction (α = 0.05; Rice, 1989).

Parameters of genetic diversity were estimated for populations with five or more genotyped
individuals and for the genetic clusters inferred by STRUCTURE. Calculation of diversity
estimates only in populations with larger sample size (≥ 10 individuals) yielded very similar
results in terms of mean genetic diversity and in the spatial interpolation analysis. We calculated
observed (H<sub>O</sub>) and expected heterozygosity (H<sub>E</sub>) using the PopGenKit R package (Rioux Paquette,
2011) in R 3.5.1 (R Core Team, 2018). Allelic richness (Ar) standardized for sample size and
rarefied private allelic richness (PAAr; calculated only at the cluster level) were calculated in HP-
RARE 1.1 (Kalinowski, 2005). Inbreeding coefficients (F<sub>IS</sub>) were estimated in FSTAT 2.9.3.2
(Goudet, 2002). We visualised geographic patterns of genetic diversity by computing a spatial
interpolation of H<sub>E</sub> and Ar values using the Inverse Distance Weighting tool implemented in
ArcGIS 10.1 (ESRI, Redlands, CA, USA).

Population structure was investigated using a Bayesian approach implemented in STRUCTURE 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). We conducted 20 independent simulations for each K value from one to 50, with 100K burn-in steps followed by 500K Markov chain Monte Carlo (MCMC) repetitions. It is highly unlikely that *C. asper* would reveal more than 50 genetic units, given that previous studies conducted at the Pyrenean scale returned a much smaller number of nuclear partitions (Milá et al., 2010; Valbuena-Ureña et al., 2018). The program was run using the admixture model with correlated allele frequencies. The analysis was conducted for the whole dataset and for each cluster separately. The optimal number of genetic clusters was determined using both the original method of Pritchard et al. (2000) and the ΔK method of Evanno, Regnaut, and Goudet (2005), as implemented in the R package pophelper (Francis, 2017). The same package was used to average replicate runs of the optimal K (Jakobsson & Rosenberg, 2007) and plot the final output. In addition, to visualise genetic divergence between populations, we constructed a neighbour-joining (NJ) tree using the program POPTREEW (Takezaki, Nei, & Tamura, 2014). We used Nei's genetic distance (D<sub>A</sub>; Nei, Tajima, & Tateno, 1983) and performed 1,000 bootstraps. Genetic relationships between STRUCTURE clusters for the optimal K were visualised by drawing a NJ tree based on net nucleotide distances (Pritchard, Wen, & Falush, 2010) using the program NEIGHBOR in the PHYLIP package 3.695 (Felsenstein, 2005).
Isolation by distance was calculated via a Mantel test (Mantel & Valand, 1970) using the R package ade4 (Dray & Dufour, 2007), to explore the relationship between genetic and geographic distances among populations. We used standardized values of F<sub>ST</sub> (F<sub>ST</sub>/(1-F<sub>ST</sub>)) and log-transformed values of geographic distance as dependent and independent variables, respectively (Rousset, 1997). Significance was estimated with 10,000 permutations. Analyses were performed between all populations and by grouping sampling localities as indicated by STRUCTURE.

The estimation of recent dispersal was conducted using a twofold approach. An assignment test was performed in GeneClass2 (Piry et al., 2004) to assign or exclude reference populations as possible origin of individuals (Paetkau, Slade, Burden, & Estoup, 2004). The test was run only for populations with 10 or more genotyped individuals (49 populations). The same program was used to detect first generation migrants, i.e. individuals born in a population other than that where they were collected. Details on parameters used in these analyses are presented in the supplement.

The sibship assignment method implemented in Colony 2.0.6.4 (Jones & Wang, 2010) was used to infer the effective size ($N_e$) of populations with more than 15 genotyped individuals (35 populations) under the hypothesis of random mating. Details on parameters used in the analysis are presented in the supplement.

We tested for sex-biased dispersal by calculating F<sub>ST</sub>, F<sub>IS</sub> and assignment values (AI<sub>C</sub>) within each sex (Goudet, Perrin, & Waser, 2002) using the hierfstat R package (Goudet & Jombart, 2015). We performed 1,000 permutations using the “two sided” alternative method (Helfer, Broquet, & Fumagalli, 2012). F<sub>ST</sub> and F<sub>IS</sub> are expected to be lower and higher for the dispersing sex compared to the philopatric sex, respectively (Goudet et al., 2002). AI<sub>C</sub> values determine the probability that an individual genotype originated from the population from which it was sampled, correcting for differences in population genetic diversity (Favre, Balloux, Goudet, & Perrin, 1997). The distribution of AI<sub>C</sub> values is centred around a mean (mAI<sub>C</sub>) of zero, with lower values expected for the dispersing sex. In contrast, the variance of AI<sub>C</sub> (vAI<sub>C</sub>) is expected to be higher for the dispersing sex.

To examine whether the genetic structure revealed by the Bayesian clustering analysis could be explained by orographic features such as tributary valleys (i.e. valleys whose brooks or rivers flow into greater ones) and ridges (i.e. a chain of mountains or hills that form a continuous elevated crest), we conducted analyses of molecular variance (AMOVA) using a nested design (Excoffier, Smouse, & Quattro, 1992). We implemented a four-level hierarchical approach and ran two separate AMOVA analyses: in the first analysis we estimated variance components among
genetic clusters identified by STRUCTURE and among tributary valleys nested within clusters; next, in the second analysis we tested for evidence of structuring among valleys and among populations within valleys, without taking genetic clusters into account (see Werth et al., 2007 for a similar approach). Further details and parameters used in the analysis are described in the supplement.

We investigated how habitat type (lakes, streams and caves) and geographic variables (latitude, longitude and altitude) explained genetic diversity estimates using multiple linear regression models. Model selection was performed in R using backward stepwise selection, where variables were dropped iteratively from the full model minimizing AIC values. Violation of the assumptions of normality, homogeneity in variance, multicollinearity and autocorrelation were checked by examining the residuals. Analyses were performed between all populations and by grouping sampling localities as indicated by STRUCTURE.

To investigate *C. asper* evolutionary history and estimate divergence times among STRUCTURE-defined genetic lineages, we employed an approximate Bayesian computation (ABC) approach, as implemented in the software DIYABC 2.1.0 (Cornuet et al., 2014). We performed the computations both combining microsatellites and mtDNA data, and separately for microsatellites to assess the impact of using different types of markers on scenario choice and posterior parameter estimation. To reduce computational demands, we selected 50 individuals from each of the five genetic groups defined by STRUCTURE. Pilot runs confirmed that varying the sample size for microsatellites (from 30 individuals per cluster to all 1,299 individuals) did not substantially affect the final outcome in terms of best supported scenario and estimated parameters (Table S2). Within each group, we selected populations representative of all habitat types, choosing among individuals with STRUCTURE ancestry coefficient $q \ge 0.9$ to exclude potentially confounding effects of contemporary gene flow (see Ortego, Noguerales, Gugger, & Sork, 2015). Following the recommendations of Cabrera and Palsbøll (2017) to improve DIYABC ability to reveal the true demographic model, we focused on simple contrasting models and reduced the number of candidate scenarios to three (Figures 2 and S1). The first type of scenario is a null model with all five lineages diverging at the same time from a common ancestor (general scenario 1). The second type is a model of initial divergence between two eastern and western ancestral lineages, keeping the eastern as ancestral, and subsequent formation of the five current genetic lineages, as suggested by Valbuena-Ureña et al. (2018) (general scenario 2). Finally, the third type is a hierarchical split model directly following results from STRUCTURE analysis, where clusters
1 and 5 were generated from cluster 3, after an initial split between clusters 2, 3 and 4 (general scenario 3). Further details on model specifications and run parameters are outlined in the supplement (Table S3).

## Results

### *Multilocus genetic diversity*

From the 258 individuals analysed for the *cyt-b* gene, we identified a total of 11 haplotypes. The haplotype network showed that adjacent haplotypes were separated by a single mutational step and confirmed the presence of two main central haplotypes separated from each other by two mutational steps (haplotype codes H5 and H9; Figure S2). The overall mean haplotype (Hd) and nucleotide (Π) diversities were 0.570 ± 0.031 and 0.003 ± 0.0002, respectively.

Regarding microsatellites, we did not find evidence of stuttering or large allele dropout. Mean null allele frequency across all loci was 0.037, ranging from 0.018 to 0.069. Global F<sub>ST</sub> values with and without correcting for null alleles were 0.377 and 0.383, respectively, and had overlapping 95% CI (0.342–0.433 for F<sub>ST</sub> using ENA and 0.350–0.443 for F<sub>ST</sub> not using ENA), indicating that the impact of null alleles is negligible. After applying the Bonferroni correction (*P* < 0.0004), significant linkage disequilibrium was found only in two populations between a total of three pairs of loci (in population NDE between locus pairs Ca1-Us3 and Us7-Ca16, and in population RVT-A between locus pair Ca22-Ca29). Significant deviations from HWE were observed in 18 (19%) localities after Bonferroni correction. 13 loci indicated significant departures from HWE in one to eight populations: Ca32, Ca25, Ca 23, Us7, Ca24 and Ca22 in one population, Us3, Ca8 and Ca1 in two populations, Ca30 in three populations and Ca16 in eight populations. However, this is probably the result of genetic structure in the populations, as most of the loci showed occasional departures from HWE in three or more populations that were not consistent across populations or loci.

We recorded variable levels of nuclear genetic diversity across the study area (Table S1). Mean values were 0.445 for H<sub>O</sub> (0.162–0.698), 0.457 for H<sub>E</sub> (0.171–0.626) and 2.659 for Ar (1.380–3.050). Westernmost populations exhibited the highest values, together with a group of central-eastern populations (Figure 3). F<sub>IS</sub> values were generally low (mean F<sub>IS</sub> = 0.069), ranging from -0.210 to 0.367.

### *Population structure analyses*
STRUCTURE analysis revealed five well-supported groups (Figure 1). Log-likelihood values showed a steady increase from K = 2 to K = 5 before slowing down and eventually reaching a plateau (Figure S3). Although ΔK values showed several peaks at different values of K, the peak at K = 5 was markedly higher and corresponded to the smallest variance. This chaotic behaviour has been reported when analysing data displaying strong isolation by distance with STRUCTURE (Ferchaud et al., 2015). Therefore, we assumed K = 5 as the clustering solution that best explained the spatial genetic structure of the species at the Pyrenean scale.

The five clusters were spatially distributed over the Pyrenean chain along a longitudinal gradient: the first cluster included the north-eastern (French) localities and four central-southern (Spanish) localities; the second cluster grouped together all Andorran localities, the south-eastern (Spanish) sites and the north-eastern population Valmanya (B10); the third cluster included the central-western localities from both sides of the Pyrenees; the fourth cluster comprised all localities at both sides of the western Pyrenees; finally, sites located on the southern (Spanish) side of the central Pyrenees in-between the first three clusters formed a fifth group (see colour codes in Figure 1: cluster 1, blue; cluster 2, light green; cluster 3, orange; cluster 4, dark green; cluster 5, pink). Relatively high levels of admixture were detected where the genetic clusters met (Figure 1). When analysing each cluster separately, further substructure emerged from clusters 1 and 2 (i.e. the easternmost clusters; Figure S4): sampling localities in cluster 1 grouped into three subclusters and those included in cluster 2 grouped into four subclusters. The NJ tree for the five clusters indicated that clusters 2 and 4, corresponding to the clusters at the eastern and western edges of the species range, respectively, were the most genetically differentiated (Figure 1). In addition, cluster 4 was the richest in terms of genetic diversity (Table 1). The NJ tree inferred from D<sub>A</sub> distances over all populations revealed the five groups identified by STRUCTURE, with geographically close populations usually grouped together (Figure 4).

A significant isolation by distance (IBD) was found between all pairs of populations (R = 0.499, *P* < 0.001; Figure S5). Similar but generally stronger IBD patterns were revealed when analysing each cluster separately (cluster 1: R = 0.469, *P* < 0.001; cluster 2: R = 0.702, *P* < 0.001; cluster 3: R = 0.687, *P* < 0.001; cluster 5: R = 0.764, *P* < 0.001; Figure S5), with the exception of cluster 4 that did not show a significant IBD signal (*P* = 0.053).

Contemporary dispersal, effective population size and sex-biased dispersal
The assignment test conducted in GeneClass2 returned an assignment rate of 82.7%, meaning that 922 individuals out of 1,115 were assigned to the localities where they were sampled (Table S4). Although the majority of misassignments were to localities belonging to the same cluster, three populations from cluster 5 and one population from cluster 3 showed ancestry to cluster 1. A total of 63 (4.9%) individuals were identified as first generation migrants: 14 and 28 individuals were selected using the L<sub>home</sub> and L<sub>home</sub>/L<sub>max</sub> approaches, respectively, and 21 were selected by both likelihood methods. Of the 63 individuals, 27 had similar migration probabilities for several localities, indicating that these samples represented individuals whose source locality could not be determined due to the presence of unsampled populations in the study area. Among the 36 migration events with estimated origin, 19 involved stream populations only, 9 involved lake populations, 7 occurred between lake and stream populations and one between cave and stream populations. In all but one instance (one individual sampled in population E2 and detected to be coming from E1, which are separated by only 1.7 km), migration was limited within groups detected with STRUCTURE and usually involved geographically close populations (Figures 5 and S6). Indeed, most individuals migrated less than 1 km (17 individuals), or between 1 and 10 km (12 individuals). However, for four individuals we found potential for recent migration between localities separated by an Euclidean distance between 24 and 33 km. Dispersal between these localities would have implied either downstream migration or migration between adjacent glacial cirques, but no data are available from some intermediate localities. The remaining putative long dispersal events were below 12 km Euclidean distance and were all amongst adjacent glacial cirques.

Colony returned low values of effective population sizes (Table S1). Values ranged from nine in the cave population Pas du Loup (B1) to 46 breeding individuals in the stream population Ruisseau de Peyrenère (E4), with a mean N<sub>e</sub> of 26.

Results from sex-biased dispersal analysis showed that F<sub>ST</sub> and F<sub>IS</sub> values were not significantly different between sexes (males: F<sub>ST</sub> = 0.377, F<sub>IS</sub> = 0.088; females: F<sub>ST</sub> = 0.367, F<sub>IS</sub> = 0.101; P<sub>Fst</sub> = 0.610, P<sub>Fis</sub> = 0.300). Similarly, there was no significant difference in either the mean or the variance of AIc between sexes (males: mAIc = 0.052, vAIc = 16.332; females: mAIc = -0.051, vAIc = 17.771; P<sub>mAIc</sub> = 0.711, P<sub>vAIc</sub> = 0.375).

Influence of orography, geography and habitat
AMOVA analyses suggested significant structure at all tested levels (Table 2). When partitioning molecular variance between genetic clusters and tributary valleys, most molecular variance was found within valleys, followed by the among clusters component. Results did not differ substantially whether including in the analysis either all valleys or only those featuring a unique genetic cluster (data not shown). Within valleys, most variation was found among individuals, as expected for polymorphic loci such as microsatellites.

At the Pyrenean scale, model selection indicated that altitude had a significant positive effect on H<sub>E</sub> and Ar, whereas longitude had a significant negative effect on Ar (Figures 6 and S7). Regarding habitat types, streams showed significantly higher levels of genetic diversity compared to lakes and caves, although this pattern was lost when performing the analysis at the genetic cluster level. Indeed, only clusters 1, 2 and 4 showed significant effects. In cluster 1, altitude was negatively associated with Ar and longitude was negatively associated with both H<sub>E</sub> and Ar, and lakes were the most diverse habitat. In cluster 2, longitude had a negative effect and latitude a positive effect on both estimates; comparison between habitats was not possible because only streams were sampled. Finally, in cluster 4, longitude and latitude were negatively associated with both estimates and streams were the most diverse habitat.

## Colonisation history

The pre-evaluation step confirmed that the chosen priors ensured a good fit between simulated and observed data sets for all tested scenarios (Figure S8). Analyses suggested highest support for scenario 3 (the multiple-refugia population model directly following Bayesian clustering analysis results) regardless of the genetic markers used (microsatellites or microsatellites + mtDNA; Figure 2). This scenario had the highest posterior probability (PP) and its 95% CI did not overlap with those for the other scenarios (Table 3). Type I and type II errors for scenario 3 were low, denoting high confidence in scenario choice (Table 3). RMedAD values were relatively small (< 0.25 in most cases), indicating precise parameter estimations (Table 4). Finally, model checking revealed that the observed dataset fell within the cloud of points of the simulated datasets obtained from the parameter posterior distribution (Figure S8).

Analyses based on either microsatellites or microsatellites + mtDNA returned similar parameter estimates (Table 4). Results suggested that peripheral genetic lineages (clusters 2 and 4), together with the central group, diverged from a common ancestor around the Last Glacial Maximum (LGM), approximately 42,000–24,000 years ago (t3). Subsequently, the central-western
lineage (cluster 3) split from the central clade ~15,000–7,500 years ago (t2), whereas the most recent divergence occurred ~12,000–5,400 years ago (t1) between the central-southern and central-eastern lineages (clusters 5 and 1; Figure 2).

Discussion

*Refugia within refugia: the Pyrenees*

Mountain systems played a crucial role in determining species diversity, and the origin of intraspecific genetic structuring has been frequently tracked back to putative glacial refugia where populations survived Quaternary ice ages (Wallis, Waters, Upton, & Craw, 2016). In Europe, the Iberian Peninsula served as one of the most important Pleistocene glacial refugia (Gómez & Lunt, 2007). The complex climatic and topographic features of this region allowed for lineage persistence in “refugia within refugia”, the Pyrenees being one of them (Abellán & Svenning, 2014; Gómez & Lunt, 2007). For this reason, the Pyrenees are considered as a biodiversity hotspot with a rich endemic flora and fauna (Wallis et al., 2016). Here, ABC-based analyses revealed that *C. asper* microsatellite lineage differentiation started either during or slightly before the LGM (~42,000–24,000 years ago) at three main focal centres (western –cluster 4–, central and eastern – cluster 2– Pyrenees) and continued within the central group through the end of the Last Glacial Period, until ~12,000–5,500 years ago (Figure 2; Table 4). Indeed, the second and third splits straddled the Pleistocene-Holocene boundary and involved the central group only, with a first divergence event consisting of the separation of the central-western lineage (cluster 3) ~15,000–7,500 years ago, followed by a split between the central Spanish and the central-eastern French lineages (clusters 5 and 1; ~12,000–5,500 years ago).

Our study describes the existence of five main genetic lineages in *C. asper*, which are distributed longitudinally along the Pyrenees. Previous studies mainly reported two or three major longitudinal splits in the Pyrenees in a number of species, such as the mountain ringlet butterfly *Erebia epiphron* (Schmitt, Hewitt, & Muller, 2006), the European beech *Fagus sylvatica* (Magri et al., 2006), the snapdragon *Antirrhinum* (Liberal et al., 2014), the rusty-leaved alpenrose *Rhododendron ferrugineum* (Charrier et al., 2014) and the ground-dwelling spider *Harpactocrates ravastellus* (Bidegaray-Batista et al., 2016). However, most of these studies either dealt with species complexes and therefore evolutionary time lags of millions of years (Bidegaray-Batista et al., 2016; Liberal et al., 2014), or did not attempt to date back the phylogeographic history of the study species across the Pyrenees (Charrier et al., 2014; Magri et al., 2006; Schmitt et al., 2006).
Other studies have only focussed on the post-glacial colonisation history (e.g. from 15,000 years ago to the present), such as in the case of the Pyrenean rock lizard *Iberolacerta bonnali* (Ferchaud et al., 2015) or the water flea *Daphnia longispina* (Ventura et al., 2014).

## Phylogeography of *C. asper*

The times of the splits approximately correspond to major cooling events in the Pyrenees. The LGM in the Pyrenees is estimated to have occurred ~22,500–18,000 years ago (González-Sampériz et al., 2006); glacial advance likely promoted species retreat to isolated refugial areas (three main refugial areas: western, central and eastern) and subsequent genetic differentiation. After the LGM, a period of increase in temperature (the Bølling-Allerød period, ~15,000–13,000 years ago) may have created favourable conditions for dispersion outside the refugia and colonization of suitable areas in the Pyrenees. This was followed by a cold period (the Younger Dryas, ~13,000–11,500 years ago; González-Sampériz et al., 2006) that likely prompted species retreat to refugial areas where further genetic differentiation was favoured (divergence of cluster 3 from the central group). The Younger Dryas marked the end of the Pleistocene and, with the beginning of the Holocene, temperatures increased again, favouring species expansion uphill and towards the central Pyrenees. An additional abrupt cooling episode took place ~8,400–8,000 years ago (8,200–yr event; Alley et al., 1997; González-Sampériz et al., 2006), which likely promoted the last split between clusters 1 and 5.

During the last glaciation, Pyrenean glaciers reached their maximum extent earlier than the LGM at > 30,000 years ago, though a later glaciar re-advance occurred during the LGM (García-Ruiz, Valero-Garcés, Martí-Bono, & González-Sampériz, 2003). During these periods, most of the Pyrenees was extensively covered with ice and likely represented an unsuitable region for *C. asper*. Although we cannot rule out that some *C. asper* populations survived glacial events in microrefugia *in situ* (e.g. in deep valleys or on southern valley slopes), optimal conditions during glacial maxima existed mostly in peripheral areas outside the mountain range, unlike other species that likely survived in nunataks along the chain (e.g. Charrier et al., 2014). We thus hypothesise that, at the time of the first split, populations took refuge in three major refugial areas (corresponding to the western, central and eastern genetic lineages) located outside the mountain range. The long branches defining these three lineages in the NJ tree and their geographic consistency support a scenario of allopatric divergence and long-term lineage persistence in separated refugia (Figure 1). After the LGM, temperatures increased and created favourable
conditions for the species to recolonise suitable habitats inside the chain. The following splits were likely prompted by cooling events occurring over shorter intervals and characterized by a lesser glacier extent (i.e. the Younger Dryas and the 8,200-yr event; González-Sampériz et al., 2006), leading to a wide availability of habitats inside the Pyrenees even during cold periods. It is reasonable to assume that *C. asper* endured these cooling periods in refugia located within the Pyrenees, where differentiation of the central group was favoured.

We would like to stress that ABC modelling has some uncertainty. Firstly, the tested models do not represent a comprehensive range of all possible scenarios, but are instead based on a selection of hypotheses that we consider are most likely to reflect our data. We focused our analysis on three simple contrasting models aimed at capturing the key demographic events, avoiding overcomplex and similar models. This approach has proven useful to increase the ability of DIYABC to reveal the true model, as well as to better estimate the error and accuracy of parameter estimates (Cabrera & Palsbøll, 2017). Secondly, ABC modelling is based on scenarios where no gene flow is permitted between populations after they initially diverge. Only single events of admixture between populations are considered, whereas recurrent gene flow due to dispersal cannot be incorporated. However, we believe that not incorporating gene flow had only a marginal effect on our ABC results, as ABC analyses run using all 1,299 individuals (and thus including admixed populations located at cluster borders) yielded parameter estimates similar to those from computations based on 50 individuals per cluster (Table S2). Thirdly, it is important to note that the time estimates presented for *C. asper* have relatively large confidence intervals, although they still embrace values broadly referred to the time of the last glaciation.

## *Mito-nuclear discordance*

Population analyses of nuclear microsatellites revealed that the Pyrenean brook newt is subdivided into five well-supported genetic groups mainly distributed along a longitudinal gradient (Figure 1), with eastern genetic groups displaying finer substructure (Figure S4). This is in agreement with previous studies investigating the nuclear genetic structure of the species (Milá et al., 2010; Valbuena-Ureña et al., 2018). However, mitochondrial DNA did not show a clear phylogeographic pattern coinciding with the five microsatellite lineages (Figure S2). Haplotype H9 partly corresponds to cluster 2 (eastern Pyrenees; but see Valbuena-Ureña et al., 2013) and haplotype H7 shows some affinity to cluster 3 (central-western Pyrenees); the remaining area is dominated by haplotype H5, which is the most widespread haplotype. The almost perfect match between ABC
analyses based on either microsatellites or microsatellites + mtDNA was possibly due to the lack
of mtDNA variation. In *C. asper*, a similar mito-nuclear discordance was detected by Milá et al.
(2010): variation at several mtDNA regions (2,040 bp) was low, whereas differentiation at AFLP
loci was high and consistent with the structure here identified with microsatellites (see also
Valbuena-Ureña et al., 2018). Milá et al. (2010) suggested that variation at AFLP loci could have
been abnormally high because of the high amount of satellite DNA in *C. asper* genome, which
possibly interfered in the amplification. However, the marked genetic structuring detected with
microsatellites, which is consistent with the genetic units revealed by AFLP, indicates that AFLP
loci variation was not an artefact but the product of real population structuring in the species.
Divergence times estimated with microsatellites approximately correspond to major cooling events
that likely impacted and shaped the genetic constitution of *C. asper*. Furthermore, the high
differentiation at AFLP and microsatellite markers is consistent with the high morphological
diversification reported among *C. asper* populations (Montori, Llorente, & García-París, 2008).
An alternative possibility is that the observed mtDNA variation could be due to female-biased
dispersal, with female-mediated gene flow and phytopatric males leading to a pattern of mito-
nuclear discordance (Prugnolle & De Meeus, 2002). However, our results do not support a sex-
biased dispersal scenario. A more plausible explanation for the observed discordance would be a
selective sweep on mtDNA, bringing haplotypes H5 and H9 close to fixation in most populations
over most of the species range (see also Valbuena-Ureña et al., 2013). Empirical evidence of
selection on mtDNA is accumulating in the literature and possible cases of selective sweep have
been reported in a number of taxa (Bazin, Glémin, & Galtier, 2006; Bensch, Irwin, Irwin, Kvist, &
Åkesson, 2006; Ferchaud et al., 2015; Rato, Carranza, Perera, Carretero, & Harris, 2010). As for
*C. asper*, a selective sweep of favourable mtDNA variants was previously suggested by Milá et al.
(2010) to explain the lack of mtDNA diversity. A selective sweep could account for the low
variation in mtDNA compared to nuclear DNA and for the geographic distribution of haplotypes.
However, further studies are needed to confirm this hypothesis.

Contemporary dispersal and influence of environmental and geographic variables

Our analyses revealed restricted contemporary gene flow and dispersal between populations of *C. asper* across the five genetic lineages (Figure 5; Table S4). This is supported by the clear pattern of isolation by distance (Figure S5) and by 19% of the observed genetic variation being explained by differences between major genetic clusters (Table 2). However, population structure analysis
revealed admixture patterns at boundaries between genetic clusters, implying potential recent gene flow across all clusters borders (Figure 1). Molecular estimates of dispersal corroborated this finding: genetic signs of contemporary dispersal, albeit weak, were detected between a number of populations located at clusters’ borders. This holds especially true for cluster 5, with three populations showing ancestry to cluster 1 (Table S4). According to ABC analyses, clusters 1 and 5 were the last to diverge and may have retained a higher degree of connectivity (Figure 2).

Moderate levels of dispersal and connectivity between habitat types were detected within genetic clusters (Figure 5; Table S4). Nevertheless, migration preferentially involved geographically close populations (0–4 km Euclidean distance; Figure S6) and it was mostly restricted within valleys. This is in agreement with Montori, Llorente, and Richter-Boix (2008), which mainly recorded short-range movements in *C. asper* using a capture-recapture framework. The short mean dispersal distances, coupled with low effective population sizes ($N_e < 50$), may explain the high levels of genetic structuring and differentiation for *C. asper* populations across the entire species range. On the other hand, our estimations suggested potential for rare long-distance dispersal (up to 33 km). This might include both movements along the stream network and overland dispersal (Grant, Nichols, Lowe, & Fagan, 2010). Some individuals could have also been carried downstream during floods (Montori et al., 2012). However, although long-distance dispersal of few individuals per population remains possible in amphibians (Cayuela et al., 2020), a plausible alternative scenario is that potential unsampled source populations located in between the study sites may have been at the origin of migrants if they shared alleles with the putative sites of origin. This is possible given the high availability of suitable habitats for *C. asper* in the study area. Nevertheless, long distance dispersal, possibly over a few successive generations (Saura, Bodin, & Fortin, 2014), is in line with our estimates of genetic diversity, as shown by most populations presenting low inbreeding coefficients (mean $F_{IS}$ = 0.069) and levels of genetic variability within the range of other urodeles and temperate amphibians (Chan & Zamudio, 2009).

The high overall $F_{ST}$ value, together with the clear pattern of isolation by distance (especially at the genetic cluster level), indicate that divergence between populations is spatially structured. The strong spatial structuring, even across contrasting habitats, suggests no support for isolation by environment (Orsini et al., 2013). Indeed, populations from different habitats clustered together in four of the five lineages, and neighbour-joining analysis showed that populations are mainly grouped by valleys rather than habitats (Figure 4). Marked genetic differentiation exists at the scale of tributary valleys, as suggested by 20.5% of the molecular variance being attributable
to differences between valleys (Table 2). Furthermore, we detected recent dispersal (as inferred by microsatellites) among populations inhabiting different habitats. In accordance with Valbuena-Ureña et al. (2018), we found evidence for a negative longitudinal and positive altitudinal gradient of genetic diversity over all *C. asper* populations, and streams showed higher values of genetic diversity compared to lakes and caves (Figures 6 and S7). This trend has been previously interpreted as evidence of preference for cooler and wetter environments, typical of the western sector of the Pyrenees and high altitudes, by *C. asper* (Valbuena-Ureña et al., 2018). However, linear models conducted at the genetic cluster level revealed contrasting patterns of genetic diversity that do not conform with the general trend. This, together with the strong isolation by distance revealed at the cluster level, suggests that the pattern detected at the Pyrenean scale is likely the result of independent drivers acting within clusters. Clusters may thus be considered as independent units as a result of independent phylogeographic histories, each being the product of separate post-glacial colonisation routes. In light of the above, isolation by colonisation remains a plausible explanation for the resulting pattern of isolation by distance (Orsini et al., 2013), but further studies focussing on local adaptation might be necessary to confirm this point (see also Oromi et al., 2018). An alternative possibility is that the contrasting patterns at the cluster level could have arisen through the combined effects of latitude, longitude and habitat type. Habitat type might have an influence on the level of genetic variation in the residing populations and the contrasting patterns among clusters could be caused by the differential availability of these habitats in different areas.

## Concluding remarks

This study highlights the importance of integrating past evolutionary processes and present-day gene flow and dispersal dynamics to shed light onto what shaped (and is currently shaping) the observed genetic composition and structure of endemic species. Here, we demonstrate that the endemic newt *C. asper* probably recolonized the Pyrenees from at least five distinct glacial refugia. Differentiation started before the LGM and continued through the end of the Last Glacial Period, leading to the formation of five well-supported genetic lineages that likely underwent separate evolutionary histories. There is currently limited gene flow between lineages, although borders represent zones of admixture resulting from postglacial recolonization of formerly glaciated areas. Within lineages, dispersal distances are relatively short, although long-distance dispersal may be accomplished by a few individuals. The incongruence between the high variation
in nuclear DNA and low variation in mtDNA could be interpreted as evidence of selective sweep
in mtDNA and underscores the importance of using a multilocus approach to achieve a complete
picture of the population structure and history of the study species. Given the age of the studied
lineages and the restricted present-day gene flow, we suggest that these broad areas should be
regarded as separate management units worthy of independent conservation consideration. At
smaller spatial scales, specific lake populations of *C. asper* have been also found to merit special
conservation focus (i.e. the paedomorphic populations described in Oromi et al., 2018).

Acknowledgements

We thank Meritxell Cases, Alba Castrillón, Eloi Cruset, Blanca Font, Ismael Jurado and Quim Pou-Rovira for field assistance. Economic support was provided by the European Commission LIFE+ project LimnoPirineus (LIFE13 NAT/ES/001210), by the Spanish Government project Funbio (RTI2018-096217-B-I00), by the Interreg POCTEFA ECTOPYR project (EFA031/15) and by the Societas Europaea Herpetologica (SEH, research grant awarded to F.L.). F.L. had a doctoral grant funded by Fundação para a Ciência e Tecnologia (FCT, grant number PD/BD/52598/2014). M.D. is research Director at Fonds de la Recherche Scientifique – FNRS.
## References

Abellán, P., & Svenning, J.-C. (2014). Refugia within refugia-patterns in endemism and genetic divergence are linked to Late Quaternary climate stability in the Iberian Peninsula. *Biological Journal of the Linnean Society*, 113(1), 13-28. doi:10.1111/bij.12309

Allentoft, M. E., Siegismund, H. R., Briggs, L., & Andersen, L. W. (2009). Microsatellite analysis of the natterjack toad (*Bufo calamita*) in Denmark: populations are islands in a fragmented landscape. *Conservation Genetics*, 10(1), 15-28. doi:10.1007/s10592-008-9510-8

Alley, R. B., Mayewski, P. A., Sowers, T., Stuiver, M., Taylor, K. C., & Clark, P. U. (1997). Holocene climatic instability: A prominent, widespread event 8200 yr ago. *Geology*, 25(6), 483-486. doi:10.1130/0091-7613(1997)025<0483:HCIAPW>2.3.CO;2

Avise, J. C. (2000). *Phylogeography: the history and formation of species*. Cambridge, MA, USA: Harvard University Press.

Bazin, E., Glémin, S., & Galtier, N. (2006). Population size does not influence mitochondrial genetic diversity in animals. *Science*, 312(5773), 570-572. doi:10.1126/science.1122033

Bensch, S., Irwin, D. E., Irwin, J. H., Kvist, L., & Åkesson, S. (2006). Conflicting patterns of mitochondrial and nuclear DNA diversity in *Phylloscopus warblers*. *Molecular Ecology*, 15(1), 161-171. doi:10.1111/j.1365-294X.2005.02766.x

Bidegaray-Batista, L., Sánchez-Gracia, A., Santulli, G., Maiorano, L., Guisan, A., Vogler, A. P., & Arnedo, M. A. (2016). Imprints of multiple glacial refugia in the Pyrenees revealed by phylogeography and palaeodistribution modelling of an endemic spider. *Molecular Ecology*, 25(9), 2046-2064. doi:10.1111/mec.13585

Bosch, J., Tejedo, M., Lecis, R., Miaud, C., Lizana, M., Edgar, P., ... Marquez, R. G., P. (2009). *Calotriton asper*. The IUCN Red List of Threatened Species 2009: e.T59448A11943040. doi:10.2305/IUCN.UK.2009.RLTS.T59448A11943040.en.

Burns, E. L., Eldridge, M. D., & Houlden, B. A. (2004). Microsatellite variation and population structure in a declining Australian Hylid *Litoria aurea*. *Molecular Ecology*, 13(7), 1745-1757. doi:10.1111/j.1365-294X.2004.02190.x

Cabrera, A. A., & Palsbøll, P. J. (2017). Inferring past demographic changes from contemporary genetic data: A simulation-based evaluation of the ABC methods implemented in DIYABC. *Molecular Ecology Resources*, 17(6), e94-e110. doi:10.1111/1755-0998.12696
Calvet, M. (2004). The Quaternary glaciation of the Pyrenees. *Developments in Quaternary Sciences*, 2, 119-128. doi:10.1016/S1571-0866(04)80062-9

Caplat, P., Edelaar, P., Dudaniec, R. Y., Green, A. J., Okamura, B., Cote, J., ... Petit, E. J. (2016). Looking beyond the mountain: dispersal barriers in a changing world. *Frontiers in Ecology and the Environment, 14*(5), 262-269. doi:10.1002/fee.1280

Carranza, S., & Amat, F. (2005). Taxonomy, biogeography and evolution of *Euproctus* (Amphibia: Salamandridae), with the resurrection of the genus *Calotriton* and the description of a new endemic species from the Iberian Peninsula. *Zoological Journal of the Linnean Society, 145*(4), 555-582. doi:10.1111/j.1096-3642.2005.00197.x

Carranza, S., Arnold, E., Mateo, J. A., & López-Jurado, L. F. (2000). Long-distance colonization and radiation in gekkonid lizards, *Tarentola* (Reptilia: Gekkonidae), revealed by mitochondrial DNA sequences. *Proceedings of the Royal Society of London. Series B: Biological Sciences, 267*(1444), 637-649. doi:10.1098/rspb.2000.1050

Cayuela, H., Valenzuela-Sanchez, A., Teulier, L., Martínez-Solano, Í., Léna, J.-P., Merilä, J., ... Schmidt, B. R. (2020). Determinants and consequences of dispersal in vertebrates with complex life cycles: a review of pond-breeding amphibians. *The Quarterly Review of Biology, 95*(1), 1-36. doi:doi.org/10.1086/707862

Chan, L. M., & Zamudio, K. R. (2009). Population differentiation of temperate amphibians in unpredictable environments. *Molecular Ecology, 18*(15), 3185-3200. doi:10.1111/j.1365-294X.2009.04273.x

Chapuis, M.-P., & Estoup, A. (2006). Microsatellite null alleles and estimation of population differentiation. *Molecular Biology and Evolution, 24*(3), 621-631. doi:10.1093/molbev/msl191

Charrier, O., Dupont, P., Pornon, A., & Escaravage, N. (2014). Microsatellite marker analysis reveals the complex phylogeographic history of *Rhododendron ferrugineum* (Ericaceae) in the Pyrenees. *PLoS One, 9*(3), e92976. doi:10.1371/journal.pone.0092976

Chiucchi, J. E., & Gibbs, H. (2010). Similarity of contemporary and historical gene flow among highly fragmented populations of an endangered rattlesnake. *Molecular Ecology, 19*(24), 5345-5358. doi:10.1111/j.1365-294X.2010.04860.x

Clergue-Gazeau, M., & Martínez-Rica, J. (1978). Les différents biotopes de l'urodèle pyrénéen, *Euproctus asper*. *Bulletin de la Société d'Histoire Naturelle de Toulouse, 114*(3-4), 461-471.
Clobert, J., Le Galliard, J. F., Cote, J., Meylan, S., & Massot, M. (2009). Informed dispersal, heterogeneity in animal dispersal syndromes and the dynamics of spatially structured populations. *Ecology Letters*, 12(3), 197-209. doi:10.1111/j.1461-0248.2008.01267.x

Cornuet, J. M., Pudlo, P., Veyssier, J., Dehne-Garcia, A., Gautier, M., Leblois, R., ... Estoup, A. (2014). DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. *Bioinformatics*, 30(8), 1187-1189. doi:10.1093/bioinformatics/btt763

Darriba, D., Taboada, G. L., Doallo, R., & Posada, D. (2012). jModelTest 2: more models, new heuristics and parallel computing. *Nature Methods*, 9(8), 772. doi:10.1038/nmeth.2109

Denoël, M., Dalleur, S., Langrand, E., Besnard, A., & Cayuela, H. (2018). Dispersal and alternative breeding site fidelity strategies in an amphibian. *Ecography*, 41(9), 1543-1555. doi:10.1111/ecog.03296

Dray, S., & Dufour, A.-B. (2007). The ade4 package: implementing the duality diagram for ecologists. *Journal of Statistical Software*, 22(4), 1-20. doi:10.18637/jss.v022.i04

Drechsler, A., Geller, D., Freund, K., Schmeller, D. S., Kuenzel, S., Rupp, O., ... Steinfartz, S. (2013). What remains from a 454 run: estimation of success rates of microsatellite loci development in selected newt species (*Calotriton asper*, *Lissotriton helveticus*, and *Triturus cristatus*) and comparison with Illumina-based approaches. *Ecology and Evolution*, 3(11), 3947-3957. doi:10.1002/ece3.764

Epps, C. W., & Keyghobadi, N. (2015). Landscape genetics in a changing world: disentangling historical and contemporary influences and inferring change. *Molecular Ecology*, 24(24), 6021-6040. doi:10.1111/mec.13454

Evanno, G., Regnaut, S., & Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. *Molecular Ecology*, 14(8), 2611-2620. doi:10.1111/j.1365-294X.2005.02553.x

Excoffier, L., Smouse, P. E., & Quattro, J. M. (1992). Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. *Genetics*, 131(2), 479-491.

Favre, L., Balloux, F., Goudet, J., & Perrin, N. (1997). Female-biased dispersal in the monogamous mammal *Crocidura russula*: evidence from field data and microsatellite patterns. *Proceedings of the Royal Society of London. Series B: Biological Sciences*, 264(1378), 127-132. doi:10.1098/rspb.1997.0019
Felsenstein, J. (2005). PHYLIP (phylogeny inference package) version 3.6. Distributed by the author. Seattle (WA): Department of Genome Sciences, University of Washington.

Ferchaud, A. L., Eudeline, R., Arnal, V., Cheylan, M., Pottier, G., Leblois, R., & Crochet, P. A. (2015). Congruent signals of population history but radically different patterns of genetic diversity between mitochondrial and nuclear markers in a mountain lizard. *Molecular Ecology, 24*(1), 192-207. doi:10.1111/mec.13011

Francis, R. M. (2017). pophelper: an R package and web app to analyse and visualize population structure. *Molecular Ecology Resources, 17*(1), 27-32. doi:10.1111/1755-0998.12509

García-Ruiz, J. M., Valero-Garcés, B. L., Martí-Bono, C., & González-Sampériz, P. (2003). Asynchroneity of maximum glacier advances in the central Spanish Pyrenees. *Journal of Quaternary Science, 18*(1), 61-72. doi:10.1002/jqs.715

Gill, D. E. (1978). The metapopulation ecology of the red-spotted newt, *Notophthalmus viridescens* (Rafinesque). *Ecological monographs*, 48(2), 145-166. doi:10.2307/2937297

Gómez, A., & Lunt, D. H. (2007). Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. In S. Weiss & N. Ferrand (Eds.), *Phylogeography of southern European refugia* (pp. 155-188). Amsterdam: Springer.

González-Sampériz, P., Valero-Garcés, B. L., Moreno, A., Jalut, G., García-Ruiz, J. M., Martí-Bono, C., ... Dedoubat, J. (2006). Climate variability in the Spanish Pyrenees during the last 30,000 yr revealed by the El Portalet sequence. *Quaternary Research, 66*(1), 38-52. doi:10.1016/j.yqres.2006.02.004

Goudet, J. (2002). FSTAT version 2.9. 3.2, a program to estimate and test gene diversities and fixation indices. Lausanne, Switzerland: Institute of Ecology. http://www2.unil.ch/popgen/softwares/fstat.htm.

Goudet, J., & Jombart, T. (2015). hierfstat: estimation and tests of hierarchical F-statistics. R package version 0.04-22.

Goudet, J., Perrin, N., & Waser, P. (2002). Tests for sex-biased dispersal using bi-parentally inherited genetic markers. *Molecular Ecology, 11*(6), 1103-1114. doi:10.1046/j.1365-294X.2002.01496.x

Grant, E. H. C., Nichols, J. D., Lowe, W. H., & Fagan, W. F. (2010). Use of multiple dispersal pathways facilitates amphibian persistence in stream networks. *Proceedings of the National Academy of Sciences, 107*(15), 6936-6940. doi:10.1073/pnas.1000266107
Helfer, V., Broquet, T., & Fumagalli, L. (2012). Sex-specific estimates of dispersal show female philopatry and male dispersal in a promiscuous amphibian, the alpine salamander (*Salamandra atra*). *Molecular Ecology*, *21*(19), 4706-4720. doi:10.1111/j.1365-294X.2012.05742.x

Hewitt, G. M. (1999). Post-glacial re-colonization of European biota. *Biological Journal of the Linnean Society*, 68(1-2), 87-112. doi:10.1006/bijl.1999.0332

Hewitt, G. M. (2000). The genetic legacy of the Quaternary ice ages. *Nature*, 405(6789), 907-913. doi:10.1038/35016000

Hewitt, G. M. (2004). Genetic consequences of climatic oscillations in the Quaternary. *Philosophical Transactions of the Royal Society of London B: Biological Sciences*, 359(1442), 183-195; discussion 195. doi:10.1098/rstb.2003.1388

Hewitt, G. M., & Butlin, R. K. (1997). Causes and consequences of population structure. In J. R. Krebs & N. Davies (Eds.), *Behavioral Ecology*, 4th edn. (pp. 350-372). Oxford: Blackwell.

Holderegger, R., & Thiel-Egenter, C. (2009). A discussion of different types of glacial refugia used in mountain biogeography and phylogeography. *Journal of Biogeography*, 36(3), 476-480. doi:10.1111/j.1365-2699.2008.02027.x

Jakobsson, M., & Rosenberg, N. A. (2007). CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. *Bioinformatics*, 23(14), 1801-1806. doi:10.1093/bioinformatics/btm233

Johnson, M. L., & Gaines, M. S. (1990). Evolution of dispersal: theoretical models and empirical tests using birds and mammals. *Annual Review of Ecology and Systematics*, 21(1), 449-480. doi:10.1146/annurev.es.21.110190.002313

Jones, O. R., & Wang, J. (2010). COLONY: a program for parentage and sibship inference from multilocus genotype data. *Molecular Ecology Resources*, 10(3), 551-555. doi:10.1111/j.1755-0998.2009.02787.x

Kalinowski, S. T. (2005). hp-rare 1.0: a computer program for performing rarefaction on measures of allelic richness. *Molecular Ecology Notes*, 5(1), 187-189. doi:10.1111/j.1471-8286.2004.00845.x

Kearse, M., Moir, R., Wilson, A., Stones-Havas, S., Cheung, M., Sturrock, S., . . . Duran, C. (2012). Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. *Bioinformatics*, 28(12), 1647-1649. doi:10.1093/bioinformatics/bts199
Kraaijeveld-Smit, F. J., Beebee, T. J., Griffiths, R. A., Moore, R. D., & Schley, L. (2005). Low gene flow but high genetic diversity in the threatened Mallorcan midwife toad *Alytes muletensis*. *Molecular Ecology*, *14*(11), 3307-3315. doi:10.1111/j.1365-294X.2005.02614.x

Kumar, S., Stecher, G., & Tamura, K. (2016). MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for bigger datasets. *Molecular Biology and Evolution*, *33*(7), 1870-1874. doi:10.1093/molbev/msw054

Li, X. Y., & Kokko, H. (2019). Sex-biased dispersal: a review of the theory. *Biological Reviews*, *94*(2), 721-736. doi:10.1111/brv.12475

Liberal, I. M., Burrus, M., Suchet, C., Thebaud, C., & Vargas, P. (2014). The evolutionary history of *Antirrhinum* in the Pyrenees inferred from phylogeographic analyses. *BMC Evolutionary Biology*, *14*(1), 146. doi:10.1186/1471-2148-14-146

[dataset] Lucati, F., Poignet, M., Miró, A., Trochet, A., Aubret, F., Barthe, L., ... Ventura, M. (2020). Data from: Multiple glacial refugia and contemporary dispersal shape the genetic structure of an endemic amphibian from the Pyrenees. Dryad Digital Repository, doi: 10.5061/dryad.5tb2rbp23

Magri, D., Vendramin, G. G., Comps, B., Dupanloup, I., Geburek, T., Gömöry, D., ... Roure, J. M. (2006). A new scenario for the Quaternary history of European beech populations: palaeobotanical evidence and genetic consequences. *New Phytologist*, *171*(1), 199-221. doi:10.1111/j.1469-8137.2006.01740.x

Mantel, N., & Valand, R. S. (1970). A technique of nonparametric multivariate analysis. *Biometrics*, *26*(3), 547-558.

Martínez-Rica, J., & Clergue-Gazeau, M. (1977). Données nouvelles sur la répartition géographique de l'espèce *Euproctus asper* Dugès, Urodèle, Salamandridae. *Bulletin de la Société d'Histoire Naturelle de Tolouse*, *113*(3-4), 318-330.

Milá, B., Carranza, S., Guillaume, O., & Clobert, J. (2010). Marked genetic structuring and extreme dispersal limitation in the Pyrenean brook newt *Calotriton asper* (Amphibia: Salamandridae) revealed by genome-wide AFLP but not mtDNA. *Molecular Ecology*, *19*(1), 108-120. doi:10.1111/j.1365-294X.2009.04441.x

Montero-Pau, J., Gómez, A., & Muñoz, J. (2008). Application of an inexpensive and high-throughput genomic DNA extraction method for the molecular ecology of
zooplanktonic diapausing eggs. *Limnology and Oceanography: Methods*, **6**(6), 218-222.
doi:10.4319/lom.2008.6.218

Montori, A., & Llorente, G. A. (2014). Tritón pirenaico–*Calotriton asper* (Dugès, 1852). In A. Salvador & I. Martínez-Solano (Eds.), *Enciclopedia Virtual de los Vertebrados* (pp. 28). Madrid: Museo Nacional de Ciencias Naturales.

Montori, A., Llorente, G. A., & García-París, M. (2008). Allozyme differentiation among populations of the Pyrenean newt *Calotriton asper* (Amphibia: Caudata) does not mirror their morphological diversification. *Zootaxa, 1945*, 39-50. doi:10.11646/zootaxa.1945.1.2

Montori, A., Llorente, G. A., & Richter-Boix, A. (2008). Habitat features affecting the small-scale distribution and longitudinal migration patterns of *Calotriton asper* in a Pre-Pyrenean population. *Amphibia-Reptilia, 29*(3), 371-381. doi:10.1163/156853808785112048

Montori, A., Richter-Boix, A., Franch, M., Santos, X., Garriga, N., & Llorente, G. A. (2012). Natural fluctuations in a stream dwelling newt as a result of extreme rainfall: a 21-year survey of a *Calotriton asper* population. *Basic and Applied Herpetology, 26*, 43-56.
doi:10.11160/bah.12001

Mouret, V., Guillaumet, A., Cheylan, M., Pottier, G., Ferchaud, A. L., & Crochet, P. A. (2011). The legacy of ice ages in mountain species: post-glacial colonization of mountain tops rather than current range fragmentation determines mitochondrial genetic diversity in an endemic Pyrenean rock lizard. *Journal of Biogeography*, **38**(9), 1717-1731.
doi:10.1111/j.1365-2699.2011.02514.x

Nei, M., Tajima, F., & Tateno, Y. (1983). Accuracy of estimated phylogenetic trees from molecular data. *Journal of Molecular Evolution*, **19**(2), 153-170. doi:10.1007/bf02300753

Nichols, R. A., & Beaumont, M. A. (1996). Is it ancient or modern history that we can read in the genes? In M. E. Hochberg, J. Clobert, & R. Barbault (Eds.), *Aspects of the Genesis and Maintenance of Biological Diversity* (pp. 69-87). Oxford: Oxford University Press.

Noguerales, V., Cordero, P. J., & Ortego, J. (2016). Hierarchical genetic structure shaped by topography in a narrow-endemic montane grasshopper. *BMC Evolutionary Biology*, **16**(1), 96. doi:10.1186/s12862-016-0663-7

Noguerales, V., Cordero, P. J., & Ortego, J. (2017). Testing the role of ancient and contemporary landscapes on structuring genetic variation in a specialist grasshopper. *Ecology and Evolution*, **7**(9), 3110-3122. doi:10.1002/ece3.2810
Oromi, N., Valbuena-Ureña, E., Soler-Membrives, A., Amat, F., Camarasa, S., Carranza, S., ... Denoël, M. (2018). Genetic structure of lake and stream populations in a Pyrenean amphibian (*Calotriton asper*) reveals evolutionary significant units associated with paedomorphosis. *Journal of Zoological Systematics and Evolutionary Research*, 57, 418-430. doi:10.1111/jzs.12250

Orsini, L., Vanoverbeke, J., Swillen, I., Mergeay, J., & De Meester, L. (2013). Drivers of population genetic differentiation in the wild: isolation by dispersal limitation, isolation by adaptation and isolation by colonization. *Molecular Ecology*, 22(24), 5983-5999.
doi:10.1111/mec.12561

Ortego, J., Noguerales, V., Gugger, P. F., & Sork, V. L. (2015). Evolutionary and demographic history of the Californian scrub white oak species complex: an integrative approach. *Molecular Ecology*, 24(24), 6188-6208. doi:10.1111/mec.13457

Paetkau, D., Slade, R., Burden, M., & Estoup, A. (2004). Genetic assignment methods for the direct, real-time estimation of migration rate: a simulation-based exploration of accuracy and power. *Molecular Ecology*, 13(1), 55-65. doi:10.1046/j.1365-294X.2004.02008.x

Piry, S., Alapetite, A., Cornuet, J. M., Paetkau, D., Baudouin, L., & Estoup, A. (2004).
GENECLASS2: a software for genetic assignment and first-generation migrant detection.
*Journal of Heredity*, 95(6), 536-539. doi:10.1093/jhered/esh074

Pritchard, J. K., Stephens, M., & Donnelly, P. (2000). Inference of population structure using multilocus genotype data. *Genetics*, *155*(2), 945-959.

Pritchard, J. K., Wen, X., & Falush, D. (2010). Documentation for STRUCTURE software, version 2.3. Department of Human Genetics University of Chicago, Department of Statistics University of Oxford.

Prugnolle, F., & De Meeus, T. (2002). Inferring sex-biased dispersal from population genetic tools: a review. *Heredity*, *88*(3), 161-165. doi:10.1038/sj.hdy.6800060

R Core Team. (2018). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. http://www.R-project.org/.

Rato, C., Carranza, S., Perera, A., Carretero, M. A., & Harris, D. J. (2010). Conflicting patterns of nucleotide diversity between mtDNA and nDNA in the Moorish gecko, *Tarentola mauritanica*. *Molecular Phylogenetics and Evolution*, *56*(3), 962-971.
doi:10.1016/j.ympev.2010.04.033

Rice, W. R. (1989). Analyzing tables of statistical tests. *Evolution*, *43*(1), 223-225.
Rioux Paquette, S. (2011). PopGenKit: useful functions for (batch) file conversion and data resampling in microsatellite datasets. R package version 1.0.

Roffler, G. H., Talbot, S. L., Luikart, G., Sage, G. K., Pilgrim, K. L., Adams, L. G., & Schwartz, M. K. (2014). Lack of sex-biased dispersal promotes fine-scale genetic structure in alpine ungulates. *Conservation Genetics, 15*(4), 837-851. doi:10.1007/s10592-014-0583-2

Ronce, O. (2007). How does it feel to be like a rolling stone? Ten questions about dispersal evolution. *Annual Review of Ecology, Evolution, and Systematics, 38*, 231-253.
doi:10.1146/annurev.ecolsys.38.091206.095611

Rousset, F. (1997). Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. *Genetics, 145*(4), 1219-1228.

Rousset, F. (2008). GENEPOP'007: a complete re-implementation of the GENEPOP software for Windows and Linux. *Molecular Ecology Resources, 8*(1), 103-106. doi:10.1111/j.1471-8286.2007.01931.x

Rozas, J., Ferrer-Mata, A., Sánchez-DelBarrio, J. C., Guirao-Rico, S., Librado, P., Ramos-Onsins, S. E., & Sánchez-Gracia, A. (2017). DnaSP 6: DNA sequence polymorphism analysis of large data sets. *Molecular Biology and Evolution, 34*(12), 3299-3302.
doi:10.1093/molbev/msx248

Saccheri, I., Kuussaari, M., Kankare, M., Vikman, P., Fortelius, W., & Hanski, I. (1998). Inbreeding and extinction in a butterfly metapopulation. *Nature, 392*(6675), 491-494.
doi:10.1038/33136

Salzburger, W., Ewing, G. B., & Von Haeseler, A. (2011). The performance of phylogenetic algorithms in estimating haplotype genealogies with migration. *Molecular Ecology, 20*(9), 1952-1963. doi:10.1111/j.1365-294X.2011.05066.x

Saura, S., Bodin, Ö., & Fortin, M. J. (2014). Stepping stones are crucial for species' long-distance dispersal and range expansion through habitat networks. *Journal of Applied Ecology, 51*(1), 171-182. doi:10.1111/1365-2664.12179

Schmitt, T. (2009). Biogeographical and evolutionary importance of the European high mountain systems. *Frontiers in Zoology, 6*, 9. doi:10.1186/1742-9994-6-9

Schmitt, T., Hewitt, G. M., & Muller, P. (2006). Disjunct distributions during glacial and interglacial periods in mountain butterflies: *Erebia epiphron* as an example. *Journal of Evolutionary Biology, 19*(1), 108-113. doi:10.1111/j.1420-9101.2005.00980.x
Smith, M. A., & Green, D. M. (2005). Dispersal and the metapopulation paradigm in amphibian ecology and conservation: are all amphibian populations metapopulations? *Ecography*, 28(1), 110-128. doi:10.1111/j.0906-7590.2005.04042.x

Smith, M. A., & Green, D. M. (2006). Sex, isolation and fidelity: unbiased long-distance dispersal in a terrestrial amphibian. *Ecography*, 29(5), 649-658. doi:10.1111/j.2006.0906-7590.04584.x

Stamatakis, A. (2006). RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. *Bioinformatics*, 22(21), 2688-2690.
doi:10.1093/bioinformatics/btl446

Taberlet, P., Fumagalli, L., Wust-Saucy, A. G., & Cosson, J. F. (1998). Comparative phylogeography and postglacial colonization routes in Europe. *Molecular Ecology*, 7(4), 453-464. doi:10.1046/j.1365-294x.1998.00289.x

Takezaki, N., Nei, M., & Tamura, K. (2014). POPTREEW: web version of POPTREE for constructing population trees from allele frequency data and computing some other quantities. *Molecular Biology and Evolution*, 31(6), 1622-1624.
doi:10.1093/molbev/msu093

Tallmon, D. A., Luikart, G., & Waples, R. S. (2004). The alluring simplicity and complex reality of genetic rescue. *Trends in Ecology & Evolution*, 19(9), 489-496.
doi:10.1016/j.tree.2004.07.003

Trochet, A., Courtois, E. A., Stevens, V. M., Baguette, M., Chaine, A., Schmeller, D. S., ...
Wiens, J. J. (2016). Evolution of sex-biased dispersal. *The Quarterly Review of Biology*,
91(3), 297-320. doi:10.1086/688097

Tucker, J. M., Allendorf, F. W., Truex, R. L., & Schwartz, M. K. (2017). Sex-biased dispersal and spatial heterogeneity affect landscape resistance to gene flow in fisher. *Ecosphere*, 8(6), e01839. doi:10.1002/ecs2.1839

Valbuena-Ureña, E., Amat, F., & Carranza, S. (2013). Integrative phylogeography of *Calotriton newts* (Amphibia, Salamandridae), with special remarks on the conservation of the endangered Montseny brook newt (*Calotriton arnoldi*). *PLoS One*, 8(6), e62542.
doi:10.1371/journal.pone.0062542

Valbuena-Ureña, E., Oromi, N., Soler-Membrives, A., Carranza, S., Amat, F., Camarasa, S., ...
Steinfartz, S. (2018). Jailed in the mountains: Genetic diversity and structure of an endemic
newt species across the Pyrenees. *PLoS One*, 13(8), e0200214.
doi:10.1371/journal.pone.0200214

Van Oosterhout, C., Hutchinson, W. F., Wills, D. P. M., & Shipley, P. (2004). MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. *Molecular Ecology Notes*, 4(3), 535-538. doi:10.1111/j.1471-8286.2004.00684.x

Ventura, M., Petrusek, A., Miró, A., Hamrová, E., Buñay, D., De Meester, L., & Mergeay, J. (2014). Local and regional founder effects in lake zooplankton persist after thousands of years despite high dispersal potential. *Molecular Ecology*, 23, 1014-1027.
doi:10.1111/mec.12656

Vos, C. C., Antonisse-De Jong, A. G., Goedhart, P. W., & Smulders, M. J. M. (2001). Genetic similarity as a measure for connectivity between fragmented populations of the moor frog (*Rana arvalis*). *Heredity*, 86, 598-608. doi:10.1046/j.1365-2540.2001.00865.x

Wallis, G. P., Waters, J. M., Upton, P., & Craw, D. (2016). Transverse alpine speciation driven by glaciation. *Trends in Ecology & Evolution*, 31(12), 916-926.
doi:10.1016/j.tree.2016.08.009

Werth, S., Gugerli, F., Holderegger, R., Wagner, H. H., Csencsics, D., & Scheidegger, C. (2007). Landscape-level gene flow in *Lobaria pulmonaria*, an epiphytic lichen. *Molecular Ecology*, 16(13), 2807-2815. doi:10.1111/j.1365-294X.2007.03344.x

Zellmer, A., & Knowles, L. L. (2009). Disentangling the effects of historic vs. contemporary landscape structure on population genetic divergence. *Molecular Ecology*, 18(17), 3593-3602. doi:10.1111/j.1365-294X.2009.04305.x
## Data Accessibility

Newly generated mtDNA sequence data were deposited in GenBank under accession numbers MT498344-MT498349. Original sequence alignments and microsatellite genotypes were deposited in Dryad (Lucati et al., 2020).

## Author Contributions

F.L., M.P., A.M., A.T. and M.V. conceived and designed the study. F.L., M.P., A.M., A.T., L.B., R.B., O.C., E.D., M.D., H.L.C., A.M.S., M.M.T., D.O'B., G.P., J.S. and J.T. collected the samples. F.L., M.P., A.T., J.C., M.R. and I.S. analysed samples and data, under the supervision of M.V.. F.L. and M.P. wrote the first draft, A.M., A.T., R.B., T.B., O.C., M.D., A.M.S., D.O'B., V.O., I.S., J.S. and M.V. improved successive versions. All authors read and approved the final manuscript.
# Tables and Figures

## Tables

<table><caption>**Table 1** Genetic diversity parameters for each genetic cluster identified by STRUCTURE analysis in *Calotriton asper*.</caption><thead><tr><th>Cluster</th><th>N</th><th>Ar</th><th>PAAr</th><th>H<sub>O</sub></th><th>H<sub>E</sub></th><th>F<sub>IS</sub></th></tr></thead><tbody><tr><td>1</td><td>470</td><td>8.120</td><td>0.530</td><td>0.484</td><td>0.647</td><td>0.253</td></tr><tr><td>2</td><td>129</td><td>7.540</td><td>0.320</td><td>0.389</td><td>0.552</td><td>0.298</td></tr><tr><td>3</td><td>160</td><td>7.680</td><td>0.220</td><td>0.369</td><td>0.626</td><td>0.414</td></tr><tr><td>4</td><td>259</td><td>10.240</td><td>1.440</td><td>0.460</td><td>0.734</td><td>0.375</td></tr><tr><td>5</td><td>281</td><td>7.690</td><td>0.410</td><td>0.422</td><td>0.633</td><td>0.335</td></tr></tbody></table>

Abbreviations: N, sample size; Ar, allelic richness standardized for sample size; PAAr, rarefied private allelic richness standardized for sample size; H<sub>O</sub>, observed heterozygosity; H<sub>E</sub>, expected heterozygosity; F<sub>IS</sub>, inbreeding coefficient.

<table><caption>**Table 2** Analysis of molecular variance (AMOVA) for *Calotriton asper* at the Pyrenean scale. Two hierarchical structures were tested: (1) among clusters identified by STRUCTURE analysis and among tributary valleys within clusters, (2) among tributary valleys and among populations within valleys.</caption><thead><tr><th rowspan="2">Source of variation</th><th rowspan="2">d.f.</th><th rowspan="2">SS</th><th colspan="2">Variance %</th><th rowspan="2">Fixation indices</th></tr><tr><th>component</th><th>Variation</th></tr></thead><tbody><tr><td colspan="5">(1)</td></tr><tr><td>Among clusters</td><td>4</td><td>2723.166</td><td>1.368</td><td>19.195</td><td>F<sub>CT</sub> = 0.192***</td></tr><tr><td>Among valleys within clusters</td><td>17</td><td>998.617</td><td>0.961</td><td>13.490</td><td>F<sub>SC</sub> = 0.167***</td></tr><tr><td>Within valleys</td><td>1928</td><td>9040.858</td><td>4.797</td><td>67.315</td><td>F<sub>ST</sub> = 0.327***</td></tr><tr><td colspan="5">(2)</td></tr><tr><td>Among valleys</td><td>20</td><td>3875.541</td><td>1.266</td><td>20.48</td><td>F<sub>CT</sub> = 0.205***</td></tr><tr><td>Among populations within valleys</td><td>48</td><td>2372.971</td><td>1.295</td><td>20.96</td><td>F<sub>SC</sub> = 0.264***</td></tr></tbody></table>
<table><tr><td>Among individuals within populations</td><td>1169</td><td>4479.044</td><td>0.212</td><td>3.43</td><td>F<sub>IS</sub> = 0.059***</td></tr><tr><td>Within individuals</td><td>1238</td><td>4218.500</td><td>3.408</td><td>55.13</td><td>F<sub>IT</sub> = 0.449***</td></tr></table>

Abbreviations: d.f., degrees of freedom; SS, sum of squares; F<sub>CT</sub>, fixation index among groups; F<sub>SC</sub>, fixation index among populations within groups; F<sub>ST</sub>, fixation index within populations; F<sub>IS</sub>, fixation index among individuals within populations; F<sub>IT</sub>, fixation index within individuals.

*** P < 0.001

<table><caption>**Table 3** Posterior probability of tested scenarios and 95% confidence intervals (CI) estimated with DIYABC analysis when considering only microsatellites and when including both mtDNA (cyt-b) and microsatellite markers. Type I and II errors for the best supported scenario are indicated. See Figure 2 for more information on the tested scenarios.</caption><thead><tr><th rowspan="2">Scenario</th><th colspan="3">Microsatellites</th><th colspan="3">Microsatellites + cyt-b</th></tr><tr><th>Posterior probability</th><th>95% CI</th><th>Type I error</th><th>Type II error</th><th>Posterior probability</th><th>95% CI</th><th>Type I error</th><th>Type II error</th></tr></thead><tbody><tr><td>1</td><td>0.002</td><td>0.002-0.003</td><td></td><td></td><td>0.007</td><td>0.006-0.008</td><td></td><td></td></tr><tr><td>2</td><td>0.002</td><td>0.002-0.003</td><td></td><td></td><td>0.010</td><td>0.009-0.012</td><td></td><td></td></tr><tr><td>3</td><td>0.996</td><td>0.995-0.996</td><td>0.034</td><td>0.041</td><td>0.983</td><td>0.980-0.985</td><td>0.039</td><td>0.029</td></tr></tbody></table>

<caption>**Table 4** Posterior parameters (median and 95% confidence intervals) and RMedAD (Relative Median Absolute Deviation) estimated with DIYABC analysis for the best supported scenario (scenario 3) when considering only microsatellites (simple sequence repeats – SSRs) and when</caption>
including both mtDNA (cyt-b) and microsatellite markers. See Figures 2 and S1 for more information on the tested scenarios.

<table><thead><tr><th rowspan="2">Parameter</th><th colspan="4">Microsatellites</th><th colspan="4">Microsatellites + cyt-b</th></tr><tr><th>Median</th><th>Q<sub>2.5</sub></th><th>Q<sub>97.5</sub></th><th>RMedAD</th><th>Median</th><th>Q<sub>2.5</sub></th><th>Q<sub>97.5</sub></th><th>RMedAD</th></tr></thead><tbody><tr><td>N<sub>1</sub></td><td>3 460</td><td>1 310</td><td>9 940</td><td>0.197</td><td>3 000</td><td>962</td><td>9 930</td><td>0.225</td></tr><tr><td>N<sub>2</sub></td><td>4 600</td><td>2 470</td><td>7 950</td><td>0.186</td><td>3 790</td><td>1 550</td><td>8 690</td><td>0.195</td></tr><tr><td>N<sub>3</sub></td><td>6 380</td><td>2 790</td><td>12 600</td><td>0.175</td><td>5 940</td><td>2 160</td><td>12 800</td><td>0.195</td></tr><tr><td>N<sub>4</sub></td><td>10 300</td><td>5 680</td><td>14 200</td><td>0.135</td><td>9 620</td><td>4 730</td><td>14 200</td><td>0.154</td></tr><tr><td>N<sub>5</sub></td><td>6 930</td><td>3 110</td><td>12 900</td><td>0.183</td><td>7 100</td><td>2 400</td><td>13 700</td><td>0.207</td></tr><tr><td>N<sub>135</sub></td><td>7 590</td><td>997</td><td>14 400</td><td>0.324</td><td>6 490</td><td>667</td><td>14 200</td><td>0.314</td></tr><tr><td>N<sub>241</sub></td><td>14 400</td><td>2 680</td><td>19 700</td><td>0.269</td><td>13 200</td><td>1 980</td><td>19 500</td><td>0.307</td></tr><tr><td>t<sub>1</sub></td><td>4 050</td><td>1 470</td><td>7 770</td><td>0.245</td><td>2 700</td><td>636</td><td>6 470</td><td>0.288</td></tr><tr><td>t<sub>2</sub></td><td>5 020</td><td>1 510</td><td>9 560</td><td>0.209</td><td>3 680</td><td>860</td><td>9 200</td><td>0.255</td></tr><tr><td>t<sub>3</sub></td><td>14 200</td><td>6 730</td><td>19 600</td><td>0.169</td><td>12 000</td><td>4 620</td><td>19 300</td><td>0.219</td></tr><tr><td>Mean μ<sub>(SSRs)</sub></td><td>1.32x10<sup>-4</sup></td><td>1.01x10<sup>-4</sup></td><td>2.45x10<sup>-4</sup></td><td>0.278</td><td>1.59x10<sup>-4</sup></td><td>1.06x10<sup>-4</sup></td><td>3.37x10<sup>-4</sup></td><td>0.253</td></tr><tr><td>Mean P<sub>(SSRs)</sub></td><td>0.229</td><td>0.122</td><td>0.300</td><td>0.188</td><td>0.195</td><td>0.110</td><td>0.291</td><td>0.180</td></tr><tr><td>Mean μ<sub>(cyt-b)</sub></td><td>-</td><td>-</td><td>-</td><td>-</td><td>1.69x10<sup>-7</sup></td><td>6.08x10<sup>-8</sup></td><td>4.00x10<sup>-7</sup></td><td>0.276</td></tr><tr><td>Mean k1<sub>(cyt-b)</sub></td><td>-</td><td>-</td><td>-</td><td>-</td><td>7.920</td><td>0.410</td><td>18.800</td><td>0.442</td></tr></tbody></table>

Abbreviations: N, effective population size for each analysed deme (1 – cluster 1; 2 – cluster 2; 3 – cluster 3; 4 – cluster 4, 5 – cluster 5; 135 – central clusters; 241 – three oldest glacial refugia: eastern, western and central); t, time of events in generations (t<sub>1</sub> – time to the most recent split; t<sub>2</sub> – time to the intermediate split; t<sub>3</sub> – time to the most ancient split); mean μ, mean mutation rate; mean P, mean coefficient P; mean k1, mean coefficient k1; Q<sub>2.5</sub>, quantile 2.5%; Q<sub>97.5</sub>, quantile 97.5%.
<figure><img src="image_2.png" /><figcaption>**Figure 1** Results of the Bayesian clustering analysis across *Calotriton asper* distribution range. Panel (a) shows the geographic distribution of the five genetic clusters identified by STRUCTURE. Sampled populations are represented by pie charts highlighting the population cluster membership obtained in STRUCTURE. Panel (b) shows STRUCTURE barplot of membership assignment for K = 5. Each individual is represented by a vertical bar corresponding to the sum of assignment probabilities to the K cluster. White lines separate populations. Panel (c) represents a neighbour-joining tree based on net nucleotide distances among clusters inferred by STRUCTURE. For population codes see Table S1.</figcaption></figure>
<figure><img src="image_3.png" /></figure>

**Figure 2** Phylogeographic scenarios tested in DIYABC during phase 2 (a). The most likely scenario, namely number 3, with the estimated time points (t1 – t3) of each split is shown in panel (b). More information on tested scenarios, estimated parameters and respective priors is given in Table 4 and S2 and in Figure S1.

<figure><img src="image_4.png" /></figure>

**Figure 3** Spatial interpolation of allelic richness (Ar; a) and expected heterozygosity (HE; b) among populations of *Calotriton asper*. Black dots denote sampling localities and black lines delimit the five genetic clusters inferred by STRUCTURE. Each cluster is identified with its corresponding number. Only populations with five or more genotyped individuals were considered in the analysis. Population codes are given in Figure 1.
<figure><img src="image_5.png" /><figcaption>**Figure 4** Neighbour-joining tree over all *Calotriton asper* populations based on D<sub>A</sub> distances. Branch colours delineate the five genetic clusters identified by STRUCTURE analysis (blue: cluster 1, light green: cluster 2, orange: cluster 3, dark green: cluster 4, pink: cluster 5), while population code colours correspond to the distinct habitat types (blue: streams, red: lakes, green: caves). See Table S1 for population codes.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5** Chord diagram tracking first generation migrants flows between *Calotriton asper* sampled populations as inferred by GeneClass2. Chord size is proportional to the number of migrants detected and arrows indicate the direction of migration. Colours delineate the five genetic clusters identified by STRUCTURE analysis (blue: cluster 1, light green: cluster 2, orange: cluster 3, dark green: cluster 4, pink: cluster 5). In the outer ring, populations belonging to the same glacial cirque or valley are connected together. Only populations where first generation migrants with known source locality were detected are shown. For population codes see Table S1.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>
**Figure 6** Partial effects of environmental (habitat type) and geographic (latitude, longitude and altitude) variables on allelic richness (Ar). (a), all populations; (b), cluster 1; (c), cluster 2; (d), cluster 4. Only variables that had a significant effect on Ar as determined by linear models selection are drawn. Latitude and longitude are in UTM coordinates and altitude is expressed in meters.