Climate Change Can Drive a Significant Loss of Suitable Habitat for Polylepis quadrijuga, a Treeline Species in the Sky Islands of the Northern Andes

It is predicted that climate change will strongly affect plant distributions in high elevation “sky islands” of tropical Andes. Polylepis forests are a dominant element of the treeline throughout the Andes Cordillera in South America. However, little is known about the climatic factors underlying the current distribution of Polylepis trees and the possible effect of global climate change. The species Polylepis quadrijuga is endemic to the Colombian Eastern Cordillera, where it plays a fundamental ecological role in high-altitude páramo-forest ecotones. We sought to evaluate the potential distribution of P. quadrijuga under future climate change scenarios using ensemble modeling approaches. We conducted a comprehensive assessment of future climatic projections deriving from 12 different general circulation models (GCMs), four Representative Concentration Pathways (R) emissions scenarios, and two different time frames (2041–2060 and 2061–2080). Additionally, based on the future projections, we evaluate the effectiveness of the National System of Protected Natural Areas of Colombia (SINAP) and Páramo Complexes of Colombia (PCC) in protecting P. quadrijuga woodlands. Here, we compiled a comprehensive set of observations of P. quadrijuga and study them in connection with climatic and topographic variables to identify environmental predictors of the species distribution, possible habitat differentiation throughout the geographic distribution of the species, and predict the effect of different climate change scenarios on the future distribution of P. quadrijuga. Our results predict a dramatic loss of suitable habitat due to climate change on this key tropical Andean treeline species. The ensemble Habitat Suitability Modeling (HSM) shows differences in suitable scores among north and south regions of the species distribution consistent with differences in topographic features throughout the available habitat of P. quadrijuga. Future projections of the HSM predicted the Páramo complex “Sumapaz-Cruz Verde” as a major area for the long-term conservation of P. quadrijuga because it provides a wide range of suitable habitats for the different evaluated climate change scenarios. We provide the first set of priority areas to perform both in situ and ex situ conservation efforts based on suitable habitat projections.

It is predicted that climate change will strongly affect plant distributions in high elevation "sky islands" of tropical Andes. Polylepis forests are a dominant element of the treeline throughout the Andes Cordillera in South America. However, little is known about the climatic factors underlying the current distribution of Polylepis trees and the possible effect of global climate change. The species Polylepis quadrijuga is endemic to the Colombian Eastern Cordillera, where it plays a fundamental ecological role in highaltitude páramo-forest ecotones. We sought to evaluate the potential distribution of P. quadrijuga under future climate change scenarios using ensemble modeling approaches. We conducted a comprehensive assessment of future climatic projections deriving from 12 different general circulation models (GCMs), four Representative Concentration Pathways (R) emissions scenarios, and two different time frames (2041-2060 and 2061-2080). Additionally, based on the future projections, we evaluate the effectiveness of the National System of Protected Natural Areas of Colombia (SINAP) and Páramo Complexes of Colombia (PCC) in protecting P. quadrijuga woodlands. Here, we compiled a comprehensive set of observations of P. quadrijuga and study them in connection with climatic and topographic variables to identify environmental predictors of the species distribution, possible habitat differentiation throughout the geographic distribution of the species, and predict the effect of different climate change scenarios on the future distribution of P. quadrijuga. Our results predict a dramatic loss of suitable habitat due to climate change on this key tropical Andean treeline species. The ensemble Habitat Suitability Modeling (HSM) shows differences in suitable scores among north and south regions of the species distribution consistent with differences in topographic features throughout the available habitat of P. quadrijuga. Future projections of the HSM predicted the Páramo complex "Sumapaz-Cruz Verde" as a major area for

INTRODUCTION
The Andean region is a South American biodiversity hotspot, home to an impressively high diversity and endemism (Myers et al., 2000;Antonelli and Sanmartín, 2011). Regrettably, tropical Andean mountain species are classified as highly vulnerable to ongoing climate change, which may lead to strong shifts in distribution ranges, populations decline, and local extinction (Mittermeier et al., 2011;Mountain Research Initiative EDW Working Group, 2015). Although future distribution predictions for tropical Andean species are often challenging and carry a high uncertainty, due to the complex topography and climate of the region, they often point toward an extreme sensitivity to climate change (Lawler et al., 2009;Tovar et al., 2013;Tejedor Garavito et al., 2015;Peyre et al., 2020). Specifically, it is hypothesized that the main factors impacting future species distribution in the Andes will be the increase of temperature, changes in water fluxes, and the upward shift of the lower snowfall limit (Vuille et al., 2008(Vuille et al., , 2015Herzog et al., 2011). Furthermore, several studies have already predicted an increase in threats and local extinction risks by the second half of the century in the Andes region (Ramirez-Villegas et al., 2014;Peyre et al., 2020;Valencia et al., 2020).
Polylepis forests are unique montane ecosystems, naturally occurring at high elevations in the Andes (Hoch and Körner, 2005;Navarro et al., 2005). These forests are generally present as sparse patches separated by vast areas of Páramo or puna grasslands (Fjeldså, 2002;Jameson and Ramsay, 2007) dominated mainly by Polylepis spp., that are endemic to South American mountain ranges, and tend to occur along the upper treeline (Simpson, 1979(Simpson, , 1986Fjeldså et al., 1996). Polylepis treeline forests are the source of essential cultural, economic, and ecosystem services, providing pivotal ecological functions such as absorbing moisture from the clouds and releasing water in springs and rivers while being optimal habitats for other threatened and endemic species (Kessler, 2006;Donald et al., 2010). These forests require special conservation efforts (Purcell and Brelsford, 2004;Pinos, 2020) as they are threatened by deforestation, fires, farming, timber extraction, and land-use change (Fjeldså et al., 1996;Etter and Villa, 2000;Kessler, 2002;Hensen et al., 2011).
Several studies already support the importance of environmental variables such as elevation, precipitation, and mean annual temperature in shaping the current distribution patterns of different Polylepis species (Zutta et al., 2012;Boza Espinoza et al., 2019). Constraints related to water availability and soil water content during the dry season, hence precipitation and low-temperature extremes define the central eco-physiologic adaptations of Polylepis species (Velez et al., 1998;Morales et al., 2004;Ramos et al., 2013). Moreover, the decline of Polylepis diversity and abundance is one of the hypothesized effects of global or regional climate change patterns on Andean ecosystems (Gareca et al., 2010;Anderson et al., 2011). In particular, recent studies predicted severe reductions of suitable habitat for Polylepis tarapacana, P. incana, and P. reticulata before the end of the twenty-first century (Ramirez-Villegas et al., 2014;Cuyckens et al., 2016). It is known that closed forest stands can maintain their own microclimate and that their irregular cohort establishment depends on the occurrence of favorable (wet) years (Rada et al., 2011). The possible importance of moisture conditions for Polylepis spp. has already been reported for northern Chile (Rundel et al., 2003). Also, it has been postulated that any increase in temperature will further restrict Polylepis australis distribution in central Argentinian Andes to the uppermost elevations (Marcora et al., 2008), therefore underlining a high sensitivity to high temperatures for the genus. Therefore, understanding the consequences of climate change on Polylepis woodlands, as one of the most sensitive Andean biomes , should be a priority for conservation biology in the Andean region, as necessary means to guide conservation policies and management strategies in the coming decades (Gosling et al., 2009;Sousa-Silva et al., 2014;Cuyckens et al., 2016). Also, it has to be taken into account that plausible niche specialization processes taking place in the different regions of the species distribution may lead to a reduced model ability to predict suitable areas under the assumption of niche equilibrium (Guisan and Thuiller, 2005;Parmesan, 2006;Hargreaves and Eckert, 2014).
In the specific case of the Colombian Andes, it is known that Polylepis quadrijuga Bitter was one of the first native species to be exploited (Rodríguez et al., 2018) and that continuous expansion of agricultural and livestock frontiers has led to increased fragmentation and isolation of P. quadrijuga woodlands. Thus, it is important to expand the understanding of the main ecological drivers of P. quadrijuga distribution, and to estimate the extent of its current distribution, while evaluating the effects of ongoing climate change on the species. Here, we selected this treeline species as a case study to address the influence of climate change on the distribution of key highland species in the northern Andean region. We apply Habitat Suitability Modeling (HSM) to infer the climate-change vulnerability of P. quadrijuga, one of the recorded species of Polylepis in Colombia (Boza Espinoza et al., 2019). Moreover, we specifically aim to (i) identify the key climate variables shaping the current geographical distribution of P. quadrijuga, (ii) estimate present and future suitable areas and, (iii) identify climatically stable areas as potential sites for the implementation of protection efforts. Additionally, we evaluate to what extent the current protected areas network and Páramo complexes of Colombia enclose the predicted future suitable areas under climate change for P. quadrijuga, and if these can be sufficient for the in situ conservation of this species. Finally, we discuss possible conservation strategies and effective management actions that can ensure the persistence of this Colombian iconic Andean treeline species, addressing the specific protected areas that could host suitable habitat for P. quadrijuga under current and future climate change conditions.

Study Species and Area
Polylepis quadrijuga (Rosaceae), locally known as "colorado" or "sietecueros" is an evergreen tree of up to 10 m height with compound leaves and a twisted trunk with continuously peeling reddish bark (Simpson, 1986). This species is easily recognized, its taxonomy resolved, and its distribution does not overlap with other Polylepis species Segovia-Salcedo et al., 2018). P. quadrijuga is wind-pollinated and the seeds are dispersed by the wind during the dry season (Simpson, 1979;van Schaik et al., 1993;Schmidt-Lebuhn et al., 2007). Particularly, P. quadrijuga is endemic to the Colombian Eastern Cordillera slopes and rocky areas that are usually immersed within the Páramo grassland matrix (Simpson, 1986;Schmidt-Lebuhn et al., 2006;Rangel-Ch and Arellano-P, 2010;Pérez-Escobar et al., 2018;Peyre et al., 2018). Also, P. quadrijuga occurs often above the timberline and forms dense stands at elevations as high as 4,000 m, usually along streams (Simpson, 1979).
The Northern Andes present a high physiographic and topographic complexity, are considered one of the most diverse areas in the world (Myers et al., 2000;Jenkins et al., 2013;Jung et al., 2020), and are characterized by a very moist climate, low thermal seasonality, and marked diurnal temperature variations (Buytaert et al., 2006). Regionally, the climate of the Colombian Eastern Cordillera ranges from humid on the eastern flank to sub-humid or slightly dry on the western flank. The precipitation regime is bimodal, with pronounced dry and wet seasons, and annual rainfall ranging from 600 to 1,300 mm/y (Etter and Villa, 2000).
Polylepis species possess physiological adaptations to low temperatures and physiological drought that allow them to grow under colder conditions than the global mean for high elevation treeline species (Fjeldså et al., 1996;Kessler et al., 2014). Key physiological traits such as high foliar concentration of flavonoids and seasonal phenology could play an important role in the adaptation of this species to the severe environmental conditions of the Andean highlands (Velez et al., 1998). Despite these appealing features, Polylepis has been little studied. Moreover, as a result of the long history of deforestation in the high Andean forests of Colombia, the remaining patches of Polylepis forests are spread out over the open modified landscape (Etter and Villa, 2000;Rangel-Ch and Arellano-P, 2010), and several studies highlight that the relict patches of Polylepis forests in Colombia should be considered as a priority for conservation programs.

Species Occurrence Data
We compiled occurrence data for P. quadrijuga from three different sources: (i) the Global Biodiversity Information Facility (GBIF), initially 104 georeferenced records, where we also checked the nomenclature for synonyms (Borja-Acosta, 2017, 2019Raz and Agudelo, 2019;García et al., 2020;Grant and Niezgoda, 2020;Marín and Moreno, 2020;Rodríguez Bolaños, 2020;Lozano Bernal, 2021;Magill et al., 2021); (ii) curated citizen science projects observations (iNaturalist, 49 records), which included geographical coordinates and a photographic record, allowing an accurate identification check (Ueda, 2021); and (iii) direct virtual exploration inside the main Páramo complexes of the Eastern Cordillera, using the Google Street View Imagery tool (Google Maps, 23 records). We surveyed all the rural roads with available street view imagery, crossing the Páramo districts "Boyacá, " "Cundinamarca, " and "Santanderes" (Sarmiento et al., 2013). Through direct inspection of hundreds of 3D georeferenced panoramas, we certified the occurrence of P. quadrijuga when we observed multi-trunk trees or shrubs showing the combination of vegetative characters as the distinctive reddish bark peeling in thin papery layers and compound leaves with grayish-green color 1 . No other group of plants could be mistakable with Polylepis in the Páramo landscape (Simpson, 1979), and P. quadrijuga is the only species of the genus recorded from the Eastern Cordillera of Colombia (Rangel-Ch and Arellano-P, 2010;Fajardo-Gutiérrez et al., 2018). Therefore we considered this methodological approach adapted from Hardion et al. (2016) to be adequate in this study.
Only records with global positioning system (GPS) coordinates and detailed localities were used, and then visually inspected for obvious geocoding errors. We included all the occurrences covering the entire known geographic distribution of the species (Figure 1 and Supplementary Table 1) except the ones belonging to the Páramo of Frontino, located in the Western Cordillera of Colombia in the Antioquia department, which were removed, as these are now believed to form a separate species, and are currently under taxonomic reevaluation (Boza and Kessler, pers. comm.). We removed records with uncertain locations and retained occurrences with a minimum distance of 1 km from each other to avoid spatial pseudo-replication, yielding a final set of 99 records: 37 from GBIF (41%), 40 from iNaturalist (40%) and 22 from Street View Imagery (22%) (Supplementary Table 1). We also evaluated potential geographical sampling bias in our species occurrence dataset due to differences in accessibility. For that, we assessed the biasing effect of roads on the sampling intensity using the Bayesian approach implemented in the R package sampbias (Zizka et al., 2020). We calculated the bias effect on the three different sources of observed records (GBIF, iNaturalist, and Google Street View Imagery) and ten simulated datasets of same size as our complete dataset under complete random sampling and other 10 highly biased simulated datasets limiting the random sampling to maximum 1 km from roadways within the species distribution range (both geographical and altitudinal). The biasing effects show that differences in the observations among the three data

Environmental Data Sets and Abiotic Predictor Variables
We gathered a total of 23 biologically relevant climatic variables: 19 bioclimatic variables derived from monthly temperature and precipitation climatology for the years 1979-2013 available at 30 arc sec resolution (1 km at the equator) were retrieved from the CHELSA dataset (Karger et al., 2017) 2 ; Slope and Aspect variables derived from a SRTM Digital Elevation Model using the raster package (Hijmans, 2015); and two topographic variables (Terrain roughness index and Topographic wetness index) at same resolution to the bioclimatic variables were obtained from the ENVIREM data set (Title and Bemmels, 2018). The bioclimatic variables have shown a high importance for defining the environmental niche of species of Polylepis (Zutta et al., 2012;Boza Espinoza et al., 2019) and are directly affected by climate change. The CHELSA climatic dataset has been proven to provide better results for topographically complex montane regions (Bobrowski and Schickhoff, 2017). Moreover, the topographyrelated variables that we incorporated, are secondary geographic predictors that are used to characterize spatial topographic diversity, terrain complexity and soil moisture patterns (Lin et al., 2006), and are very sensitive to local differences in elevation, being able to capture the presence of steep gradients and high spatial heterogeneity (Różycka et al., 2017). 2 http://chelsa-climate.org/bioclim/ We cropped and aligned the explanatory environmental layers to an area encompassing 0 • -12.5 • N and 79.5 • -69.1 • W, using the raster package (Hijmans, 2015) in R. This geographic extent includes the full known extant range of P. quadrijuga. In order to avoid multicollinearity-related noise, we performed variable selection using the function corselect of the R package fuzzySim (Barbosa, 2015). This function calculates Pairwise Pearson correlations among all explanatory variables and then, among each pair of highly correlated variables (
Following Barve et al. (2011), we sampled the pseudoabsences within the accessible area for the species. We outlined the area including the Colombian Eastern Cordillera above 1,500 m of elevation and outside a 5 km buffer around the presence records (Supplementary Figure 3). Within this area, we randomly selected 10,000 pseudo-absences. Prevalence was maintained at 0.5, allowing the sum of presence-weights to be equal to the sum of pseudo-absence weights in the model calibration process (Barbet-Massin et al., 2012). We ran ten cross-validation replicates where presence records were split into sets of 75% for training and 25% for testing the models. For each cross-validation replicate we performed five runs with different pseudo-absence sets, therefore completing a total of 400 models. To assess the predictive performance of the models, we used both the threshold-independent area under the receiver operating characteristic curve statistic-AUC (Phillips et al., 2006) and the threshold-dependent true skill statistic-TSS (Allouche et al., 2006). Consensus models were obtained using a TSS-weighted average method to account for the predictive power of each model. Models with low predictive power (TSS < 0.7) were discarded.
As an additional characterization of the multi-dimensional occupied environmental space of P. quadrijuga, we extracted selected predictor values in correspondence of our species occurrences and performed a principal components analysis (PCA), with the function prcomp in R (R Core Team, 2020).

Evaluation of Future Range and Elevation Shifts and Representation in the Protected Areas System
We assessed the possible impacts of climate change on the distribution of P. quadrijuga through the projection of the HSMs to two time periods (2041-2060 and 2061-2080). For future projection we retrieved the projected bioclimatic variables from the CHELSA dataset, according to four greenhouse gas concentration scenarios, including low (RCP 2.6), moderate (RCP 4.5), high (RCP 6.0), and strong (RCP 8.5) increases in global radiative forcing (IPCC, 2013). These scenarios express the intensity of climatic changes by representative concentration pathways, using predicted values for radiative forcing for 2100, in comparison to the preindustrial era (Harris et al., 2014). The RCP2.6 scenario assumes CO 2 concentrations of 450 ppm in 2100, that global mean temperatures will increase 0.2-1.8 • C and has no strict reference in previous (Fourth Report) IPCC guidelines. RCP4.5 assumes 650 ppm CO 2 and 1.0-2.6 • C in 2100 and refers to the previous guideline scenario B1. RCP6.0 assumes 850 ppm CO 2 and 1.3-3.2 • C in 2100 and refers to the previous guideline scenario B2. RCP8.5 assumes 1350 ppm CO 2 and 2.6-4.8 • C in 2100 and refers to the A1F1 scenario of previous IPCC guidelines (Harris et al., 2014).
In addition, to consider the inherent uncertainty derived from each different general circulation model, we performed a comprehensive assessment including twelve different CMIP5 General Circulation Models (GCMs) [CESM1(CAM5), CSIRO-Mk3-6-0, FIO-ESM, GFDL-ESM2G, GFDL-ESM2M, GISS-E2-H, IPSL-CM5A-MR, MIROC-ESM-CHEM, MIROC5, MRI-CGCM3, NorESM1-M, and CCSM4]. For the resulting future climate projections, we generated consensus distribution maps for each time period and RCP scenarios, outlining the suitable areas predicted by at least three binary-transformed ensemble projections for the different GCMs. Future bioclimatic variables for all the GCMs, RCPs and time periods were retrieved from the CHELSA database. The topographic variables were kept invariant for all future projections.
We intersected the current and future HSMs projections to estimate the sensitivity to climate change at the species level, and assess the distribution range stability for all future RCPs scenarios and time periods. To this aim, we calculated the percentage of predicted area that remains stable, as well as range expansions or contractions in relation to the current projected area. Values of less than 100% indicate a decrease in environmental suitability over time (i.e., a high vulnerability of the species to climate change), whereas values greater than 100% indicate an increased extent of suitable habitat (i.e., a low vulnerability of the species to climate change).
In order to identify the specific Protected Areas (PAs) that could be relevant for the protection of P. quadrijuga in the future, we measured the overlapping area between the different projections of suitable habitat and the areas of 10 out of 36 Páramo complexes of Colombia (PCC), belonging to the Páramo districts "Boyacá, " "Cundinamarca, " and "Santanderes" (Sarmiento et al., 2013). Colombia has defined 17 Páramo districts and 36 Páramo complexes, by grouping natural areas that are divided by geographical features or human transformed landscapes that constitute special management and protection zones (Sarmiento et al., 2013). We also calculated the overlap between the suitable habitat projections and 36 protected areas of the National System of Protected Areas of Colombia (SINAP, 2020), including five National Natural Parks 3 . We used the PCC and the SINAP maps and calculated the portion of the suitable area for P. quadrijuga within these protected areas under the current conditions and future scenarios (RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5), for the periods 2041-2060 and 2061-2080.

Model Performance and Current Suitable Area Distribution
The projection of the ensemble-HSM to current climate conditions predicted a suitable area of 4444.9 km 2 , retrieving the highest and spatially continuous values of suitability throughout the southern range of the species distribution (SCV and CHG Table 1 and Figure 1), and less suitable and  Figure 1). Less than a half of the predicted suitable area (46.1%, 1813.4 km 2 ) was covered by the protected areas network of the SINAP, but most of the suitable areas (95.8%, 3765.3 km 2 ) were located inside of the PCC boundaries, with only 4.2% was located outside of the PCC protection (163.4 km 2 ) (Figure 1 and Table 1). The algorithms returning higher model performance metrics were GBM, GLM, Maxent, and RF and, the HSM returned overall high model evaluation metrics, with average AUC = 0.96 ± 0.04 and TSS = 0.9 ± 0.07 (Supplementary Figure 4), underscoring a good model performance in discriminating the species' fundamental climatic niche.

Importance of Environmental Predictors
The PCA identified several important climatic variables that highlight differences in niche breadth through the species distribution. The analysis revealed differences in niche occupancy with a strong geographic signal represented by the different Páramos complexes (Figures 1, 2). The first axis of PCA (35% of the total variation) synthesized a distribution gradient between the eastern and western flanks of the Eastern Cordillera, with the variables Mean Diurnal Range and Precipitation of Warmest Quarter showing the highest loadings, as well as Temperature Seasonality and Topographic wetness index, but with an opposite direction (Figure 2). The second axis (ca. 21% of the variation) with Precipitation of Driest Quarter showing the highest loading, highlighted a north-south gradient, with the exception of the localities from the SCV Páramo that had a wide niche breadth along this axis (Figure 2). We found that Mean Temperature of Warmest Quarter (Bio 10) was the variable with consistently higher individual influence (mean = 77% ± 11.3) for all modeling techniques applied. Isothermality (Bio 3, mean = 51% ± 19.1) and Precipitation of Driest Quarter (Bio 17, mean = 45% ± 25.4) had also relatively high importance values (Supplementary Figure 5). Topographic variables appeared to be less informative in predicting the species occurrences but showed more relevance for differentiating the different geographical groups within the species (Figure 2 and Supplementary Figure 5).

Predicted Suitable Distribution of P. quadrijuga Under Climate Change
The potential geographic range projected for P. quadrijuga under different climate change scenarios presented strong reductions in suitable areas (Figure 3), and also in areas encompassed by PAs ( Table 2). Our results predicted a substantial decline of the suitable habitat under all climate change scenarios by 2040s, the fragmentation of currently suitable areas, and the loss of most of the suitable habitat patches in the north of the Eastern Cordillera (Páramo districts "Boyacá" and "Santanderes"). In fact, the predicted future distribution of the suitable habitat for P. quadrijuga shows that only the population of SCV, would be able to persist under all the four different RCP scenarios (Figure 3), while for all the other areas, a strong reduction and fragmentation (or even disappearance) of suitable areas is predicted. Moreover, the HSMs predicted a potential loss of suitable area especially in the northern part of the Eastern Cordillera, even for the most optimistic scenario (RCP 2.6 for 2041-2060; Figure 3).
For all RCPs and for both the evaluated time periods, a decrease of more than 80% of the suitable habitat for P. quadrijuga is evident, and underscores a general scenario of dramatic range reduction and loss. The effect of global climate change on the potential distribution of P. quadrijuga is particularly drastic under the most pessimistic scenario (RCP 8.5), leaving only small relictual patches of the current suitable area in the northern species range (SNC and PSB Páramo complexes; Figure 1), for the 2061-2080 period. Interestingly, for the 2041-2060 period, the HSM predicted a higher loss of suitable area under the RCP 4.5 scenario than under all the other RCPs, followed by the RCP 8.5 (Figure 3). The potential loss of climatically suitable areas ranges from 76.4% (RCP 2.6-2041-2060) to 98.5% (RCP 4.5 2041(RCP 4.5 -2060, even for the more climatically stable regions (i.e., SCV Páramo complex; Figure 1).
Additionally, the ensemble projection of the HSMs to future climate conditions predicted a fast contraction of the altitudinal range around medium-high elevation values, together with more fragmentation of suitable areas and even local population extinction in several of the northernmost localities (Figures 3, 4).

Distribution of P. quadrijuga in the Protected Areas Network
Within the predicted current potential distribution of P. quadrijuga (4444.9 km 2 ), we identified about 2036.9 km 2 (45.8%) overlaps SINAP network areas. Moreover, 4168.8 km 2 (93.8% of the current suitable area) overlays Páramo complexes areas ( Table 1). From the 99 presence records, 41 occur inside the protected areas of the SINAP. The Natural Parks Sumapaz and Pan de Azúcar el Consuelo are the protected areas with the highest number of records (16 and 11 occurrence points, respectively; Supplementary Table 2). Among the National Natural Parks (PNN) we found species records in the PNN Chingaza, PNN Sierra Nevada del Cocuy, and PNN Sumapaz. With regard to Regional Natural Parks (PNR), P. quadrijuga occurrences were recorded inside the PNR Pan de Azúcar el Consuelo, PNR Mutiscua Pamplona, and PNR Páramo de Santurbán.
The future projections of P. quadrijuga predicted a strong decline of protected areas covering current highly to moderately suitable habitats. The most dramatic shift of suitable areas was observed for the period 2041-2060, compared to the current distribution, affecting the northernmost portion of the species distribution. Although the species is of high conservation interest, the known occurrences located inside protected areas in Colombia are limited to only 41% of the reported occurrences (Supplementary Table 2).

DISCUSSION
This study reinforces the growing evidence of significant range shifts and drastic reduction of high Andean plant species driven by global climate change. The predicted potential future distribution of the Colombian Eastern Cordillera highland endemic Polylepis quadrijuga, highlights a strong decline in its suitable habitat. Thus, the implementation of management and conservation strategies becomes necessary as a basis for ensuring P. quadrijuga persistence and the preservation of the entire ecosystems that support this species. To this aim, our study identifies the areas of highest vulnerability for P. quadrijuga under future climatic conditions and provides useful information about stable suitable areas, serving as a baseline for decisionmakers to implement effective conservation and management strategies to preserve Polylepis forest in Colombia. Lastly, we underline the importance of Protected Natural Areas in conservation studies and planning in the higher Andes, and their relevance for mountain species whose distribution is often naturally fragmented and more susceptible to habitat loss.

Current Suitable Habitat and Climatic Predictors for the Distribution of P. quadrijuga
We present the most detailed model to date regarding the potential distribution of P. quadrijuga, which shows four main areas with high habitat suitability values (Figure 1). Compared to a previous assessment (Fajardo-Gutiérrez et al., 2018), we included an updated dataset of herbarium records and verified records from citizen science projects and street imagery sources, enhancing the model ability to characterize the species niche (Wisz et al., 2008). In addition, we applied a multiple algorithm ensemble modeling approach using a climatic dataset suited for montane regions, overall increasing the robustness in the identification of the climatic predictors of P. quadrijuga distribution and the spatial prediction of suitable areas.
The HSM predicts large suitable areas for P. quadrijuga in the western slopes of the Eastern Cordillera of Colombia, specifically in the SCV and CHG Páramo complexes, located in the southern range of the species, but lower and more fragmented suitable areas in the northern Páramos (JSB and ALM; Figure 1). In fact, our occurrences dataset shows a northward decrease in records, which may suggest a sampling bias problem (Table 1). However, we interpret this observed pattern of decreasing suitability toward the north as a natural response to the change of available niches driven by three interrelated factors: slope, topographic complexity, and soil water content. SCV and CHG Páramos host larger and topographically less complex valleys, allowing the establishment of bigger populations of P. quadrijuga. At higher latitudes, the suitable climatic belt for the species presents a more complex topography, with steep slopes and lower water availability. This restricts the species establishment to smaller forests on selected-facing slopes that provide a sheltering effect and enable the species to avoid long-term freezing, which limits tree growth, as well as seedling dehydration and photoinhibition caused by strong solar radiation (Toivonen et al., 2018).
Interestingly, the PCA grouped the populations occurring in the Páramos complexes of the larger Eastern Cordillera highland plateau (i.e., CHG and SCV) along similar values of the variables Topographic wetness index, Bio 4 and Bio 6, and separated them from the populations of the northeastern corner of the Eastern Cordillera (i.e., GLR, JSB, PSB, SNC), that were instead more related with variables such as Slope and Bio 3 (PC1, Figure 2), showing a temperature and topographic gradient. The PC2 highlighted a precipitation gradient spanning from the CHG Páramo (with high values of Bio 14, Bio 16, and Bio 19), to PSB, JSB, and SNC, that showed low values of those variables, and that in other words classify as atmospherically dry Páramos. Two slightly different trends were observed regarding the environmental niche of P. quadrijuga, one for the 35 records found in GLR, which forms part of the Boyacá Páramo district, and another for the 42 records of the SCV Páramo complex in the Cundinamarca district (Table 1 and Figure 2), where the variables Bio 2 and Bio 4 play an important explicative role in the PC1. The GLR Páramo complex had a higher temperature diurnal range (Bio 2) compared to SCV, and SCV had a higher temperature seasonality (Bio 4) compared to GLR (Figure 2).
An alternative explanation of the lower amount of occurrences of P. quadrijuga in the northernmost part of the species distribution (Boyacá and Santanderes Páramo districts) is that these areas, where environmental conditions are suitable to host high andean forest and Páramo communities, had undergone severe deforestation through the past decades, thus resulting in the impoverishment of the local flora. In fact, 22-47% of the land belonging to the ALM, GLR, GUE, IGM, JSB, and PSB Páramo complexes is highly transformed, and as of 2011, these hosted complessively 168 km 2 of legal mining titles (Sarmiento et al., 2013). Oppositely, SCV and CHG Páramo complexes are also protected areas of the SINAP, hosting only 10 and 6% of transformed area, respectively, and 3556.4 km 2 of natural ecosystems, mostly Páramos and high andean forests (Sarmiento et al., 2013). Also, the northern east corner of the Eastern  Values indicate the mean (sd) proportion of the stable, loss, or gain area over the projections from the 12 Global Circulation Models compared to the current conditions. Stable suitable habitats (Stable), predicted loss of suitable habitat (Loss) and predicted gained suitable habitat (Gain).
Cordillera displays a remarkably high topographic complexity, due to the presence of numerous high peaks surrounding the deep Chicamocha inter-Andean canyon (Peyre et al., 2019;Valencia et al., 2020). As a consequence, we found that these Páramo complexes (SNC, GLR, ALM, and PSB) host smaller suitable areas for P. quadrijuga, as already pointed out for the Espeletia complex (Valencia et al., 2020).
Polylepis trees are exposed to the harsh climatic conditions that characterize highland areas. Constraints related to water availability, hence precipitation, soil water content during the dry season, and low-temperature extremes define the main ecophysiologic adaptations of Polylepis species (Velez et al., 1998;Morales et al., 2004;Ramos et al., 2013). Consistently, our results showed that temperature-related variables targeting extreme diurnal/annual ranges are among the better climatic predictors of P. quadrijuga occurrence (Supplementary Figure 5). These results are in agreement with previously published studies focusing on other Polylepis species showing that relatively warm diurnal and freezing night temperatures are expected to limit the distribution in Polylepis species (Purcell and Brelsford, 2004;Toivonen et al., 2014;Cuyckens et al., 2016;Zutta and Rundel, 2017). Field observations reports support that P. quadrijuga reproduces sexually inside the forest fragments, while in edge zones, clonal reproduction becomes the primary natural propagation strategy. This change relates to the scarce success of seedling establishment in ecotonal zones, where the open vegetation does not protect nor create favorable conditions for seedling survival.
Several studies predict a higher sensitivity to climate change of highland species and mainly relate it to upward distribution shifts along the altitudinal ranges, thus leading to a contraction of the available suitable area Morueta-Holme et al., 2015). In the case of P. quadrijuga, our results predict significant reductions in the projected area and an altitudinal range shrinkage to the uppermost edge as a possible consequence of climate change (Figures 3, 4). Our projections show that only 11-15% of the suitable habitat of P. quadrijuga will be able to persist under the different 2041-2060 RCP scenarios and that only 8-17% will persist under the 2061-2080 RCP scenarios. These results are compatible with the projection of a severe reduction of P. tarapacana in the South American Altiplano (Cuyckens et al., 2016) and with a climate vulnerability assessment of the Espeletia complex, which shows that several Páramos of the Eastern Cordillera are highly vulnerable to a changing climate (Valencia et al., 2020). Overall, our results agree with the prediction that ongoing climate change could compromise the long−term persistence of the genus Polylepis (Cuyckens et al., 2016;Pinos, 2020).
Warming temperatures and changes in precipitation patterns are causing many high Andean biomes to retreat upslope and suffer changes in their composition as higher-altitude environments undergo faster temperature changes than loweraltitude environments (Costion et al., 2015;Mountain Research Initiative EDW Working Group, 2015). In particular, for some species that will need to track the isotherms throughout the Tropical Andes, a general upslope migration of 600 m is predicted to happen by the year 2100 (Anderson et al., 2011). However, forests' range shift dynamics in response to climate change are a complex process, and a possible Andean treeline upward displacement to areas occupied by Páramos can be limited by many factors, including high radiation levels, freezing events, and human disturbances (Rehm and Feeley, 2015a,b).
Even if P. quadrijuga forests' extent will not decrease, their distribution is unlikely to increase at higher elevations. For example, some Polylepis woodlands from the Cordillera de Vilcanota (Perú) can be highly stable (Jameson and Ramsay, 2007;Feeley, 2013, 2015b). However, since our result predicted an altitudinal range shrink under climate change conditions, the species reproduction needs to be considered because treeline stands may show slow and species-specific responses (Camarero and Gutiérrez, 2004;Dullinger et al., 2004).
Interestingly, Toivonen et al. (2018) found that the topographic heterogeneity of the high-Andes could help P. subsericans forests to persist regionally, facilitating these treeline species by providing dispersal and establishment microsites, achievable by local relocation.

Importance of Protected Areas Network to the Conservation of P. quadrijuga
We observed that PAs covered almost 75% of all remaining P. quadrijuga forests in Colombia. However, in all the examined PAs, the total suitable area for this species tends to be reduced under future climate change scenarios, especially toward the northern Páramos complex (Figure 3). In this context, our results serve as a baseline to adequately protect Polylepis woodlands, while highlighting relevant areas to expand the current PAs network to include suitable areas for P. quadrijuga (e.g., SCV, GLR, SNC, and GUE; Figure 1). The future HSMs projections define some important climatically stable areas where P. quadrijuga could be able to persist. However, these stable areas are located principally inside the SCV and CHG Páramos, where the model predicts a lower decrease in suitable areas. Nonetheless, remaining Polylepis patches might also undergo further fragmentation and occupy sub-optimal conditions at the periphery of their future ranges. A severe concern underlines the areas located at the northeast edge of the Eastern Cordillera, where a dramatic decrease or even a complete loss of suitable habitat is predicted ( Table 2).
Our study reiterated the value and importance of an adequate representation and protection of biological entities in the Páramo complexes network in Colombia. According to our predictions, the SCV area supports stable conditions within the intervened local landscape despite the regional biotic and abiotic change. Interestingly, Zutta and Rundel (2017) predicted refugial areas along the Andes for many Polylepisassociated species indicating that temporally stable areas will be crucial for conservation during anthropogenic-induced climate change, and that their protection and management may reflect the persistence of Polylepis woodlands in the future. Correspondingly, our findings suggest that conservation efforts targeting P. quadrijuga should be directed first toward the large National Natural Parks (Sumapaz, Chingaza, Sierra Nevada del Cocuy; Supplementary Table 2).
According to our future projections, the current PAs network will protect ∼36-48% of the extent of the predicted P. quadrijuga suitable area. Additionally, the few areas that could maintain part of its future suitable habitat (CE-CM-SCV) are currently strongly threatened by human activities. For example, only 142,112 of the 333,420 ha of the Sumapaz Páramo are currently protected as part of the Sumapaz Natural National Park (IAvH, 2012), and fracking projects currently threaten the area. Thus, for P. quadrijuga, habitat protection and connectivity restoration should become a central objective of conservation strategies. In this sense, management plans of already established protected areas such as SCV and CHG Páramos, which our study reports as climatically suitable, should receive special attention to preserve Polylepis woodlands. We recommend additional protection efforts in the north of the species distribution, where areas under protection are currently scarce.
Here we pinpoint that several areas within SCV, CHG, and SNC Páramos could be considered as stable suitable areas for P. quadrijuga, since our models predicted their persistence under some of the evaluated future time frames and for most of the evaluated RCPs. Interestingly, these Páramos coincide with identified high diversity and Pleistocene refuge areas (Cleef, 2008;Zutta and Rundel, 2017). These Páramos were already identified as future persistence areas for different highland plant species (Ramirez-Villegas et al., 2014;Valencia et al., 2020). Thus, we assume that our results and their conservation implications can be of high relevance for protecting the highland Andean biodiversity in the future. However, it is very worrying that Páramo ecosystems have long been at risk due to agriculture, mining and that now they are susceptible to climate change. For example, for the CHG complex, it was found that increasing temperatures and changing precipitation regimes will reduce the 39-52% of the currently suitable area for the Páramo ecosystem during the dry season, and 13-34% during the wet season (Cresso et al., 2020).
We suggest that five of the PAs that showed at least a minimal area of persistence in some of the future scenarios and time frames, should be initially targeted for management: "Parque Nacional Natural Chingaza, " "Parque Nacional Natural Sumapaz, " "Parque Nacional Natural Sierra Nevada del Cocuy, " "Parque Nacional Natural Pisba, " "Pan de Azúcar-El Consuelo, " and "Páramo de Guantivá y La Rusia" regional parks. Prioritizing and expanding small protected areas in the correspondent mountain ranges and PCCs will be essential for protecting this species in situ and ex situ, and pivotal to safeguard it from further habitat loss. Interestingly, our model projections to future scenarios highlighted several areas of climatic stability or even suitable habitat expansion in many protected areas that currently host no records of P. quadrijuga (e.g., Pisba and Páramo de las Oseras; Supplementary Table 2), or small private reserves (e.g., La Reserva, Villarica; Supplementary Table 2). These can be effective for ex situ conservation measures for the scattered mature P. quadrijuga individuals and seedlings, and can be particularly relevant to carry out seed germination or seedling transplant experiments in the field.

Implications and Recommendations for the Conservation of Polylepis Woodlands Under Climate Change
Currently, few sustainable management strategies have been proposed for Polylepis woodlands, and where implemented, they do not match the widespread management practices in Colombian rural areas (Gareca et al., 2010;Vásquez et al., 2014;Boza Espinoza et al., 2019). Since many Polylepis woodlands are located in densely populated areas and are under considerable pressure from human disturbance, it is urgent to identify hotspots of vulnerability and resilience and define management interventions. We also emphasize the need to maximize protection efforts targeted around treeline ecotones biodiversity in the Andes, as these are essential providers of ecosystem services and home to a diverse array of unique wildlife (Fjeldså, 1993;Arnal et al., 2014;Sevillano-Ríos et al., 2018).
We consider P. quadrijuga as an appropriate umbrella species to monitor the upward shift of treeline Andean ecosystems. However, it is essential to gather and consider more information about its interactions with other species for the selection of appropriate sites for protection. For example, high habitat specialization with strong confinement inside Polylepis forests was observed for some specialist bird species (Fjeldså, 1993;Cahill and Matthysen, 2007;Gareca et al., 2010;Meneses and Herrera, 2013). In agreement with Sevillano-Ríos and Rodewald (2017), we suggest protecting large remnants at lower elevations and maintaining all relictual patches of Polylepis forests irrespective of size become two key strategies constituting the cornerstones of future conservation efforts. Also, protected areas in the higher Andes must incorporate the study of potential species distributions and species responses to local climate changes.
Although it is known that the presence of a species inside PAs offers no guarantees for effective species long-term conservation (Hoffmann et al., 2019), protection efforts are the most successful tools to maximize landscape connectivity of Polylepis forest. Simultaneously, identifying areas of historical and potential future refugia or climatically stable areas is important to direct the focus of conservation efforts (Keppel et al., 2015). Moreover, studies predicting species' future distribution considering the idiosyncratic species responses (i.e., phenotypic acclimation, rapid evolution, ecophysiological studies of selected species, and local adaptation patterns) are necessary to set the methodological baseline for the management of Polylepis woodlands inside PAs. Some examples include Ramos et al. (2013), showing that the physiological performance of P. quadrijuga was significantly affected by landscape fragmentation, and Velez et al. (1998) supporting that P. quadrijuga trees develop a high protection strategy against UV-B radiation, which tends to increase with altitude. However, more detailed studies are still needed to clarify relevant aspects of species biology and deliver comprehensive eco−evolutionary data in the context of global climate change. Toivonen et al. (2018) suggested that a better understanding of the environmental niche preferences of Polylepis species is crucial to identifying potential climatic resilience areas for the Polylepis forests. Accordingly, we show that HSM could be crucial for possible conservation strategies to address changing environmental conditions. However, the incorporation of species−specific dispersal constraints and the influence of local topographic and edaphic factors are essential additions to projections of plant species distribution models in the future (Rull and Vegas-Vilarrúbia, 2017). Moreover, since treeline species responses to climate changes are species-specific, climate, landuse, fragmentation, and soil degradation are all critical factors determining whether an upwards migration is possible. For this reason, intensive fieldwork might be relevant to monitor plant altitudinal migration and could increase the knowledge of the potential responses of threatened treeline species in the higher Andes.
Under climate change conditions, plants face various barriers to their dispersal that may hamper gene-flow and recruitment, leading to extinction. Although Polylepis species are windpollinated, palynological and population genetics analyses show that pollen dispersal is limited to relatively short distances (Schmidt-Lebuhn et al., 2007), underlining that forest fragmentation can lead to an increase of genetic drift, inbreeding, and a consequent loss of adaptive capacity in local populations. On average, for the northern Andean region, is expected an increase in annual precipitation by values as high as 300 mm/year (Buytaert et al., 2011). This issue is crucial for biological processes such as Polylepis species dispersion. Also, the dispersion ability of tropical treeline species could be limited by the complex topography of the Andes and the harsh microclimatic conditions near the soil surface, which may be one of the reasons limiting seedlings recruitment . We still have much to learn about the genetic mechanisms behind climate change adaptations of high-Andean treeline species, and the recognition of regional patterns of genetic diversity of Polylepis populations could offer the opportunity not only to preserve ecosystems with high biological diversity but also to understand and protect evolutionary processes from climate change (i.e., Schierenbeck, 2017).
In this work, we incorporated data deriving from citizen science projects to increase the number of occurrences of P. quadrijuga. We observed that integrating data from alternative databases could substantially increase sample size and could open up new opportunities to identify declining populations earlier to evaluate the effectiveness of policy interventions. We remark that citizen science has the potential to generate valuable biological observations suitable for species distribution modeling (Tiago et al., 2017;Roy-Dufresne et al., 2019;Petrovan et al., 2020) and ecological restoration projects (Edwards et al., 2018;Callaghan et al., 2019). Citizen science projects also provide the users with informal learning experiences that can contribute to the knowledge on ecological traits, seed collection, plant germination, and propagation of species of the higher Andes, which are often neglected in the scientific literature. However, our results should be interpreted with caution, given that the distribution of utilized species records is heavily biased toward the PAs and the central part of Colombian Eastern Cordillera. Probably, due also to the fact that patches of Polylepis are challenging to identify using current air photographs and satellite images (Aucca and Ramsay, 2005), most of our occurrences of P. quadrijuga are located in PAs situated and along roads, whereas large areas in the northeast and southeast of the Eastern Cordillera remain poorly sampled or not sampled at all.
Our model is subjected to limitations, since we did not account for physiological processes or any information about the species' ecology and natural history (e.g., interspecific competition, migration rate), and recent human activity and disturbance (e.g., land use change, fires). Moreover, we could not incorporate data on microclimatic conditions into the modeling, due to the lack of available high resolution data for our study area. Nevertheless, an improved understanding of P. quadrijuga distribution and the possible future response to climate change represents a precondition for more precise modeling efforts of potential range shifts of treeline species in the Andes. Our results provide a better comprehension of the role of environmental constraints in shaping treeline species potential ranges. Specifically, we highlight a "Change in the climatically suitable area" and a "Change in average climatic suitability in already occupied cells" for P. quadrijuga, which have been defined as two quantifying components of risk for woody species under climate change (Ohlemüller et al., 2006). For this reason, based on our results, we suggest implementing an evaluation criterion to categorize the conservation status of P. quadrijuga, targeting a key treeline forming species on the northern Andes, a region that is still being heavily transformed and fragmented by agricultural and mining activities. Immediate implementation of in situ conservation programs could be crucial to ensure the persistence of the species. Monitoring the upward distributional shift would help us to detect the predicted decline in Polylepis populations over time. Lastly, the creation of a germplasm bank from specimens from all the different localities would be fundamental.

CONCLUSION
According to our future projections, the threat status of P. quadrijuga is alarming and needs attention due to the predicted considerable reduction in climatically suitable habitat for the species. The detailed information provided here could represent a valuable tool to guide the future establishment of new and efficient conservation efforts. Moreover, our findings can help restore Polylepis forests and preserve them through in situ and ex situ strategies, planned according to the predicted future stable climatic areas. Medium-and long-term conservation strategies for P. quadrijuga woodlands should be focused on actions to prevent deforestation of the extant patches and to implement ecological restoration measures in the Cruz Verde-Sumapaz and Páramo complex, which is the only predicted future suitable area for six out of eight of the future scenarios employed here. Our results can be used to define high-priority areas in the Andean region for conservation and management plans for Polylepis forest and can offer a baseline for similar analyses of other endangered and threatened treeline species in the region. Finally, we reinforce the need to strengthen management strategies inside those protected areas that comprise the remaining patches of P. quadrijuga forest in Colombia.

DATA AVAILABILITY STATEMENT
The original occurrence data used in the study is included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
GS-A and LC-V designed the study and conducted the formal analysis. FF-G and MC performed data curation. GS-A, FF-G, and MC prepared the figures, tables, and Supplementary Data. LC-V and MC wrote the first draft of the manuscript. All authors interpreted the results, contributed substantially to revisions and worked together to define the structure and content of the manuscript, and gave final approval for publication.