Modeled Shifts in Polylepis Species Ranges in the Andes from the Last Glacial Maximum to the Present

Polylepis woodlands, the dominant high-elevation woodland species of the Andes of South America, are an increasingly important focus of conservation and restoration efforts as a buffer to the regional effects of climate change. However, the natural extent of these woodlands before the arrival of human populations is still debated. One significant approach to this question is an assessment of the change in woodland extent from a hypothetical peak at the time of the Last Glacial Maximum (LGM) to the present where distributions have been altered by both Holocene climate oscillations and anthropogenic pressures of pre-Colombian and modern communities. LGM and present distributions for 21 Polylepis species were modeled using Maxent with environmental data obtained from the WorldClim database. Overall, potential woodland extent is 36% smaller today than at LGM, however a few species have experienced a projected increase in potential range of 180%. This has occurred at the interface of the southern Amazonian Basin with the Altiplano where Polylepis species richness is highest. Bioclimatically stable areas for each species averaged 20 ± 4% of the modeled range and mostly occurred in disjunct pockets from central Peru to northern Argentina and Chile.


Introduction
The high-elevation Polylepis (Rosaceae) woodlands of the Andean Cordillera of South America are an ecologically significant but highly vulnerable ecosystem.Known to support many endemic species, these species form an integral part of cloud forests and treeline landscapes at high-elevation tropical and subtropical Andean ecosystems [1][2][3].Their presence in often otherwise treeless alpine landscapes make Polylepis woodlands a driver of biodiversity by providing sheltered, richly-structured woodlands for a variety of animal and plant species [3,4].A number of authors have hypothesized that Polylepis woodlands have played a significant influence as a refuge for relict plant and animal species during climate change [3,5].Moreover, Polylepis woodlands have played an important role in the life of indigenous human communities for thousands of years [6].
There has long been a major controversy as to the original natural extent of Polylepis woodlands and whether the current highly fragmented distribution of these stands are the result of widespread pressure from thousands of years of human populations in the Andes or instead abiotic limiting factors [2].It has long been hypothesized that these woodlands are not climatically controlled in distribution and that they once covered major areas of the Andean highlands [7].An alternative hypothesis is that the disjunct ranges of Polylepis woodlands represent favorable microsites for survival in the harsh growing conditions of the high Andes [8,9].There is weak evidence for species-specific leaf traits related to climatic niche optima as evidence of ecological specialization but these data have not resolved the question [10].
Understanding the response of mountain ecosystems to climatic change is essential for both conservation efforts and to mitigate any adverse effects to local human communities.One avenue to better understand these changes is through the investigation of previous global climate change and the response of near modern ecosystems [11].Various future emission scenarios are similar in magnitude to climate changes that have occurred within the Quaternary period (the past 2.59 million years) and have had profound impacts on modern-day biota distribution, evolution and extinction [12].
A notably important global event was the last glacial maximum (LGM); the most recent interval in Earth's history when global ice sheets and mountain glaciers reached their fullest extent about 21,000 years ago (21 ka).Changes in climate were in response to forcing influenced from decreases in Northern Hemisphere summer insolation, lowered tropical Pacific sea surface temperatures of as much as 4 • C, and lowered atmospheric CO 2 [13].In addition, the local last glacial maximum (LLGM), while generally considered to be 21 ka, may have occurred in portions of the tropical Andes as early as 34 ka [14].Whatever the date, these glacial conditions contributed to substantial ecological change, including changes in plant community composition and lowering of the Andean treeline [15,16].
Species distribution modeling (SDM) with current and paleoclimate regimes can be an invaluable tool to predict changes in distribution in response to climate change and identify potential stable (refugial) areas [17].Also known as hindcasting, SDMs that use paleoclimate envelopes are increasingly being used to address questions regarding the past distribution of different taxa, ecological niche conservation, and changes within biodiversity hotspots [18].Hindcasting SDMs applied to relict woodlands within the Mediterranean Basin [19] and the spatially restricted hotspot of the Brazilian Atlantic rain forest [17] predicted variation and persistence in species distribution, and patterns of biodiversity following LGM.This approach can add a spatial perspective on the contraction, expansion, and degrees of continuity for spatially restricted forest habitats.
Hindcasting the currently fragmented Polylepis woodlands of the tropical and subtropical Andes addresses several important questions regarding past response of high Andean vegetation to climate change.Although Quaternary records indicate range fluctuations along an elevational gradient [20][21][22], the absence of species level paleoecological data limits our understanding of Polylepis species response to changes in glacial extent and climate.This important Andean ecosystem with its high number of endemic species is a current keystone habitat for biodiversity in an often treeless alpine landscape and may have had a similar function in the past.Hindcasting presents an approach to predict refugial areas and the extent of contraction and expansion of Polylepis species woodlands at the LGM, and thereby stable areas of species endemism and high species richness in the high-elevation Andes.
Our study has three primary objectives.The first objective was to use current distribution of Polylepis species to predict distribution during the LGM, assuming that species climatic niches have not been changed.We use these data to estimate the change of each species distributional area from current patterns of geographic occurrence along the length of the Andes.Our second objective was to use the predicted fluctuations in distribution of each species to identify changes in potentially suitable areas for several Polylepis species and specifically within biogeographic zones identified using major Andean hydrological basins (see methods).Finally, we analyzed the change in Polylepis species distribution to identify refugial areas that have remained spatially stable since LGM.This approach is designed to improve the understanding and conservation of this significant component of a biodiversity hotspot that has had major spatial responses to climate change.

Study Species
The approximately 27 species of the genus Polylepis (Rosaceae) are woody shrubs or trees growing to 10 m in height and characterized by multi-layered red exfoliating bark, small pinnately-lobed leaves [1].Species grow in a variety of high elevation habitats from cloud forests to tree line woodlands to isolated woodlands above treeline in the tropical and subtropical Andes.They range in distribution from western Venezuela to northern Argentina and Chile, and from the cold and wet eastern flank of the Andes to high elevations along the much drier western slopes.Known as queñua by Andean peoples, Polylepis populations often form highly fragmented evergreen woodlands at elevations above 3000 m among high-elevation grassland ecosystems known as páramo from Venezuela to Ecuador, puna from central Peru to Argentina, and jalca in the transition zones between these two ecosystems in northern Peru.
Present georeferenced point localities for 21 species of Polylepis were obtained from an extensive review of online herbarium databases (Association of Andean Ecosystems-ECOAN, TROPICOS, GBIF), published research conducted from April 2006 to January 2008, and by field surveys by BRZ in March 2003 (Table S1).The number of georeferenced points obtained ranged from a pre-set minimum of 20 to more than 100 per species and screened to include georeferenced points obtained from recent surveys (within the last 15 years).To help verify geographic accuracy, points were overlaid on images from Google Earth (version 4.2; Google, California, CA, USA), compared to descriptions from the source (e.g., country, department, elevation, proximity to landmark, etc.), and visibility in moderate and high resolution images which allowed verification of locations.Erroneous and duplicate presence records were removed.

Climate Data
Present climate data was obtained from 19 bioclimatic metrics from Worldclim (version 1.4) [23].Often used in SDMs, these metrics are derived from monthly temperature and precipitation records over a 50-year period (1950-2000) and gridded to 5-km spatial resolution.The metrics include 11 temperature and eight precipitation metrics, which represent annual means, seasonality, and extreme or limiting environmental factors.After data reduction to minimize covariance, a total of 11 bioclimatic layers were used for this analysis.Covariance between the bioclimatic metrics was estimated using cross-correlation matrixes based on Polylepis point localities and 1000 random points drawn from South America.This reduced set of metrics also correspond to environmental factors affecting Polylepis distribution in published research [3,10].We included annual mean temperature, mean diurnal temperature range, temperature seasonality, maximum temperature of warmest month, minimum temperature of coldest month, annual precipitation, precipitation seasonality, precipitation of wettest quarter, precipitation of driest quarter, precipitation of warmest quarter, and precipitation of coldest quarter.Minimum temperatures and minimum amounts of precipitation are generally expected to limit woodland distribution in elevation and latitude based on ecophysiological limitations.Maximum temperature is known to restrict woodland growth in high elevations.Distinct bioclimatic conditions in the Andes and annual rainfall tolerance are expected to limit Polylepis distribution in particular [3].
LGM estimates of these climate layers were also obtained from Worldclim and are based on the National Center for Atmospheric Research (NCAR) Community Climate System Model (CCSM3) for the Last Glacial Maximum (LGM) following the protocols established by the Paleoclimate Modeling Intercomparison Project (PMIP2) [24,25].CCSM3 is a coupled atmosphere-ocean-sea ice-land surface General Circulation Model (GCM) with a horizontal atmosphere and land grid spacing of approximately 311 km [24].Modeled forcing for the LGM resulted from large changes in greenhouse gases, aerosols, ice sheets, sea level, land cover and plant functional types.Past climate change projections (Figure 1) were statistically downscaled to a 5-km spatial resolution and made available through Wordclim.

Species Distribution Modeling
The SDM algorithm, Maxent (version 3.2.1),was used for this study.Maxent uses the known presence of a species, along with predictor variables, to predict distribution over a study area [26].Maxent runs on presence-only point occurrences and has been shown to have a high predictive performance with relatively few point localities [27].These characteristic makes it especially suitable for predicting the distribution of Polylepis species since fragmented woodlands and the limited availability of accurate geographic data in several countries may reflect anthropogenic impacts rather than environmental limits [2].
We selected Maxent default settings for regularization and feature classes for all model runs, which includes linear, quadratic, threshold and hinge features [28].The default setting of regularization allows Maxent to select an amount of regularization that is appropriate for the number of point localities and types of features used.The set of features represent environmental factors and interactions between predictor variables that constrain the geographic distribution of the species being modeled [26].We used the default prevalence value of 0.5, which refers to the proportion of occupied locations by the species being half of all possible locations in the study area.Although this may not be appropriate for rare species, the literature indicates that the geographic distribution of the 21 Polylepis species in this study are not currently exceptionally restricted [3,4,6].Maxent models (from version 3.2.1)produce a logistic output format in the form of an image, with pixel values ranging from 0 to 1, which represent an estimate of probability of presence [28].Pixels with a value nearer to 1 are sites predicted to be the most suitable for the species of interest [28].In this study, pixels with values above or equal to 0.5 are identified as areas with a greater probability of presence, than lower thresholds, and are expected to contain typical bioclimatic conditions of ecological niche of a species.For analysis of spatial accuracy of each species' model prediction, test localities were generated by making 10 partitions of each species' point localities and randomly selecting 70% of points for model training and 30% for model testing [26].
Maxent first predicted Polylepis species distribution based on present climate, represented by the subset of 11 bioclimatic layers (Figure S1).The relative contribution (%) of each climate layer to the final predictive distribution model for each species was recorded (Table S2).Then, the model was projected on the CCSM3 past climate layers, which reflect estimates of LGM temperature and precipitation.The difference between past and present predicted species distribution indicated areas of decreasing or increasing probability of occurrence, evident with negative or positive pixel values, respectively.Areas of present and past Polylepis species richness were identified by summing predicted species distributions as performed in Gap analysis [29] and biodiversity studies [30,31].Only pixels with value at or above a 0.5 logistic output were considered.The resulting analysis provided the area of species distribution at LGM, the present area, and the area of stable distribution from LGM to the present.Also predicted was the elevational range for each species at LGM and present.The difference between past and present Polylepis species richness indicated areas of decreasing or increasing richness as a result of climate change.Historically stable areas (refugia) were defined as those pixels for which presence was predicted in each species model and two climate scenarios.
The spatial accuracy of model predictions of individual species was evaluated using threshold-dependent and threshold-independent tests.We used the extrinsic omission rate, or the fraction of test localities that fall outside the predicted area, and the proportional area, or fraction of all pixels predicted be suitable by the model, both at a fixed 10% cumulative probability threshold as two threshold-dependent tests.A one-tailed binomial test was used to determine if the models had a lower omission rate than a random prediction.
The two threshold-independent tests used to evaluate overall model performance were the area under the receiver operating characteristic (ROC) curve (AUC) of the test localities and the test gain.The ROC curve is obtained by plotting model sensitivity (omission rate) against 1-specificity (commission rate) and measures a model's ability to correctly predict presence and absence.Although Maxent uses presence-only data, "pseudo-absences" are generated by randomly selecting a set number of background pixels.For our study, 10,000 background pixels were randomly chosen as pseudo-absences for each model run.Therefore, the AUC statistic can be interpreted as the probability that a presence site (pixel) is ranked above a random background pixel [26].An AUC value for a specific model scenario will range from 0.5 (random) to 1.0 (perfect discrimination).A one-tailed Wilcoxon rank sum test was used to determine if the AUC of a prediction was significantly better than random.AUC data are presented as means ± S.E.(n = 10 models).The test gain is the average log probability of the presence sites used in the model [32] and has been used as an additional measure of overall model performance [33,34].The uniform distribution of a model has a gain of 0, whereas the model test gain value of x indicates that the average likelihood of a test presence site is exp(x) times greater than a random background pixel [32].Statistical analysis was carried out with SPSS 13.0 (SPSS Inc., Chicago, IL, USA).

Biogeographic Zones
In order to divide the Andean cordillera into biogeographic zones, major hydrological basins were selected from the South American HYDRO1k data [35].HYDRO1k is a digital elevation model (DEM) derived dataset that includes global and continental topographical data at 1-km spatial resolution.Major basins include San Juan, Orinoco, Amazon, Altiplano, Argentina, and Paraná (Figure 2).Several basins were further separated into subdivisions to analyze specific high-elevation biogeographic zones where Polylepis woodlands are often present (Figure 2b).
The San Juan Basin, forming the western slopes of the Andes with Pacific drainage, extends from northern Colombia and Venezuela to northern Chile.Within this region, high tropical montane vegetation above the continuous timberline from Venezuela to Ecuador is considered páramo and has a generally cold and humid climate compared to the analogous puna of the typically more xeric and lower humidity highlands of central Peru to northern Argentina [36].The transition zone, or jalca, between these two ecosystems in northern Peru is often delineated by the Huancabamba Depression to the north and the Cordillera Blanca and Cordillera Huayhuash to the south.Therefore for this study, the basin is divided by the Amotape-Huancabamba Zone in northern Peru [37] and the southern extent of the Cordillera Blanca and Cordillera Huayhuash in central Peru.The resulting biogeographic zones are identified as Northwest Andes, West Jalca, and Central Andes.Also, the Amazon Basin, which borders much of the San Juan Basin to the east, was subdivided into North Amazon, East Jalca, and South Amazon regions (Figure 2).

Polylepis Species Distribution
For all model runs, the extrinsic omission rates at the 10% cumulative threshold were statistically significant compared to a random prediction (p < 0.001, one-tailed binomial test).The low omission rates indicated that only a small percentage of test points fell outside of the model area predicted to be suitable.The AUC values for 21 Polylepis species distribution predictions were statistically significant (p < 0.001, one-tailed Wilcoxon rank sum test), indicating that all model scenarios and data partitions were significantly better than random.The average test AUC value, based on 10 random partitions, for predictive models was 0.995 ± 0.003 and the average training AUC value of 0.993 ± 0.002 which indicates that predictive models were successful in discriminating suitable from unsuitable habitats (Table S3).The model test gains across all species indicated that the average likelihood of a test presence site was at least 90 times greater than a random background pixel, which would also reflect the restrictive climatic envelope of Polylepis.
An overall contraction of Polylepis woodland distribution was predicted to have occurred since the LGM based on present climate and CCSM estimations of LGM climate.Overall Polylepis woodland distribution is 36% smaller today than at LGM, with about 249,000 km 2 at that time compared to only 158,800 km 2 today over its full range from western Venezuela to northern Argentina and Chile.Predicted suitable area for Polylepis species richness along the total length of the Andes reached as many as ten species in a 25 km 2 pixel (Figure 1).
The geographic distribution of most Polylepis species was predicted to have contracted or expanded in range by only a small amount since the LGM although in many case by a small amount (Table 1).However two species had major losses in their area of occurrence.P. australis and P. tarapacana are predicted to have contracted to nearly one fourth of their former extent, and the rare P. rugulosa shows a more than a 90% drop in area.These three species are in the south-central Andes near the southern limit of the genus.Mean elevational occurrence increased under present climate conditions by about 1450 m for P. australis and 348 m for P. tarapacana.Polylpeis weberbaueri from Ecuador and northern Peru is predicted to have lost more than half of its LGM area.At the other extreme, two species showed an increase in area of more than 50%.Polylepis besseri distribution in the Altiplano region increased from LGM nearly doubled in predicted area from about 17,180 to31,440 km 2 .Mean elevation of this species moved upslope by 615 m, with a minimum elevation of 1270 m at LGM to a minimum of 2220 m with present climate Polylepis pepei from southern Peru and Bolivia increased its range by 56% (Table 1).
For the fifteen species with contracting distribution, the mean (± S.E.) stable area of coverage since LGM was 4785 ± 1110 km 2 or 20 ± 4% (Table 1).Among the species with the largest contraction in range, P. australis, P. rugulosa, and P. tarapacana from the south-central Andes had among the smallest area predicted to have remained stable to present climate with 4%, 2%, and 0.2% from their maximum extent at LGM, respectively (Table 1).The mean stable area for species with predicted expansion was 3103 ± 289 km 2 or 20 ± 2% of their past range.Polylepis besseri, with a large expansion since the LGM, had a stable area of 18% of its predicted extent at LGM.

Biogeographic Zones
Polylepis species richness per 25 km 2 pixel was predicted to decrease within most biogeographic zones since the LGM, except for woodlands in the South Amazon (Table 2, Figure 3).In particular, the North Amazon Basin, which contains significant portions of the northern Andes, was predicted to decrease in the number of Polylepis species within 9960 km 2 , or about 62% of the total woodland area (Figure 3a).Similarly, a decrease in richness occurred in 67% of the woodlands predicted within the East Jalca Basin.Polylepis species richness was predicted to have remained relatively steady in the West Jalca Basin with areas that increased in richness or remain stable over time to be 93% and 87% of the 1975 km 2 extent at LGM, respectively (Figure 3b).The South Amazon Basin was the only biogeographical zone predicted to have greater woodlands area (20,775 km 2 ) with an increase in richness compared to areas that would decrease or remain stable in richness (Figure 3c).This region was predicted to contain the highest number of Polylepis species and pixel areas of increasing richness under both LGM and present climate conditiond.In contrast, the Argentina biogeographic zone was predicted to only contain P. australis at LGM and a mere 325 km 2 range would remain today.Each biogeographic zone was predicted to have changes in Polylepis species richness following changes in mean annual temperature and total annual precipitation (Table 2).As the mean elevation of species rich areas moved upslope in the North and South Amazon zones, areas that lost species had a greater difference between LGM and present annual mean temperature than other areas, with −1.5 and −2.9 • C, respectively.This trend also occurred on the western slopes of the Northwestern Andes.However, both the West and East Jalca zones had the opposite trend with predicted areas of increasing richness possessing a greater difference in mean annual temperature, at or above −2.3• C, than stable or declining richness areas.
The greatest fractional change in annual precipitation also occurred in predicted areas of increasing richness in the North Amazon, East Jalca, and South Amazon zones, relative to areas of decreasing richness.Species-rich areas today, at higher mean elevations in these zones, tended to be wetter than lower elevation areas that lost Polylepis species richness after the LGM.Only the North Andes and Orinoco zones, at the northern woodland extent in the Andes, and the Paraná, at the southern extent, had either predicted stable or increasing richness in areas currently drier than in the LGM.
Stable areas of species richness tended to have intermediate mean elevation, area, and changes in temperature and precipitation relative to areas with significant changes in species richness.Across the western tropical Andean slopes, in the North Andes, West Jalca and Central Andes biogeographic zones, the mean elevation of stable areas was 3759 m with a mean temperature (± S.E.) of −1.9 ± 0.2 • C. The eastern slope of the tropical Andes, which drains into the Amazon Basin and divided into North Amazon, East Jalca, and South Amazon zones, had similar stable areas with a lower mean elevation of 3688 m and a similar annual mean temperature of −2.0 ± 0.4 • C. Mean fractional change in precipitation indicates drier conditions from LGM to present for both slopes of the tropical Andes with −11.8 ± 6.7% on eastern slopes and −6.8 ± 2.1% on western slopes.The magnitude of climatic change in climatically stable relative to non-stable areas depended on the particular zone.In both the East and West Jalca zones, the mean temperature response was <1 • C and the fractional change in precipitation was no more than 5% in either species richness trend.The Orinoco Basin also had <1 • C mean temperature difference but as much as 24% fractional change in precipitation in stable areas relative to areas with decreasing richness.

Discussion
A growing number of paleoecological records have clarified the changing composition and distribution of vegetation belts along the latitudinal length and elevational extent of the Andes during glacial/interglacial episodes.Fossil pollen in the northern Andes suggests that the species-diverse lower vegetation belts moved upslope under warmer interglacial climatic conditions and species-poor high elevation vegetation moved downslope during colder glacial episodes [35].Portions of the Southern Hemisphere tropical Andes show little fluctuations in plant composition, but contain significant changes in the dominant plant type that would have created plant communities with no modern Andean analogue [22].Nonetheless, certain vegetation belts such as high elevation páramo grasslands in the northern Andes and puna in the central Andes, and woodlands dominated by Polylepis (Rosaceae) were consistently present throughout Quaternary records [20,36].Unfortunately, palynological analysis often lack spatial precision since fossil pollen collection sites are often far apart and rely on factors such as wind direction and lake size, resulting in regional rather than local proxies [38].The identification of Polylepis paleo-woodlands is difficult since fossil pollen is nearly identical to the closely related genus Acaena (Rosaceae) [2].Although pollen at elevations above 4800 m is more likely to be Polylepis rather than the mostly cloud forest restricted Acaena, paleo-woodlands are often identified as Polylepis/Acaena or Polylepis type [22,38].
The use of species distribution models with climate scenarios for the present and LGM predicted a broad but variable increase in mean elevation of occurrence and contraction of Polylepis woodland distribution since the LGM.Extensive paleoecological records in the central Andes indicate an upslope movement of taxa by more than 1000 m in response to warming since the last glacial period [15].Polylepis woodlands were predicted to reach their maximum extent during the LGM as colder conditions would have forced them to occupy lower elevations [21] that have greater topographical area.
At high elevations, moraine sites within the Cordillera Blanca of Peru contain remains of Polylepis trees killed by glacial advances up to 400 years BP [39].Paleoecological evidence suggests climate conditions, including increased wetness and insolation, in the Altiplano 26-14.9ka may have been beneficial to the expansion of Polylepis populations [21].Warmer and drier conditions during the Late Pleniglacial (21-14 ka) in the northern Andes widened vegetation belts and upper limits increased in elevation as glaciers retreated [20].Pollen record provide evidence of the development of Polylepis forests in the southern part of the so-called Andean Depression of southern Ecuador, reflecting warmer and drier climatic conditions during early and mid-Holocene [40].
The range expansion of several Polylepis species in the Altiplano region since LGM reflects the unusual topography of the central Andes.Most predictive models assume range contraction of species distribution with increased elevation following climate warming, since land area typically shrinks with increasing elevation [41].The topography of the northern and southern Andes follow this pattern, however the central Andes contain the Altiplano, an extensive high-elevation plateau that occurs above 3600 m in elevation in western Bolivia and southern Peru, extending into northern Chile and Argentina.The 250 km-wide plateau will potentially allow more land area following upslope range shift of species in response to climate change [42].Therefore, previous fluctuations of Polylepis woodlands extent and complex regional topography suggest the potential non-uniform response of species along the length of the Andes following climate change.
The lower limits of Polylepis woodland distribution during and after the LGM are rather unclear.Although fire is expected to be among the most relevant factors limiting the present formation of Polylepis woodland edges at lower elevations [43] the lower elevational limit tends to transition to other forest types [2].Paleoecological data suggests an approximate lower limit of 1400 m at the time of LGM for Polylepis as the dominant component of upper montane forests in the northern Andes [20]; however our study predicts niche presence at elevations below 1000 m in the Argentina, Altiplano, and Paraná basins at LGM.Nonetheless, these lower elevational niches were predicted to have been lost in the transition to present climate conditions and resulted in decrease in Polylepis species richness in these regions.
On a species level, the magnitude of predicted contraction of distribution range varied by species, with several Polylepis species expanding their range.Fourteen species were predicted to contract their ranges after the LGM, particularly the southernmost limits of P. australis and P. tarapacana and P. rugulosa in the western arid regions.The six species predicted to have expanded their range mostly occur along the southeastern Andes and Altiplano where the central Andean plateau widens to 400 km and an average elevation of 4000 m [44], compared to the narrower 50 km wide and 3000 m elevation sections in Venezuela [45].In addition, predicted changes in Polylepis species richness is supported by the Quaternary pollen diversity record of species-rich vegetation altitudinal shifts with past climate change [35].Our study also indicated that climatic changes correlated with shifts in species richness varied by biogeographic zone with marked difference between northern and southern latitudes and east and west slopes of the Andean cordillera.Although precipitation and cold night temperatures may limit Andean vegetation [20], the different predicted changes in mean elevation and species richness highlight the complexity of the Andean topography and climate regimes.
Along with dynamic changes in habitat, the predicted existence and extent of stable areas through Holocene climate fluctuations indicate that conserved areas of woodland continuity served to function as refugia to maintain species richness and narrow endemism.Major refugial areas were predicted to occur in the eastern Andean slopes of the southern Amazon Basin, Altiplano, and Paraná.The predicted lower limit of refugia, near 2000 m in elevation for most of the tropical Andes, coincides with paleoecological data indicating the lower limits of the Polylepis vegetation belt from the LGM to present in the northern Andes [20].However, the extent of the refugial belt was predicted to vary by biogeographic zone and far lower limits were predicted to occur in the Altiplano and southern ranges of the Andes.
Although Quaternary records indicate the widespread occurrence of Polylepis woodlands [21], other studies have proposed that Polylepis woodlands did not form permanent continuous woodlands before human arrival [9].The predicted stable areas in this study suggest the existence of conserved Polylepis species niches with climate change and the potential for continuous woodland associations in regions of high Polylepis species richness.The aggregated pattern of Polylepis-specialist birds in the Cordillera Blanca and Cusco region of Peru suggests isolation from other forest types during Pleistocene glacial periods [5].These stable areas also indicate key refugial zones for many associated endemic species following climate change and may help understand the mechanisms underlying local endemism and genetic diversity of flora and fauna [46].These regions appear to support multiple core species populations, which include a variety of Polylepis-associated bird and mammal species, including the pampas cat (Leopardus colocolo) [7,47].
The use of bioclimatic models and SDM algorithms comes with many significant advantages and inherent limits that should be addressed to properly put predictive models into context.Generally, climatic models based on historic climate station data are required to extrapolate in time and space due to the lack of continuous data in various locations causing significant gaps in information [23].The confidence in these temperature data are reasonable due to the relationship between elevation and temperature, but precipitation numbers are often questioned in mountainous regions due to the low density of climate stations [23].Therefore, caution must be taken on how distribution models are constructed and interpreted.Nonetheless bioclimatic models, as the model used in this study, can provide a useful approximation of species response to changes in climate using the central biogeographic premise that climate has a dominant control over the natural distribution of a species [48,49].Biotic interactions, a species' dispersal ability, genetic variation, and anthropogenic activity are major factors to species distribution that are not incorporated directly in SDM frameworks and have been repeatedly criticized due to this inherent limitation [50,51].
Indeed, much of Polylepis current distribution may reflect extensive anthropogenic pressures resulting in fragmented populations in different regions, reduced quantity of geographic point localities and an altered realized niche.The debate includes if the present fragmentation largely occurred before or after Spanish colonization of the region and if the current location of populations reflect land use rather than natural environmental extent.Regardless, the representativeness of the current distribution directly affects present and past modeling with the potential of under estimating distribution in both scenarios.We find a few models of our study to have disjunctive potential habitats, which would be difficult to support but are a result of these distant and isolated populations (Table 1).Nevertheless, the use of a representative sample size or a diverse set of point localities should provide enough information to achieve a statistically robust prediction a species' geographic range through spatial modeling.It is also notable to address issues in species identification, since isolated Polylepis species populations have been recorded far from any neighboring populations.As with many taxa, there are uncertainty in present species identification and taxonomic concepts and the Polylepis genus is no exception.There has been extensive hybridization in various landscapes, including from intentional introductions [52].Care must be taken in including the most reasonable point locations from digital taxonomic databases and those made available by local experts.
One final characteristic of SDMs worth discussing is the spatial resolution of the environmental layers and models.The spatial resolution of present and past climate variables at 25 km 2 does not capture all the variation that can occur within this pixel area.Andean topography is complex and fine-scale distribution patterns of Polylepis woodlands can be confined to slopes and habitats defined by microclimates and land use pressures, which would not be accurately represented within the pixel.This limits our approach to answer the larger question of the fate of these woodlands after civilizations established themselves in the Andes.Although the lack of detailed resolution is a common limitation of geospatial modeling, our objective was to identify the fluctuations in suitable areas across the region.
The large scale patterns across the biogeographic basins can tell us more on distribution patterns and limitations than small-scale patterns associated with local conditions.
SDMs have the potential to increase our understanding species response to past climate fluctuations and therefore understanding future distribution dynamics with anthropogenic climate change.Paleoecological data has provided key insights into the extent of Polylepis woodlands, but without species level identification as we have shown, the generalized response of Polylepis species distribution to temperature and precipitation will remain unclear.SDMs can provide a spatial perspective on species response and species richness that is directly applicable to current conservation efforts to preserve Polylepis woodlands.Temporally stable areas, as a refuge for many Polylepis associated species, will be key factor for conservation during anthropogenic climate forcing since their protection and management may reflect the persistence of Polylepis woodlands in the future.

Conclusions
The fragmented distribution of Polylepis woodlands present important questions regarding past response of high Andean vegetation to climate change.Although paleoeclimatic records from the Quaternary indicate that major fluctuations in elevational changes in climate have occurred, the absence of species level paleoecological data have limited our understanding of Polylepis species response to changes in glacial extent and climate.Hindcasting presents an approach to predict refugial areas and the extent of contraction and expansion of Polylepis species woodlands since the LGM, and thereby stable areas of species endemism and high species richness in the high-elevation Andes.The extent of Polylepis woodland today is 36% smaller today than at LGM, however, a few species have experienced a projected increase in range of 180%.Much of those changes have occurred at the interface of the southern Amazonian Basin with the Altiplano where Polylepis species richness is highest.Bioclimatically stable areas for each species from the LGM to the present averaged 20 ± 4% of the modeled range and mostly occurred in disjunct pockets from central Peru to northern Argentina and Chile.Polylepis woodlands with their high number of endemic species are keystone habitats for biodiversity in an often treeless alpine landscape and may have played a similar function in the past.

Supplementary Materials:
The following are available online at www.mdpi.com/1999-4907/8/7/232/s1, Figure S1: Potential present species distribution for each Polylepis species in this study, Table S1: Geographic location and source information for the presence data of 21 species of Polylepis, Table S2: Relative contributions (%) of environmental variables to Maxent models for Polylepis species, Table S3: Means of proportional predicted area (Area) and test omission rates (OR) at the 10% cumulative threshold, area under the receiver operating characteristic curve (AUC) and test gains for all models of Polylepis species.

Figure 1 .
Figure 1.Past (a), present (b), and overall change (c) of Polylepis species richness predicted with Community Climate System Model (CCSM) global circulation model.Overlaying the ≥50% probability of each species predicted distribution was used to estimate richness.

Figure 2 .
Figure 2. (a) Present elevation and (b) hydrological basins and subdivisions of the tropical and subtropical Andean cordillera.The San Juan Basin is divided into the subdivisions of North Andes (IA), West Jalca (IB), Central Andes (IC), and the Orinoco Basin (II).The Amazon Basin is divided into the subdivisions of North Amazon (IIIA), East Jalca (IIIB), and South Amazon (IIIC).The southern basins include Argentina (IV), Altiplano (V), and Paraná (VI).Present elevation greater or equal to 1000 m is indicated in gray scale.The annual mean temperature difference (c) and fractional change in precipitation (d) from present is given using CCSM simulation of LGM climate.

Figure 3 .
Figure 3. Regional change of Polylepis species richness for (a) North Andes (IA), Orinoco (II), and North Amazon (IIIA); (b) West and East Jalca (IB, IIIB); and (c) the southern zones of Central Andes (IC), South Amazon (IIIC), Altiplano (V), and Paraná (VI).Refugial areas are highlighted in blue.Present elevation greater or equal to 1000 m is indicated in gray scale.

Table 1 .
Area and middle elevation of Polylepis species distribution predicted with the CCSM global circulation model of past climate and Worldclim present climate model.

Table 2 .
Changes in predicted Polylepis species richness per 25 km 2 pixel since Last Glacial Maximum (LGM) within basins and biogeographic zones of tropical and subtropical Andes.
See methods for the calculation of mean temperature and precipitation changes from present at LGM. Predicted ranges of elevation are shown in parenthesis.