A cc ep te d A rt ic le 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.13275 This article is protected by copyright. All rights reserved. Received Date : 12-Mar-2015 Revised Date : 09-Jun-2015 Accepted Date : 12-Jun-2015 Article type : Original Article Testing the role of ecology and life history in structuring genetic variation across a landscape: A trait-based phylogeographic approach. Andrea Paz1, *, Roberto Ibáñez2, 3, Karen R. Lips2, 4, Andrew J. Crawford1, 2, 3 1 Departamento de Ciencias Biológicas, Universidad de los Andes, A.A. 4976, Bogotá, Colombia 2 Smithsonian Tropical Research Institute, Apartado 0843–03092, Panamá, Republic of Panama 3 Círculo Herpetológico de Panamá, Apartado 0824-00122, Panamá, Republic of Panama 4 Department of Biology, University of Maryland, College Park MD 20742-4415 * Corresponding author Email addresses: paz.andreita@gmail.com (AP), ibanezr@si.edu (RI), klips@umd.edu (KRL), andrew@dna.ac (AJC) Keywords: Dispersal ability, DNA barcoding, Ecological niche modeling, Hierarchical approximate Bayesian computation, Landscape resistance, phylogeography. Running title: Trait-based phylogeography of Panamanian frogs Type of article: Original Article Online Material: Table S1, Table S2, Table S3, Table S4 and Figure S1, Figure S2 Figure S3, Figure S4 and Figure S5 Correspondence: Andrea Paz, Facultad de Administración Universidad de los Andes Bogotá, Colombia, South America E-mail:paz.andreita@gmail.com Telephone: +57 1 3394949 x4893 A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Abstract Hypotheses to explain phylogeographic structure traditionally invoke geographic features, but often fail to provide a general explanation for spatial patterns of genetic variation. Organism’s intrinsic characteristics might play more important roles than landscape features in determining phylogeographic structure. We developed a novel comparative approach to explore the role of ecological and life-history variables in determining spatial genetic variation and tested it on frog communities in Panama. We quantified spatial genetic variation within 31 anuran species based on mitochondrial DNA sequences, for which hierarchical approximate Bayesian computation analyses rejected simultaneous divergence over a common landscape. Regressing ecological variables on genetic divergence allowed us to test the importance of individual variables revealing that body size, current landscape resistance, geographic range, biogeographic origin, and reproductive mode were significant predictors of spatial genetic variation. Our results support the idea that phylogeographic structure represents the outcome of an interaction between organisms and environment, and suggest a conceptual integration we refer to as ecophylogeography. Keywords: Dispersal ability, DNA barcoding, Ecological niche modeling, Hierarchical approximate Bayesian computation, Landscape resistance, Phylogeography. Introduction A key first step in the generation of species diversity is limiting gene flow between populations, whether by natural selection or by simple physical isolation due to allopatry (Mayr A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. 1942). Understanding the factors that determine species’ distributions therefore becomes a central goal in evolutionary biology (Gaston 2000; Wiens 2004). Historical processes that may have influenced the history of species’ distributions may be inferred from shared geographical or temporal patterns of genetic subdivision across co-occurring species (Avise 1992; Bermingham & Martin 1998). When phylogeographic structure is shared among multiple, co-distributed species, the most parsimonious explanation is that they have responded to the same geographical or environmental events (Bermingham & Avise 1986; Zink 1996; Avise 2000; Arbogast & Kenagy 2001). Comparative phylogeographic studies often support the hypothesis that geographic features such as mountains and rivers serve as vicariant barriers driving allopatric divergence and creating possible geminate species (Avise et al. 1987; Marko 2002). By reducing the movement of individuals, barriers reduce gene flow between populations and increase population structure (Wright 1951). Various hypotheses regarding types of geographic barriers have been tested. The riverine barrier hypothesis for example states that rivers are important in separating populations (Wallace 1852). A number of studies have been conducted that support (Ayres & Clutton-Brock 1992; Bates et al. 2004; Fouquet et al. 2012; Ribas et al. 2012) or reject (Gascon et al. 2000; Gehring et al. 2012) this hypothesis for different groups of animals. The Pleistocene refugia hypothesis states that during the cooler and drier period of each Pleistocene climatic cycle suitable habitat for wet-forest species was broken into isolated blocks (Haffer 1969). While the refugia hypothesis has been supported in a few cases (Waltari et al. 2007; Byrne 2008), it has been widely rejected (Lessa et al. 2003), especially in cases where lineages are older than the A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Pleistocene (Rull 2008). Finally, while Pleistocene climate change has provided a general explanation for the phylogeographic patterning of many North American and European taxa (Hewitt 2000; Johnson & Cicero 2004), such geographic mechanisms of lineage diversification seem to fail as a general explanation across the tropical realm (Antonelli et al. 2010). Given the lack of generality of geographic explanations of phylogeographic structure, we hypothesize that variation among co-distributed species in phylogeographic histories could be predicted from variation in key ecological traits. Dispersal limitation is caused by the interaction between an organism and its environment, whose outcome is determined as much by the extrinsic environmental conditions as by the intrinsic biological properties of the organism (Bernardo et al. 2007). The physical height of a potential montane barrier, for example, may not limit species’ distributions so much as the variation in moisture or temperature along an environmental gradient (Janzen 1967). Thus we can think of the mountain range as a physiological rather than a physical barrier, suggesting that ecological rather than geographical variables may be more predictive of dispersal potential and thus phylogeographic structure (Crawford et al. 2007; Burney & Brumfield 2009; Fine et al. 2013; Wang et al. 2013). Amphibians are thought to be physiologically intolerant of novel environments due to their permeable skin (Navas & Otani 2007) as well as evolutionarily conservative in the breadth of their environmental niche (Wiens et al. 2006). Amphibians show substantial variation in their life history characteristics (Crump 1974; Wells 2007), however, and we might expect them to respond differentially to variation in environmental factors. Studies of Central American frogs A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. have shown that phylogeographic histories of different species show few common spatial or temporal patterns across a shared landscape, leaving classical vicariant hypotheses with little predictive power (Weigt et al. 2005; Crawford et al. 2007; Wang et al. 2008; Robertson et al. 2009). These observations suggest that amphibians may be an ideal group to evaluate whether intrinsic variables can explain the observed divergence where geographic variables have thus far failed. Here we examine three frog communities in Panama using multiple co-occurring taxa to evaluate whether variation among species in intraspecific genetic divergence across a common landscape may be predicted by ecological variables potentially related to a frog’s ability to disperse, survive and reproduce in a heterogeneous landscape. Each of the three communities represents a mid-elevation montane site separated by stretches of lower elevation habitat. We evaluate four hypotheses of increasing complexity to explain average pairwise genetic divergence between localities, our surrogate for phylogeographic structure. First, the simplest hypothesis assumes that physical landscape features alone are responsible for the separation of populations, thus a single divergence event (i.e., similar levels of genetic divergence) should apply to all species across the shared landscape. Second, we ask whether variation in the altitudinal distribution of species might predict dispersal ability across the lowlands, and therefore cause variation in genetic divergence and divergence times among co-distributed species. In other words, is it more difficult for highland species to disperse through lowland habitats? In this analysis, however, elevation may be just a proxy for environmental A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. factors (Janzen 1967), thus we also sought to quantify variation among species in the environmental suitability of a common landscape. Thus, our third hypothesis is that interspecific variation in genetic divergence is predicted by variation in abiotic niche requirements of species, as quantified by species distribution models (Peterson 2001; Soberón & Nakamura 2009). Rather than looking at single species across a diverse landscape, as commonly done in landscape genetics (Storfer et al. 2007), here we study a diversity of species across a single landscape in a comparative phylogeographic approach. Finally, we ask whether life history traits themselves, such as developmental mode and body size, limit dispersal ability and thereby predict variation among species in pairwise genetic divergence between populations (Wollenberg Valero 2015). This latter approach is rarely attempted, but could become more common as phylogeographic studies are eventually applied to whole communities, e.g., through DNA barcoding campaigns (e.g., Murphy et al. 2013), thus providing genetic data on a diversity of co-distributed species. METHODS Geographic and taxonomic sampling Sampling was performed in three Panamanian national parks: G. D. Omar Torrijos H. National Park located in the province of Coclé, Chagres National Park located in the provinces of Panamá and Colón, and Cana field station in the Darién National Park located in the Darién Province, referred to here as El Copé, Brewster and Cana, respectively (Fig. 1). We selected 31 species of frogs present in at least two of the three localities and with at least two individuals sampled per locality for a total of 481 individuals, 179 from El Copé, 190 from Brewster and 112 from Cana (Table S1). Sampling at each site varied in elevation: El Copé was sampled at 700 – A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. 750 m, Brewster was sampled at 135 – 810 m and Cana was sampled from 462 – 1320. Samples were collected prior to epizootic outbreaks and population declines at each site (Crawford et al. 2010). Molecular data collection Molecular data for this study consisted of DNA sequences from two mitochondrial gene fragments: cytochrome oxidase subunit 1 (COI) and the 16S ribosomal RNA gene. Molecular data from El Copé were published previously (Crawford et al. 2010), and sequence data for the other two localities were obtained with the same primers and protocols. Full length COI sequences were 658 base pairs (bp) while full length 16S fragments ranged from 491 to 684 bp with a species mean of 561 bp. Alignment of conspecific COI and 16S sequences was unambiguous and could be accomplished by eye. Gapped sites were excluded from the analyses. Specimen vouchers and GenBank numbers are provided in Table S2. Hierarchical approximate Bayesian computation analysis One particularly effective method to detect correspondence in the temporal dimension between co-occurring taxa is hierarchical approximate Bayesian computation (HABC). HABC includes data from multiple loci and species in a single analysis within a coalescent framework to estimate interspecific congruence in divergence times (Hickerson et al. 2006a, 2010; Hickerson & Meyer 2008). Using coalescent simulations, this method accounts for uncertainty in current and ancestral population sizes and the stochastic nature of the mutational and coalescent A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. processes in estimating genetic divergence among populations (Knowles & Maddison 2002; Hickerson et al. 2006b; Nielsen & Beaumont 2009). To test for simultaneous divergence we performed HABC analyses with the program MTML-msBayes (Huang et al. 2011), following recent recommendations for implementation and analyses using this software (Oaks et al. 2012; Hickerson et al. 2014). Thirty species were analyzed in the comparison between El Copé and Brewster, and fifteen species were compared between Cana and Brewster (Fig. 1). All pairwise comparisons involved conspecific samples except for four lineages treated as conspecific due to their low genetic divergence or lack of taxonomic resolution: the three species of Atelopus, one endemic to each site, Pristimantis museosus versus its close relative, P. aff. museosus, Craugastor opimus versus C. megacephalus, and the as yet unnamed Diasporus spp. Conspecific samples included a minimum of two samples per population, as required for estimation of genetic divergences (Hickerson et al. 2007). MTML-msBayes is run as a 3-step process. First, several summary statistics including the net average pairwise differences (πnet, Nei & Li 1979), average pairwise differences between populations (πB, Takahata & Nei 1985), average pairwise differences within populations (πS, Takahata & Nei 1985), and Tajima’s D (DT, Tajima 1989) are calculated from the data (Hickerson et al. 2006a). Second, the program performs coalescent simulations of a dataset of the same size as the observed dataset, drawing randomly parameter values from a prior distribution and calculating summary statistics for each of the simulated datasets. Finally, summary statistics of the simulated datasets are compared to the observed estimates, and those parameter values that produced statistics differing from the observed data by less than a value ε are accepted A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. (Hickerson et al. 2007; Huang et al. 2011). The values obtained from the acceptance/rejection step were subjected to a local linear regression to improve estimation. We used the average pairwise differences between each population pair (πB) as our summary statistic. This statistic is easily computed and biologically relevant even in species with low sample sizes per population (Hickerson et al. 2007). πB has been shown to perform well in the face of small sample sizes, and rejecting the null hypothesis of simultaneous divergence becomes more difficult using a single summary statistic. Thus, relying on only πB as our summary statistic makes our test more conservative (Hickerson et al. 2007). To evaluate relative posterior support for number of divergence intervals we used Bayes Factors (BF, Hickerson & Meyer 2008). We used the BF results to constrain the number of intervals to their most likely values and re-ran the HABC analysis to obtain the most likely number of species in each divergence interval and the species membership within each of these temporal intervals. Ecological niche modeling We sought a composite measure of species-specific habitat suitability, as data are not available for all life-history variables of potential interest (e.g., resting metabolic rate or moisture uptake efficiency). We used ecological niche models (ENM, Peterson 2001) as a proxy for the abiotic environmental requirements of each species, and used these niche models to quantify variation among species in habitat suitability across a common landscape, as measured by ‘landscape resistance’ (see below). Algorithms for generating models from presence-only data result in a continuous prediction of the species occurrence probability for a given a set of environmental conditions (Soberón & Peterson 2005; Elith et al. 2006; Elith & Leathwick 2009). A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Raw data for modeling consisted of a set of georeferenced localities obtained from the Global Biodiversity Information Facility (http://www.gbif.org/), HerpNET (http://www.herpnet.org/), Instituto de Ciencias Naturales de la Universidad Nacional de Colombia (http://www.biovirtual.unal.edu.co/ICN/), the literature (e.g., Duellman 2001) and our own fieldwork. All previously georeferenced localities were visually screened for obvious coordinate errors and localities without coordinates were georeferenced by the authors. Number of locality records for each species ranged from 15 for Sachatamia ilex to 283 for Engystomops pustulosus, with only four species having fewer than 25 locality records (Table S1). We generated one ENM for each of 30 species (Pristimantis museosus was not used, as we could obtain only four locality records) using the maximum entropy algorithm implemented in Maxent 3.3.3 using presence and randomly selected background points (Phillips et al. 2006; Phillips & Dudík 2008). All models were generated for an area ranging from southern Honduras to northwestern Colombia. In all cases dataset was randomly divided into two sub-datasets, a training group containing 80% of the points and a test group containing 20% of the points for model validation. Although more locality points are preferred, Maxent models perform better than other algorithms with as few as 5 locality points (Hernandez et al. 2006; Pearson et al. 2007). Climatic data consisted of the 19 bioclimatic variables available in the WorldClim database with a 30 arc second resolution (Hijmans et al. 2005). These variables are derived from three climate parameters, maximum temperature, minimum temperature and rainfall, and are highly correlated (Xu & Hutchinson 2010). We used the receiver operating characteristic or area under the curve (AUC, Phillips et al. 2006) and the Akaike information criterion, AIC (Akaike A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. 1973), to compare the performance of models based on all 19 variables versus a reduced 7- variable model. AUC is expected to favor models with more parameters while AIC penalizes the use of unnecessary or unidentifiable parameters (Warren & Seifert 2011). The 7-variable model was obtained by first performing a two-tailed Spearman correlation analysis between all pairwise combinations of the 19 variables based on their values at 1000 points selected randomly from across the landscape. Variables with a correlation coefficient ≥ 0.75 were considered redundant and we removed the variable that least likely limited geographic distributions (e.g., removing mean values rather than maxima or minima; Wiens et al. 2006). Both the AUC and AIC criteria suggested that the full 19-variable model performed better (Fig. S1), thus we performed all further analyses on the 19-variable niche models. All 30 resulting ENMs yielded AUC scores above 0.7 and were considered good predictors of species distributions (Elith et al. 2006). Palaeo-climatic as well as current climatic conditions may predict the observed genetic divergence among conspecific populations. ENMs based on current distributional records can be projected onto spatial hypotheses of environmental conditions during the last glacial maximum (LGM) or other geological time periods (Martínez-Meyer et al. 2004; Peterson et al. 2004). We therefore projected or ‘back-cast’ species distribution models to the LGM, 21,000 years before present under two alternative historical scenarios generated using different General atmospheric Circulation Models (GCM) available in the WorldClim database (Hijmans et al. 2005): the Community Climate System Model (CCSM, Collins et al. 2006) and the Model for Interdisciplinary Research on Climate (MIROC, K-1-Model-Developers 2004). A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Elevation and geographic range Using georeferenced localities for each species we extracted elevation data from the shuttle radar topography mission (SRTM) digital elevation model (DEM) with a 90 m resolution (Jarvis et al. 2008). We calculated mean elevation and altitudinal range based on the central 90% of observations (i.e., we discarded the extreme upper and lower 5% of elevation records) in order to minimize the influence of outliers that may potentially represent sink populations atypical of the species in question (Table S3). Geographic range for all 31 species was computed from the species range shapefiles distributed by the Global Amphibian Assessment program (IUCN 2011). Landscape resistance When a species confronts a given landscape, the strength of the potential barrier will depend on the species’ intrinsic ability to move through the environment and withstand the prevailing abiotic conditions (Janzen 1967; Navas 2002). If an environment is suitable for a given species then the passage between two points will be easier and, ignoring biotic factors such as potential competitors, demographic connectivity across the landscape should be higher. If abiotic climatic requirements were relevant to explaining variation among species in phylogeographic structure, then we would expect the degree of genetic isolation of conspecific populations to co- vary inversely with degree of habitat suitability among species. To measure the effect of the landscape as a barrier to gene flow one can use an approach based on circuit theory. An ENM assigns to each pixel a probability of detection of a given A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. species. If we assume that probability of detection is a measure of habitat suitability, then a circuit theory-based approach estimates the effective resistance (or unsuitability) of the landscape to the passage of an organism between two focal regions (McRae 2006; McRae et al. 2008). Although habitat unsuitability is not necessarily equivalent to resistance to movement, this is a reasonable assumption for frogs with low dispersal capabilities and a high dependency on their local environment. This measure of resistance takes into account many possible pathways between two localities and can be interpreted as a measure of the effective environmental distance between foci for a given organism (McRae et al. 2008). The analysis was performed using the software, Circuitscape 3.5.4 (McRae 2006; McRae et al. 2008; Shah & McRae 2008) and comparisons were done between two pairs of foci: El Copé versus Brewster and Brewster versus Cana (Fig. 1, Table S3). Life history variables The ability of organisms to move across a landscape, survive and successfully reproduce may be determined by a variety of life-history traits. In the case of amphibians, cohorts of individuals may need to move long distances while resisting unfavorably dry conditions, to locate potential breeding sites, and maximize the survival of eggs and juvenile stages. We selected three variables available in the literature that potentially influence the ability of amphibians to disperse, survive and reproduce in a marginal environment: reproductive mode, body size and breeding habitat (Table 1). We predicted that species with larval stages may have lower dispersal ability given their higher dependence on water bodies relative to direct-developing frogs, and that larger bodied species will have higher dispersal ability than smaller ones (Austin et al. 2004; Bocxlaer A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. et al. 2010). On the other hand, pond breeders may be expected to move among ponds thus having high migration rates relative to terrestrial breeders that may show higher territorial behavior and thus lower migration (Zeisset & Beebee 2008), although we are not aware of any previous studies that have made an explicit comparison across a common landscape. Breeding habitat and reproductive mode were treated as a single variable, however, since they were tightly correlated within our dataset (i.e., all direct-developers are terrestrial breeders). In addition to the above life-history variables, we included the biogeographic origin of each lineage (North or South; Table 1), as recent arrivals on the landscape might be expected to show less genetic structure than older lineages (Cody et al. 2010). Geological and paleontological reconstructions suggest that northern lineages could have been present on the Isthmus of Panama well before southern lineages (Whitmore & Stewart 1965; Kirby & MacFadden 2005). Ecological and life history determinants of genetic diversity We used a set of nine predictive variables and one response variable, with continuous variables transformed to conform to normality assumptions (Table 1). Based on biological criteria we defined 28 a priori general lineal models (GLM) to explain the observed range of genetic divergences among species (Table 2). Models varied in the number and combination of predictor variables and in the presence of interactions between variables. Linear models were evaluated using the software R 2.14.1 (R Core Team 2013). To evaluate which variables were supported by the data and might be more relevant in determining genetic divergence we performed an AIC model selection procedure among all models (Table 2) using the base R package. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Results Hierarchical approximate Bayesian computation Although sample size (n) varied among species and among sites, n did not significantly predict average pairwise distance between sites (πB, P-value = 0.48, r2 = 0.011). For pairwise comparisons, El Copé vs. Brewster and Brewster vs. Cana the posterior probability for the hypothesis of one simultaneous divergence time (Ψ=1) was nearly 0 in both cases and was thus rejected (πB ranged from 0.001438—0.152888 and from 0.005482—0.111562 respectively for the two comparisons), implying that the intraspecific divergence times of the various species were spread among multiple temporal intervals (Fig. 2). For the comparison of El Copé vs. Brewster, Bayes factor (BF) analysis showed highest support for three divergence time intervals among 30 taxa (BF = 109, Fig. S2a). For the Cana vs. Brewster comparison among 15 taxa, four divergence intervals received the highest support (BF = 3.4, Fig. S2b). To determine the specific taxa grouped in each divergence time interval, we used the BF results to constrain the number of intervals to their most likely values and re-ran the analysis (Table S4, Hickerson & Meyer 2008). The timing and specific membership of each interval are given in Figs. S3 and S4 and Table S4). ENMs and landscape resistance The effective environmental resistance was estimated for each species using the present ENM and two projections to the LGM (Fig. S5, CCSM and MIROC, see Methods). Resistance for both past models is positively and significantly correlated with each other (R = 0.62, P-value = 1.04×10-5). Present resistance is positively and significantly correlated with past resistance A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. both with the MIROC (R = 0.66, P-value = 1.13×10-6) and the CCSM models (R = 0.55, P-value = 1.50×10-4). For those species present at all three sites, present and past landscape resistance was generally higher between Cana and Brewster than between El Copé and Brewster (Fig. 3), suggesting the environment in eastern Panama has been more inhospitable for amphibians than central Panama (Fig. 1). Ecological and life history determinants of genetic diversity Our model selection procedure identified a clearly preferred (ΔAICc = 6.43) best-fit model containing five variables: maximum male body size, geographic range, landscape resistance in the present, reproductive mode and biogeographic origin and an interaction term between biogeographic origin and reproductive mode (Table 2). This model explains a large portion of the variation (R2 = 0.64, Adjusted R2=0.57). The second best model lacked the interaction term and included mean elevation. Models with fewer variables but that included past climate (MIROC and CCSM), mean elevation or elevational range tended to perform poorly, suggesting that these parameters may be less informative in predicting genetic divergence across the landscape. Body size and reproductive mode tended to occur across the best-performing models (Table 2). Discussion Our results show that, indeed, variation in life-history traits may explain variation in phylogeographic structure among species. HABC analyses show that among-species variation in A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. intraspecific divergence could not be explained by a common barrier or by coalescent stochasticity, as among species the observed range of genetic divergences between conspecific population pairs was too great (Fig. 2). If the genetic divergence between populations was strictly stochastic in nature, or caused by independent biogeographic histories of species, we would not expect a relationship between ecological variables and genetic divergence. Instead, ecological traits including body size and reproductive mode explain a large portion (R2 = 0.64, Adjusted R2 = 0.57) of the genetic variation among our three frog communities. The significant association of genetic divergence and landscape resistance as estimated from niche modeling also suggests that ENM successfully identifies environmental variation important to the evolution of diverse species of amphibians (Glor & Warren 2011). Thus, by effectively holding geography constant and sampling many taxa from the same geographic localities, we have successfully integrated ecological and life-history variables into a model explaining population divergence (Johnson et al. 2009; Buckley 2009). Only one variable showed a phylogenetic signal thus phylogenetic correlations did not affect the present analysis, but may need to be addressed in future studies. Subsequent studies should include additional life-history traits, such as metabolic rate or desiccation resistance; as such data are not yet available for most species studied here. Our hope is that this study of Panamanian frogs may provide a general approach of wide applicability, as obvious geographic features have failed to constitute a general explanation for observed divergences among allopatric populations, especially in the tropics (Antonelli et al. 2010). Vicariant obstacles such as mountains or rivers certainly divide many phylogeographic clades and sister species (Avise et al. 1987; Ribas et al. 2012), yet there may be a great many pairs of allopatric clades for which there are no apparent barriers separating them. Furthermore, A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. many physical barriers such as mountains do not impede gene flow directly, rather they promote environmental and habitat gradients that may be inhospitable to the dispersing organism (Janzen 1967). The dichotomy between physical and environmental barriers may be a false one, at least in some cases. Recent research has focused not on the type of barrier but on the type of mechanism impeding gene flow, i.e., isolation by distance (IBD) or by physical impassibility versus isolation by environment (IBE) or ‘immigrant unviability,’ in which natural selection acts against migrating individuals (Nosil et al. 2005; Bradburd et al. 2013; Wang 2013). The IBE model highlights why co-distributed species may vary widely in their phylogeographic structure, while our comparative approach using general linear modeling allows one to evaluate specific ecological or physiological traits that may contribute to divergence. Our approach to investigating the effect of ecology on phylogeography, which we refer to as ‘ecophylogeography,’ follows Pabijan et al. (2012), with some modifications. We employed an HABC analysis to first determine whether the variation in mtDNA divergence cannot be explained simply by coalescent stochasticity in the context of isolation by distance or isolation by physical barrier (Hickerson et al. 2007). In addition to life-history variables, we incorporated environmental heterogeneity (IBE) as an explanatory variable through ENM and ‘isolation by resistance’ analyses (McRae 2006). Our dataset also covered a larger phylogenetic diversity, including species from several taxonomic families, which may explain why we recovered a greater number of significant explanatory variables in addition to body size (Pabijan et al. 2012; Wollenberg Valero 2015). A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Our approach does not supplant other ecophylogeographic inference methods but does offer some advantages. While landscape genetic studies frequently uncover an association between population genetic structure and environmental heterogeneity (Storfer et al. 2007; Richards-Zawacki 2009), our comparative phylogeographic approach allows one to search for common patterns across species and quantify to what extent phylogeographic patterns may be predictable from a knowledge of organismal biology. Previous studies have shown the global importance of dispersal potential in explaining genetic divergence particularly for marine and freshwater organisms as well as some terrestrial invertebrates (Bohonak 1999), however the importance of continuous variation in specific ecological and life-history traits has not been evaluated. Using a phylogeographic approach we can infer the action of evolutionary processes acting over longer time scales than those revealed by landscape genetics (Wang 2010; Lawson 2013). By including explicit life-history characteristics in a general linear modeling exercise, our approached evaluated the contribution of specific organismal traits to genetic divergence (Pabijan et al. 2012) coupled with a contribution of IBE (Wang 2013). We estimated phylogeographic structure as genetic divergence at mitochondrial DNA (mtDNA). While variation in mtDNA may not be representative of variation in the nuclear genome (Fay & Wu 1999; Bazin et al. 2006; Toews & Brelsford 2012), this standardized marker is widely used in comparative phylogeographic studies (Carnaval et al. 2009; Wang et al. 2013) due to its particular evolutionary genetic properties. High mutations rates of anuran mtDNA allow one to estimate recent divergence times among species (Crawford 2003). Smaller effective population sizes (Ne) permit a rapid achievement of reciprocal monophyly between populations (Hudson & Coyne 2002), as observed here. The HABC approach takes into account the A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. stochasticity inherent in single-locus coalescence, and moreover, estimation of population divergence time E(τ) does not improve appreciably until eight nuclear loci are included (Huang et al. 2011). While mtDNA datasets are still the norm for comparative ecophylogeographic analyses (Pabijan et al. 2012; Wang et al. 2013), genomic-scale phylogeographic datasets may soon be routinely available (McCormack et al. 2013). Although ecology was a main predictor of genetic divergence, our model also included an effect of biogeographic origin of species. Species of northern origin showed significantly less genetic divergence across the same landscape than southern lineages, which is the opposite of what we predicted given that southern lineages would be more recent arrivals and would therefore have had less time to build up spatial genetic structure (Weir et al. 2009). Alternatively, many southern lineages likely arrived in Panama perhaps across an early Isthmian connection (Montes et al. 2012) around eight million years ago (Pinto-Sánchez et al. 2012), which may have been sufficient time to erase the potential genetic imprint of recent colonization. We found that species ranges predicted genetic divergence, thus revealing that species with wider geographic ranges have less genetic divergence between populations (Duminil et al. 2007). Wider geographic ranges are indicative of habitat generalists that by definition can tolerate a wider range of environments and should experience lower IBE and increased dispersal. Although this pattern was not detected in Malagasy frogs (Pabijan et al. 2012), perhaps because that study covered less phylogenetic diversity, this relationship should be relatively straight- forward to test in other species. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. The relevance of life history variables in explaining phylogeographic structure has not been widely tested. In plants, life-history traits such as breeding system, life form and seed mass are known to influence population structure (Hamrick & Godt 1996; Duminil et al. 2007). For marine invertebrates, developmental mode is thought to be important in limiting dispersal and promoting population structure (Arndt & Smith 1998; Collin 2001). Habitat preference in birds and more generally their ability to move through the landscape has been shown to predict genetic differentiation and eventually speciation (Burney & Brumfield 2009; Smith et al. 2014). We found two life history characteristics of anurans that are important in predicting genetic divergence. Larger animals showed less genetic divergence among populations, implying they have higher dispersal capacity than smaller species (Pabijan et al. 2012; Wollenberg Valero 2015). Larger frogs may be better dispersers simply because of scale (they can move faster than smaller frogs) or may be better survivors, e.g., better able to resist desiccation. Reproductive mode was also a significant predictor of genetic divergence, but the effect was in the opposite direction. We predicted that direct developers, being independent of bodies of water for reproduction, would have greater dispersal capabilities, yet it was the toads and treefrogs with larval stages, which showed less genetic divergence and, by implication, greater dispersal potential. Higher phylogeographic structure could be due to low dispersal associated with site faithfulness and territorial behavior in terrestrial frogs (Stewart & Rand 1991; Ryan et al. 2008; Fouquet et al. 2012). If true, then behavioral ecology may be another variable that could be taken into account in ecophylogeographic studies, though this may be difficult using the comparative approach advocated here. Alternatively, the direct developers may simply show greater A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. vulnerability to desiccation and greater niche conservatism, even though they do not need bodies of water to reproduce. Using our comparative multivariate approach we were able to study both intrinsic (life- history) and extrinsic factors (environmental resistance obtained from ENMs). We recognize, however, that IBE can be interpreted as the interaction between an organism’s intrinsic characteristics and the abiotic landscape it confronts. Because we looked at variation in genetic divergence across a common landscape, in effect we held the landscape constant to reveal variation among species. Thus, we interpret the landscape resistance variable more as a composite property of the organism. The differing responses of organisms to a common landscape has been suggested by a few studies using pairs of species (Goldberg & Waits 2010; Richardson 2012), but rarely is such a high number of multiple co-distributed species synthesized in one analysis. To successfully delimit intrinsic versus extrinsic factors and quantify their relative roles in determining how populations interact with landscape characteristics, biologists may need to incorporate direct studies of organismal physiology, ecomorphology and even behavior (Janzen 1967; Navas 2002; Bernardo et al. 2007). Global DNA barcoding initiatives have begun to provide genetic data on whole communities and regional biotas (Kerr et al. 2007; Dick & Kress 2009; Murphy et al. 2013). Such massive datasets should provide opportunities to implement the present comparative approached over much larger regions and more diverse species. These data will lend themselves naturally to a phylogeographic approach (Hajibabaei et al. 2007). Looking beyond the ENM-based approach (Richards et al. 2007), comparative phylogeography should benefit greatly from the direct integration of organismal biology, which A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. may interface with new approaches to studying adaptive genetic variation in a phylogeographic context (Hickerson et al. 2010; Günther & Coop 2013). Acknowledgments Specimens were collected under Panama’s Autoridad Nacional del Ambiente permit number SE/A-37-07 to R.I.D. Field work was funded by National Geographic CRE grant 8133- 06 to A.J.C., K.R.L., and R.I.D. DNA sequence data were kindly obtained by A. C. Driskell at the Smithsonian Institution’s Laboratories of Analytical Biology. We are grateful to C. Ortiz and J. Velásquez for help with ENM analyses. For valuable comments on this work we are thankful to M. Hickerson, C. D. Cadena, N. R. Pinto, L. S. Barrientos, C. E. Guarnizo, L. F. Dueñas, K. C. Wollenberg Valero. To T. Vines, L. Rissler and two anonymous reviewers for feedback via Axios Review and to three anonymous reviewers from this journal for valuable feedback. References Akaike H (1973) Information Theory and an Extension of the Maximum Likelihood Principle. In: International Symposium on Information Theory, 2 nd, Tsahkadsor, Armenian SSR , pp. 267–281. Antonelli A, Quijada-Mascareñas A, Crawford AJ et al. (2010) Molecular studies and phylogeography of Amazonian tetrapods and their relation to geological and climatic models. In: Amazonia, Landscape and Species Evolution: a Look into the Past (eds Hoorn C, Wesselingh F), pp. 386–404. Wiley-Blackwell, Oxford, UK. Arbogast BS, Kenagy GJ (2001) Comparative phylogeography as an integrative approach to historical biogeography. Journal of Biogeography, 28, 819–825. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Arndt A, Smith MJ (1998) Genetic diversity and population structure in two species of sea cucumber: differing patterns according to mode of development. Molecular Ecology, 7, 1053–1064. Austin JD, Lougheed SC, Boag PT (2004) Discordant temporal and geographic patterns in maternal lineages of eastern north American frogs, Rana catesbeiana (Ranidae) and Pseudacris crucifer (Hylidae). Molecular Phylogenetics and Evolution, 32, 799–816. Avise JC (1992) Molecular population structure and the biogeographic history of a regional fauna: a case history with lessons for conservation biology. Oikos, 63, 62–76. Avise JC (2000) Phylogeography: The History and Formation of Species. Harvard University Press, Cambridge, Massachusetts. Avise JC, Arnold J, Ball RM et al. (1987) Intraspecific phylogeography: The mitochondrial DNA bridge between population genetics and systematics. Annual Review of Ecology and Systematics, 18, 489–522. Ayres JM, Clutton-Brock TH (1992) River boundaries and species range size in Amazonian primates. The American Naturalist, 140, 531–537. Bates JM, Haffer J, Grismer E (2004) Avian mitochondrial DNA sequence divergence across a headwater stream of the Rio Tapajós, a major Amazonian river. Journal of Ornithology, 145, 199–205. Bazin E, Glemin S, Galtier N (2006) Population size does not influence mitochondrial genetic diversity in animals. Science, 312, 570–572. Bermingham E, Avise JC (1986) Molecular Zoogeography of freshwater fishes in the southeastern United States. Genetics, 113, 939–965. Bermingham E, Martin AP (1998) Comparative mtDNA phylogeography of neotropical freshwater fishes: testing shared history to infer the evolutionary landscape of lower Central America. Molecular Ecology, 7, 499–517. Bernardo J, Ossola RJ, Spotila J, Crandall KA (2007) Interspecies physiological variation as a tool for cross-species assessments of global warming-induced endangerment: validation of an intrinsic determinant of macroecological and phylogeographic structure. Biology Letters, 3, 695–699. Bocxlaer I Van, Loader SP, Roelants K et al. (2010) Gradual adaptation toward a range- expansion phenotype initiated the global radiation of toads. Science, 327, 679–682. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Bohonak AJ (1999) Dispersal, gene flow, and population structure. The quarterly review of biology, 74, 21–45. Bradburd GS, Ralph PL, Coop GM (2013) Disentangling the effects of geographic and ecological isolation on genetic differentiation. Evolution, 67, 3258–3273. Buckley D (2009) Toward an organismal, integrative, and iterative phylogeography. BioEssays, 31, 784–793. Burney CW, Brumfield RT (2009) Ecology predicts levels of genetic differentiation in neotropical birds. The American naturalist, 174, 358–368. Byrne M (2008) Evidence for multiple refugia at different time scales during Pleistocene climatic oscillations in southern Australia inferred from phylogeography. Quaternary Science Reviews, 27, 2576–2585. Carnaval AC, Hickerson MJ, Haddad CFB, Rodrigues MT, Moritz C (2009) Stability predicts genetic diversity in the Brazilian Atlantic forest hotspot. Science, 323, 785–789. Cody S, Richardson JE, Rull V, Ellis C, Pennington RT (2010) The Great American Biotic Interchange revisited. Ecography, 33, 326–332. Collin R (2001) The effects of mode of development on phylogeography and populations structure of North Atlantic Crepidula (Gastropoda: Calyptraeidae). Molecular Ecology, 10, 2249–2262. Collins WD, Bitz CM, Blackmon ML et al. (2006) The Community Climate System Model Version 3 (CCSM3). Journal of Climate, 19, 2122–2143. Crawford AJ (2003) Relative rates of nucleotide substitution in frogs. Journal of Molecular Evolution, 57, 636–641. Crawford AJ, Bermingham E, Polania C (2007) The role of tropical dry forest as a long-term barrier to dispersal: a comparative phylogeographical analysis of dry forest tolerant and intolerant frogs. Molecular Ecology, 16, 4789–4807. Crawford AJ, Lips KR, Bermingham E (2010) Epidemic disease decimates amphibian abundance, species diversity, and evolutionary history in the highlands of central Panama. Proceedings of the National Academy of Sciences USA, 107, 13777–13782. Crump ML (1974) Reproductive strategies in a tropical anuran community. Miscellaneous Publications - Museum of Natural History, University of Kansas, 61, 1–67. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Dick CW, Kress WJ (2009) Dissecting tropical plant diversity with forest plots and a molecular toolkit. BioScience, 59, 745–755. Duellman WE (2001) The Hylid Frogs of Middle America. Society for the Study of Amphibians and Reptiles ; University of Saint Louis in cooperation with the Natural History Museum of the University of Kansas. Duminil J, Fineschi S, Hampe A et al. (2007) Can population genetic structure be predicted from life-history traits? The American Naturalist, 169, 662–672. Elith J, Graham CH, Anderson RP et al. (2006) Novel methods improve prediction of species’ distributions from occurrence data. Ecography, 29, 129–151. Elith J, Leathwick JR (2009) Species distribution models: Ecological explanation and prediction across space and time. Annual Review of Ecology, Evolution, and Systematics, 40, 677–697. Fay JC, Wu C-I (1999) A human population bottleneck can account for the discordance between patterns of mitochondrial versus nuclear DNA variation. Molecular Biology and Evolution, 16, 1003–1005. Fine PVA, Zapata F, Daly DC et al. (2013) The importance of environmental heterogeneity and spatial distance in generating phylogeographic structure in edaphic specialist and generalist tree species of Protium (Burseraceae) across the Amazon Basin. Journal of Biogeography, 40, 646–661. Fouquet A, Ledoux J-B, Dubut V, Noonan BP, Scotti I (2012) The interplay of dispersal limitation, rivers, and historical events shapes the genetic structure of an Amazonian frog. Biological Journal of the Linnean Society, 356–373. Gascon C, Malcolm JR, Patton JL et al. (2000) Riverine barriers and the geographic distribution of Amazonian species. Proceedings of the National Academy of Sciences of the United States of America, 97, 13672–13677. Gaston KJ (2000) Global patterns in biodiversity. Nature, 405, 220–227. Gehring PS, Pabijan M, Randrianirina JE, Glaw F, Vences M (2012) The influence of riverine barriers on phylogeographic patterns of Malagasy reed frogs (Heterixalus). Molecular Phylogenetics and Evolution, 64, 618–632. Glor RE, Warren D (2011) Testing ecological explanations for biogeographic boundaries. Evolution, 65, 673–683. Goldberg CS, Waits LP (2010) Comparative landscape genetics of two pond-breeding amphibian species in a highly modified agricultural landscape. Molecular ecology, 19, 3650–63. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Günther T, Coop G (2013) Robust identification of local adaptation from allele frequencies. Genetics, 195, 205–220. Haffer J (1969) Speciation in amazonian forest birds. Science, 165, 131–137. Hajibabaei M, Singer GAC, Hebert PDN, Hickey DA (2007) DNA barcoding: how it complements taxonomy, molecular phylogenetics and population genetics. Trends in Genetics, 23, 167–172. Hamrick JL, Godt MJW (1996) Effects of Life History Traits on Genetic Diversity in Plant Species. Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences, 351, 1291–1298. Hernandez PA, Graham CH, Master LL, Albert DL (2006) The effect of sample size and species characteristics on performance of different species distribution modeling methods. Ecography, 29, 773–785. Hewitt GM (2000) The genetic legacy of the Quarternary ice ages. Nature, 405, 907–913. Hickerson MJ, Carstens BC, Cavender-Bares J et al. (2010) Phylogeography’s past, present, and future: 10 years after Avise 2000. Molecular Phylogenetics and Evolution, 54, 291–301. Hickerson MJ, Dolman G, Moritz C (2006a) Comparative phylogeographic summary statistics for testing simultaneous vicariance. Molecular Ecology, 15, 209–223. Hickerson M, Meyer C (2008) Testing comparative phylogeographic models of marine vicariance and dispersal using a hierarchical Bayesian approach. BMC Evolutionary Biology, 8, 322. Hickerson MJ, Stahl EA, Lessios HA (2006b) Test for simultaneous divergence using approximate Bayesian computation. Evolution, 60, 2435–2453. Hickerson MJ, Stahl E, Takebayashi N (2007) msBayes: Pipeline for testing comparative phylogeographic histories using hierarchical approximate Bayesian computation. BMC Bioinformatics, 8, 268. Hickerson MJ, Stone GN, Lohse K et al. (2014) Recommendations for using MsBayes to incorporate uncertainty in selecting an ABC model prior: A response to Oaks et al. Evolution, 68, 284–294. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A (2005) Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology, 25, 1965–1978. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Huang W, Takebayashi N, Qi Y, Hickerson M (2011) MTML-msBayes: Approximate Bayesian comparative phylogeographic inference from multiple taxa and multiple loci with rate heterogeneity. BMC Bioinformatics, 12, 1. Hudson RR, Coyne JA (2002) Mathematical consequences of the genealogical species concept. Evolution, 56, 1557–1565. IUCN (2011) IUCN Red List of Threatened Species. Version 2011.2. Janzen DH (1967) Why mountain passes are higher in the tropics. The American Naturalist, 101, 233. Jarvis A, Reuter HI, Nelson A, Guevara E (2008) Hole-filled SRTM for the globe Version 4, available from the CGIAR-CSI SRTM 90m Database (http://srtm.csi.cgiar.org). Johnson NK, Cicero C (2004) New mitochondrial DNA data affirm the importance of pleistocene speciation in North American birds. Evolution, 58, 1122–1130. Johnson JB, Peat SM, Adams BJ (2009) Where’s the ecology in molecular ecology? Oikos, 118, 1601–1609. K-1-Model-Developers (2004) K-1 Coupled GCM (Miroc Description)1. Kerr KCR, Stoeckle MY, Dove CJ et al. (2007) Comprehensive DNA barcode coverage of North American birds. Molecular Ecology Notes, 7, 535–543. Kirby MX, MacFadden B (2005) Was southern Central America an archipelago or a peninsula in the middle Miocene? A test using land-mammal body size. Palaeogeography, Palaeoclimatology, Palaeoecology, 228, 193–202. Knowles LL, Maddison WP (2002) Statistical phylogeography. Molecular Ecology, 11, 2623– 2635. Lawson LP (2013) Diversification in a biodiversity hot spot: landscape correlates of phylogeographic patterns in the African spotted reed frog. Molecular Ecology, 22, 1947– 1960. Lessa EP, Cook JA, Patton JL (2003) Genetic footprints of demographic expansion in North America, but not Amazonia, during the Late Quaternary. Proceedings of the National Academy of Sciences of the United States of America, 100, 10331–10334. Marko PB (2002) Fossil calibration of molecular clocks and the divergence times of geminate species pairs separated by the Isthmus of Panama. Molecular Biology and Evolution, 19, 2005–2021. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Martínez-Meyer E, Townsend Peterson A, Hargrove WW (2004) Ecological niches as stable distributional constraints on mammal species, with implications for Pleistocene extinctions and climate change projections for biodiversity. Global Ecology and Biogeography, 13, 305–314. Mayr E (1942) Systematics and the origin of species: from the viewpoint of a zoologist. Columbia University Press, New York, New York, USA. McCormack JE, Hird SM, Zellmer AJ, Carstens BC, Brumfield RT (2013) Applications of next- generation sequencing to phylogeography and phylogenetics. Molecular Phylogenetics and Evolution, 66, 526–538. McRae BH (2006) Isolation by resistance. Evolution, 60, 1551–1561. McRae BH, Dickson BG, Keitt TH, Shah VB (2008) Using circuit theory to model connectivity in ecology, evolution, and conservation. Ecology, 89, 2712–2724. Montes C, Cardona a., McFadden R et al. (2012) Evidence for middle Eocene and younger land emergence in central Panama: Implications for Isthmus closure. Geological Society of America Bulletin, 124, 780–799. Murphy RW, Crawford AJ, Bauer AM et al. (2013) Cold Code: the global initiative to DNA barcode amphibians and nonavian reptiles. Molecular Ecology Resources, 13, 161–167. Navas CA (2002) Herpetological diversity along Andean elevational gradients: links with physiological ecology and evolutionary physiology. Comparative Biochemistry and Physiology - Part A: Molecular & Integrative Physiology, 133, 469–485. Navas CA, Otani L (2007) Physiology, environmental change, and anuran conservation. Phyllomedusa, 6, 83–103. Nei M, Li W-H (1979) Mathematical model for studying variation in terms of restriction endonucleases. Proceedings of the National Academy of Sciences USA, 76, 5269–5273. Nielsen R, Beaumont MA (2009) Statistical inference in phylogeography. Molecular Ecology, 18, 1034. Nosil P, Vines TH, Funk DJ (2005) Reproductive isolation caused by natural selection against immigrants from divergent habitats. Evolution, 59, 705–719. Oaks J, Sukumaran J, Esselstyn J et al. (2012) Evidence for climate-driven diversification? A caution for interpreting ABC inferences of simultaneous historical events. Evolution, 67, 991–1010. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Pabijan M, Wollenberg KC, Vences M (2012) Small body size increases the regional differentiation of populations of tropical mantellid frogs (Anura: Mantellidae). Journal of Evolutionary Biology, 25, 2310–2324. Pearson RG, Raxworthy CJ, Nakamura M, Peterson AT (2007) Predicting species distributions from small numbers of occurrence records: a test case using cryptic geckos in Madagascar. Journal of Biogeography, 34, 102–117. Peterson AT (2001) Predicting species´ geographic distributions based on ecological niche modeling. The Condor, 103, 599–605. Peterson AT, Martínez-Meyer E, González-Salazar C (2004) Reconstructing the Pleistocene geography of the Aphelocoma jays (Corvidae). Diversity and Distributions, 10, 237–246. Phillips SJ, Anderson RP, Schapire RE (2006) Maximum entropy modeling of species geographic distributions. Ecological Modelling, 231–259. Phillips SJ, Dudík M (2008) Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography, 31, 161–175. Pinto-Sánchez NR, Ibáñez R, Madriñán S et al. (2012) The Great American Biotic Interchange in frogs: Multiple and early colonization of Central America by the South American genus Pristimantis (Anura: Craugastoridae). Molecular Phylogenetics and Evolution, 62, 954–972. R Core Team (2013) R: A language and environment for statistical computing. Ribas CC, Aleixo A, Nogueira ACR, Miyaki CY, Cracraft J (2012) A palaeobiogeographic model for biotic diversification within Amazonia over the past three million years. Proceedings of the Royal Society B: Biological Sciences, 279, 681–689. Richards CL, Carstens BC, Knowles LL (2007) Distribution modelling and statistical phylogeography: an integrative framework for generating and testing alternative biogeographical hypotheses. Journal of Biogeography, 34, 1833–1845. Richards-Zawacki CL (2009) Effects of slope and riparian habitat connectivity on gene flow in an endangered Panamanian frog, Atelopus varius. Diversity and Distributions, 15, 796–806. Richardson JL (2012) Divergent landscape effects on population connectivity in two co-occurring amphibian species. Molecular ecology, 21, 4437–51. Robertson JM, Duryea MC, Zamudio KR (2009) Discordant patterns of evolutionary differentiation in two Neotropical treefrogs. Molecular Ecology, 18, 1375–1395. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Rull V (2008) Speciation timing and neotropical biodiversity: the Tertiary–Quaternary debate in the light of molecular phylogenetic evidence. Molecular Ecology, 17, 2722–2729. Ryan MJ, Lips KR, Eichholz MW (2008) Decline and extirpation of an endangered Panamanian stream frog population (Craugastor punctariolus) due to an outbreak of chytridiomycosis. Biological Conservation, 141, 1636–1647. Shah VB, McRae B (2008) Circuitscape: A Tool for Landscape Ecology. In: Proceedings of the 7th Python in Science Conference (eds Varoquaux G, Millman J), pp. 62–65. Pasadena, CA USA. Smith BT, McCormack JE, Cuervo AM et al. (2014) The drivers of tropical speciation. Nature. Soberón J, Nakamura M (2009) Niches and distributional areas: Concepts, methods, and assumptions. Proceedings of the National Academy of Sciences, 106, 19644–19650. Soberón J, Peterson AT (2005) Interpretation of models of fundamental ecological niches and species’ distributional areas. Biodiversity Informatics, 2, 1–10. Stewart MM, Rand AS (1991) Vocalizations and the defense of retreat sites by male and female frogs, Eleutherodactylus coqui. Copeia, 4, 1013–1024. Storfer A, Murphy MA, Evans JS, Goldberg CS, Robinson S (2007) Putting the landscape in landscape genetics. Heredity, 98, 128. Tajima F (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics, 123, 585–595. Takahata N, Nei M (1985) Gene genealogy and variance of interpopulational nucleotide differences. Genetics, 110, 325–344. Toews DPL, Brelsford A (2012) The biogeography of mitochondrial and nuclear discordance in animals. Molecular Ecology, 21, 3907–3930. Wallace AR (1852) On the monkeys of the Amazon. Proceedings of the Zoological Society of London, 20, 107–110. Waltari E, Hijmans RJ, Peterson AT et al. (2007) Locating pleistocene refugia: comparing phylogeographic and ecological niche model predictions. PLoS ONE, 2, e563. Wang IJ (2010) Recognizing the temporal distinctions between landscape genetics and phylogeography. Molecular Ecology, 19, 2605–2608. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Wang IJ (2013) Examining the full effects of landscape heterogeneity on spatial genetic variation: a multiple matrix regression approach for quantifying geographic and ecological isolation. Evolution, 67, 3403–3411. Wang IJ, Crawford AJ, Bermingham E (2008) Phylogeography of the Pygmy Rain Frog (Pristimantis ridens) across the lowland wet forests of isthmian Central America. Molecular Phylogenetics and Evolution, 47, 992–1004. Wang IJ, Glor RE, Losos JB (2013) Quantifying the roles of ecology and geography in spatial genetic divergence (F Adler, Ed,). Ecology Letters, 16, 175–182. Warren DL, Seifert SN (2011) Ecological niche modeling in Maxent : the importance of model complexity and the performance of model selection criteria. Ecological Applications, 21, 335–342. Weigt LA, Crawford AJ, Rand AS, Ryan MJ (2005) Biogeography of the túngara frog, Physalaemus pustulosus: a molecular perspective. Molecular Ecology, 14, 3857–3876. Weir JT, Bermingham E, Schluter D (2009) The Great American Biotic Interchange in birds. Proceedings of the National Academy of Sciences, 105, 21737–21742. Wells KD (2007) The Ecology and Behavior of Amphibians. The University of Chicago Press, Chicago. Whitmore FC, Stewart RH (1965) Miocene mammals and Central American seaways: Fauna of the Canal Zone indicates separation of Central and South America during most of the Tertiary. Science, 148, 180–185. Wiens JJ (2004) Speciation and ecology revisited: Phylogenetic niche conservatism and the origin of species. Evolution, 58, 193–197. Wiens J, Graham C, Moen D, Smith S, Reeder T (2006) Evolutionary and Ecological Causes of the Latitudinal Diversity Gradient in Hylid Frogs: Treefrog Trees Unearth the Roots of High Tropical Diversity. The American Naturalist, 168, 579–596. Wollenberg Valero KC (2015) Evolutionary and Population Genetics Evidence for an intrinsic factor promoting landscape divergence in Madagascan leaf-litter frogs. Frontiers in Genetics, 6. Wright S (1951) The genetical structure of populations. Annals of Eugenics, 15, 323–354. Xu T, Hutchinson M (2010) ANUCLIM User Guide, Version 6.1. Australian National University, Canberra. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Zeisset I, Beebee TJC (2008) Amphibian phylogeography: a model for understanding historical aspects of species distributions. Heredity, 101, 109–119. Zink RM (1996) Comparative phylogeography in North American birds. Evolution, 50, 308–317. Data Accessibility DNA sequences: Genbank accession numbers for new sequences of 16S are KR863128 - KR863371 and KR911916-KR911918, and for new sequences of COI are: KR862873 - KR863113 and KR911913-KR911915. Genbank accession numbers for previously published sequences are provided in Table S2. Localities and climatic variables used for constructing environmental niche models, DNA sequence alignments per species per locality, R code for generating GLMs and database used in all analyses are available in the Dryad repository (DOI: 10.5061/dryad.t430k). Author Contributions AP and AJC designed the study. RI, KRL and AJC collected data. AP and AJC analyzed data and wrote the manuscript with vital input from RI and KRL. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Table 1. Summary of nine independent variables used in general linear modeling to predict the dependent variable, average between-population pairwise genetic distance (πB) in 31 species of frogs present in at least two of the three sampling sites referred to as El Copé, Brewster and Cana. Prior to statistical analyses (Table 2), six variables were transformed to meet assumptions of normality. Since reproductive mode and breeding habitat were tightly correlated in our dataset (all terrestrial breeders had direct development) we decided to use only reproductive mode. Variable Type Units Transformation Resistance present Continuous Ohm (Ω) None Resistance at LGM* (CCSM)† Continuous Ohm (Ω) Logarithmic Resistance at LGM* (MIROC)† Continuous Ohm (Ω) Logarithmic Body size (Maximum snout-vent-length) Continuous Millimeters (mm) Logarithmic Mean elevation Continuous Meters (m) Logarithmic Elevation range Continuous Meters (m) None Geographic range Continuous Square Kilometers (km2) Logarithmic Reproductive mode Categorical Direct/Indirect None Origin Categorical North/South None πB Continuous Average between- population pairwise genetic distance Square root *Last Glacial Maximum at 21,000 years ago. † CCSM and MIROC are two different global circulation models used to generate predictions of past climate (see Methods). A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Table 2. 28 a priori models hypothesized to explain pairwise between-population genetic distance (πB) of 31 species of frogs across a common landscape, their biological interpretation, and resulting corrected Akaike Information Criterion (AICc) scores. Models involved combinations of nine independent variables (Table 1). Following the notation in the statistical package, R, a + indicates a linear combination of variables, while * between two variables indicates the same as + as well as an interaction between terms. Note, in R, AICc values are often negative, with the smallest AICc value still indicating the preferred model. CCSM and MIROC are two spatial models of environmental variation during the Last Glacial Maximum (LGM). See text for details. Model Biological Significance AICc ∆AICc Body size + Present landscape resistance + Geographic range + Origin * Reproductive mode Large body size, low present environmental resistance, and generalist life history predict low πb with contributions from biogeographic origin and its interaction with reproductive mode. -111.58 0.00 Present landscape resistance + Mean elevation + Geographic range + Body size + Reproductive mode + Origin Full model, minus past resistance variable. -105.15 6.43 Present landscape resistance + Landscape resistance at LGM (MIROC)+ Mean elevation + Geographic range + Body Ssize + Reproductive mode + Origin Full model, assuming past environment based on the MIROC model. -103.61 7.96 Present landscape resistance + Landscape resistance at LGM (CCSM) + Mean elevation + Geographic range + Body size + Reproductive mode + Origin Full model, assuming past environment based on the CCSM model. -101.94 9.64 Origin * Reproductive mode Southern frogs show distinctive life histories although there is a potential bias with 6 of 7orthern species being Craugastor. -96.72 14.86 Reproductive mode + Body size Intrinsic life-history traits alone determine πB. -91.96 19.62 A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Reproductive mode Higher dependence on water may drive species to move less and stay close to water bodies, leading to higher πB. -91.04 20.53 Origin * Geographic range * Body size Biogeographic factor interacts with range size and body size. -90.00 21.57 Geographic range Species with wider geographic ranges are habitat generalists that move more and thus have lower πB. -89.37 22.21 Body size Large species disperse more and thus should have lower πB. -89.33 22.25 Reproductive mode * Body size * Present landscape resistance life-history variables interact with each other and with landscape resistance. -88.65 22.93 Present landscape resistance + Geographic range + Elevation range Combined habitat generalist traits plus landscape resistance. -87.89 23.69 Body size + Present landscape resistance + Geographic range Larger species have less environmental restrictions and habitat generalists have more dispersal capacity. -87.74 23.84 Body size + Present landscape resistance Larger species have higher dispersal ability, and environmental resistance also plays a role. -87.70 23.87 Origin Biogeography as sole driver. -87.38 24.20 Present landscape resistance Present resistance alone explains πB. The interaction between the organism and its environment predicts genetic isolation. -86.27 25.30 Geographic range + Elevation range Species with wider geographic ranges and wider elevation ranges will have wider tolerances and thus be able to disperse more. -86.25 25.33 Body size + Present landscape resistance + Elevation range Larger species have less environmental restrictions and habitat generalists, as measured by elevation range, have more dispersal capacity. -86.16 25.42 Body size + Present landscape resistance + Geographic range + Origin Larger species have less environmental restrictions and habitat generalists have more dispersal capacity, while biogeographic origin also plays a role. -85.14 26.44 Mean elevation + Present landscape resistance Elevation and climatic resistance together predict isolation and divergence. -84.49 27.09 Landscape resistance at LGM (CCSM) Environmental resistance at the LGM (based on CCSM) is sufficient to explain πB. -84.11 27.47 A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Present landscape resistance + Landscape resistance at LGM (CCSM) Present and LGM (CCSM) environments together predict genetic divergence. -84.04 27.53 Present landscape resistance +Landscape resistance at LGM (MIROC) Present and LGM (MIROC) environments together predict genetic divergence. -83.90 27.68 Landscape resistance at LGM (MIROC) Environmental resistance at the LGM (based on MIROC) explains πB. -83.29 28.29 Body size + Present landscape resistance + Geographic range + Mean elevation * Origin Large body size, low present environmental resistance, and generalist life history predict low πB with contributions from biogeographic origin and its interaction with elevation. -83.28 28.30 Mean elevation * Present landscape resistance Higher elevation species may experience stronger environmental barriers across lowlands. -82.02 29.55 Mean elevation Elevation alone may explain genetic divergence, although elevation may be a proxy for environmental factors. -81.89 29.69 Mean elevation * Elevation range Narrow elevation range predicts larger πB at high elevations, but not at low elevations. -77.47 34.10 A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. Figure Legends Figure 1. Map of central and eastern Panama showing the three sites where all frog specimens were obtained: G. D. Omar Torrijos H. National Park (in black), Chagres National Park (horizontal lines) and Cana field station in the Darién National Park (vertical lines), referred to as El Copé, Brewster and Cana, respectively. White points within each site indicate exact localities where tissue samples were collected for genetic analyses. Figure 2. Average between-population pairwise genetic distance (πB) computed for each of 31 species of frog based on the concatenated mitochondrial gene sequences, COI and 16S. Data points consist of 30 intraspecific comparisons between El Copé and Brewster (squares) and 15 intraspecific comparisons between Brewster vs. Cana (triangles). Complete generic names and numeric πB estimates are given in Table S3. Figure 3. Landscape resistance in log(ohms) for each of 30 frog species as measured between El Copé and Brewster (black bars), and between Cana and Brewster (white bars), based on present climatic conditions. Resistance was obtained from environmental niche models (see Methods) using the software CircuitScape 3.5.4. Resistance could not be estimated for Pristimantis museosus as only four locality records were available for this species. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved. A cc ep te d A rt ic le This article is protected by copyright. All rights reserved.