Distribution shifts in habitat suitability and hotspot refugia of Andean tree species from the last glacial maximum to the Anthropocene

ABSTRACT

Forecasting the effects of climate change on the distribution of Andean trees (Polylepis, Rosaceae) is important to understand how species respond to climate variability and to assess their resilience to the ongoing climate crisis. Here, paleodistribution modelling is used to assess distribution shifts of 17 Polylepis species during the Last Glacial Maximum (LGM), mid-Holocene (MH), and the Anthropocene in the central Andes. Species distribution models (SDMs) were computed by combining presence records and current climatic conditions using Maxent and projected onto three climatic scenarios for the LGM (~22,000 yr BP), the MH (~6,000 yr BP), and the Anthropocene (1,970–2,000). Subsequently, current refugia hotspots were identified by intersecting suitable habitat over the LGM, MH, and anthropogenic conditions for all the studied species. SDMs for the LGM and MH showed a contraction of climate suitable habitat for most of the species in comparison to the Anthropocene. Four current refugia hotspots were identified: central Cordillera of Peru, Lake Titicaca basin, western Cordillera of Bolivia, and northern Chile. In general, SDMs were consistent with patterns proposed with pollen records, and partially with available phylogeographic studies. Current hotspots are highly important areas for the conservation of Polylepis and associated biota. This study assists in understanding distribution shifts over millennia of Andean tree species in response to climate change and identifying key conservation areas for the delineation of future conservation strategies.

RESUMEN

La predicción de los efectos del cambio climático en la distribución de los árboles andinos es importante para comprender cómo responden las especies a la variabilidad climática y evaluar su resiliencia a la crisis climática en curso. Aquí, modelamiento de paleodistribución se utiliza para evaluar cambios en la distribución de 17 especies de Polylepis (Rosaceae) durante el Último Máximo Glacial (UMG), el Holoceno medio (HM) y el Antropoceno en los Andes centrales. Los modelos de distribución de especies (MDS) se calcularon combinando registros de presencia y condiciones climáticas actuales utilizando Maxent, y fueron proyectados a tres escenarios climáticos para el LGM (~ 22,000 años BP), el HM (~ 6,000 años BP) y el Antropoceno (1,970-2,000). Subsecuentemente, se identificó refugios climáticos intersectando el hábitat durante el LGM, HM y las condiciones actuales para todas las especies estudiadas. Los SDMs para LGM y HM mostraron una reducción del hábitat climático para la mayoría de las especies en comparación con el Antropoceno. Se identificaron cuatro zonas de refugio actuales: la Cordillera central del Perú, la cuenca del Lago Titicaca, la Cordillera Occidental de Bolivia y el norte de Chile. En general, los SDMs fueron consistentes con los patrones propuestos con registros de polen, y parcialmente con los estudios filogeográficos disponibles. Los refugios climáticos actuales son áreas muy importantes para la conservación de los bosques de Polylepis y la biota asociada. Este estudio ayuda a comprender los cambios de distribución durante milenios de las especies de árboles andinas en respuesta al cambio climático y a identificar áreas prioritarias para la conservación para la delimitación de estrategias de conservación futuras.

KEYWORDS:

• Climate refugia
• habitat suitability models
• niche modeling
• Polylepis

Introduction

The distribution of species has changed over time, contracting or expanding due to climate variability and human influence. As the cold conditions that prevailed during glacial cycles began to recede after the Last Glacial Maximum, major reorganizations of ecosystems produced biodiversity distribution patterns worldwide that involved changes in the abundance of species and large-scale migrations of taxa [1–3]. Nevertheless, a few important regions served as climate refugia and remained relatively stable, which helped to sustain certain specialized species and enable their persistence over time [4,5]. Understanding how species distributions changed in the past has important implications when defining conservation strategies, as it can provide useful knowledge on how species responded to extreme climatic events and to identify key conservation areas particularly for highly sensitive and endemic species.

Ongoing global climate change is presently impacting plant composition and ecosystem functions of Andean ecosystems. Recent research demonstrates that Andean ecosystems are changing as a consequence of increasingly warmer and drier conditions [6–9]. For instance, in the central Andes, between 1965 and 2012, maximum temperature increased at a rate of 0.15°C decade⁻¹ [10] and between 1985 and 2004, precipitation decreased while droughts episodes increased in frequency and intensity [11–13]. Global circulation models project a much more intense warming in the Andean highlands, with a potential increase in temperature ranging from 2.5°C to 5.9°C and a 20% decline in precipitation between 2070 and 2099, with potentially catastrophic consequences in ecosystem composition, structure, and the services they provide [14]. Hence, identifying key areas for the conservation of Andean biota is important for maintaining biodiversity in the future.

Polylepis (Rosaceae) includes 20–45 species of small trees and shrubs that inhabit the Andean highlands from Venezuela to Chile and Argentina [15,16]. Polylepis forms small and fragmented woodlands, where the genus usually constitutes the dominant or exclusive arboreal component [17]. In the Central Andes, Polylepis is found in a diversity of high-altitude vegetation formations, including the Puna, humid montane forests, inter-Andean dry forests, and the Tucuman-Bolivian forest [18–20]. These woodlands support a particular biota, including several endemic species that rely on these forests for foraging, nesting, and breeding [21–23]. Despite their ecological importance and high-conservation priority, Polylepis woodlands have decreased dramatically over the last decades due to the land-use intensification and global climate change [24–29].

Paleoecological studies suggest that Polylepis distribution and abundance have fluctuated over time due to climate variability and human disturbance [30,31]. During the Last Glacial Maximum (~25,000–21,000 years BP), ice and cold weather restricted the presence of Polylepis above 4000 m asl and, as a result, woodlands retracted to protected sites at lower elevations [30,32,33]. Subsequently, during the glacial–interglacial transition, conditions became warmer and more humid, which enable upslope migration of Polylepis and its expansion over the Altiplano. However, during the mid-Holocene (9,500–5,000 years BP), aridity prevailed, lowering lake levels and Polylepis distributional range contracted [30,33]. When mesic conditions returned between 5,000 and 4,000 yr BP, Polylepis was not able to expand as anthropogenic disturbance, particularly wood extraction and burning intensified, especially after agropastoralism became the main subsistence strategy throughout much of the Andes around 3,500 years ago [34–36]. After the European conquest, forest loss and habitat degradation increased due to intensification of economic activities, such as farming, grazing, and mining as well as an increase of firewood demand [37,38]. During the last several decades, loss of Polylepis woodlands has continued due to unsustainable land-use activities and ongoing global climate change [24–29].

Species distribution models (SDM) have become a common tool to reconstruct species responses to environmental changes, identifying climate refugia and predicting how species will respond to climate change [5,25,39]. SDM project suitable areas for a given species occurrence based the relationship between occurrence records and environmental characteristics that define the environmental niche [40]. These models can be used to identify patterns of modern distribution of a species based on suitable areas and use it to project changes in the species distribution under past or future climatic scenarios, assuming a constant biotic response or niche stability [5,41–44].

In this study, we seek to complement previous paleoecological reconstructions and identify possible areas of climate refugia by estimating the potential distribution of 17 species of Polylepis during the Last Glacial Maximum (LGM, ~22,000 yr BP), mid-Holocene (MH, ~6,000 yr BP) and the Anthropocene (1,970–2,000). We use species distribution models (SDM) to examine the geographic occurrence of suitable climate under paleoclimate scenarios and modern conditions based on present occurrence records of the species. This approach is based on the assumption of climate niche conservatism, meaning that past changes in the distribution of Andean trees were mainly the product of migration or local extinction in response to changing climatic conditions.

The main goal is to identify areas of long-term climatic stability in the central Andes (Peru, Bolivia, northern Chile, and northern Argentina) and to identify regions of high conservation value for Polylepis woodlands. We specifically aim to: a) define the environmental ecological niche for each species, b) assess if species distribution during the LGM, middle Holocene, and Anthropocene overlaps, and whether the range size shifts over time, and c) identify the current refugia, i.e. regions climatically suitable from the Last Glacial Maximum to the Anthropocene.

Materials and methods

Study area and species

A total of 17 species of Polylepis were selected that are distributed in the Central Andes. These Polylepis species were specifically selected because they occur in dry (i.e. Puna, inter-Andean dry forests and Bolivian-Tucuman forest) or in humid mountain regions (i.e. humid Puna and Yungas mountain forest). The Andean dry Puna (3,500–4,100 m) is characterized by low temperatures (2–7°C), scarce precipitation (250–400 mm), and high solar radiation [45]. In contrast, the humid Puna is distributed at elevation between 3,700 and 4,000 m asl (up to 4,100 m) and the weather is cold (7–9°C) and relatively humid (500–1,600 mm) [45]. The inter-Andean dry forest is found between 500 and 3,300 m asl, in regions with a (12–16°C) and dry (500–700 mm) weather [45]. The Bolivian-Tucuman forest (18–28° S) is distributed between 800 and 3,900 m asl and in areas with a relatively warm (5–23°C) and humid (700–2,000 mm) weather [45]. The humid mountain forests or Yungas is distributed between 1,000 and 4,200 m asl, the weather is warm (7–24°C) and humid (1,500 – >6,000 mm) with high cloud cover [45].

Species occurrences

We obtained occurrence records of the selected Polylepis species in the Central Andes from the Global Biodiversity Information Facility (https://www.gbif.org/) and herbarium collections (Missouri Botanical Garden, Herbario Nacional de Bolivia, and Herbario del Sur). All herbarium specimens were checked for accurate species identification with the support of specialists. Additional occurrences were obtained in the published literature, field reports, and field surveys conducted by the authors. To control for the quality of the occurrence data, we reduced spatial uncertainty by removing non-georeferenced records, outliers, and duplicate records for a given species. To account for the occurrence of several woodlands of the same species at a finer spatial scale, we imported presence records into QGIS (version 3.8 [46]) and kept only one record within 50 m of Euclidean distance. Locality data varied between 14 and 190 unique point localities per species.

Climatic variables and ecological niche models

To generate SDM, we employed 19 bioclimatic variables available at WorldClim (https://www.worldclim.org/). Current climatic variables (1,970–2,000) at 2.5-arc minute resolution were used to produce the current distribution models. To generate past distribution models during the Last Glacial Maximum (~21,000 years BP) and the middle Holocene (~6,000 years BP), three different circulation models were selected: The community Climate System Model (CCSM4), the Model for Interdisciplinary Research on Climate (MIROC-ESM), and the Max-Planck-Institute fur Meteorologie (MPI-ESM-P) with a spatial resolution of 2.5 arc minutes. As the study region, we delimited a rectangle that surrounded the extent of the Central Andes (latitude −80 to −70° S) and clip all environmental layers. This area was used to model the distribution of all the species. We acknowledge that it includes large regions that species do not inhabit because of limitations to dispersal; thus, potential distribution maps should be taken into consideration.

To model the potential distributions of species, we employed Maxent v. 3.3.3k [47]. Maxent utilizes a maximum entropy algorithm to model species distribution with presence-only data (Phillips et al., 2006). Given that a high correlation is expected between many climatic variables, we conducted an initial run in Maxent using all 19 bioclimatic variables for current conditions for each species. The selection of environmental variables was based on a stepwise removal of environmental variables. This technique assumes that the best performing combination of environmental variables is unknown and it aims to approximate this optimum through a stepwise removal of the least contributing variables [48]. We employed a MaxEent jackknife test to identify which bioclimatic variables contribute the least to the overall model for each species. The variables contributing with less than 0.01 were removed from the list of environmental variables, and this process was repeated for all the 17 Polylepis species [48,49]. The value of variable contribution within each species was determined heuristically via a jackknife approach in the Maxent algorithm, and it fluctuates slightly between each species. We set up Maxent to select randomly 75% occurrence localities for model training and the rest 25% localities for model testing, and 1,000 iterations with a convergence threshold of 0.00001 [50,51]. We subsequently generated modeled potential distribution models using current, LGM, and MH scenarios for three circulation models (CCSM, MIROC, and MPI-ESM-P). We used the receiver operating characteristic curve analysis (AUC) and receiver operating characteristic (ROC) to assess the predictive capacity of the models [52]. ROC analysis is a method designed to evaluate the specificity (absence of commission) and sensitivity (absence of omission error) at a diagnostic test [53]. The AUC provides a threshold-independent measure of model performance as compared with that of null expectations [51].

To identify if species distributions changed in the past, an occurrence binary prediction was generated for each species at each time period. The distribution threshold was set at the minimum probability area containing all the training localities (Phillips et al., 2006). Maxent produces a modeled species distribution with relative probabilities of occurrence between 0 and 100, however the minimum probability threshold allows for the identification of a standardized percentage of the distribution area, thus allowing for comparisons of different models (e.g. Anthropocene, MH, and LGM) [47]. We combined maps to identify grid cells classified as suitable across all scenarios, producing “agreement maps”, representing a conservative approach to account for variation between climate scenarios. Agreement maps were used to identify presumed species-specific refugia, by intersecting cells classified as suitable over consecutive times. “Current refugia” distribution models were constructed by identifying areas considered suitable during the LGM, MH, and Anthropocene (Supplementary information 1) [54].

Results

The distribution of sampling points in the environmental space (Figure 1) shows that several species share a similar climatic niche. For example, several species (P. rugulosa, P. subsericans, and P. tarapacana) currently distributed in the semiarid Puna occur under the same climatic conditions but P. tomentella occupies warmer areas. The species found commonly in the humid puna and inter-Andean dry forests (P. besseeri, P. subtusalbida, and P. triacontandra) also occupy areas with similar climates, except for P. incarum that occupies areas comparatively drier. The species that inhabit the humid montane forests (P. flavipila, P. lanata, P. neglecta, P. pauta, P. pepei, P. racemosa, and P. weberbaueri) showed a high overlap and little disaggregation in the climate space. The only two species distributed in the Tucuman-Bolivian forest, P. crista-galli and P. hieronymi, also exhibit high similarity in their climate niche.

Figure 1. Climatic envelope of 17 species of Polylepis in the central Andes as a function of: a) mean annual precipitation vs. mean annual temperature, b) maximum temperature of warmest month vs. minimum temperature of coldest quarter, and c) precipitation of driest quarter vs. precipitation of wettest quarter in the central Andes used for species distribution models

Across the 17 species, the average cross-validated test AUC was 0.98. AUC values ranged from 0.93 in P. racemosa and 0.998 in P. rugulosa and P. subsericans, which indicates good high classifier performance. The most important predictor variable for each SDM varied across species (Table 1). Nevertheless, temperature seasonality was one of the main predictors for 37.5% of the species, including P. besseri, P. neglecta, P. pepei, P. subsericans, P. substusalbida, and P. tarapacana. Minimum temperature of coldest month for 25% of the species, such as P. crista-galli, P. hieronymi, P. lanata, P. pauta, and P. subtusalbida. Some variables were important predictors for a very few species. For example, P. rugulosa distribution was mainly explained by the mean temperature of the coldest quarter and P. flavipila distribution was mostly associated with precipitation of coldest quarter.

Table 1. Most explanatory variables for model fit, according to the jackknife test of Maxent for 17 Polylepis species in the central Andes

Species	AUC	Bio1	Bio2	Bio3	Bio4	Bio5	Bio6	Bio7	Bio8	Bio9	Bio10	Bio11	Bio12	Bio13	Bio14	Bio15	Bio16	Bio17	Bio18	Bio19
Polylepis besseri	0.991				47.7	0.8	1.9	0.1		35.9			0.3	2		1.4	7.2	0.4	2.2
Polylepis crista-galli	0.996		0.1	6.9	0.1		74.9		0.4				1.1		14.6	0.1		0.6	1	0.1
Polylepis flavipila	0.993			4.8	5.1	0.2			1.3								2.1			86.5
Polylepis hieronymi	0.996		0.1	1.8	2.9	20.5	41.4		9.3	0.5				4.7	0.3			0.4	18
Polylepis incarum	0.993			0.8	9.9	0.4			71.5	0.9	3.5			0.2	0.1	1.1	10.4		0.9	0.4
Polylepis lanata	0.995		0.1		11.1		75.3			0.9			0.8	0.5	0.7		4.7		5.4	0.5
Polylepis neglecta	0.986			2.5	67.5	14.4	3.2		0.5					10.6		0.5		0.7
Polylepis pauta	0.988			0.1	30.7	0.4	42.8						0.3	3.9	0.5	8.8	1.2	2	0.2	9.1
Polylepis pepei	0.993			10.4	50.1		1.8	1		29.7			0.3	3.8		0.3	0.4	0.4	0.7	0.9
Polylepis racemosa	0.933			5.6	15.5		17.7	2		1.3	0.2	0.2	0.1	17.8	3	7.1	16.7		3.3
Polylepis rugulosa	0.998		0.3	2.7	9.9							78.6				7.7			0.1	0.7
Polylepis subsericans	0.998		0.9		61.3		0.3			10.1				0.5	2.2	0.8			23.3	0.6
Polyelpis subtusalbida	0.995			7.8	21.9		48.2			1.5			2.5	8.3	2.3				7.5
Polylepis tarapacana	0.989	0.3		0.3	16.6	0.3		0.5	65.9	0.5			0.3	4.9		3.5	0.3	5.8	0.6
Polyelpis tomentella	0.986	0.1	0.2	14.4	56	3.1	1.4		0.3	0.7	0.2		0.4	5.6	9.1	6.2	0.6		1.1	0.7
Polylepis triacontandra	0.996			25.2	59	4.2				1.1				0.9		3.4	0.7			5.5
Polylepis webarbeueri	0.993			9.8		2.1		3	23.5		32.6	2.8			4.4		11.2	5.3

Paleoclimate models projected different suitable habitats for the studied Polylepis species. During the LGM, the CCSM4, MPI-ESM-P, and MIROC-ESM models projected a consistent increase in suitable area for P. hieronymi, P. pepei, P. tarapacana, and P. tomentella and a decrease in suitable area for P. crista-galli, P. flavipila, P. pauta, and P. subsericans (Figure 2). SDMs showed inconsistent outcomes (i.e. contractions under some climate scenarios and expansions under other climate scenarios) for some species, including P. lanata, P. neglecta, P. racemosa, P. subtusalbida, and P. weberbaueri. The SDMs of P. besseri project a climatically suitable area on the region occupied by Lake Titicaca, we considered that this is an over-projection and that it is unlikely that the species colonized this region in the past.

Figure 2. Potential distribution maps for 17 species of Polylepis in the central Andes generated with Maxent under the Last glacial Maximum (~22,000 yr BP) on the left, the mid-Holocene (~6,000 yr BP) center, and current conditions (1,970–2,000) on the right. Areas where probability of occurrence was predicted by only one climate scenario are indicated with the lightest color (1), areas predicted by two climate scenarios are indicated with an intermediate color (2), while areas where probability of occurrences were predicted by three climate scenarios are depicted with the darkest color (3)

Figure 2. (Continued)

Figure 2. (Continued)

During the middle Holocene, SDMs projected a relative increase in suitable area compared to the LGM suitable area for P. incarum, P. triacontandra, P. lanata, and P. racemosa (Figure 2). In contrast, SDMs showed a decrease in suitable area for P. neglecta, P. rugulosa, and P. weberbaueri. SDMs generated for P. pauta, P. subsericans, P. tomentella, and P. triacontandra showed inconsistent projections.

Agreement maps showed that suitable habitat spanned a smaller spatial extent in both previous time periods compared to current conditions for all the species (Table 2). During the LGM, the area of suitable habitat was on average 67% smaller that the warmer modern conditions, with the exception of P. hieronymi, P. rugulosa, P. tarapacana, and P. triacontandra which showed a larger suitable area in the past (range: 141–327%). During the MH, despite conditions became warmer and drier, the suitable area increased for 76% of the species, although for most of the species it was still smaller than the current suitable area. Exceptions were P. incarum, P. lanata, P. racemosa, and P. triacontandra, for which suitable areas were larger than modern scenarios, between 134% and 1158%.

Table 2. Predicted distribution change (%) under three different climate scenarios for 17 Polylepis species in the Last Glacial Maximum (22,000 yr BP) and the middle Holocene (6,000 yr BP) relative to the present-day distribution area

Species	Last Glacial Maximum				Middle Holocene
	CCSM4	MPI-ESM-P	MIROC-ESM	Consensus	CCSM4	MPI-ESM-P	MIROC-ESM	Consensus
Polylepis besseri	86.16	133.94	432.38	43.08	167.89	205.22	441.25	98.43
Polylepis crista-galli	27.42	40.32	25.00	16.53	78.63	14.52	1.61	1.21
Polylepis flavipila	51.63	14.18	69.89	4.25	18.93	65.47	54.96	13.59
Polylepis hieronymi	120.23	262.99	219.08	221.84	5.06	2.30	0.23	0.23
Polylepis incarum	136.18	72.36	167.34	35.68	1327.14	1753.27	2108.54	1158.29
Polylepis lanata	100.00	89.20	501.70	38.64	442.61	346.59	854.55	245.45
Polylepis neglecta	35.31	70.39	160.31	10.31	19.08	35.75	0.00	17.32
Polylepis pauta	48.12	139.38	53.59	6.75	127.22	81.95	40.05	16.25
Polylepis pepei	101.28	1019.87	268.59	10.26	294.87	1028.21	1026.28	85.90
Polylepis racemosa	100.37	74.21	29.16	30.09	256.07	344.67	236.07	134.77
Polylepis rugulosa	260.44	262.22	327.11	327.11	7.56	0.00	10.67	18.22
Polylepis subsericans	55.12	4.43	9.70	1.94	210.25	119.11	62.33	58.45
Polyelpis subtusalbida	5.34	53.74	146.26	3.91	270.11	352.67	1246.98	5.69
Polylepis tarapacana	232.93	192.02	296.61	165.87	195.81	134.73	128.34	75.65
Polyelpis tomentella	387.33	675.21	269.42	78.51	130.03	97.52	54.55	8.54
Polylepis triacontandra	235.38	246.19	44.96	141.52	408.35	397.05	1295.58	238.08
Polylepis webarbeueri	101.10	17.28	1927.21	11.76	77.57	69.12	58.82	17.65

The current refugia model (the intersection between LGM, MH, and modern suitable habitat) showed a similar distribution of the centers of climate refugia for Polylepis species (Figure 3). The areas identified were the central Cordillera of Peru, the Lake Titicaca basin, the western Cordillera of Bolivia, and the semiarid Andes of Bolivia, Chile, and Argentina.

Figure 3. Centers of current refugia for 17 species of Polylepis in central Andes, based on habitat suitability generated with Maxent. Current refugia are defined as the intersection of areas classified as suitable habitat over the Last Glacial Maximum, mid Holocene and Anthropocene. Highlighted areas show climate refugia for a) central Cordillera of Peru, b) Lake Titicaca basin, c) western Cordillera of Bolivia, and d) semiarid Andes of Chile . Areas where probability of occurrence was predicted by only one climate scenario are indicated with the lightest color (1), areas predicted by two climate scenarios are indicated with an intermediate color (2), while areas where probability of occurrences were predicted by three climate scenarios are depicted with the darkest color (3)

Discussion

The Last Glacial Maximum (~22,000 yr BP) was characterized by cooler conditions (5–8°C colder than modern conditions) that restricted vegetation cover and resulted in the contraction and downslope migration of tree species [30,55,56]. SDMs suggest that the suitable habitat for most of the studied Polylepis species was between 22% and 98% smaller during the LGM than under current conditions. However, our analysis does not take into account the area extent covered by glaciers, and it is feasible that ice cover further restricted the presence of Polylepis at higher elevations (>4,000 m asl), particularly P. tarapacana an P. incarum as they currently inhabit such elevations. Paleovegetation reconstructions have shown that at elevations above 4,000 m asl, the occurrence of Polylepis was scarce due to the presence of glaciers and cold conditions during the Last Glacial Maximum (20,000 and 17,000 BP) [30–32]. Palynological studies have also shown that Polylepis remained abundant between 2,800 and 3,400 m asl during the LGM [31]. This niche could have been occupied by a diversity of species that currently inhabit lower elevations, such as P. pauta, P. pepei, P. incarum, and P. racemosa.

In a few species, paleoclimate models projected an increase in suitable area during the LGM. These species currently inhabit a diversity of ecosystems, as cold-adapted P. tarapacana in the arid Puna. Palynological studies have suggested that P. tarapacana woodlands persisted near Salar de Uyuni during the LGM [57] but were possibly absent near Nevado Sajama [58]. Furthermore, phylogeographic studies of P. tarapacana in the semiarid Central Andes found low levels of genetic differentiation between populations [59,60]. The low levels of genetic differentiation could be either the product of extensive genetic flow among populations or larger but disjunct woodlands during the late Pleistocene.

During the mid-Holocene, temperatures increased substantially (1–3°C warmer than today) and precipitation decreased considerably, generating arid conditions in the central Andes [56]. Species distribution models suggest that suitable habitat expanded from the LGM to the mid-Holocene for most of the species, in particular, for P. hieronymi, P. rugulosa, P. tarapacana, and P. triacontandra. These results agree with fossil records that have shown an increase in Polylepis abundance as well as its expansion during the MH (~9,500–5,000 yr BP) in several areas of the central Andes, including Lake Titicaca, Lake Miski, Chocheras, Huamanmarca, and Chochos [30–33,56].

It is likely that increasing fire activity may have further resulted in the contraction of suitable habitat of Polylepis during the MH. Evidence has shown that fire regime also increased in the central Andes during this period [55,56]. Fire is an important modulatory mechanism for the distribution of Polylepis, as it promotes patchiness by increasing adult mortality and reducing seedling recruitment [24,61]. However, the effect of fire on Polylepis cover is heterogeneous across the Andes. For example, in Khomer Kocha (central Bolivia) and Pacucha (Peru), a marked decline in Polylepis during the MH was associated with increasing fire activity [30–32,55,62]. By contrast, despite fire was a persistent element of the landscape in Huamanmarca and Miski (southern central Peru), Polylepis remained as a relative dominant component of the vegetation during periods of increased burning between ~7,000 and 2,500 yr BP [33]. Topographic relief and water balance have been attributed to play an important role in the permanence of Polylepis woodlands in these sites, as an irregular terrain likely acted as a physical barrier by protecting plants and generating humid microenvironmental conditions [29,30,33].

It is also worth mentioning that isolated and small populations of Polylepis species likely survived outside of the modelled range boundaries likely favored by specific and fluctuating microclimatic and microenvironmental conditions [31,33,63]. For example, Polylepis forests likely thrived during the middle Holocene in the high-rugosity catchments of Huamanmarka and Miski in the Peruvian southern cordillera [33]. Nevertheless, despite the importance of micro-refugia for Polylepis permanence, the spatial resolution used to run the simulations cannot be used to identify them, highlighting the potential that using increasingly finer-scale datasets can have for improving the fidelity and precision of the results.

Current refugia were defined as areas that have remained suitable over millennia and across the LGM, MH, and Anthropocene. The modeled results suggest that climate refugia overlap for several species of Polylepis and we were able to identify four refugia hotspots located in the Central Cordillera of Peru, the Lake Titicaca basin, the western Cordillera of Bolivia, and the semiarid Andes of northern Chile. These putative climate refugia hotspots cover a diversity of ecosystems, including the semiarid Puna, the meso-Andean valleys, the humid montane forests, the Inter-Andean dry forests, and the Tucuman-Bolivian forest. Pollen records provide evidence of the existence of Polylepis woodlands in several locations corresponding to putative refugia hotspots [32,55,56,58,62,64–67]. The agreement between past and present-day distributions for several Polylepis species could be elucidated by relative climate stability in the refuge areas. Local conditions seem to have buffered the extreme effects of the climate variability during the Late Pleistocene and Holocene, contributing to the survival of Polylepis species.

The identification of current refugia hotspots has important conservation implications for the conservation of Polylepis woodlands in the central Andes. These hotspots also include areas that have been previously identified as high conservation priority areas due to high species richness and endemism of plants and animals, including Polylepis-associated fauna [21,68–70]. The SDMs also support the importance of the southern-central Andes as a high conservation area [71]. The Bolivian inter-Andean dry forests and the Tucuman-Bolivian forest are highlighted as regions that merit further investigation efforts and conservation attention.

The current refugia hotspots could constitute stable areas for the permanence and survival of Polylepis species. These areas may harbor Polylepis woodlands with high levels of genetic diversity in comparison to areas non-climatically stable. At the moment, few regional phylogenetic studies of Polylepis have been conducted, and it is unknown if this pattern hold for the other species. Results have shown, for example, that the semiarid northern Andes of Chile constitutes an area of relatively high-genetic diversity for P. tarapacana [59]. Multi-species genetic diversity studies between putative refugia and colonization areas could help to further understand the complex biogeographic history of Polylepis in the Andes.

This study provides evidence that SDM can be used as a complementary approach to paleoecological studies, offering spatially explicit projections that can be used to draw hypothesis about past and future geographic distributions of Andean species. The development of SDMs under different scenarios could greatly be benefited by taken into consideration the distributions of glaciers, paleolakes, fire occurrences, and archeological data to model the past role of environmental variability and humans on the distribution of Andean trees. Furthermore, the integration of ecological-niche and palaeoecological studies constitutes a useful approach for identifying critical areas in which to prioritize research and conservation efforts. This is particularly important considering the challenges of designing long-term conservation strategies under the ongoing climate crisis.

Author contributions

AID conceived and designed the study. AID and JMC analyzed the data and wrote the manuscript.

Acknowledgments

We thank the Missouri Botanical Garden, Herbario Nacional de Bolivia, and the Herbario del Sur for providing access to plant collections. Stephan Beck, Rosa Isela Meneses, Martha Serrano, and Isabel Gomez for the logistic support. Sarah Ivory provided valuable comments and reviewed preliminary versions of the manuscript. This work was supported by the International Foundation of Science (IFS) under Grant D-4839-2.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the International Science Foundation [D-4839-2].