Testing effects of Pleistocene climate change on the altitudinal and horizontal distributions of frogs from the Colombian Andes: a species distribution modeling approach

. Recent climatic models suggest the late Pleistocene was colder and had different precipitation regimes from the present. If this climatic shift occurred more rapidly than species could adapt, species likely shifted their ranges as populations moved in concert with suitable environmental conditions. We examined changes in altitudinal and horizontal distribution in response to past climate change of amphibian species from different elevational zones and habitat requirements in the Eastern Cordillera of the Colombian Andes. We used environmental information and species occurrence data to model the distribution of 14 amphibian species (seven highland and seven lowland) which we projected to the Last Glacial Maximum (LGM) using two past climatic reconstructions. For these 14 species, we studied the predicted elevational and areal shifts. In agreement with palynological-derived models for Andean flora, we predicted that the elevation of montane amphibians shifted downwards increasing their total altitudinal range. We did not detect any evidence of compression related to drier lowlands. In some cases, the wider distribution areas of high-elevation amphibians during the LGM overlapped with contemporary distributions implying that these areas are present-day refugia for some species. Lowland species showed little or no elevational changes across time, but their areal distributions changed depending on habitat requirements. Four lowland frog species occurring in present-day xeric environments showed substantial range expansion across the lowlands during the LGM, while two species occurring in humid habitats likely expanded their ranges since the LGM. Since the LGM ended, ranges of mid- to high-elevation species shrank and shifted to higher elevation. Lowland species in xeric or open habitats have also experienced shrinking ranges, with some evidence that they have been moving upwards. Thus, low- and high-elevation species may be at risk under predicted anthropogenic climate change. Our results generate spatial hypotheses about amphibian responses to climate change that can be tested with phylogeographic data.


Introduction
Climatic fluctuations can have dramatic impacts on species' distribution and survival. Populations must either adapt to a changing environment or modify their distributions to track their preferred environmental conditions (Gavin et al. 2014).
During the Late Pleistocene Last Glacial Maximum (LGM; ~21 000 years ago) global temperatures dropped dramatically, often accompanied by changes in precipitation regimes (Haffer 1969). The LGM had major impacts on many species that shifted, Frontiers of Biogeography 2019, 11.1, e37055 reduced, or expanded their distributions relative to current ranges (e.g., Hewitt 2000).
In the face of global climate fluctuations, some areas may have remained stable such that populations may have persisted without the need to move or adapt to new conditions. These temporally stable and geographically isolated pockets of suitable habitat are known as 'refugia' (Haffer 1969, Gavin et al. 2014). The Pleistocene refugia hypothesis states that reduced precipitation during glacial maxima reduced the extent of lowland wet-forest environments, thus restricting wet-forest species to refugia isolated by a matrix of open, drier habitat (Haffer 1969). During interglacial times, such as the present day, these wet-forest species are expected to have larger ranges relative to their distribution during glacial periods (Waltari et al. 2007). On the other hand, the LGM may have represented a time of geographic expansion for species associated with open or xeric habitats, such as seasonally dry tropical forests (SDTF), with the present-day distribution of xeric-habitat species representing refugia (Pennington et al. 2000, Gür 2013). While some evidence exists for expansion of xeric-habitat species during the LGM (Quijada-Mascarenas et al. 2007), palynological and molecular evidence in support of the predictions of Pleistocene refugia for wet-forest species is equivocal (Colinvaux 1996, Lessa et al. 2003. A global drop in temperature of 6-8 ˚C degrees is estimated for the LGM (Caballero et al. 2010). The highlands likely experienced significantly more cooling than the lowlands, generating even steeper gradients in temperatures across elevations (Van der Hammen 1974). Such changes might have had differential consequences for species in different latitudes and at different points along altitudinal gradients. This drop in temperature likely forced high-elevation species downward, as corroborated by several examples in temperate-zone species. A study based on the fossil record of the Great Basin of North America showed that the altitudinal extent of some montane mammals shifted downward during the LGM compared to the present (Waltari and Guralnick 2009). Similarly, a study of high altitude populations of the American Pika showed genetic connectivity during the LGM, suggesting they shifted to lower elevations during colder conditions (Galbreath et al. 2009).
In Neotropical regions, however, temperature drops were potentially more heterogeneous and accompanied by variable changes in precipitation (Ramírez-Barahona and Eguiarte 2013). The tropical Andes experienced glaciation during the LGM that ranged from limited to extensive across South America (Smith et al. 2005). During this period, vegetation belts at high-to mid-elevations descended by 1000 m or more Van der Hammen 1993, 2004). In response, the associated fauna and flora presumably migrated downslope, although humid micro-refugia likely persisted at high altitudes in northern South America even while the coastal lowlands were very dry (González et al. 2008). Thus, present-day montane populations may have been 'compressed' between declining temperatures at higher elevation and increasing aridity in lower elevations (Van der Hammen 1974, Ramírez-Barahona andEguiarte 2013).
Amphibians are highly sensitive to extreme environmental changes in temperature and humidity (Navas and Otani 2007) and appear to have high niche conservatism with climate limiting dispersal for many tropical species (Wiens et al. 2006). Furthermore, previous studies suggest some Neotropical lowland species persisted in Pleistocene refugia (Carnaval et al. 2009). Amphibians, therefore, offer an excellent model system for evaluating the differential effects of late Quaternary climatic changes on the distributions of species across an altitudinal gradient. We studied 14 amphibians in the Eastern Andes of Colombia, an area with high amphibian diversity (Bernal andLynch 2008, Guarnizo et al. 2015). We selected species ranging from high and mid-elevations (seven species total) to lowlands (seven species with a range of habitat preferences, from wet to xeric environments). We constructed species distribution models trained on species' occurrence data and projected these models to the LGM to estimate the distributional and altitudinal changes experienced by these species, assuming niche conservatism (Wiens et al. 2009). We hypothesize that there was a downward movement in mean elevation for montane species during glacial periods relative to the current interglacial period that would have been restricted to middle elevations by the aridification of the lowlands (Hooghiemstra et al. 2006). We further hypothesize that in the lowlands, wet-forest species would experience range contractions during the LGM, while dry forest species would undergo range expansions, potentially expanding their elevational range upward during the LGM to find wetter conditions. To test these hypotheses, we evaluated species' predicted elevational shifts and changes in total geographic area between the LGM and the present.

Specimen and environmental data collection
To predict how frog species with contrasting elevational ranges and habitat preferences may have differentially responded to Pleistocene climate changes, we focused on 14 species from the Eastern Cordillera of the Colombian Andes that had at least 7 locality records to proceed with the generation of distribution models (Hernandez et al. 2006, Pearson et al. 2007). We identified seven highland and mid-elevation species (hereafter referred to as highland or montane) in three families whose altitudinal distributions reach above 2000 m but do not reach 0 m (mean elevation above 1000 m, Table 1), plus seven lowland species in four families whose lower elevation limits reach sea level and with mean elevations below 1000 m (Table 1). Among the seven montane species, we included two members of the treefrog family, Hylidae (Dendropsophus molitor and Hyloscirtus bogotensis), three in the direct-developing frog family, Craugastoridae (Pristimantis bogotensis, P. nervicus, and P. savagei), one species from the family of cryptic forest frogs, Aromobatidae (Rheobates palmatus), and one glass

Present-MIROC
High to mid-elevation Five out of seven lowland species are representatives of two of the previously mentioned families, Hylidae (D. microcephalus, Boana xerophylla, Smilisca phaoeta, and Scinax ruber), and Centrolenidae (Rulyrana flavopunctata). The last two lowland species belong to the diverse family, Leptodactylidae (Leptodactylus fuscus and Engystomops pustulosus). Several of our 14 species have been subject to recent taxonomic changes and analyses of cryptic diversity. We recognize the recent separation of the northern South American Boana xerophylla from its southern Brazilian counterpart, B. crepitans (Orrico et al. 2017). On the other hand, it has been suggested that the name Dendropsophus molitor contains two species, D. molitor and D. luddeckei (Guarnizo et al. 2012), and Pristimantis bogotensis reportedly contains two species, P. bogotensis and P. lynchi (Duellman and Simmons 1977). However, our own unpublished data suggest these are genetically almost indistinguishable and live in very similar high-elevation environments, so we decided to model these nominal species pairs together. Locality records for each species were obtained from electronic databases, museum records, and the literature. We used three online databases: Global Biodiversity Information Facility (http://www.gbif.org), HerpNET (http://www.herpnet.org), and the online database of the Instituto de Ciencias Naturales de la Universidad Nacional de Colombia (ICN, www.biovirtual.unal.edu.co). Furthermore, Colombian museums that provided additional locality data for this study include the Museo de Historia Natural ANDES at the Universidad de Los Andes, Bogotá, the ICN (additional data not available online) and the Museo de Historia Natural UIS at the Universidad Industrial de Santander. We filtered our database to remove duplicated records. Localities that did not have GPS coordinates were georeferenced using gazetteers and Google Earth (ver. 6.2.1.6014). We acknowledge there is a potential for taxonomic errors in databases, and we eliminated points that were obviously outside the known distributions for these species including, for example, points in the wrong cordillera or obviously in the wrong elevation. However, we did not use automatic elimination based on IUCN distributions since they are known to underestimate the distributions of species in the area. To reduce bias due to oversampling, for each species we thinned our presence datasets using a 5 km distance with the 'spThin' package for R using 5 replicates (Aiello-Lammens et al. 2015).
Environmental data consisted of 19 bioclimatic variables from the WorldClim database at 30" resolution (~1 km) (Hijmans et al. 2005). Because of substantial uncertainty in estimating environmental conditions during the LGM, we used reconstructions based on two contrasting, independent general circulation models (GCM): the Community Climate System Model (CCSM) and the Model for Interdisciplinary Research on Climate (MIROC). The basic components of the two models are atmosphere, ocean, sea ice, and land surface, while MIROC also accounts for rivers (Hasumi andEmori 2004, Collins et al. 2006). In the tropics, different GCMs differ mainly in their precipitation predictions, whereas in temperate zones they differ in their temperature estimates (Varela et al. 2015). To be conservative, we preferred to use both GCMs because there is not enough evidence to select one GCM over the other.

Species distribution modeling
We generated species distribution models by associating georeferenced species presence data with environmental data using Maxent (ver. 3.3.3k; Phillips and Dudík 2008). Since we only had presence-data available we preferred to use a presence-background method, Maxent, which unlike other methods makes no use of the controversial pseudo-absences but rather samples points randomly from a given area (which could include the presence points). Furthermore, it can produce accurate models even with small numbers of presence-only data; we used as a minimum 7 unique, spatially separated, filtered occurrences to create species' distribution models (Table S1, Hernandez et al. 2006, Pearson et al. 2007). Models were trained using 10,000 background points within a minimum convex polygon of the presence localities, buffered by 5 degrees. For species with sample sizes of > 26 occurrences, we used the random k-fold cross-validation for model training with k of 5 (Radosavljevic and Anderson 2014). That is, we divided each dataset into 5 groups of equal sizes and constructed the models using k-1 groups, withholding one for evaluation; the process is repeated until all folds are withheld once. For species with < 26 localities (Pearson et al. 2007), we trained the models using the jackknife method, a special case of the random k-fold where k is the total number of observations (Shcheglovitova and Anderson 2013). To avoid over-fitted models we evaluated several combinations of settings for model construction using the ¨ENMeval¨ R package (Muscarella et al. 2014). Specifically, we modified the transformation of the predictor variables (feature classes) and the penalty for the inclusion of parameters (regularization multiplier). The more feature classes were used and the smaller the regularization parameter, the more complex and overfitted was the resulting model (Elith et al. 2011, Muscarella et al. 2014. We used the following combinations of feature classes: 'L', 'LQ', 'H', 'LQH', 'LQHP', 'LQHPT', in combination with regularization multipliers from 1 to 4, in 0.5 increments (Muscarella et al. 2014). We used the corrected Akaike Information Criterion (AICc) to choose the most appropriate model for each species.
To obtain an estimate of each species' potential distribution in the present, all models were projected to the same area in northwestern South America (Fig. 1). To estimate the past potential distribution of each species at the LGM, we assumed two GCM's widely used in the literature that have the most extreme differences in projected climates. CCSM and MIROC have been shown to differ widely in both their modeled climates in the past (at the individual variable level) and in the resulting SDM projections (Guevara et al. 2018a(Guevara et al. , 2018b. In this study, we intended to investigate the potential distributional shifts of amphibian species in the past. Thus, we wanted to include contrasting climate models to encompass as much variability as possible in the projections. These differences are particularly dramatic regarding precipitation estimates (Varela et al. 2015), which are the most contentious in these past projections in the Neotropics and are particularly relevant to our hypotheses regarding the responses of species to drier conditions in the lowlands. A recently published study on transferability shows that the Maxent algorithm is in fact liberal in its predictions, thus making it a good choice for transferring when data on the fundamental niche are incomplete (Qiao et al. 2018). Maxent was shown to be less conservative than other algorithms such as presence-only or presence-absence algorithms (e.g., ENFA and GLMs). In the present study, then, Maxent is more likely to be catching the extremes of the distributional shifts.
For all species, we compared the results from the two GCMs using Schoener's D (Rödder and Engler 2011), a measure of niche overlap well suited to compare predicted distributions derived from SDMs, as implemented in ENMeval. Schoener's D values range from 0 to 1, where 1 indicates complete overlap (identical distribution models) and 0 indicates no overlap. All models and past projections were converted from continuous to binary predictions using the species-specific 10 th percentile training presence threshold (Pearson et al. 2007). With this threshold, suitability of presence localities used for building each species model were sorted, and after excluding the lowest 10% (accounting for potential sink populations and errors), the lowest suitability value was used as our threshold (i.e., all cells with estimated suitability values equal or bigger that that are considered as presence).

Elevation
We used the resulting binary SDMs and projections to estimate current and past elevational ranges. To estimate change in elevation we first resampled the present model to match the coarser resolution of the GCM's. From within the full projection area of the models for each species and all three maps (present, CCSM, and MIROC), we extracted elevation values for pixels in suitable habitat areas (binary prediction of 1) from an altitude layer obtained from the WorldClim database (Hijmans et al. 2005). All elevation value extractions were done using the raster package for R (Hijmans and van Etten 2013).
We evaluated whether changes in elevational distribution occurred between the two time periods: the present versus LGM. We used a Kolmogorov-Smirnov test of normality of the distribution of elevation values among the random points in suitable habitat areas (see above) from the present and past distributions of each species. All variables rejected normality even when variously transformed (P-value < 0.05). Therefore, we used the non-parametric Kruskal-Wallis test and Tukey-like post-hoc tests as implemented in the 'pgirmess' (Giraudoux et al. 2018) package for R (R Core Team 2015) to examine the potential difference in elevation between the present and either of the two LGM projections for each species, and for the pooled groups (highland versus lowland species). For each test we used the time period as a predictor variable with three categories (present, CCSM, and MIROC) and responses were either elevation of each individual species (14 different tests), elevation of all highland species together (1 test) or elevation of all lowland species together (1 test).

Changes in distribution area and refugia
We calculated the area of potential distribution for each species in the present and for each past projection as the area in square kilometers by converting the model output to a binary prediction of suitable areas (using the 10 th percentile threshold, see above) in ArcMap 10.1. SDMs were projected onto the Magna Colombia Bogota projection and the total area of presence was calculated using the Zonal geometry tool in ArcMap 10.1. We then computed the changes in distribution area between the present SDM and each past projection ( Fig. S1 and S2). Finally, we tested for the significance of the effect of time period and elevational preference on the area of potential distribution for each GCM (CCSM and MIROC) separately. For each GCM used for projection, we used a two-way ANOVA with time period (present or LGM) and elevational preference (highland or lowland) as predictors and area of potential distribution of each species (14 total) as response variable.
Refugial areas are defined as areas of environmental stability where species are predicted to have persisted across time in the same geographic space. To infer refugial areas, we intersected past and present distributions of each species for each LGM projection independently. To evaluate whether individual areas of stability might be part of shared refugia for montane species in the present, or shared refugia for lowland species at the LGM, we also intersected refugia separately for all montane species and for all lowland species.

Species distributions
Occurrence data ranged from 7 to 271 records per species after thinning, with a median of 62.5 records (Table 1). Various combinations of parameters were chosen for different species and all models had AUC values > 0.52 for test points and > 0.74 for training points, with omission rates ranging from 0.008 to 0.25 for each species (Table S1). Species' distribution maps showed lowland species generally with wider distributions than montane species (Fig. S1 and Fig. S2, respectively), e.g., D. molitor versus L. fuscus (Fig. 1a, b). Comparisons of projected past distributions using two different GCMs (CCSM vs. MIROC) yielded Schoener's D values ranging from 0.60 (P. bogotensis) to 0.97 (Ru. flavopunctata) with a mean value of 0.80 (Table 1), indicating that results based on the two GCM's were quite comparable though not identical.

Changes in altitudinal distribution
When comparing the altitudinal distributions of species between the past [LGM at 21 thousand years ago (kyr)] and present, we observed that all 14 species showed small but significant changes in their predicted altitudinal distributions between the present relative to one or both LGM projections ( Table 1). Only one species, the lowland En. pustulosus, did not show differences between the present and the MIROC projection. We present the montane D. molitor and the lowland L. fuscus elevational profiles as two examples of our results (Fig. 2). When grouping elevational data by species' altitudinal distribution (7 highland species together and 7 lowland species together), we found significant change in median altitude by comparing the LGM versus the present in highland and lowland species, shifting downwards at the LGM [Kruskal-Wallis tests with P-value < 2.2 × 10 -16 for both highland and lowland species and post-hoc tests with all three comparisons significantly different (CCSM vs. MIROC, CCSM vs. present, and MIROC vs. present), Fig. 3]. Montane species overall showed a > 10-fold greater shift downward during the LGM (a drop of 877 m or 255 m in median elevation for CCSM and MIROC respectively) relative to lowland species (change of 21 m or 26 m in median elevation for CCSM and MIROC respectively, Fig. 3).
Six out of seven montane species showed a decrease in mean elevation at the LGM under both the CCSM and the MIROC models (Kruskal Wallis test with P-value < 2.2 × 10 -16 , and post hoc tests comparing present vs. CCSM and present vs. MIROC were also different; Fig. S3, Table 1, Table S2). The one exception was P. savagei, which showed significant differences and a median upwards movement during the LGM, irrespective of the GCM (Fig. S3f). When lower elevations were accessible, five of the seven lowland species shifted their ranges slightly but significantly downwards (median change of 34 m and 43 m for CCSM and . To obtain elevation values, we generated 10,000 random points in the entire study area and extracted elevation values from points in areas of predicted suitable habitat (binary prediction of 1 with a 10 th percentile training presence threshold). Elevation was obtained from a digital elevation model provided by the WorldClim database (Hijmans et al. 2005). Similar graphs for the other 12 species are in the supplemental material, with highland species shown in Fig. S3 and lowland species in Fig. S4. masl = meters above sea level. Figure 3. Elevational ranges predicted for species in each of the two elevational bands, high-to mid-elevation species (n=7) versus low elevation species (n=7), in two time periods, the present and the last glacial maximum (LGM). Predicted LGM climate conditions were obtained from two general circulation models: the Community Climate System Model (CCSM) and the Model for Interdisciplinary Research on Climate (MIROC). To obtain elevation values, we generated 10,000 random points in the whole study area, and elevation values were extracted from values for points in suitable habitat areas (binary prediction of 1) from an altitude layer obtained from the WorldClim database (Hijmans et al. 2005). For both elevational bands there were significant differences in elevation between time periods, irrespective of the LGM projection used (df=2, P-value < 2.2 × 10 -16, for both high and low elevation species with χ 2 = 9105.3 and χ 2 = 1614.8 respectively). masl = meters above sea level. MIROC, respectively) during the LGM, irrespective of the GCM used. This pattern appeared just for the CCSM projection for a sixth species, Sm. phaeota (Table S2). Contrary to our prediction that drier climates during the LGM would force some species upward towards greater precipitation, our results show upward shifts in elevation for only two species, D. microcephalus under both GCMs and Sm. phaeota under the MIROC projection (Table 1).

Changes in horizontal distribution
According to our SDM analyses, ten out of 14 species had more extensive geographic ranges at the LGM than in the present, irrespective of the GCM used ( Table 1). Of the four species showing a different pattern, two lowland species, Sm. phaeota and D. microcephalus, have larger distributional areas now than during the LGM. Two highland species, Rh. palmatus and P. savagei, are predicted to have a more extensive geographic area in the present according to the MIROC projection but not according to the CCSM projection (Table 1). Across all species, the average geographic range size of species in our two elevational categories (highland vs lowland) were significantly different (P-value = 5.81 × 10 -6 , Fig. 4) but mean range size was not significantly different comparing LGM to present day (P-value = 0.37, Fig. 4).
We computed areas of stability across time for each species under both projections to the LGM separately. Refugia are smaller distributions contained within larger distributions in either the present or the LGM, thus while refugial areas are typically thought of as smaller distributions for humid-habitat species during the LGM, they can also correspond to the present distribution of xeric habitat species. Intersecting areas of past and current distribution for each montane species yielded extensive, continuously distributed areas of stable habitat for individual species (Fig. S5). However, intersecting refugia for all montane species did not result in any common areas of stability (Fig. S7). We show refugia inferred for D. molitor as an example for highland species (Fig. 5a). For lowland species, we observed two patterns regarding the distribution and size of refugia among the seven species. Five species had widespread refugia while two species, S. phaeota and D. microcephalus, had smaller refugia (Fig. S6). . Difference in area obtained from binary predictions of suitable area derived from species distribution models (using a 10 th percentile training presence threshold) for present and two past projections to the last glacial maximum (CCSM and MIROC). a) Difference in area between the present and CCSM projection in square kilometers for montane species (n=7) and lowland species (n=7). b) Difference in area between present and MIROC projection in square kilometers of montane (n=7) and lowland species (n=7). Two-way ANOVA showed an effect of elevational category but no effect of time period or the interaction between elevation and time in predicting changes in area of distribution (P-value = 5.81 × 10 -6 , 0.367 and 0.671 respectively with F values of 28.196, 1.030 and 0.404, respectively, and 1, 2 and 2 degrees of freedom, respectively).
We show the areas of stability for L. fuscus (Fig. 5b) as an example of a widespread lowland refugia.

Discussion
Climatic cycles of the Quaternary had a profound effect on biodiversity (e.g., Bennett 2004). For individual species, such effects could range from extinction to adaptation or simply tracking suitable habitat by shifting ranges (Gavin et al. 2014). Species' tracking of preferred climatic conditions during these cycles could have led to shifting elevational ranges, dynamic geographic distributions, and potential areas of long-term persistence. Our results support the existence of climatically suitable areas in the LGM for all 14 species studied, but the size and location of such areas varied in relation to altitudinal preferences and different habitat requirements. We are currently in an interglacial period, and we used past data from the LGM to better understand how species distributions shifted in response to these two contrasting climates that repeated throughout the Pleistocene climatic cycling.
Palynological evidence shows that mid-elevation Andean forests shifted 700 m downslope in their upper limit but their lower limit remained constant; whereas upper montane forests shifted downwards around 1000 m in both their upper and lower limits (Hooghiemstra and Van der Hammen 2004). The Andean forests and the Sub-Andean vegetation belt in particular experienced altitudinal narrowing during the LGM as a result of drier lowlands preventing mid-elevation plant species from moving further downwards (Hooghiemstra et al. 2006). Our results suggest that montane species shifted their mean elevation downwards during the LGM by either 829 m or 352 m (median change was 877m or 255 m), for CCSM and MIROC models, respectively (Table 1), with the former estimate being consistent with the idea that fauna, such as frogs, followed the mid-elevation forest downwards. Contrary to our To obtain stability areas, species distribution models and both projections to the last glacial maximum were made into binary maps using the 10 th percentile training presence threshold (see Methods). Stable areas were those in which the species is predicted to be present during both the present and the past projection. expectation based on palynological models, however, we found that almost all highland species expanded their elevational ranges during the LGM even reaching near 0 m elevations, and thus were apparently not limited by drier lowlands (Fig S3). This observed mismatch in downward shifts between climate projections and palynological evidence raises a question of whether in this case these highland species might have been constrained by the absence of forest, rather than by drier climate itself, to not go lower.
We are unaware of other studies evaluating past altitudinal shifts in temperate or tropical amphibians, but similar downward shifts were also previously reported for temperate zone mammals of the Great Basin (Waltari and Guralnick 2009) and in the tropical cloud forests of Central America (Ramírez-Barahona and Eguiarte 2014). Besides downward elevational shifts, projections to the LGM for six out of seven montane species showed an expansion in distributional area during the LGM (for R. palmatus this was true only under the CCSM). Range expansion during the LGM has been previously reported for other montane species in both temperate and tropical zones (Guralnick 2007, Carnaval et al. 2014) and is consistent with palynological evidence of highland forests expanding. Our dataset includes three species of direct developers, all of which are in the high-elevation category. Among the 14 species studied, these should be more water-independent species and thus potentially less affected by variation in precipitation. However, our results here do not support any difference in response to climate change for direct-versus indirect-developing frogs (though this could be due to the low sample size of direct developers). For one species, the direct-developing P. savagei, our projections show opposite results: a range contraction and, depending on the GCM used, an upward elevational shift (with MIROC) or no shift (with CCSM) during the LGM. This species, whose mean elevation today is 1200 m (the lowest of all our high-elevation species), might indeed have a pattern consistent with palynological evidence of contraction of lower montane forests being prevented from going downwards and potentially pushed higher by the drier lowlands, thus contracting its range during the LGM.
During glacial periods, the coastal lowlands of northern South America were drier with more open habitat (González et al. 2008), leading to the expectation that the range of lowland frogs associated with moist environments would have been reduced and driven to higher elevation to track wetter climate, while lowland frogs associated with xeric habitats could have expanded their distributions. Contrary to our prediction, our models suggest five out of seven lowland species stayed in the lowlands and slightly shifted their lower limit downward. These results point to low temperatures being so limiting for these species that even species associated with wet areas were unable to go upwards towards more humid places. Two lowland species closely associated with either wet forest (Sm. phaeota) or wetter portions in more xeric areas (D. microcephalus) shifted upwards as expected given their water requirements. Lowland species showed contrasting patterns regarding the temporal change in their geographical extent and related to their different ecological requirements. The two species associated with moist environments had refugia during the LGM (Figs. S2c and 2Sg,respectively). This is consistent with previous results from wetter lowlands in South America, suggesting that species associated with moister environments persisted in refugia of wet forest (Haffer 1969, Carnaval et al. 2009). On the other hand, species not closely associated with wet environments, either associated with xeric habitats or generalists, (i.e., L. fuscus, En. pustulosus, Sc. ruber, and B. xerophylla) expanded their ranges during the LGM (Fig. S1, Table 1). The present distribution of these species is within their inferred distributions during the LGM; thus, their current distribution may be regarded as a present-day dry-forest refuge. However, one lowland species associated with wet forests, the glass frog Ru. flavopunctata, apparently expanded its distributional area during the LGM rather than reducing it. Such a pattern could be potentially due to this species' preference for slightly colder temperatures or could be an artifact of the fewer occurrence points available for modelling, reflecting the currently restricted distribution of this species. For lowland species, it appears that differences among species with respect to their association with water determined their distributional response to Late Quaternary climate change, with more xeric-associated species benefiting more from the drier conditions of the LGM.
Our results generate hypotheses regarding changes in species' horizontal and altitudinal ranges across time. The two GCM's used (CCSM and MIROC) have different assumptions and contrasting estimates of precipitation in the tropics during the LGM (Varela et al. 2015). However, inferred historical distributions were robust to the GCM used, as the high niche overlap measured in geographic space indicates (mean Schoener's D = 0.80; Rödder and Engler 2011). This robustness to GCM with different precipitation estimates, even when modelling species with distributions closely linked to water availability, gives confidence to our results and to the phylogeographic hypotheses generated from them. Although partly concordant with hypotheses generated from palynological data, our results suggest species moved beyond the predicted ranges of wet forests during the LGM. Our results prompt hypotheses and explicit predictions about historical demography that can be tested using genetic data. For example, the hypothesized synchronicity of demographic responses among species associated with similar elevations and environments can be tested (Chan et al. 2014). Under that hypothesis, we predict that the demography of highland species as well as lowland xeric-habitat species should all show signals of expansion during glacial maxima, while lowland species associated with wet environments should show a corresponding demographic contraction during these colder and drier periods. Furthermore, we hypothesize highland species will show signatures of range contractions over the past 20 kyr except for Rh. palmatus that will show signatures of expansion. We also hypothesize that lowland populations associated with more xeric habitats will show population genetic signatures of historical demographic contraction since the LGM, whereas species associated with wet forests will show signatures of post-LGM expansions (concordant with range expansions seen in the species distribution models). Finally, we hypothesize that populations of lowland species associated with wet forest and living in areas of LGM refugia will have higher genetic diversity than populations living outside those areas (Fig. S5, Fig. S6).
Understanding how species respond and adapt to changing environmental conditions is of tremendous relevance for predicting future biotic responses and informing plans to mitigate the further loss of biodiversity (Jezkova et al. 2011). Current climate change is dramatic, with annual mean temperatures exceeding records set during the past millennia, and it is clearly a result of human activity (IPCC 2014). Climate projections agree on a continued warming during the 21 st Century that might well reach 3ºC or more in the northern Andes, along with complex changes in precipitation that will likely include decreased rainfall for the coastal Caribbean of Colombia coupled with increased rainfall elsewhere in the country (Buytaert et al. 2011). However, responses to environmental variation are species-specific . Therefore, species' unique natural history and habitat requirements must be considered when predicting the impacts of climate change. The impact of climate change on amphibians is far from understood (Winter et al. 2016), but some studies suggest that under current climate warming some Andean frogs have recently increased their upper elevation limit, along with their pathogens (Bustamante et al. 2005, Seimon et al. 2007). Our results show that there is a differential response to climate change according to species altitudinal ranges and habitat associations (i.e., xeric vs. humid environments). The lowland biotic attrition hypothesis states that lowland systems will lose species unable to tolerate warming (Colwell et al. 2008). In accordance with this hypothesis, we found that high-elevation species responded to the global warming trend from the LGM to the present by shifting their elevations upwards (leaving the lowlands), whereas lowland species varied their horizontal distributions, which could lead to local extinction if no suitable area remains. Potential future impacts should be evaluated differently for highland and lowland species and for wet-versus xeric environment species because changes in temperature might have bigger impacts for high-elevation species while precipitation changes impact lowland species further.
Santander, Bucaramanga, Colombia. We would like to thank the @CrawLab at the Universidad de los Andes, especially Lucas Barrientos and Carolina Ortiz, for their help in the data collection and analyses. Many thanks to the Anderson Lab at the City College of New York for valuable insights into this manuscript.

Supplementary Materials
The following materials are available as part of the online article from https://escholarship.org/uc/fb Table S1. Summary data for selected species distribution models for each species. Table S2. Results of the Tukey-like post-hoc tests to examine the potential difference in elevation between the present and either of the two LGM projections for each species. Table S3. Original presence points used for analyses for all 14 studied species.
Figures S1-S2. Present and past distributions of high and mid elevation frog species. Figures S3-S4. Histograms of the frequency of each elevation in the species distribution models and both past projections to the last glacial maximum, the Community Climate System Model (CCSM) and the Model for Interdisciplinary Research on Climate (MIROC). Figure S5. Maps representing areas of stability for montane species. Figure S6. Maps representing areas of stability for lowland species. Figure S7. Maps representing the combined refugia for all highland and lowland species using the CCSM and MIROC projections. Figure S8. Maxent response curves for all 14 species modelled in this study.