Molecular Ecology (2013) 22, 3261–3278

doi: 10.1111/mec.12310

Tracking climate change in a dispersal-limited species: reduced spatial and genetic connectivity in a montane salamander
, * ¶ J . L . P A R R A , † G . P A R R A - O L E A ‡ and K . R . Z A M U D I O * G . V E L O - A N T ON *Department of Ecology and Evolutionary Biology, Cornell University, Ithaca, NY 14853, USA, † Facultad de Ciencias Exactas y Naturales, Instituto de Biolog
ıa, Universidad de Antioquia, Medell
ın, Colombia, ‡Departamento de Zoologı´a, Instituto de Biolog
ıa, Universidad Nacional Aut
onoma de M
exico, Distrito Federal 04510, M
exico

Abstract Tropical montane taxa are often locally adapted to very specific climatic conditions, contributing to their lower dispersal potential across complex landscapes. Climate and landscape features in montane regions affect population genetic structure in predictable ways, yet few empirical studies quantify the effects of both factors in shaping genetic structure of montane-adapted taxa. Here, we considered temporal and spatial variability in climate to explain contemporary genetic differentiation between populations of the montane salamander, Pseudoeurycea leprosa. Specifically, we used ecological niche modelling (ENM) and measured spatial connectivity and gene flow (using both mtDNA and microsatellite markers) across extant populations of P. leprosa in the Trans-Mexican Volcanic Belt (TVB). Our results indicate significant spatial and genetic isolation among populations, but we cannot distinguish between isolation by distance over time or current landscape barriers as mechanisms shaping population genetic divergences. Combining ecological niche modelling, spatial connectivity analyses, and historical and contemporary genetic signatures from different classes of genetic markers allows for inference of historical evolutionary processes and predictions of the impacts future climate change will have on the genetic diversity of montane taxa with low dispersal rates. Pseudoeurycea leprosa is one montane species among many endemic to this region and thus is a case study for the continued persistence of spatially and genetically isolated populations in the highly biodiverse TVB of central Mexico. Keywords: climate change, connectivity, dispersal, genetic structure, landscape genetics, Pseudoeurycea leprosa, Trans-Mexican Volcanic Belt Received 18 October 2012; revision received 18 February 2013; accepted 26 February 2013

Introduction Variation in climate and landscape features define the degree of gene flow among natural populations, shape species’ distributions (Slatkin 1987), and in some cases, threaten current biodiversity (Pimm 2008). Climate and landscape play a role in population divergence;

Correspondence: Guillermo Velo-Ant
on, Fax: + 351 252661780; E-mail: [email protected] ¶Present address: CIBIO/InBIO - Centro de Investigacß~ao em Biodiversidade e Recursos Gen
eticos da Universidade do Porto, Instituto de Ci^encias Agr
arias de Vair~ao, R. Padre Armando Quintas, 11, 4485-661, Vair~ao, Portugal © 2013 John Wiley & Sons Ltd

however, species-specific life history and demography ultimately determine the degree of differentiation and spatial structure among populations (Mart
ınez-Solano et al. 2007; Bell et al. 2010; Safner et al. 2011). Multiple studies in comparative phylogeography show that repeated historical cycles in climate forced species to retract or spread and that species responded differently to these climatic changes: low-elevation temperate species retreated and lost habitat during glacial periods (e.g. Taberlet et al. 1998; Avise 2000), whereas highelevation species expanded their distributional ranges and increased gene flow among populations with the expansion of suitable habitats (e.g. Hewitt 2000; Galbreath et al. 2009; Parra-Olea et al. 2012; Devitt et al.


ET AL. 3262 G . V E L O - A N T ON 2013). Currently, rapid anthropogenic climate change is causing shifts in species distributions globally and has the potential to increase their vulnerability due to abrupt changes in ecological conditions (Root et al. 2003; Parmesan 2006; Sherry et al. 2007; Bellard et al. 2012; Devictor et al. 2012). Therefore, understanding how species-specific traits or attributes limit dispersal and contribute to genetic population structure is important for predicting how species will respond to future anthropogenic changes in climate and landscape. How species respond to climate fluctuations will vary depending on the direction and extent of climate change relative to their evolved ecological niche, and whether that niche is evolutionarily conserved or not. High biodiversity in montane regions has been linked to niche conservatism of montane biota (Kozak & Wiens 2010), which can retain niche characteristics over hundreds or thousands of years (Peterson 2011). In addition to the potential for adaptation to new environments, individual dispersal capacity also determines how species will react to the challenge of shifting optimal environments (Massot et al. 2008). Species that are highly specialized to specific climatic zones and have low vagility will be less capable of shifting in geographical distribution in response to climate change and habitat loss. Thus, if they are to persist, these species either must adapt to variable environments or survive only in those areas where climate remains stable through time (Hugall et al. 2002; Visser 2008). Tropical montane species are particularly sensitive to climate fluctuations because their narrow thermal tolerances typically restrict their ability to persist outside of optimal habitats to which they are adapted (Janzen 1967; Parmesan 2006; Colwell et al. 2008; Bell et al. 2010). Global warming forces montane species towards higher elevations, in some cases causing drastic changes in the extent of available habitat, and increasing their probability of extinction (Williams et al. 2003; Thomas et al. 2004). Montane amphibians in particular are often dependent on specific moisture and temperature regimes, and thus, global warming is one of the most important threats to this group (Kiesecker et al. 2001). Obviously, future changes in species’ distributions do not depend only on climate; landscape features such as topographic complexity and anthropogenic land-use can influence potential dispersal pathways (Feeley & Silman 2010; Schloss et al. 2012). As with climate, landscape changes can potentially make some regions uninhabitable or unreachable, creating patchiness in the distribution of optimal habitats. These discontinuities in distributions then lead to a reduction in gene flow among fragments (Fischer & Lindenmayer 2007), genetic isolation, potential decreases in neutral and functional genetic diversity (Reed & Frankham 2003)

and reduced colonization success after local extinctions (With & King 1999; Ferraz et al. 2007). Here, we quantify the influence of both climate and landscape characteristics on the distribution, and population connectivity of the plethodontid salamander Pseudoeurycea leprosa, a montane species endemic to oak and pine forests between 2000 and 3500 m in the TransMexican Volcanic Belt (TVB) of central Mexico (Fig. 1). The TVB is a biodiversity hotspot characterized by relatively young volcanic activity (maximum age estimated at 3.7 Mya; Ferrari et al. 2000) that gave rise to geographically distinct mountain chains and promoted diversification in TVB taxa (e.g. Bryson & Riddle 2012; Bryson et al. 2011; Parra-Olea et al. 2012). The high topographic complexity of this region makes the TVB an excellent study area to track the effects of climate change on distributions and genetic connectivity among populations after the last glacial maximum (LGM). Pseudoeurycea leprosa is particularly well suited for this study because it occurs only at high elevations (Parra-Olea et al. 2012) and belongs to a phylogenetic group characterized by low rates of dispersal (Mart
ınez-Solano et al. 2007). Our earlier mtDNA phylogeographic study of this species showed that isolation in multiple refugia in north-eastern and eastern mountain ranges played an important role in structuring historical population genetic diversity in the TVB, followed by population expansion and genetic divergence of populations in the central and western part of the species’ range (Parra-Olea et al. 2012). These data indicate that colder periods during the LGM may have promoted increased connectivity among mountain ranges permitting expansion of P. leprosa into the westernmost sections of its distribution. Landscape genetic studies integrate geospatial and genetic statistics to analyse factors that alter the distribution of genetic variation within and among populations (Wang et al. 2008; Goldberg & Waits 2010; Murphy et al. 2010; Row et al. 2010; Chan et al. 2011). Our goal is to understand how climate and landscape influence the contemporary genetic structure of natural populations by testing the following hypotheses: (i) the spatial connectivity among P. leprosa populations decreased since the LGM as populations shifted their distribution upslope in response to warmer climates; (ii) the reduction in spatial connectivity due to historical climate warming and habitat fragmentation, together with limited dispersal ability, led to reduced gene flow in these terrestrial salamanders. We integrate genetic signatures from historical (mtDNA) and contemporary (microsatellites) markers, niche modelling, and spatial connectivity analyses to investigate how climate, topography and land cover explain patterns of genetic structure in P. leprosa across the TVB. We use our findings to © 2013 John Wiley & Sons Ltd

T R A C K I N G C L I M A T E C H A N G E I N A M O N T A N E S A L A M A N D E R 3263

TVB

4 2

Perote

11 10

3

1

5

7 8

19

Tlatlauquitepec

12

Sierra de las Cruces

Nevado de Toluca

TVB

20

22

15 14

9

13

N

21 0

16

Malinche

6

18 Orizaba

Popocatepetl Iztaccihuatl (Sierra Nevada)

Central clade

0–1000 1000–1500 1500–2000 2000–2500 2500–3000 3000–3500 3500–4000 4000–4500 4500–5500

19. Vigas

21. G. Ortega

17. Texmalaquilla 18. Xometla

23. Tres Mogotes

20. Tlatlauquitepec

22. Tlaxco

16. Malinche E

15. MAlinche N

14. Malinche W

13. Malinche S

7. Popocatepetl 8. Calpan 9. Atzompa 10. Llano Gradce 11. Río Frío 12. Nanacamilpa

3. D. Leones 4. Ajusco 5. Santa Marta 6. Zempoala

2. Texcalyacac

1. Nevado

50

Elevation (m)

17

23

1–16

25

km

Eastern clades 17–18

19, 20 23 20, 22

19, 21

Fig. 1 Upper panel: sampling localities throughout the entire known distributional range of Pseudeurycea leprosa. Lower panel: Bayesian phylogenetic tree for P. leprosa populations (adapted from Parra-Olea et al. 2012). The Central Trans-Mexican Volcanic Belt (TVB) clade includes recently diverged populations following a rapid population expansion from eastern source populations.

predict the potential distribution of suitable environment for P. leprosa in the future (2050) and discuss the implications of reduced spatial and genetic connectivity driven by rapid climate warming for the highly diverse montane fauna of the TVB.

Materials and methods Genetic data collection and analyses of population genetic diversity We obtained tissue samples from 351 Pseudoeurycea leprosa from 23 localities (mean sample size = 15; Fig. 1, Tables 1 and S1, Supporting information). We extracted genomic DNA from tissue samples using QIAGEN extraction kits (Qiagen Inc., Valencia, CA, USA) for use as templates in PCR amplification of 11 microsatellite loci (Pld005, Plt009, Plt28, Plt39, Plt42, Plt45, Plt64, Plt66, Plt87, Plt107, Plt109; Velo-Ant
on et al. 2009). PCR amplifications were performed in 15 lL reaction volumes including 1 lL of DNA template (~100 ng/lL), 0.5 U Taq polymerase (Applied Biosystems), 1 9 PCR buffer with MgCl2, 0.4 mM dNTPs, and 0.5 lM forward (5′-labelled with VIC, PET, 6-FAM, NED) and reverse primers. PCR cycling © 2013 John Wiley & Sons Ltd

conditions consisted of an initial denaturation at 94° (5 min), followed by 30 cycles of 94° (1 min), annealing at 49–62 °C (1 min) and extension at 72° (1 min), followed by a final elongation at 72° (10 min). Loci were multiplexed and fragment sizes determined using LIZ500 or LIZ-600 standards. Results were scored and binned using GENEMAPPER v3.7 (Applied Biosystems). We used MICRO-CHECKER v2.2.3 to test for null alleles and/or large allele dropout (Van Oosterhaut et al. 2004). We tested for significant deviations from Hardy– Weinberg equilibrium and the presence of linkage disequilibrium in ARLEQUIN v3.11 (Excoffier et al. 2005), using a Markov chain method with 100 000 dememorization steps, 1000 batches of 10 000 iterations per batch and sequential Bonferroni correction (Rice 1989). We used GENEALEX v6 (Peakall & Smouse 2006) to calculate genetic diversity values (NA, average number of alleles; HE, expected heterozygosity; AR, allelic richness; P-AR, private allelic richness). Populations with less than eight samples were excluded from these analyses (Table S2, Supporting information). To correct for sample size differences, we used the program HP-RARE v1.0 (Kalinowski 2005), for rarefied estimates of allelic richness and private allelic richness across populations.


ET AL. 3264 G . V E L O - A N T ON Table 1 Pseudoeurycea leprosa populations from the Trans-Mexican Volcanic Belt (TVB) included in this study. For each locality, we list the sample size (N), mountain chain, region within the TVB, and GPS coordinates of collection sites. Population ID numbers correspond to those in Figs 1 and 2 ID

Locality

N

Mountain chain

Region

Latitude

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

Nevado de Toluca Texcalyacac Desierto de los Leones Ajusco Santa Marta Zempoala Popocatepelt Calpan Atzompa Llano Grande Rio Frio Nanacamilpa MalincheS MalincheW MalincheN MalincheE Texmalaquilla Xometla Vigas Tlatlauquitepec Tlaxco Tres Mogotes Gonz
alez Ortega

12 29 7 8 3 5 30 34 43 15 27 14 8 17 5 13 8 16 48 4 1 1 3

Nevado de Toluca Sierra de las Cruces Sierra de las Cruces Sierra de las Cruces Sierra de las Cruces Sierra de las Cruces Sierra Nevada Sierra Nevada Sierra Nevada Sierra Nevada (Popocatepetl-Iztaccihuatl) Sierra Nevada (Popocatepetl-Iztaccihuatl) Sierra Nevada (Popocatepetl-Iztaccihuatl) Malinche Malinche Malinche Malinche Pico de Orizaba Pico de Orizaba Cofre de Perote TVB TVB Tres Mogotes Cofre de Perote

Central West Central West Central West Central West Central West Central West Central Central Central Central Central Central Central East Central East Central East Central East Southeast Southeast Northeast Northeast Northeast Southeast Northeast

19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 18 18 19 19 19 18 19

Population genetic structure We inferred population structure based on Bayesian clustering of multilocus genotypes in STRUCTURE v2.2 (Pritchard et al. 2000), with 20 runs for numbers of clusters, K, ranging from 1 to 20. Twenty was the total number of populations (those with large sample sizes) included in the analysis, and therefore, this maximum number allows for the possibility that all samples populations belong to independent demes. We used a burn-in period of 500 000 and a run length of 5 000 000 iterations, with a correlated allele frequencies admixture model and without prior information on population origin. The number of populations (K) with the highest likelihood was estimated by the highest value of log probability data L(K) (Falush et al. 2003) and DK (Evanno et al. 2005). We used CLUMPP v1.1.2 (Jakobsson & Rosenberg 2007) to average the outcomes from the STRUCTURE runs with the highest L(K) value. We estimated two indices of genetic differentiation among populations: pairwise FST values in ARLEQUIN v3.11 (Excoffier et al. 2005) and the harmonic mean of Jost’s DEST (Jost 2008) in SMOGD v1.2.5 (Crawford 2010). For both estimates, we excluded populations with low sample size (Santa Marta, Gonz
alez Ortega, Tlatlauquitepec, Tlaxco and Tres Mogotes for DEST; Tlaxco and Tres Mogotes for FST; Table 2). For both indices, we

11 07 16 10 03 03 04 07 10 20 21 28 11 15 16 13 56 58 37 42 41 39 21

37 14 00 58 41 12 22 53 50 20 58 49 14 28 37 48 31 30 51 47 45 28 02

Longitude N N N N N N N N N N N N N N N N N N N N N N N

99 99 99 99 99 99 98 98 98 98 98 98 98 98 98 97 97 97 97 97 98 97 97

50 30 18 18 19 18 42 35 33 43 41 35 01 05 02 58 17 11 05 32 04 29 14

53 00 02 00 17 35 42 30 35 14 41 46 19 42 40 30 24 26 26 24 33 24 49

W W W W W W W W W W W W W W W W W W W W W W W

tested significance of pairwise comparisons using an exact test, bootstrapping with 1000 iterations, and P-values corrected by sequential Bonferroni (Rice 1989).

Migration rates To estimate migration rates for P. leprosa among mountain ranges, we used a coalescent-based Monte Carlo Markov Chain (MCMC) method implemented in MIGRATE 3.0 (Beerli & Felsenstein 1999). We compared our results to previously published mitochondrial results for the same individuals (Parra-Olea et al. 2012) and used the two independent data sets to estimate historical (mtDNA) and contemporary (microsatellites) migration rates among mountain chains. We compared groups by randomly selecting 15–20 samples per mountain range (the number of samples recommended by the authors, http://popgen.sc.fsu.edu/migratedoc.pdf). In the few cases, where populations within one mountain range showed high divergence, we subsampled populations independently, so that our inferred immigration rates would not be biased by pooling mixed genetic demes. We analysed the mtDNA and microsatellite data sets under F84 and Brownian approximation models, respectively. To optimize parameter estimates, we selected 15 representative populations that correspond to different genetic demes (Table 3), excluding © 2013 John Wiley & Sons Ltd

T R A C K I N G C L I M A T E C H A N G E I N A M O N T A N E S A L A M A N D E R 3265 localities with fewer than three individuals (Tlaxco and Tres Mogotes). We performed five independent runs using one long chain with a run of 1 9 107 recorded genealogies sampled every 1000 steps and a burn-in of 1 9 106. MIGRATE infers approximate confidence intervals (CIs) around the maximum probable estimate (MPE) for each parameter. Parameter convergence was verified by examining the stationarity in parameter trends over the length of the chains and the Effective Sample Sizes (ESS) obtained for each migration rate using TRACER v1.5; ESS values greater than 300 indicated that the sampled trees were not correlated and thus represented independent samples. To calculate more recent migration rates (over the last few generations), we estimated m (proportion of migrants) among populations of P. leprosa using BAYESASS v1.3 (Wilson & Rannala 2003). The Markov Chain Monte Carlo (MCMC) run included 5 9 106 iterations sampled every 2000 steps, with the first 106 iterations discarded as burn-in. We repeated runs five times to check for convergence of the chains and compared migration rate means and 95% confidence intervals with those obtained when there is no information in the data (significant estimates result when the 95% confidence intervals do not overlap). Based on preliminary runs using different delta values (i.e. maximum parameter change per iteration), we used delta values of 0.25, 0.30 and 0.30 for allele frequency, migration and inbreeding parameter estimates, respectively.

Ecological niche modelling To identify the potential distribution of P. leprosa at different time periods, we first developed an ecological niche model (ENM) based on current climatic data using the maximum entropy algorithm implemented in MAXENT v3.3.3 (Phillips et al. 2006) and geo-referenced points from museum specimens and our own collections (73 unique localities at 1 km2 spatial resolution). We set our study region to include all present populations of P. leprosa and surrounding mountainous areas where the species could potentially persist (Fig. 3). We first tested for correlation among the 17 bioclimatic layers (WorldClim http://www.worldclim.org/) using R v2.9.2 and chose uncorrelated climatic variables (<0.7 Pearson correlation coefficient: annual mean temperature, temperature seasonality, temperature annual range, precipitation seasonality and precipitation of driest quarter) to model the current distribution at a scale of 1 km2 resolution. To identify stable climatic areas that possibly served as montane refugia, we projected the current ENM to two paleoclimate scenarios from the last glacial maximum (LGM, 21 000 years B.P.) characterized by cooler and drier climates relative to current conditions: © 2013 John Wiley & Sons Ltd

CCSM3 (Community Climate System Model, available at http://www.ccsm.ucar.edu), and Model of Interdisciplinary Research on Climate (MIROC, available at http://www.ccsr.utokyo.ac.jp/kyosei/hasumi/MIROC/ tech-repo.pdf). Finally, we projected the current distribution model to one future scenario (2050; IPPC 3rd Assessment data). These three sets of climate surfaces were all downscaled using the same current layers (1950–2000) offered by the WorldClim database. We checked the current model accuracy by using a random sample of 25% of the data set as training data and the remaining 75% as testing data for model validation. The resulting model was evaluated with the ‘area under the curve’ (AUC) of the Receiver Operating Characteristics Curve (ROC) (Fielding & Bell 1997; Manel et al. 2001). The AUC ranges from 0.5 (random accuracy) to a maximum value of 1.0 (perfect discrimination). In addition, we developed two independent ENMs for the present period based only on land cover and topographic complexity surfaces using MAXENT (under the same parameters used for the climatic data) to estimate resistance distances among populations. We used a 1 km2 resolution digital land cover map of Mexico (Instituto Nacional de Estad
ıstica y Geograf
ıa, INEGI, Mexico) that contains 50 detailed land classes (e.g. urban areas, forests, agriculture and grassland). To quantify topographic complexity, we used a 250 m resolution digital elevation model from the Shuttle Radar Mission and for each nonoverlapping 1 km2 pixel, we calculated the standard deviation in elevation using the Block Statistics function in the Neighbourhood toolbox in ARCMAP 9.3. The distribution models based on climatic, land cover and topographic complexity were used to generate spatial resistance surfaces (an index of spatial connectivity) to investigate the role of these factors in shaping contemporary dispersal and potential gene flow in P. leprosa.

Spatial connectivity and spatial-genetic correlations Each ENM map was used as input for the program CIR3.1 (McRae 2006) to generate spatial resistance distances to dispersal among sites based on habitat suitability maps (conductance grid). CIRCUITSCAPE calculates pairwise resistances to gene flow among populations based on all possible paths, not just the least cost path, thus better explaining the movement of genes among widely separated regions over many generations (McRae & Beier 2007; McRae et al. 2008). The input for CIRCUITSCAPE is a raster data set (habitat map) in which each cell is assigned a conductance value corresponding to the probability of the study organism moving through the habitat type encoded by the cell. We chose a conductance grid in which higher cell values denote CUITSCAPE


ET AL. 3266 G . V E L O - A N T ON greater ease of movement and applied a connection scheme that allowed gene flow among the four nearest cells. We used the above ENM rasters as input maps to quantify pairwise resistance distances between P. leprosa populations across its distribution and within populations of the Central TVB clade. To evaluate the influence of specific climate or landscape factors on P. leprosa connectivity and genetic structure, we established 12 models that differ in their combination of landscape features (Table 4). First, we calculated pairwise landscape resistances on a ‘flat’ resistance matrix (1) (all cells have equal conductance) and estimated isolation by landscape resistance (IBR) for comparison with other models that include climate or landscape variables. Although we expect the IBR to yield similar results as isolation by Euclidean distance (IBD), the IBR model accounts for a finite size of the landscape and thus is more appropriate for comparison with competing models (Lee-Yaw et al. 2009). Then, we generated pairwise resistance distances based on habitat suitability maps obtained for climatic (2, LGM-CCSA; 3, LGM-MIROC and 4, Present), land cover (5) and topographic complexity projections (6) in MAXENT. Each raster contains habitat suitability scores derived from the appropriate ENM and exported as ASCII grid format by ArcGIS. We standardized the raw values obtained for the above ENMs before using them in CIRCUITSCAPE. To standardize raster maps, we first calculated the mean (l) and standard deviation (r) for each raster using ARCMAP and estimated Z-scores for each data point [Z = (X l)/r, where X are the values prior to normalization]. We then summed normalized values to generate mixed models (7, 8, 9 and 10). In addition, we generated two stability resistance models (11 and 12). We multiplied the continuous ENMs from both climatic times (LGM and Present) to generate stability maps for P. leprosa. Our measure of stability represents the probability of a lineage persisting from the LGM to the present in a single pixel. This measure assumes no dispersal and does not capture directional change in habitat suitability through time (i.e. whether it was more suitable in the past than in the present). We performed a separate CIRCUITSCAPE analysis for each raster, generating independent pairwise resistances distances as output for each model. To investigate the relative importance of climate and landscape features in shaping population genetic structure in P. leprosa, we conducted a series of Mantel tests (Mantel 1967) examining significant associations between microsatellite pairwise genetic distances (FST and DEST) and increasingly more complex models of pairwise landscape resistance distances obtained in CIRCUITSCAPE (Table 4). We used partial Mantel tests (Smouse et al. 1986) to examine the association between

pairwise genetic and resistance distances while controlling for Euclidean distance. We estimated Mantel and partial Mantel correlation coefficients (r) using ZT software (Bonnet & Van de Peer 2002), with 10 000 matrix permutations to estimate significance and 95% confidence intervals. Comparison of the r values allowed for the assessment of relative model performance. In addition, because our earlier studies indicated early isolation in north-eastern and eastern mountain ranges and posterior population expansions into the central and western ranges of the species’ range (Parra-Olea et al. 2012), we examined the relationships between genetic structure and connectivity among populations within the most recently diverged clade (Central TVB clade, Fig. 1). To estimate differences in spatial connectivity among populations over time, we compared the distribution of spatial resistances across the study area obtained from CIRCUITSCAPE for all combinations of populations using boxplots for past, present and future climatic projections.

Results Genetic data collection and analyses of population genetic diversity Two of the 11 microsatellites (Plt028 and Plt066) were excluded from all downstream analyses because they showed evidence of null alleles and HW disequilibrium in all populations. For the remaining 9 loci, significant departures from Hardy–Weinberg equilibrium occurred in 12 of 207 tests after sequential Bonferroni correction (adjusted P-value 0.00027). Nonetheless, these 12 deviations did not show any pattern across populations, loci or in directionality, and thus, the loci were retained in the final data set. No linkage disequilibrium was detected among loci. A summary of genetic diversity values (HO; HE; NA) is presented in Table S2 (Supporting information).

Population genetic structure Pseudoeurycea leprosa showed evidence of significant range wide population genetic structure. STRUCTURE inferred seven genetic demes (K = 7) that are geographically associated to mountain chains, with the largest genetic deme distributed among multiple central TVB mountains (Fig. 2). Among the older and more diverse eastern populations (Parra-Olea et al. 2012), we identified two main genetic demes with limited genetic admixture: a south-east genetic deme (Texmalaquilla and Xometla) that also includes a few animals genotyped from eastern populations (Tlaxco, Gonz
alez © 2013 John Wiley & Sons Ltd

T R A C K I N G C L I M A T E C H A N G E I N A M O N T A N E S A L A M A N D E R 3267 Ortega and Tres Mogotes); and an isolated deme found only in two north-east populations (Vigas and Tlatlauquitepec). Within the recently expanded Central TVB clade, we detected five genetic demes: two isolated genetic demes restricted to the western populations of Nevado de Toluca and Texcalyacac, a third widespread genetic deme that groups populations of Sierra de las Cruces and southern populations of Sierra Nevada (Popocatepetl-Iztaccihuatl), a fourth genetic deme clustering northern populations of Sierra Nevada (Llano Grande, Rio Frio and Nanacamilpa) and fifth genetic demes in an isolated central mountain (Malinche). We found evidence of limited admixture between the two eastern genetic demes, and between a south-east population (Xometla) and Central TVB populations from Malinche, a pattern that was observed in our earlier mtDNA study (Parra-Olea et al. 2012). Across the Central TVB clade, admixture is clearly evident between northern and southern populations of Sierra Nevada, and surprisingly, most individuals from the easternmost population at Malinche (Malinche E) are assigned to the widespread TVB genetic deme distributed along populations of Sierra de las Cruces and Sierra Nevada. Pairwise FST values between populations ranged from 0.01 (within Malinche populations) and 0.40 [between northeast (Tlatlauquitepec) and central-western (Ajusco) populations]. Similar genetic differentiation is observed for DEST values, ranging from 0.01 (within Malinche populations) and 0.85 [between north-east (Vigas) and central-western (Ajusco) populations]. All pairwise comparisons were significant after sequential Bonferroni correction (Table 2).

Migration rates Because the Texcalyacac population is located between two mountain chains and belongs to an independent genetic deme (Fig 2), we considered it an independent group in our migration analyses. The MIGRATE analyses inferred significant rates of dispersal from both mitochondrial and nuclear markers (Table 3). The most probable estimates (MPE) of migration using mitochondrial data indicate migration among several montane groups within the Central TVB (between Sierra de las Cruces and Texcalyacac, Popo-Izta, and Malinche; between Malinche and Popo-Izta and Nevado) and eastern clades (Orizaba/Perote); however, we found no evidence of migration between Central TVB and eastern clades (Table 3). In contrast, contemporary migration rates estimated with microsatellites show only two high migration rates between montane chains (Popo-Izta/ Malinche; Sierra de las Cruces/Popo-Izta). Although we obtained high ESS values for all migration rate estimates, most MPE 95% confidence intervals are broad © 2013 John Wiley & Sons Ltd

and overlap; thus, the estimates should be interpreted with caution (Table 3). BAYESASS analyses only detected significant contemporary migration (m) within mountain chains: Popocatepetl (28% immigrant ancestry from Atzompa); Calpan (19% immigrant ancestry from Atzompa); Llano Grande (21% immigrant ancestry from R
ıo Fr
ıo) and Malinche N (13% immigrant ancestry from Malinche W).

Ecological niche modelling and spatial connectivity The variable that most influenced the distribution of P. leprosa in all climatic models was annual mean temperature (Table S1, Supporting information). All models had high values of AUC (>0.90, Table S1, Supporting information) indicating overall adequate performance. Predictive maps are shown above a minimum probability threshold of 10% in Fig. 3. The MAXENT model for P. leprosa correctly predicted species occurrence under current conditions: a fragmented distribution with isolated high-elevation TVB areas (Fig. 3). The model also predicted the occurrence of the species in nearby montane regions, which could represent overprediction of the model, extinct populations or populations not yet detected in fieldwork (Fig. 3). Both the CCSA and MIROC projections for the cold period of the LGM predicted a wider distribution than observed currently for P. leprosa, although MIROC projection was more geographically restricted than the CCSA projection. These results indicate that P. leprosa populations during the LGM had a continuous distribution across the TVB and that populations were more connected than they are now (Fig. 3). The projection for the near future showed a severe and dramatic restriction of distribution for P. leprosa, with clear signs of fragmentation and low probability of suitable areas. While climatic spatial resistances surfaces showed connectivity corridors between neighbour mountain chains, resistance surfaces based on land cover and topographic complexity showed connectivity only within mountain chains (Fig. 4). Furthermore, pairwise resistance matrices obtained from CIRCUITSCAPE showed a continuous increase in spatial resistance between populations over time (Fig. 5) and thus a decrease in connectivity among populations.

Spatial-genetic correlations Overall, our results show significant but weak positive correlations between both population genetic distance indices (FST and DEST) and spatial resistance distances among populations over time (last glacial maximum and present time), regions of stability (potential refugia since the LGM) and landscape features (topographic


ET AL. 3268 G . V E L O - A N T ON Table 2 Genetic differentiation between populations. Pairwise FST (below the diagonal) and DEST values (above the diagonal). All P-values were significant (P < 0.01) after sequential Bonferroni correction. DEST is the harmonic mean across loci as calculated in SMOGD (Crawford 2010). NA represents the absence of analyses for DEST estimates due to the low samples size

FST/ DEST

Nevado

Texcalyacac

Leones

Ajusco

Santa Marta

Zempoala

Popocatepetl

Calpan

Atzompa

Llano Grande

Nevado Texcalyacac Leones Ajusco Santa Marta Zempoala Popocatepetl Calpan Atzompa Llano Grande R
ıo Fr
ıo Nanacamilpa Malinche S Malinche W Malinche N Malinche E Texmalaquilla Xometla Gonz
alez Ortega Vigas Tlatlauquitepec

* 0.18 0.20 0.28 0.22 0.21 0.14 0.16 0.19 0.19 0.26 0.21 0.25 0.27 0.21 0.29 0.20 0.17 0.32 0.29 0.36

0.40 * 0.12 0.12 0.19 0.17 0.08 0.07 0.10 0.16 0.20 0.16 0.19 0.19 0.18 0.19 0.18 0.14 0.32 0.25 0.34

0.38 0.18 * 0.16 0.17 0.13 0.07 0.06 0.10 0.14 0.17 0.16 0.22 0.21 0.18 0.22 0.14 0.13 0.29 0.22 0.39

0.43 0.25 0.26 * 0.28 0.17 0.11 0.12 0.12 0.22 0.27 0.19 0.24 0.25 0.27 0.17 0.23 0.17 0.42 0.30 0.46

NA NA NA NA * 0.15 0.08 0.08 0.07 0.07 0.14 0.17 0.18 0.19 0.20 0.23 0.22 0.12 0.37 0.23 0.44

0.34 0.26 0.33 0.25 NA * 0.06 0.08 0.07 0.11 0.16 0.11 0.12 0.16 0.18 0.12 0.16 0.12 0.30 0.22 0.38

0.33 0.25 0.20 0.27 NA 0.22 * 0.02 0.02 0.06 0.11 0.08 0.11 0.11 0.09 0.11 0.10 0.09 0.24 0.18 0.28

0.39 0.17 0.11 0.41 NA 0.03 0.06 * 0.04 0.06 0.12 0.09 0.12 0.11 0.09 0.10 0.13 0.09 0.25 0.18 0.28

0.50 0.29 0.29 0.29 NA 0.20 0.08 0.14 * 0.04 0.10 0.09 0.13 0.14 0.13 0.11 0.13 0.12 0.27 0.20 0.30

0.35 0.33 0.18 0.51 NA 0.26 0.10 0.13 0.07 * 0.03 0.14 0.12 0.13 0.10 0.17 0.18 0.13 0.29 0.21 0.31

complexity and land cover) (r = 0.25–0.60; P < 0.01; Table 4). Mantel test analyses showed significant patterns of IBR. Likewise, spatial resistance based on climatic and landscape features showed similar values and did not outperform the IBR model. Several models combining resistance surfaces from climatic (past and present times) and landscape features (land cover and topographic complexity), and stability models from climatic data, showed similar results and did not achieve higher associations between spatial resistance and genetic differentiation than those based on single features or IBR. The combined resistance model (land cover + topographic complexity) showed the highest association (r = 0.58–0.60; P < 0.001). Because the IBD pattern found in our data set could be influenced by deeper evolutionary divergences within P. leprosa (Parra-Olea et al. 2012) and limit detection of climatic and landscape effects on genetic structure, we also analysed the Central TVB populations alone (Central TVB mtDNA clade). For the central TVB populations, we observed a slight improvement in the correlation between genetic structure and spatial resistance for most of the models with pairwise FST, but not DEST. Partial Mantel tests were all nonsignificant (results not shown), and thus, none of the resistance surfaces based on climatic and landscape features explained further variation in genetic distance for Central TVB populations after controlling for Euclidean distance.

Discussion Montane regions are biodiversity hotspots that harbour a high number of endemic and often highly differentiated species (Myers et al. 2000). Tropical mountain species are highly susceptible to climate fluctuations and the predicted increases in temperature variance will disproportionately increase their probability of extinction (Williams et al. 2003, 2007; Deutsch et al. 2008; Huey et al. 2009; Cadena et al. 2012). In this study, the combination of ENM, spatial connectivity analyses and genetic markers with historical and contemporary signatures, enhanced our understanding of spatial connectivity loss, the main mechanism that shaped both the current distribution and genetic structure in Pseudoeurycea leprosa. Our results also highlight the susceptibility of tropical montane species with limited dispersal to ongoing global warming.

Spatial and genetic structure in Pseudoeurycea leprosa after the LGM In our previous mtDNA phylogeography of P. leprosa (Parra-Olea et al. 2012), we concluded that this species originated in the eastern TVB, and populations expanded from multiple refugia in north-eastern and eastern mountains to the central and western ranges during the Pliocene and Pleistocene (Fig. 1). While the © 2013 John Wiley & Sons Ltd

T R A C K I N G C L I M A T E C H A N G E I N A M O N T A N E S A L A M A N D E R 3269

R
ıo Fr
ıo

Nanacamilpa

Malinche S

Malinche W

Malinche N

Malinche E

Texmalaquilla

Xometla

Gonz
alez Ortega

Vigas

Tlatlauquitepec

0.49 0.38 0.20 0.45 NA 0.23 0.27 0.30 0.25 0.06 * 0.16 0.18 0.19 0.17 0.25 0.20 0.19 0.36 0.27 0.38

0.36 0.41 0.31 0.29 NA 0.19 0.17 0.21 0.17 0.37 0.35 * 0.22 0.24 0.22 0.22 0.19 0.17 0.34 0.24 0.37

0.64 0.35 0.44 0.28 NA 0.15 0.23 0.25 0.23 0.19 0.32 0.34 * 0.01 0.01 0.18 0.24 0.11 0.33 0.26 0.38

0.69 0.43 0.37 0.39 NA 0.19 0.25 0.28 0.34 0.25 0.29 0.53 0.02 * 0.02 0.17 0.24 0.13 0.34 0.26 0.36

0.39 0.44 0.39 0.25 NA 0.15 0.05 0.05 0.18 0.12 0.21 0.33 0.01 0.01 * 0.20 0.20 0.07 0.32 0.25 0.38

0.57 0.35 0.30 0.28 NA 0.22 0.24 0.20 0.18 0.24 0.48 0.36 0.28 0.23 0.17 * 0.23 0.14 0.34 0.26 0.40

0.46 0.52 0.22 0.35 NA 0.26 0.23 0.41 0.30 0.41 0.28 0.35 0.78 0.60 0.53 0.56 * 0.06 0.18 0.22 0.30

0.71 0.48 0.47 0.52 NA 0.43 0.29 0.27 0.43 0.37 0.45 0.54 0.34 0.33 0.23 0.47 0.31 * 0.09 0.17 0.20

NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA * 0.20 0.34

0.74 0.71 0.55 0.85 NA 0.60 0.54 0.51 0.61 0.54 0.72 0.57 0.69 0.78 0.47 0.60 0.57 0.53 NA * 0.20

NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA *

Table 3 Asymmetric immigration rates inferred in MIGRATE: maximum probable estimate (MPE) and 95% confidence intervals (CI) (in parentheses). Above the diagonal, numbers refer to migration rates into the montane groups listed at the top of the table. Below the diagonal, numbers are migration rates into the montane groups at the left. (A) Historical asymmetric migration rates inferred in MIGRATE for Pseudoeurycea leprosa using mtDNA. (B) Contemporary asymmetric migration rates inferred in MIGRATE using microsatellites. Estimated MPE values higher than 200 are in bold to show the main immigration rates obtained for both data sets

(A) Immigration Nevado Texcalyacac S. Cruces Popo-Izta Malinche Perote Orizaba (B) Immigration Nevado Texcalyacac S. Cruces Popo-Izta Malinche Perote Orizaba

Nevado

Texcalyacac

S. Cruces

Popo-Izta

Malinche

* 2.5 137.5 2.5 387.5 2.5 2.5

(0–140) (0–370) (0–325) (0–500) (0–315) (0–300)

2.5 * 42.5 77.5 2.5 2.5 42.5

2.5 982.5 * 732.5 202.5 2.5 17.5

2.5 2.5 2.5 * 602.5 2.5 2.5

2.5 32.5 2.5 942.5 * 2.5 2.5

* 17.5 37.5 42.5 37.5 37.5 22.5

(5–25) (20–45) (30–55) (20–45) (20–45) (5–30)

17.5 * 27.5 57.5 12.5 22.5 2.5

(0–265) (0–265) (0–310) (0–295) (0–300) (25–65) (5–25) (15–35) (45–70) (0–20) (10–30) (0–15)

37.5 42.5 * 982.5 37.5 12.5 37.5

(0–385) (750–990) (700–760) (110–330) (0–330) (0–340) (20–45) (30–55) (940–995) (20–45) (0–20) (20–45)

dynamic climate and volcanic changes during the PlioPleistocene explain the intraspecific diversification in P. leprosa and population expansion into the Central © 2013 John Wiley & Sons Ltd

2.5 12.5 12.5 * 27.5 17.5 2.5

(0–40) (0–30) (0–30) (430–755) (0–20) (0–25) (0–10) (0–20) (0–20) (15–35) (5–25) (0–10)

17.5 2.5 17.5 382.5 * 22.5 22.5

(0–225) (0–180) (0–210) (635–995) (0–125) (0–240) (5–25) (0–15) (5–25) (360–430) (10–30) (10–30)

Perote

2.5 2.5 2.5 2.5 2.5 * 2.5

(0–30) (0–30) (0–35) (0–45) (0–25)

2.5 2.5 2.5 22.5 2.5 * 2.5

(0–10) (0–10) (0–10) (10–30) (0–10)

(0–50)

(0–10)

Orizaba

7.5 2.5 2.5 2.5 72.5 292.5 * 2.5 2.5 7.5 27.5 2.5 22.5 *

(0–245) (0–205) (0–245) (0–255) (0–325) (260–380)

(0–10) (0–10) (0–15) (0–30) (0–10) (10–30)

TVB (Parra-Olea et al. 2012), we have found evidence that later climatic changes explain the more recent isolation of Central TVB populations.


ET AL. 3270 G . V E L O - A N T ON

TVB

22

TVB

20 19

12

N

11 10

3 1

15 9

4 2

5

21 25

0

16

50

km

8

7

6

14 13

18

Elevation (m)

17

0–1000 1000–1500 1500–2000 2000–2500 2500–3000 3000–3500 3500–4000 4000–4500 4500–5500

23

23. Tres Mogotes

21. G. Ortega 22. Tlaxco

20. Tlatlauquitepec

19. Vigas

18. Xometla

16. Malinche E

17. Texmalaquilla

15. MAlinche N

14. Malinche W

13. Malinche S

12. Nanacamilpa

11. Río Frío

10. Llano Grande

9. Atzompa

8. Calpan

7. Popocatepetl

3. D. Leones 4. Ajusco 5. Santa Marta 6. Zempoala

2. Texcalyacac

1. Nevado

GD

Fig. 2 Genetic structure in Pseudoeurycea leprosa. The seven genetic demes inferred in STRUCTURE are colour coded according to the key in the lower right hand corner. Individual assignment to the seven genetic demes is shown for each population below the map. Pie charts represent the average population assignment of P. leprosa individuals to genetic groups for each locality. Localities with less than four genotyped individuals (Gonz
alez Ortega, Tlaxco and Tres Mogotes) are not included in pie charts on the map but are included on the genetic assignment figure.

Table 4 Results of Mantel tests (r) testing the relationship between pairwise genetic structure (FST; DEST) and models of spatial resistance between Pseudoeurycea leprosa populations for the complete data set (All) and a subset of populations from the Central mtDNA clade [Central Trans-Mexican Volcanic Belt (TVB)]. Highest r values are represented in bold. All correlations were significant after Bonferroni correction (P < 0.01) DEST

FST Variables

Type of model

All

Central TVB

All

Central TVB

1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12.

Euclidean resistance Climatic resistance

0.40 0.28 0.31 0.34 0.39 0.39 0.37 0.38 0.39 0.39 0.25 0.40

0.40 0.41 0.40 0.40 0.47 0.40 0.41 0.45 0.45 0.46 0.41 0.46

0.50 0.40 0.41 0.41 0.55 0.53 0.46 0.49 0.60 0.51 0.34 0.41

0.43 0.39 0.42 0.39 0.53 0.44 0.42 0.48 0.58 0.47 0.35 0.40

IBR CCSA MIROC Present Landcover Topo complex Current + Topo Current + Land Land + Topo Current + Land + Topo CCSA 9 current MIROC 9 current

Landscape resistance Combined resistance

Stability resistance

Our analyses show a recent reduction in the range of P. leprosa and reduced spatial connectivity after the LGM. Our ecological niche models of historic (LGM) distributions show a large range and lower landscape resistance among populations, even at lower elevations, indicating that climatic conditions across the TVB were more suitable to the species during the LGM.

Subsequent climatic warming led to upward shifts in distribution and population structure. A continuous distribution and increased intermountain gene flow in P. leprosa during the LGM fits with the pattern observed in several montane species (Wiens 2004; Galbreath et al. 2009), but it is in contrast with other tropical montane specialists that persisted in isolated refugia over time © 2013 John Wiley & Sons Ltd

T R A C K I N G C L I M A T E C H A N G E I N A M O N T A N E S A L A M A N D E R 3271

TVB

LGM_CCSA

LGM_MIROC

Probability of occurence (%) 90 – 100 80 – 90 70 – 80 60 – 70 50 – 60 40 – 50 30 – 40 20 – 30 10 – 20 0 – 10

Present

Future (2050)

Fig. 3 Potential distribution of Pseudoeurycea leprosa identified using ecological niche modelling (ENMs) under historical [cooler, last glacial maximum (LGM) (21 000 years B.P.), CCSA and MIROC], present, and future (warmer, 2050) climatic times. Geo-referenced points from museum specimens and our own collections (73 unique localities at 1 km2 spatial resolution) were used to model P. leprosa distribution on an area that include all present populations of P. leprosa and surrounding mountainous areas where the species could potentially persist. Green dots on the maps represent the sampled localities used in this study to assess the genetic structure of the species.

(Garc
ıa-Par
ıs et al. 2000; Bowie et al. 2006; Cadena et al. 2007; Bell et al. 2010; Bryson & Riddle 2012). The high levels of connectivity observed among P. leprosa populations during LGM is likely due to the highly elevated plateau along the central and eastern TVB (c. 1500– 2000 m; Fig. 1). Cooler climates during the late Pleistocene glacial advances shifted pine-oak forests downward by 730–930 m (V
azquez-Selem & Heine 2004) forming dispersal corridors and increasing gene flow among mountain ranges on the TVB plateau. Our results show different genetic patterns in the younger Central TVB populations and the older source populations in the north-eastern mountain ranges, further supporting the recency of climate-induced population isolation in our focal species. The north-eastern and south-eastern populations continue to be isolated (as they were in the past, grouping in distinct mtDNA lineages). In contrast, based on microsatellite data, the most recently diverged Central TVB mtDNA clade is subdivided into five geographically structured genetic demes (Fig. 2). We found signals of admixture between some isolated mountain chains (Malinche and Xometla) that are most likely explained by ancestral polymorphism maintained in those populations, although colonization from neighbouring areas with ‘central’ genes that either have not been detected during fieldwork or have become extinct cannot be ruled out. A third line of © 2013 John Wiley & Sons Ltd

evidence for post-LGM genetic isolation across TVB populations is based on our estimates of migration. The genetic markers we used have very different, yet complementary, mutation rates and thus infer older (mtDNA) and more contemporary (microsatellites) migration events in P. leprosa. Thus, the population expansion and high gene flow from east to west inferred from the mtDNA data set (Parra-Olea et al. 2012) reflect historical migration across TVB populations, likely when the range of P. leprosa was continuous during colder climates. The lower migration estimates found with rapidly evolving genetic markers (microsatellites) indicate a reduction in gene flow among populations as habitat was progressively reduced and fragmented during warmer climates. Recent migration estimates from BAYESASS (within the last few generations) also show lower gene flow than migration estimates from MIGRATE, further corroborating the pattern of reduced gene flow among extant populations. Although all our analyses support continuous loss of connectivity among TVB mountain ranges since the LGM, we note that our analyses potentially violate assumptions for parameter estimation. MIGRATE assumes that no unsampled population are exchanging genes with sampled populations and that population size and migration rate do not change over time, both conditions that are unlikely to be met in our study.


ET AL. 3272 G . V E L O - A N T ON A CCSA (LGM)

B MIROC (LGM)

C Present

D 2050

E Topography

F Landcover

4 2 0 –2 –4

Pairwise resistance

6

Fig. 4 Connectivity maps among populations (black points) for Pseudeurycea leprosa. Colder (blue) colours indicate areas with higher current density; areas where connectivity is most tenuous are shown in warmer colours. Note that the resolution in present and future connectivity maps is higher than for the last glacial maximum (LGM) maps.

LGM_CCSA

LGM_MIROC

Present

Future (2050)

Fig. 5 Past and present spatial connectivity rates in Pseudoeurycea leprosa from three time periods [last glacial maximum (LGM); Present; Future]. Spatial resistance among populations was estimated using CIRCUITSCAPE and based on habitat suitability from ecological niche modelling (ENM) reconstructions.

We compared our results with other microsatellitebased landscape genetic studies of recently diverged salamander populations. A recent study of the fire salamander, Salamandra salamandra (Velo-Ant
on et al. 2012)

compared seven populations separated by 30 km, which have been isolated c. 9000 years. Our FST results in P. leprosa are similar to those found in S. salamandra (FST range: 0.06–0.27 and mean pairwise FST = 0.15), despite large differences in population size (higher in S. salamandra) and geographic distances (larger in P. leprosa) between the two species. Other fine-scale studies on salamanders using microsatellite markers showed lower FST values in continuous forest populations (Plethodon cinereus, FST = 0.019, Cabe et al. 2007) and similar (Phaeognathus hubrichti, FST = 0.15, Apodaca et al. 2012) or slightly higher FST values (Ambystoma tigrinum, FST = 0.24, Spear et al. 2005; Ambystoma macrodactylum, FST = 0.27, Savage et al. 2010) in highly fragmented landscapes. Finally, a recent study evaluated the effect of climate-induced habitat shifts on population genetic structure in the Large-blotched Ensatina (Ensatina eschscholtzii klauberi), a plethodontid salamander endemic to mid and high-elevation conifer forest in the Transverse and Peninsular Ranges of southern California and northern Baja California. (Devitt et al. 2013). As in our study, the authors combined ENMs from the © 2013 John Wiley & Sons Ltd

T R A C K I N G C L I M A T E C H A N G E I N A M O N T A N E S A L A M A N D E R 3273 present and two late Quaternary time periods with mitochondrial DNA and microsatellite data and concluded that recent climatic change has played a role in shaping patterns of population persistence and connectivity. Thus, recent population differentiation, coinciding with reduction in habitat spatial connectivity after the LGM, may be a more general pattern for low-vagility montane-adapted species. One difficulty of reconstructing general patterns in historical demography is the potential mismatch in temporal scale, when genetic markers of choice do not reflect processes that are temporally concordant with landscape or climatic changes (Zellmer & Knowles 2009; Anderson et al. 2010). Likewise, a spatial scale mismatch between genetic marker resolution and the geographic scale of sampling (in our case across distances much larger than the dispersal capacity of individual salamanders) can also reduce the signal of important abiotic factors in landscape genetic analyses. These issues of temporal and spatial scale will be exacerbated in cases of low-vagility species in very complex landscapes and are an ongoing subject of study in landscape genetics (Zellmer & Knowles 2009; Anderson et al. 2010; Cushman & Landguth 2010; Landguth et al. 2010). It is in these cases that the integration of past and present environmental data will be most useful in informing the potential causal link between genetic and environmental signals. In addition to climate warming, anthropogenic landscape alteration is a major threat to Mexican highland forests and the fauna that depends on them. Deforestation in the TVB is prevalent (Masera et al. 1997) and despite established protected areas, large-scale changes in vegetation and extinction are expected throughout the region (Villers-Ruiz & Trejo-V
azquez 1998). Pseudoeurycea leprosa is currently associated with pine, oak and fir forests that occur above 2000 m, and thus, the potential impacts of climate change cannot be easily disassociated from the direct negative anthropogenic effects on forest ecosystems (Villers-Ruiz & TrejoV
azquez 1998; Parra-Olea et al. 2005). We used our genetic data to examine the relative roles of these different factors by comparing 12 landscape models that combined climate, topography and land cover in estimates of conductance; Table 4). Isolation by geographic resistance (IBR) plays a large role in shaping population divergences in P. leprosa (Table 4) and our models incorporating landscape resistance based on climatic reconstructions, land cover and topographic complexity did not further explain population genetic variance beyond that already explained by IBR. However, landscape genetic analyses are often confounded by ‘noisy’ data (Storfer et al. 2007). Possible reasons for this nonsignificant result include the fact that multiple © 2013 John Wiley & Sons Ltd

demographic events, such as bottlenecks and genetic drift, have contributed to genetic structure over the history of these populations, thus increasing the historical ‘noise’ in the data and obscuring the signal of more current habitat changes. Alternatively, anthropogenic habitat change in this region is relatively recent and may not yet have left a population genetic signature beyond that which was already present due to historical evolutionary processes (Sumner et al. 2004 Dixo et al. 2009). However, reduction in genetic diversity, or changes in haplotypic diversity due to genetic drift can occur rapidly if migration is curtailed and population sizes are small (Lacy 1987; Peakall & Lindenmayer 2006), and this is most likely to happen in species with low dispersal capacity and low habitat availability (Luoy et al. 2007). Finally, The lack of significance in partial Mantel tests might be due to the fact that resistances calculated based on climatic and landscape features are highly correlated with Euclidean distances (IBD), a common finding in partial Mantel test analyses controlling for Euclidean distances (e.g. Klug et al. 2011). We also note that partial Mantel test can have inflated type I error with sample sizes less than ~ 50 (Raufaste & Rousset 2001; Legendre & Fortin 2010).

Constraints in population connectivity under scenarios of future climate change The ability of species to keep pace with modern climate change will depend on both the speed of the climatic changes and species-specific dispersal abilities (Parmesan & Yohe 2003; Pearson 2006; Moritz et al. 2008; Lawler et al. 2009; Loarie et al. 2009; Tingley et al. 2009; Schloss et al. 2012). Species with inherently low dispersal rates (such as plethodontid salamanders) are less likely to cope with temperature increases, will accumulate ‘climatic debts’ (Devictor et al. 2012; Schloss et al. 2012) and may only be able to track climatic shifts in habitats when environmental changes are gradual. In addition, recent anthropogenic landscape barriers (e.g. habitat fragmentation) also impede movement by species with low vagility, and render habitat-tracking and long-term population persistence unlikely (Travis 2003; Thomas et al. 2004). Understanding how species responded to past climate change can help us predict whether or not a species will cope with future global warming. We infer that large niche shifts in P. leprosa did not happen after the LGM because ecological niche characteristics are typically conserved over short to moderate time-spans (Peterson 2011) and niche conservatism is common among montane salamanders (Kozak & Wiens 2010). Our migration analyses indicate that P. leprosa tracked the historic displacement of habitats since the LGM; thus, climate warming has been a major


ET AL. 3274 G . V E L O - A N T ON driver in shaping the historical distribution of this species. Therefore, we expect that future increased temperatures will cause upwards shifts and/or local extinction of isolated populations, further reducing connectivity and gene flow among populations. Indeed, our projected distribution for P. leprosa in 2050 (Fig. 3) shows dramatic habitat loss, low geographic connectivity and very low probability of occurrence in most parts of its distribution. These projections do not bode well for P. leprosa; however, it is possible the models overestimate the impact of increased temperatures by not taking into account environmental variables at finer spatial scales (i.e. microclimate) (Pearson & Dawson 2003). Unfortunately, there are no historical fine-scale environmental data that would allow us to relate the distribution of microhabitats with population persistence or genetic structure. The role of microclimate and microhabitats in shaping spatial genetic structure will be of increasing importance in predicting species responses to future climate change. Whether montane populations persist will also depend on the nature, magnitude and rate of environmental shifts. Pseudoeurycea leprosa has several specific attributes that increase their vulnerability to new ecological conditions: first, tropical species typically have more restricted physiological tolerances and narrower climatic niches than temperate species, and thus, global warming is expected to exacerbate the decline of ectothermic animals in the tropics (Deutsch et al. 2008; Tewksbury et al. 2008; Rovito et al. 2009; Duarte et al. 2011). In addition, P. leprosa populations occur only in high elevation tropical montane regions and thus have evolved narrow thermal tolerances (Janzen 1967), and thus, this species will not easily acclimate or adapt to large temperature fluctuations (Ghalambor et al. 2006; Sheldon et al. 2011). Second, plethodontid salamanders are distinguished by the lack of lungs, conducting respiration through the skin and thus are prone to negative effects of rising temperatures and atmospheric CO2 (Bernardo & Spotila 2006). Third, P. leprosa populations are isolated (and will be even more so in the future) and thus are prone to genetic drift, loss of genetic diversity, inbreeding depression and the loss of evolutionary potential to adapt to new environmental conditions (Lande 1988; Frankham et al. 2002). Thus, maintenance of gene flow among natural populations is essential to the maintenance of genetic variation necessary for evolutionary tracking of environmental change. Fourth, the low dispersal ability of P. leprosa in combination with its highly fragmented distribution will also directly impede this species from tracking habitats. In addition to species-specific limitations, other factors we do not consider here could also influence the persistence of species over time, such as changes in biotic

interactions (e.g. diet, competitors, predators and species interactions). Parra-Olea et al. (2005) used GARP niche modelling to compare the predicted distribution of P. leprosa and P. cephalica in 2050. Currently, both species occur at the same elevations, but the overall distribution of P. cephalica is much broader than that of P. leprosa. That study showed a more dramatic predicted range reduction in P. leprosa with c. 75% of its current distributional area lost by the year 2050. Our study corroborates those results, and in addition, predicts the continuous reduction in gene flow among current TVB populations and the future loss of spatial connectivity between potential source regions, underscoring the high risks of local extinctions in P. leprosa. Management efforts for this and other TVB taxa need to preserve current populations and continuous habitat that might be optimal in the future. Dispersal corridors between fragmented populations in species with limited dispersal capabilities are crucial for genetic interchange and to preclude the loss of adaptive genetic diversity (Massot et al. 2008). However, the highly fragmented distribution of habitat predicted for the future will preclude effective corridors throughout most of the species’ range, and thus, active management, such as habitat restoration and translocations, should be considered and studied in detail as potential conservation actions to guarantee persistence of this and other highland-adapted species in the TVB.

Conclusions Our study integrated paleodistribution modelling, spatial connectivity analyses and multilocus genetic data (mtDNA and microsatellites), to evaluate factors affecting the distribution of genetic diversity in a montane landscape, in the past and in the future. Our integrative efforts revealed some shortcomings of landscape genetic approaches. First, species-specific characteristics are critical for accurate predictions, yet the particular mechanisms (e.g. physiological tolerance, habitat specialization) limiting species’ distributions are rarely known. Researchers will certainly benefit from ongoing improvements in ENMs, which incorporate species adaptive potential by including physiological parameters in addition to environmental data (Kearney & Porter 2009). Likewise, efforts to use finer-scale landscape and microclimatic variables to predict distributions (Pearson & Dawson 2003) will, in turn, better explain the effect of spatial features on population genetic structure, especially for cryptic and secretive organisms that are often specialized to particular microhabitats. Our study underscores that identifying biotic and abiotic features that enabled historical species persistence in © 2013 John Wiley & Sons Ltd

T R A C K I N G C L I M A T E C H A N G E I N A M O N T A N E S A L A M A N D E R 3275 space over time, improves our predictions of future climate change impacts (Williams et al. 2008; Bell et al. 2010), particularly for montane tropical species (Williams et al. 2007). Increasing our resolution and integrating mechanisms for species responses to environmental change are important goals for the field of landscape genetics in these times of rapid anthropogenic global change.

Acknowledgements Molecular data were collected in the Evolutionary Genetics Core Facility and the Core Laboratories Genotyping and Sequencing Facility at Cornell. Analyses benefited from resources provided by the Computational Biology Service Unit at Cornell University, a facility partially funded by Microsoft Corporation. GVA was supported by a postdoctoral fellowship from the Spanish Ministerio de Ciencia e Innovaci
on (Ref: 2008-0890) and a post-doctoral fellowship from Fundacß~ ao para a Ci^encia e Tecnologia (Ref: SFRH/BPD/74834/2010; FCT, Portugal). Fieldwork and laboratory efforts were funded by grants from PAPIIT-UNAM (number IN212111 to G. P.-O.) and the National Science Foundation (to KZ). We thank Dr. Stephen Spear and three anonymous referees for constructive comments on earlier versions of the manuscript.

References Anderson CD, Epperson BK, Fortin M-J et al. (2010) Considering spatial and temporal scale in landscape-genetic studies of gene flow. Molecular Ecology, 19, 3576–3591. Apodaca JJ, Rissler LJ, Godwin JC (2012) Population structure and gene flow in a heavily disturbed habitat: implications for the management of the imperiled Red Hills salamander (Phaeognathus hubrichti). Conservation Genetics, 13, 913–923. Avise JC (2000) Phylogeography: The History and Formation of Species. Harvard University Press, Cambridge, Massachusetts. Beerli P, Felsenstein J (1999) Maximum-likelihood estimation of migration rates and effective population numbers in two populations using a coalescent approach. Genetics, 152, 763–773. Bell RC, Parra JL, Tonione M, Hoskin CJ (2010) Patterns of persistence and isolation indicate resilience to climate change in montane rainforest lizards. Molecular Ecology, 19, 2531–2544. Bellard C, Bertelsmeier C, Leadley P, Thuiller W, Courchamp F (2012) Effects of climate change on the future of biodiversity. Ecology Letters, 15, 365–377. Bernardo J, Spotila J (2006) Physiological constraints on organismal response to global warming; mechanistic insights from clinally varying populations and implications for assessing endangerment. Biology Letters, 2, 135–139. Bonnet E, Van de Peer Y (2002) ZT: a software tool for simple and partial Mantel tests. Journal of Statistical Software, 7, 1–12. Bowie RCK, Fjelds a J, Hackett SJ, Bates JM, Crowe TM (2006) Coalescent models reveal the relative roles of ancestral polymorphism, vicariance, and dispersal in shaping phylogeographical structure of an African montane forest robin. Molecular Phylogenetics and Evolution, 38, 171–188. Bryson RW Jr, Riddle BR (2012) Tracing the origins of widespread highland species: a case of Neogene diversification

© 2013 John Wiley & Sons Ltd

across the Mexican sierras in an endemic lizard. Biological Journal of the Linnean Society, 105, 382–394. Bryson RW Jr, Murphy RW, Graham MR, Lathrop A, LazcanoVillareal D (2011) Ephemeral Pleistocene woodlands connect the dots for highland rattlesnakes of the Crotalus intermedius group. Journal of Biogeography, 38, 2299–2310. Cabe PR, Page RB, Hanlon TJ, Aldrich ME, Connors L, Marsh DM (2007) Fine-scale population differentiation and gene flow in a terrestrial salamander (Plethodon cinereus) living in a continuous habitat. Heredity, 98, 53–60. Cadena CD, Klicka J, Ricklefs RE (2007) Evolutionary differentiation in the Neotropical montane region: molecular phylogenetics and phylogeography of Buarremon brush-finches (Aves, Emberizidae). Molecular Phylogenetics and Evolution, 44, 993–1016. Cadena CD, Kozak KH, G
omez JP et al. (2012) Latitude, elevational climatic zonation, and speciation in New World vertebrates. Proceedings of the Royal Society B: Biological Sciences, 279, 194–201. Chan LM, Brown JL, Yoder AD (2011) Integrating statistical genetic and geospatial methods brings new power to phylogeography. Molecular Phylogenetics and Evolution, 59, 523–537. Colwell RK, Brehm G, Cardel
us CL, Gilman A, Longino JT (2008) Global warming, elevational range shifts, and lowland biotic attrition in the wet tropics. Science, 322, 258–261. Crawford NG (2010) SMOGD: software for the measurement of genetic diversity. Molecular Ecology Resources, 10, 556–557. Cushman SA, Landguth EL (2010) Scale dependent inference in landscape genetics. Landscape Ecology, 25, 967–979. Deutsch CA, Tewksbury JJ, Huey RB et al. (2008) Impacts of climate warming on terrestrial ectotherms across latitude. Proceedings of the National Academy of Sciences, 105, 6668–6672. Devictor V, van Swaay C, Brereton T et al. (2012) Differences in the climatic debts of birds and butterflies at a continental scale. Nature Climate Change, 2, 121–124. Devitt T, Cameron Devitt S, Hollingsworth B, McGuire J, Moritz C (2013) Montane refugia predict population genetic structure in the Large-blotched Ensatina Salamander. Molecular Ecology, 22, 1650–1665. Dixo M, Metzger JP, Morgante JS, Zamudio KR (2009) Habitat fragmentation reduces genetic diversity and connectivity among toad populations in the Brazilian Atlantic Coastal Forest. Biological Conservation, 142, 1560–1569. Duarte H, Tejedo M, Katzenberger M et al. (2011) Can amphibians take the heat? Vulnerability to climate warming in subtropical and temperate larval amphibian communities. Global Change Biology, 18, 412–421. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology, 14, 2611–2620. Excoffier L, Laval G, Schneider S (2005) Arlequin ver. 3.0: an integrated software package for population genetics data analysis. Evolutionary Bioinformatics Online, 1, 47–50. Falush D, Stephens M, Pritchard JK (2003) Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics, 164, 1567–1587. Feeley KJ, Silman MR (2010) Land-use and climate change effects on population size and extinction risk of Andean plants. Global Change Biology, 16, 3215–3222. Ferrari L, Conticelli S, Vaggelli G, Petrone CM, Manetti P (2000) Late Miocene volcanism and intra-arctectonics during


ET AL. 3276 G . V E L O - A N T ON the early development of the Trans-Mexican Volcanic Belt. Tectonophysics, 318, 161–185. Ferraz G, Nichols JD, Hines JE, Stouffer PC, Bierregaard RO Jr, Lovejoy TE (2007) A large-scale deforestation experiment: effects of patch area and isolation on Amazon birds. Science, 315, 238–241. Fielding AH, Bell JF (1997) A review of methods for the measurement of prediction errors in conservation presence/ absence models. Environmental Conservation, 24, 38–49. Fischer J, Lindenmayer DB (2007) Landscape modification and habitat fragmentation: a synthesis. Global Ecology and Biogeography, 16, 265–280. Frankham R, Briscoe DA, Ballou JD (2002) Introduction to Conservation Genetics. Cambridge University Press, New York. Galbreath KE, Hafner DJ, Zamudio K (2009) When cold is better: climate-driven elevation shifts yield complex patterns of diversification and demography in an alpine specialist (American pika, Ochotonaprinceps). Evolution, 63, 2848–2863. Garc
ıa-Par
ıs M, Good DA, Parra-Olea G, Wake DB (2000) Biodiversity of Costa Rican salamanders: implications of high levels of genetic differentiation and phylogeographic structure for species formation. Proceedings of the National Academy of Sciences, 97, 1640–1647. Ghalambor CK, Huey RB, Martin PR, Tewksbury JJ, Wang G (2006) Are mountain passes higher in the tropics? Janzen′s hypothesis revisited. Integrative and Comparative Biology, 46, 5–17. Goldberg CS, Waits LP (2010) Comparative landscape genetics of two pond-breeding amphibian species in a highly modified agricultural landscape. Molecular Ecology, 19, 3650– 3663. Hewitt GM (2000) The genetic legacy of the Quaternary ice ages. Nature, 405, 907–913. Huey R, Deutsch CA, Tewksbury JJ et al. (2009) Why tropical forest lizards are vulnerable to climate warming. Proceedings of the Royal Society B: Biological Sciences, 276, 1939–1948. Hugall A, Moritz C, Moussalli A, Stanisic A (2002) Reconciling paleodistribution models and comparative phylogeography in the Wet Tropics rainforest land snail Gnarosophia bellendenkerensis (Brazier 1875). Proceedings of the National Academy of Sciences, 99, 6112–6117. Jakobsson M, Rosenberg NA (2007) CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics, 23, 1801–1806. Janzen DH (1967) Why mountain passes are higher in the tropics. American Naturalist, 101, 233–249. Jost L (2008) GST and its relatives do not measure differentiation. Molecular Ecology, 17, 4015–4026. Kalinowski ST (2005) HP-Rare: a computer program for performing rarefaction on measures of allelic diversity. Molecular Ecology Notes, 5, 187–189. Kearney M, Porter W (2009) Mechanistic niche modelling: combining physiological and spatial data to predict species’ ranges. Ecology Letters, 12, 334–350. Kiesecker JM, Blaustein AR, Belden LK (2001) Complex causes of amphibian population declines. Nature, 410, 681–684. Klug PE, Wisely SM, With KA (2011) Population genetic structure and landscape connectivity of the Eastern Yellowbelly Racer (Coluber constrictor flaviventris) in the tallgrass prairie of northeastern Kansas. Landscape Ecology, 26, 281–294.

Kozak KH, Wiens JJ (2010) Niche conservatism drives elevational diversity patterns in Appalachian salamanders. American Naturalist, 176, 40–54. Lacy RC (1987) Loss of genetic diversity from managed populations: interacting effects of drift, mutation, immigration, selection and population subdivision. Conservation Biology, 1, 143–158. Lande R (1988) Genetics and demography in biological conservation. Science, 241, 1455–1460. Landguth E, Cushman S, Schwartz MK, McKelvey K, Murphy MA, Luikart G (2010) Quantifying the lag time to detect barriers in landscape genetics. Molecular Ecology, 19, 4179– 4191. Lawler JJ, Shafer SL, White D et al. (2009) Projected climate induced faunal change in the Western Hemisphere. Ecology, 90, 588–597. Lee-Yaw JA, Davidson A, McRae BH, Green DM (2009) Do landscape processes predict phylogeographic patterns in the wood frog? Molecular Ecology, 18, 1863–1874. Legendre P, Fortin MJ (2010) Comparison of the Mantel test and alternative approaches for detecting complex multivariate relationships in the spatial analysis of genetic data. Molecular Ecology Resources, 10, 831–844. Loarie SR, Duffy PB, Hamilton H, Asner GP, Field CB, Ackerly DD (2009) The velocity of climate change. Nature, 462, 1052– 1055. Luoy D, Habel JC, Schmitt T, Assmann T, Meyer M, Mu ller P (2007) Strongly diverging population genetic patterns of three skipper species: the role of habitat fragmentation and dispersal ability. Conservation Genetics, 8, 671–681. Manel S, Williams HC, Ormerod SJ (2001) Evaluating presenceabsence models in ecology: the need to account for prevalence. Journal of Applied Ecology, 38, 921–931. Mantel N (1967) The detection of disease clustering and a generalized regression approach. Cancer Research, 27, 209–220. Mart
ınez-Solano I, Jockusch EL, Wake DB (2007) Extreme population subdivision throughout a continuous range: phylogeography of Batrachoseps attenuatus (Caudata: Plethodontidae) in western North America. Molecular Ecology, 16, 4335–4355. Masera OR, Ordo~ nez MJ, Dirzo R (1997) Carbon emissions from Mexican forests: current situation and long-term scenarios. Climatic Change, 35, 265–295. Massot M, Clobert J, Ferriere R (2008) Climate warming, dispersal inhibition and extinction risk. Global Change Biology, 14, 461–469. McRae BH (2006) Isolation by resistance. Evolution, 60, 1551– 1561. McRae BH, Beier P (2007) Circuit theory predicts gene flow in plant and animal populations. Proceedings of the National Academy of Sciences, 104, 19885–19890. McRae BH, Dickson BG, Keitt TH, Shah VB (2008) Using circuit theory to model connectivity in ecology and conservation. Ecology, 10, 2712–2724. Moritz C, Patton JL, Conroy CJ, Parra JL, White GC, Beissinger SR (2008) Impact of a century of climate change on smallmammal communities in Yosemite National Park, USA. Science, 322, 261–264. Murphy M, Dezzani R, Pilliod D, Storfer A (2010) Landscape genetics of high mountain frog metapopulations. Molecular Ecology, 19, 3634–3649.

© 2013 John Wiley & Sons Ltd

T R A C K I N G C L I M A T E C H A N G E I N A M O N T A N E S A L A M A N D E R 3277 Myers N, Mittermeier RA, Mittermeier CG, da Fonseca GA, Kent J (2000) Biodiversity hotspots for conservation priorities. Nature, 403, 853–858. Parmesan C (2006) Ecological and evolutionary responses to recent climate change. Annual Review of Ecology, Evolution, and Systematics, 37, 637–669. Parmesan C, Yohe G (2003) A globally coherent fingerprint of climate change impacts across natural systems. Nature, 421, 37–42. Parra-Olea G, Mart
ınez-Meyer E, Ponce P
erez, de Le
on G (2005) Forecasting climate change effects on salamander distribution in the highlands of Central Mexico. Biotropica, 37, 202–208. Parra-Olea G, Windfield JC, Velo-Ant
on G, Zamudio KR (2012) Isolation and habitat refugia promote rapid diversification in a montane tropical salamander. Journal of Biogeography, 39, 353–370. Peakall R, Lindenmayer D (2006) Genetic insights into population recovery following experimental perturbation in a fragmented landscape. Biological Conservation, 132, 520–532. Peakall R, Smouse PE (2006) GenAlEx 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes, 6, 288–295. Pearson RG (2006) Climate change and the migration capacity of species. Trends in Ecology and Evolution, 21, 111–113. Pearson RG, Dawson TP (2003) Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful? Global Ecology and Biogeography, 12, 361–371. Peterson AT (2011) Ecological niche conservatism: a time-structured review of evidence. Journal of Biogeography, 38, 817–827. Phillips SJ, Anderson RP, Schapire RE (2006) Maximum entropy modeling of species geographic distributions. Ecological Modelling, 190, 231–259. Pimm SL (2008) Biodiversity: climate change or habitat loss which will kill more species? Current Biology, 18, 117–119. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics, 155, 945–959. Raufaste N, Rousset F (2001) Are partial Mantel tests adequate? Evolution, 55, 1703–1705. Reed DH, Frankham R (2003) Correlation between fitness and genetic diversity. Conservation Biology, 17, 230–237. Rice WR (1989) Analyzing tables of statistical tests. Evolution, 43, 223–225. Root TL, Price JT, Hall KR, Schneider SH, Rosenzweig C, Pounds JA (2003) Fingerprints of global warming on wild animals and plants. Nature, 421, 57–60. Rovito SM, Parra-Olea G, Vasquez-Almazan CR, Papenfuss TJ, Wake DB (2009) Dramatic declines in neotropical salamander populations are an important part of the global amphibian crisis. Proceedings of the National Academy of Sciences, 106, 3231–3236. Row JR, Blouin-Demers G, Lougheed SC (2010) Habitat distribution influences dispersal and fine-scale genetic population structure of eastern foxsnakes (Mintonius gloydi) across a fragmented landscape. Molecular Ecology, 19, 5157–5171. Safner T, Miaud C, Gaggiotti O et al. (2011) Combining demography and genetic analysis to assess the population structure of an amphibian in a human-dominated landscape. Conservation Genetics, 12, 161–173.

© 2013 John Wiley & Sons Ltd

Savage WK, Fremier AK, Shaffer HB (2010) Landscape genetics of alpine Sierra Nevada salamanders reveal extreme population subdivision in space and time. Molecular Ecology, 19, 3301–3314. Schloss CA, Nu~ nez TA, Lawler JJ (2012) Dispersal will limit ability of mammals to track climate change in the Western Hemisphere. Proceedings of the National Academy of Sciences, 109, 8606–8611. Sheldon KS, Yang S, Tewksbury JJ (2011) Climate change and community disassembly: impacts of warming on tropical and temperate montane community structure. Ecology Letters, 14, 1191–1200. Sherry RA, Zhou X, Gu S et al. (2007) Divergence of reproductive phenology under climate warming. Proceedings of the National Academy of Sciences, 104, 198–202. Slatkin M (1987) Gene flow and the geographic structure of natural populations. Science, 236, 787–792. Smouse PE, Long JC, Sokal RR (1986) Regression and correlation extensions of the Mantel Test of matrix correspondence. Systematic Zoology, 35, 627–632. Spear SF, Peterson CR, Matocq M, Storfer A (2005) Landscape genetics of the blotched tiger salamander, Ambystoma tigrinum melanostictum. Molecular Ecology, 14, 2553–2564. Storfer A, Murphy MA, Evans JS et al. (2007) Putting the ‘‘landscape’’ in landscape genetics. Heredity, 98, 128–142. Sumner J, Jessop T, Paetkau D, Moritz C (2004) Limited effect of anthropogenic habitat fragmentation on molecular diversity in a rain forest skink, Gnypetoscincus queenslandiae. Molecular Ecology, 13, 254–269. Taberlet P, Fumagalli L, Wust-Saucy A-G, Cosson J-F (1998) Comparative phylogeography and postglacial colonization routes in Europe. Molecular Ecology, 7, 453–464. Tewksbury JJ, Huey RB, Deutsch CA (2008) Ecology. Putting the heat on tropical animals. Science, 320, 1296–1297. Thomas CD, Cameron A, Green RE et al. (2004) Extinction risk from climate change. Nature, 427, 145–148. Tingley MW, Monahan WB, Beissinger SR, Moritz C (2009) Birds track their Grinnellian niche through a century of climate change. Proceedings of the National Academy of Sciences, 106, 19637–19643. Travis JMJ (2003) Climate change and habitat destruction: a deadly anthropogenic cocktail. Proceedings of the Royal Society B: Biological Sciences, 270, 467–473. Van Oosterhaut C, Hutchinson WF, Wills DPM, Shipley P (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes, 4, 535–538. V
azquez-Selem L, Heine K (2004) Late Quaternary glaciation of Mexico. In: Quaternary Glaciations: Extent and Chronology, part III (eds Ehlers J & Gibbars PL), pp. 233–242. Elsevier, Amsterdam. Velo-Ant
on G, Windfield JC, Zamudio K, Parra-Olea G (2009) Microsatellite markers for Pseudoeurycea leprosa, a plethodontid salamander endemic to the Transmexican Neovolcanic Belt. Conservation Genetics Resources, 1, 5–7. Velo-Ant
on G, Zamudio KR, Cordero-Rivera A (2012) Genetic drift and rapid evolution of viviparity in insular fire salamanders (Salamandra salamandra). Heredity, 108, 410–418. Villers-Ruiz L, Trejo-V
azquez I (1998) Climate change on Mexican forests and natural protected areas. Global Environmental Change, 8, 141–157.


ET AL. 3278 G . V E L O - A N T ON Visser ME (2008) Keeping up with a warming world; Assessing the rate of adaptation to climate change. Proceedings of the Royal Society B: Biological Sciences, 275, 649–659. Wang YH, Yang KC, Bridgman CL, Lin LK (2008) Habitat suitability modelling to correlate gene flow with landscape connectivity. Landscape Ecology, 23, 989–1000. Wiens JJ (2004) Speciation and ecology revisited: phylogenetic niche conservatism and the origin of species. Evolution, 58, 193–197. Williams SE, Bolitho EE, Fox S (2003) Climate change in Australian tropical rainforests: an impending environmental catastrophe. Proceedings of the Royal Society B: Biological Sciences, 270, 1887–1892. Williams JW, Jackson ST, Kutzbacht JE (2007) Projected distributions of novel and disappearing climates by 2100 AD. Proceedings of the National Academy of Sciences, 104, 5738–5742. Williams SE, Shoo LP, Isaac JL, Hoffmann AA, Langham G (2008) Towards an integrated framework for assessing the vulnerability of species to climate change. Plos Biology, 6, 2621–2626. Wilson GA, Rannala B (2003) Bayesian inference of recent migration rates using multilocus genotypes. Genetics, 163, 1177–1191. With KA, King AW (1999) Extinction thresholds for species in fractal landscapes. Conservation Biology, 13, 314–326. Zellmer AJ, Knowles L (2009) Disentangling the effects of historic vs. contemporary landscape structure on population genetic divergence. Molecular Ecology, 18, 3593–3602.

question with contributions from coauthors, collected genetic data, analysed and interpreted results and led the writing of the paper with important contributions from all coauthors; G.P.-O. funded the fieldwork and laboratory data collection; K.Z. provided important feedback that improved analyses and the manuscript; J.P. provided significant intellectual contributions to spatial analyses.

Data accessibility Microsatellite genotypes and input files (geographic coordinates for localities and conductance grids) used to construct species distribution models and resistance matrices have been deposited in DRYAD (doi:10.5061/ dryad.tm093). DNA sequences were obtained from Parra-Olea et al. 2012 and were deposited in GenBank (accessions GQ468424–GQ468493).

Supporting information Additional supporting information may be found in the online version of this article. Table S1 Sample sizes and genetic diversity values for each population.

This collaboration stems from our combined interest in processes generating diversity in tropical amphibians and their conservation. G.V.-A. developed the research

Table S2 Bioclimatic variables ranked according to their overall model contribution for Pseudoeurycea leprosa populations.

© 2013 John Wiley & Sons Ltd