Assessing the impact of climate change on threatened endemic vascular plants of Argentina

Biodiversity loss due to climate change is one of the most critical current environmental problems. Global warming is causing substantial species-range shifts and local extirpations, especially for species with restricted distribution ranges. Studies of impact of the climate change on species ranges and environmental suitability have become a fundamental tool for evaluating conservation strategies. However, one important limitation of these approaches is that only species with an adequate number of spatially distinct occurrence records can be modelled, generally excluding threatened rare species from the analyses, a situation referred to as the ‘rare species modelling paradox’. To overcome this limitation, we analysed the effect of climate change on the richness of threatened endemic plants of Argentina employing a macroecological modelling approach, using three different modelling techniques (generalized linear mixed models – GLMM, generalized additive models – GAM, and boosted regression trees – BRT), four general circulation models, two representative CO2 concentration pathways (RCPs), and two time periods (2050 and 2070). We identified grid cells with the greatest decline in numbers of threatened endemics, determined species composition in these cells and characterized their vulnerability using three indices. A loss of species richness is observed in ca 83% of the cells, and both protected areas and hotspots of threatened species show significant decrease in future species richness. We identified 32 most affected cells under future climatic projections, including a total of 370 threatened endemic species and exhibiting high beta diversity values (high dissimilarity) among most of the cells. Cells with the highest vulnerability were located along the Central Andes of northwestern Argentina, along the Southern Andean Yungas, High Monte and Central Andean Puna ecoregions, and including a total of 118 threatened endemics (15% of those registered for Argentina) with greater representation of Asteraceae, Apocynaceae, Cactaceae and Iridaceae. However, coverage of protected area network is less than 5% for each of these cells. Our results highlight the urgent need for both in situ and ex situ conservation policies and strategies for the vascular flora of Argentina.


Introduction
Anthropogenic climate change has accelerated and become more intense in recent years. The IPCC's sixth assessment report (AR6) on the impacts of global warming on the earth indicates an increase of ca 1.5°C above pre-industrial levels . In addition, global temperatures will continue to increase until at least mid-century under all emissions scenarios considered, requiring immediate, rapid and largescale reduction in greenhouse gas emissions to limit global warming to 1.5-2°C by 2081-2100 (IPCC 2021). Climate change, changes in land use and direct exploitation of organisms have been reported as the three main sources of threat to biodiversity (IPBES 2019). Likewise, biodiversity loss is one of the most critical current environmental problems (Dirzo and Raven 2003;Sharrock 2012;Ceballos et al. 2015) on earth. Global warming and habitat destruction are causing significant species-range shifts, local contractions and extirpations, which can be important drivers of biodiversity loss (Chen et al. 2011;Urban 2015), especially for species with restricted ranges that are exposed to high levels of vulnerability (Gaston 1994;Urban 2015). Due to climate change and human activities, species extinction rates are currently 1,000 times higher than extinctions due to natural causes, and are expected to increase ten-fold during this century as a result of the accumulation of human impacts on natural ecosystems (De Vos et al. 2015).
One of the main impacts of climate change on biodiversity is the expected change in the geographical distribution of species (Chen et al. 2011;Garcia et al. 2014). Increased temperatures associated to global warming would be expected to cause poleward or upward shifts of species distribution ranges (Parmesan 2019), and about half of species for which long-term data exist (e.g. mountain plants, trees, birds, butterflies, fish), exhibited significant changes in their distributions over the past 20 to 140 years (Parmesan and Yohe 2003;Root et al. 2003;Rosenzweig et al. 2008). In this way, current high rates of climate change are causing range loss and displacement in many species (Lenoir and Svenning 2015;Pecl et al. 2017;Steinbauer et al. 2018), with an observed lag in population dynamic responses to climate trends (Dullinger et al. 2012). When analysing responses of biodiversity to recent climate change, it was found that restricted-range species are the first group to respond with the highest observed extinction rates (Parmesan 2006). The size of the occupied area [e.g. extent of occurrence (EOO) and area of occupancy (AOO) from IUCN assessments - IUCN 2012IUCN , 2019 are among the most relevant variables associated with high extinction risk (Pearson et al. 2014). In addition, rare and narrowly distributed species often have unique trait combinations, which underlines the importance of these species for conservation (Mouillot et al. 2013). Therefore, endemic and rare species with restricted ranges and small population sizes are highly sensitive to the current threats (Lavergne et al. 2005;Van Calster et al. 2008), and require an assessment of their vulnerability for the development of conservation strategies because their extirpation contributes disproportionately to the current extinction crisis.
Approaches using species distribution models (SDMs) are widely used to study the impact of climate change on species ranges and environmental suitability (Guisan et al. 2013). However, an important limitation of SDMs is that only species with an adequate number of spatially distinct occurrence records can be modelled (Stockwell and Peterson 2002;van Proosdij et al. 2016), which is uncommon in the case of threatened endemic taxa (Platts et al. 2014). This implies that although rare and narrowranging species are the most in need of models that support their conservation, in most cases they are excluded from the analyses due to the low number of occurrences, a situation referred to as the 'rare species modelling paradox' (Lomba et al. 2010). Several approaches have been proposed to overcome this limitation, as the ensemble of small models (ESM; Lomba et al. 2010, Breiner et al. 2015, joint species distribution modelling (JSDM; Pollock et al. 2014;Ovaskainen et al. 2017), climate-niche factor analysis (cENFA; Rinnan and Lawler 2019), and hybrid mechanistic or population models, which can incorporate environmental requirements, but also the dispersion of species and local interactions (Zurell et al. 2016;Soriano-Redondo et al. 2019). However, a minimum of 5-10 different occurrences (separated for example by at least 1 km) is still necessary to account for cell resolution and spatial autocorrelation (van Proosdij et al. 2016). In the case of threatened endemic species, this condition is often not achieved, as many species are only known from the type collection or two to three specimen records.
In Argentina, of the 1,683 endemic vascular plants (February 2021), around one-third of species are represented by less than five specimens (529, 36%), most of them preliminarily categorized as threatened (66% of the 800 threatened endemics; Salariato et al. 2021). Most threatened plants in Argentina grow in valleys and highlands of the Central and Southern Andes ( Fig. 1; delimitation of the Andes sensu Luebert and Weigend 2014), along the Central Andean Puna, Central Andean Dry Puna, Southern Andean Yungas, High Monte and Southern Andean steppe ecoregions (Olson et al. 2001;Fig. S1 in the Electronic supplementary material). These mountain ecosystems are particularly sensitive to climate change (Hughes 2000;Dullinger et al. 2012;Guisan et al. 2019), since mountainous regions appear to be warming faster than other regions (Beniston et al. 1997;Nogués-Bravo et al. 2007;Guisan et al. 2019), including the Andes of the southern cone (Chevallier et al. 2011;Barros et al. 2000Barros et al. , 2015Rivera et al. 2017). Thus, predictions of future climate change for mountain ranges around the world show an average temperature change of 2-3°C by 2070 and 3-5°C by the end of the century (Nogués-Bravo et al. 2007). Although this lack of information on the geographic distribution of rare / threatened plants (i.e. the Wallacean shortfall - Whittaker et al. 2005;Bini et al. 2006;Oliveira et al. 2016) should be mitigated by intensive field explorations, the risk of extinction due to climate change, coupled with habitat loss and degradation, represent an imminent major threat to many species. Therefore, it is of vital importance to promptly generate preliminary studies to identify vulnerable species and areas, in order to guide policies and conservation efforts (Root and Schneider 2006;Stanton et al. 2015).
An option for analysing the impact of climate change on rare threatened species is to use macroecological modelling (MEM; 'assemble first, predict later'; Ferrier and Guisan 2006;Guisan and Rahbek 2011;Dubuis et al. 2011). The MEM approach models species richness (SR), directly relating this variable (the number of species within a geographic unit) to values of environmental variables that characterize the same unit (Ferrier and Guisan 2006;Dubuis et al. 2011). The main limitation of the MEM approach remains that it does not provide information on the identity of the species occurring at a given location and thus cannot predict changes in species composition, unlike other methods such as S-SDM, SESAM (Guisan and Rahbek 2011), or J-SDM (Ovaskainen et al. 2017). However, MEM models can accurately predict diversity after accounting for the effects of spatial autocorrelation (Algar et al. 2009;Dubuis et al. 2011;D'Amen et al. 2018) and, allows us to incorporate rare species with few occurrences, as is the case for most threatened species. In this way, general patterns of threat can be detected, and areas sensitive to species loss (in terms of species richness reduction) due to climate change can be identified (Newbold et al. 2009;Biber et al. 2020).
In this work, we analyse the effect of climate change on the richness of threatened endemic plants of Argentina. In particular, using macroecological modelling with different modelling techniques, general circulation models (GCMs), representative concentration pathways (RCPs) and time periods, we first address the following major questions: (1) Is the area in Argentina that experience a loss of species under future climate projections greater than the area that do not lose species or that increase their count? and (2) Is there a significant decrease in the number of species within protected areas and hotspots of threatened species (cells with large numbers of threatened endemic species) in future models? In a second step, we characterized the areas in Argentina with the greatest decrease in richness of threatened endemics by their ecoregions and habitats, the species they include (according to our occurrence data), and by three complementary indices that consider the loss of species due to climate change (I a ), the threat categories of the species (I b ) and the level of human pressure (I c ). For these cells, we also estimated the percentage of their areas included in the protected areas (PA) network. The list of species registered in these most vulnerable cells is also supplied for application in both in situ and ex situ conservation projects. This work seeks to contribute to the knowledge of the general patterns of threat for Argentinean vascular flora and to be useful as a proxy to determine future conservation priorities in the country.

Species and climate data
For our study we included all accepted Argentinean vascular plant species (Tracheophyta) endemic to   Salariato et al. (2021). Although categorization under threat does not directly imply the status of rare species and/ or restricted distribution, 66% of the 800 threatened endemics are represented by less than five specimens (529 spp.; Fig. S2), which implies that their distribution cannot be adequately estimated by species distribution modelling. We included 3,354 occurrence points from mainland Argentina (i.e. excluding South Georgia, South Sandwich and the Malvinas) representing 785 threatened endemic species, while fifteen species were excluded because of the lack of adequate data for accurate georeferencing), consulting for this purpose published floristic treatments (see Zuloaga and Belgrano 2015;Zuloaga et al. 2019), the Argentina online flora and the Documenta Florae Australis database (DFA 2017) (http:// www. darwin. edu. ar/ iris). The DFA database includes records from more than half a million specimens deposited in national and international herbaria, representing the most up-todate and comprehensive database of the Argentinean Flora (Zuloaga and Belgrano 2015;Zuloaga et al. 2019).
The geographic coverage of our dataset for threatened endemics is particularly sparse in some regions of Argentina, partly due to collection biases (Fig. S1c). Thus, species may incorrectly appear to be absent from a given location simply because that location has not been sampled (Ferrier 2002). Although analyses at fine scale may have the advantage of estimating patterns of species richness at higher resolution and may reveal new associations between richness and different environmental predictors (since both climatic variables and species richness can vary considerably within a small area), to reduce the problem of false-absences we used coarse grid-cells that are large enough so that each cell has a high probability of containing enough data (Linder 2001;Newbold et al. 2009). We conducted preliminary assessments using two cell sizes: 0.5 arc-degrees (approx. 55 × 55 km) and 0.25 arc-degrees (approx. 27 × 27 km). We compared the performance of these two resolutions on models with the selected climatic variables (see below) by estimating the C-index (Briscoe et al. 2021) through ten cross-validations using the package 'Hmisc' v. 4.6.0 (Harrell 2021) implemented in R v. 4.0.4 (R Core Team 2021). Since models with 0.5 arc-degrees cells showed higher values of the C-index (Fig. S3), this cell size was selected for all subsequent analyses. Occurrence data are included in the Electronic supplementary material (file ESM_1.txt) and also deposited in the Figshare data repository (https:// doi. org/ 10. 6084/ m9. figsh are. 15180 195).
As current environmental predictors we used data on current climatic conditions  included in the nineteen bioclimatic variables from the CHELSA v. 1.2 climatic dataset (Karger et al. 2017; https:// doi. org/ 10. 5061/ dryad. kd1d4). All layers were clipped to the Argentina shapefile and resampled from 30 arc second resolution to 0.5 arc degrees using bilinear interpolation. For future climate data we used GCM compareR (Fajardo et al. 2020) to select general circulation models (GCMs) that represent all possible climate change scenarios for the study area. With this tool, estimating difference between values for pixels in each GCM and the average among the projections from all GCMs selected, we selected from the range of CMIP5 (Coupled Model Intercomparison Project Phase 5) climate model projections that represent possible alternative future conditions. We selected four GCMs (BCC-CSM1-1, IPSL-CM5A-LR, HadGEM2-AO and MRI-CGCM3), two representative concentration pathways (RCPs) of a moderate (RCP 4.5) and a severe (RCP 8.5) greenhouse gas emission scenario, and two time periods: 2050 (average for 2041-2060) and 2070 (average for 2061-2080). Projections were downloaded from CHELSA at 30 arc second resolution and resampled to 0.5 arc degrees with the same method as for current climate data. Data manipulation and visualization were done using the R packages 'raster' v. 3.4.5 (Hijmans 2020), 'rasterVis' v. 0.50 (Perpiñán and Hijmans 2021) 'sp' v. 1.4.5 (Bivand et al. 2013) and 'maptools' v. 1.0.2 (Pebesma and Bivand 2005;Bivand and Lewin-Koh 2020).

Macrocological modelling approach
To study the impact of climate change on the richness of threatened endemic species, we conducted MEMs with species richness (SR) of threatened endemics as a response variable and values of climatic variables that characterize the same geographic unit as predictor variables (Ferrier and Guisan 2006;Dubuis et al. 2011;D'Amen et al. 2018). We selected three different modelling techniques: generalized linear mixed models (GLMM), generalized additive models (GAM) and boosted regression trees (BRT), representing different approaches to modelling species richness (parametric, semi-parametric and machine learning, respectively), allowing us to evaluate uncertainties arising from the choice of methods when combining them in ensemble models (Araújo and New 2007). We first performed a preliminary variable selection using the nineteen bioclimatic variables from CHELSA and evaluating univariate variable importance with the AIC. We fitted generalized linear models (GLM) separately for each predictor to rank the variables according to their importance, and only retaining six best-ranked variables with Spearman's correlation < 0.70 (Dormann 2013) with the R package 'corrplot' v. 0.84 (Wei and Simko 2017;Fig. S4 in the Electronic supplementary material). Then, we evaluated the importance of the six pre-selected variables (bio2, bio3, bio5, bio9, bio13, bio15) in a full GLM model using bidirectional stepwise AIC model selection and considering linear and quadratic terms for each predictor (using the 'poly' function in R). Because fitting too many independent variables in GLMs and GLMMs can lead to overfitting and the selection of meaningless variables in the final model (Wintle et al. 2005), we decided to retain only three bioclimatic variables that, when combined, contributed the most (in terms of AIC values) to the model for our response variable (richness). Subsequently, we verified that the variance inflation factor (VIF) for these three selected predictors was < 10 with the R package 'usdm' v. 1.1.18 (Naimi et al. 2014). The selected variables were maximum temperature in the warmest month (BIO5), precipitation of the wettest month (BIO13) and precipitation seasonality (BIO15; Fig. S4), and besides they proved to be important predictors of the current richness of threatened endemic species in Argentina, they also reflect the boundaries of the climatic niche (Grinnellian niche) of the species for both temperature and precipitation, thus quantifying more directly the impact of climate change on species distribution (La Sorte and Jetz 2012). Because spatial autocorrelation was detected in residuals of the selected model (global Moran's I = 0.053, P < 0.001) with the R packages 'ape' v. 5.4.1 (Paradis and Schliep 2019) and 'pgirmess' v. 1.7.0 (Giraudoux 2021), we incorporated a spatial random effects including the latitude and longitude of the centre of each sampled cell. Generalized linear mixed models were conducted with the R package 'spaMM' v. 3.7.34 (Rousset and Ferdy 2014) with Poisson error, link log, linear and quadratic terms for predictors previously selected by bidirectional stepwise AIC (BIO5: linear and quadratic terms, BIO13: only linear term, BIO15: only quadratic term), and a spatially correlated random effect assuming a Matérn covariance function for latitude and longitude. Generalized additive models were performed with the 'mgcv' package v. 1.8.33 (Wood 2011(Wood , 2017 with Poisson error, link log, thin plate regression splines for BIO5, BIO13, BIO15, and a tensor product of a spatial smoother for latitude and longitude. Boosted regression trees were conducted with the climatic variables (BIO5 + BIO13 + BIO15) and the spatial structure (lat + long) in the R packages 'gbm' v. 2.1.8 (Greenwell et al. 2020) and 'dismo' v. 1.3.3 (Hijmans et al. 2020), using a Poisson distribution, tree complexity = 2, maximum number of trees = 50,000, learning rate = 0.001 and bag fraction = 0.5. In all analyses, cells with richness of zero were excluded for the model fit (leaving 434 cells with richness values > 0), since cells with a recorded richness of zero may occur simply because they have not been sampled and results could be biased by the inclusion of false zero values (Newbold et al. 2009). Performance of the final models was evaluated using tenfold cross-validation and calculating for each test set the root-mean-square error (RMSE, Potts and Elith 2006) between fitted vs observed values, and comparing the RMSE values between modelling techniques (GLMM, GAM, and BRT) using the Kruskal-Wallis test. Additionally, we also evaluate the performance of the different modelling techniques calculating for each test set the intercept (ideally close to zero), the slope (ideally close to one) and Pearson's correlation coefficient.
The future richness of threatened endemics was then estimated by projecting fitted models from the three techniques (GLMM, GAM, BRT) to the four different global circulation models (BCC-CSM1-1, IPSL-CM5A-LR, HadGEM2-AO, and MRI-CGCM3) for each representative concentration pathway (RCPs 4.5 and 8.5) and time period (2050 and 2070). First, we evaluated the environmental analogy between current and future climate data using the multivariate environmental similarity surface (MESS) analysis (Elith et al. 2010) implemented in the R package 'dismo'. To obtain final MEM projection for each RCP and time period, we combined individual models into a final ensemble model (Araújo and New 2007), first by weighted-averaging the individual models of the different techniques but the same GCM proportionally to their RMSE scores in a similar way to the Akaike weights (Wagenmakers and Farrell 2004): and then averaging (unweighted average) the assembled models corresponding to different GMCs.

Climate change impact
To assess the impact of climate change on the richness of threatened endemics, we first calculated the percentage of cells in the ensemble future projections that gained or lost richness compared to the estimated current values. Next, we analysed how climate change affects the richness of threatened endemics: (1) within the network of protected areas (PAs) of Argentina, (2) on hotspots of threatened endemics. To study the impact of climate change on the richness within PAs, shapefiles of protected areas of Argentina were downloaded from the World Database of Protected Areas (UNEP-WCMC and IUCN 2021, www. prote ctedp lanet. net) considering all the IUCN categories of PAs available within Argentina. Then, we used Wilcoxon signed-rank test and the Kolmogorov-Smirnov test in R to assess for significant differences between current and projected future richness of cells within PAs. Alternatively, to explore the effect of climate change on the richness of the hotspots, we first defined hotspots of threatened species as grid cells with threatened species richness at or above the 95th percentile. We used 95% as the threshold because previous analyses have shown that the richest 95-97.5% cells can effectively represent a substantial proportion of species (Prendergast et al. 1993, Myers et al. 2000, Orme et al. 2005. Then, we registered percentage of cells within hotspots that gained and lost richness in ensemble future projections, and compared richness of hotspots for present and projected future using the Wilcoxon signed-rank and the Kolmogorov-Smirnov tests. In a second step, we explored the effect of climate change identifying cells with the greatest richness decrease. For this purpose, we delimited for each of the four assembled models (RCPs 4.5 and 8.5 for 2050 and 2070) cells with a decrease in richness at or below the 2.5th percentile, and then we retained only those cells recovered in at least three of the four models (75% consensus cells). For these consensus cells, we first identified ecoregions of Argentina (sensu Olson et al. 2001) using the shapefile available at https:// www. world wildl ife. org/ publi catio ns/ terre strial-ecore gions-of-the-world, and registered the elevation (mean and maximum) and the terrain ruggedness (difference between the maximum and the minimum value of a cell and its surrounding cells) from the GTOPO30 digital elevation model (USGS 2021; https:// www. usgs. gov/ cente rs/ eros/ scien ce/ usgs-erosarchi ve-digit al-eleva tion-global-30-arc-second-eleva tion-gtopo 30) at 30 arc-second resolution. Then, we studied the presence and composition of threatened species in these cells using the specimen occurrence data. For the 32 selected cells we determined the level of similarity in terms of species composition, estimating beta diversity among cells with the R package 'betapart' v. 1.5.4 (Baselga et al. 2021), using the Sørensen dissimilarity index (βsor) which includes both dissimilarity by turnover (βsim) (species replacement between cells) and nestedness (species loss from cell to cell) (βsor = βsim + βsim ) (Baselga 2012). Next, we explored vulnerability of these cells using three complementary indices. The first index (I a ) considers the loss of species due to climate change in each cell and was estimated as the cumulative decrease in the richness of threatened species over the four future models (∆richness 4.5-2050 + ∆richness 4.5-2070 + ∆richness 8.5-2050 + ∆richness 8.5-2070 ). The second index (I b ) quantifies the level of vulnerability of each cell in relation to the number and threat categories and was calculated as the number of threatened endemic species observed in each cell weighted by their IUCN category (evaluation from Salariato et al. 2021), with VU = 1, EN = 2, CR = 3; list of threatened endemics with their categories are supplied in the electronic supplementary material as file ESM_2.xlsx. The third index (I c ) measures the level of human pressure for each cell using data from the global map of the cumulative human footprint on the environment for 2009 (Venter et al. 2016(Venter et al. , 2018 https:// sedac. ciesin. colum bia. edu/ data/ set/ wilda reas-v3-2009-human-footp rint), downloaded at 1 km and resampled to our grid resolution (0.5 arc-degrees) using bilinear interpolation. Each index was rescaled so that the maximum value was 1 (dividing each cell value by the maximum value obtained among the 32 cells) and used to calculate a fourth index: which represents the global risk for threatened endemics in each cell considering these three sources of threat (climate change, IUCN threat categories, and human activity), where values close to 1 indicate higher risk. For these cells, the percentage of their areas included within the PA network was also estimated using the R package 'sf' 1.0.4 (Pebesma 2018).

Results
Macroecological models for richness of threatened endemics obtained with GLMM showed the lowest RMSE values (Fig. 2, Fig. S5 and Table S1 in the Electronic supplementary material). However, the Kruskall-Wallis test did not show significant differences for the RMSE of the different modelling techniques (P = 0.533). All models are consistent in recovering the highest current richness of threatened endemics along the Central Andes and the northern portion of the Southern Andes (Fig. 2, Fig. S6 in the Electronic supplementary material), with secondary richness areas located at 'Sierras Pampeanas' mountain system in central Argentina and forests at the 'Sierra de Itacuara' and the southernmost portion of the 'Sierra Central', in northeastern Argentina. Alternatively, the 'Ventania' mountain system also appears as a secondary richness zone mainly under the GLMM and ensemble models (Fig. 2, S6). Comparisons between current and future climate data using MESS analysis revealed higher climatic analogy primarily for regions located in southwestern Argentina, while novel climatic conditions were identified mainly in northern, northeastern and central Argentina, with the greatest dissimilarity compared to current conditions for projections under RCP 8.5 (Fig. S7). All projections using the individual GCMs show a main trend of richness loss, less severe with the MRI-CGCM3 model and more severe with the IPSL-CM5ALR model (Table S2, Fig. S8).
Overall, the decline in richness is accentuated for RCP 8.5, becoming the most severe for the future scenario RCP 8.5-2070 (Fig. 3, Table 1). Based on ensemble models, a loss of richness of threatened endemics is observed for most cells (approx. 83%, Table 1, Fig. 3) under both RCPs (4.5-8.5) and time periods . Similarly, approximately 80% of the cells within the PA network exhibit a decrease in future species richness (Table 2), and both Kolmogorov-Smirnov and Wilcoxon signedrank tests are significant (P < 0.001 for all future scenarios; Fig. 3, S9). Hotspots of threatened endemics (Fig. 2, 3) included cells with six or more threatened endemics (median = 9 spp., percentiles 2.5-97.5% = 6-39 spp., max = 47 spp.). Future richness projections show that between 79% and 85% of hotspot cells decrease their richness of threatened endemics (Table 2, Figs. 3, S9), while only 15-21% increase it. The largest decreases are observed in hotspot located along the Central Andes at the transition zone between the Puna, Prepuna and Yungas biogeographical provinces (Fig. 3). In addition, both Kolmogorov-Smirnov and Wilcoxon signed-rank tests show a significant decrease in future richness for hotspot cells of threatened endemics (P < 0.001 for all future scenarios; Figs. 3, S9).
Selecting cells with the greatest decrease in richness (≤ 2.5% percentile) in at least three of the four assembled climate models (RCPs 4.5 and 8.5 for 2050 and 2070), we identified 32 consensus cells most affected in future projections grouped into nine main clumps (Fig. 4a). These cells included 370 species (47% of all 785 threatened endemic species analysed here) mainly from Asteraceae (84 spp.), Fig. 2 Current species richness (SR) for threatened endemic species in Argentina and estimates from different modelling techniques: a -empirical (observed) species richness (SR) for threatened endemic species used to fit the models (response variable); b -boxplots (median, 25%-75% percentiles, and min-max) for RMSE values obtained using ten-fold cross validation and different modelling techniques: generalized linear mixed model (glmm), generalized additive model (gam), boosted regression trees (brt), and an ensemble model (ensemble) using weighted-average over the three above techniques; c -predicted SR using an ensemble model of weighted-average over the different modelling techniques (glmm, gam, brt); d -weighted standard deviation of SR displaying dispersion of richness for each cell associated with the ensemble model; and e -coefficient of variation (CV) for SR and each cell associated with the ensemble model.  Fig. S1). The βsor index recovered high values of beta diversity (high dissimilarity) among most cells (Fig. 4c). Values of I a (loss of species due to climate change), I b (IUCN threat categories for species) and I c (human footprint) indices for consensus cells showed correlation values < 0.7, with I a and I b being the most correlated (Fig. S10). When these three indices were combined to calculate the I d index (which represents the risk of each consensus cell considering these three sources of threat), cells with the higher values were located along the Central Andes of northwestern Argentina ( Fig. 5; Table S3) in the Southern Andean Yungas, High Monte, and Central Andean Puna ecoregions. Cells 3, 10 and 12 (cell numbers are shown in Fig. 4a), which exhibited the highest I d values, included 47 (5.9% of total threatened endemics), 54 (6.9%) and 44 (5.6%) species, respectively. Alternatively, other cells with high I d values not included along the Central Andes were recovered for the Alto Paraná Atlantic forest ecoregion of eastern Argentina (cell 18; Fig. 5, Table S3) and the 'Sierras Pampeanas' mountain system along the Dry Chaco ecoregion at central Argentina (cell 23), exhibiting also these cells high values for I c (human footprint). Likewise, cells with higher values of I d also present high values for elevation/ roughness (Fig. S11), with the exception of the cell 18 located in the Alto Paraná Atlantic forest. The   4 (−15.5, 4.5) relationship between elevation / roughness was stronger for index I b , showing the positive association of these topographic variables with the number of threatened species (Spearman's ρ I b ~ roughness = 0.48; Fig. 4b, S11).
Finally, we recovered that twelve most affected consensus cells were not associated with any PA (0%), and only eight of the 32 consensus cells had coverage greater than 5% (Fig. 6, Table S3). In particular, the three cells with the highest I d values (cells 3, 10, and 12) exhibited a PA coverage lower than 5% (1.6%, 0%, and 4.7%, respectively; Fig. 6, Table S3).

Discussion
Natural habitats are changing at an accelerating rate (Dirzo and Raven 2003;Chen et al. 2011;Urban 2015;De Vos et al. 2015) and, given that resources to address biodiversity loss are limited, it is necessary to identify priority sites for conservation. The inclusion of rare and threatened species in modelling studies has traditionally been complicated by the scarcity of distribution data Lomba et al. 2010;Platts et al. 2014). The low representativeness and geographic precision of rare and threatened species imposes important limitations on the study of the  c Fig. 4 Most affected cells under future climatic projections: a -Cells with the greatest decrease in richness (≤ 2.5% percentile) in at least three of the four assembled climate models (RCPs 4.5 and 8.5 for 2050 and 2070). Each cell was identified with a number for reference and discussion in the text. Total number of species corresponds to each of the nine main clumps of cells (the number of species for each cell is given in Table S3); b -terrain roughness vs. number of observed spe-cies for each most affected cells, regression line corresponds to smooth local regression (loess), shaded area shows 0.95 confidence interval; c -heatmap plot showing values of Beta diversity among most affected cells based in the specimen occurrence data and the βsor index. Numbers correspond to cell numbers and coloured boxes represent the nine main clumps where cells are included. association between geographic distribution and environmental niche (Breiner et al. 2015). However, rare and threatened species are probably the most vulnerable to environmental change and most in need of such studies, leading to the 'rare species modelling paradox' (Lomba et al. 2010). In this study, using macroecological modelling together with present and future climate data, we have highlighted that climate change has an overall negative impact on the richness of threatened endemics of the Argentine vascular flora, including regions within the PA network and hotspots for those species. The possibility of directly including these taxa in the analyses allows us to take the first steps to characterize their potential threat to climate change, since threatened endemics are highly sensitive to global warming due to their small distribution range and low abundances (Gaston 1994;Lavergne et al. 2005;Parmesan 2006;Pearson et al. 2014).
Human-induced climate change is already affecting many weather and climate conditions in all regions of the world. Overall, with 1.5°C of global warming, there will be a global increase in heat waves, with longer hot seasons and shorter cold seasons, as well as changes in precipitation patterns leading to floods and droughts, which will change the distribution ranges of species and may lead to the extinction of the most sensitive ones. Compared to 1850-1900, the global surface temperature averaged over 2081-2100 is most likely to increase between 1.0°C and 1.8°C in the very low greenhouse gas (GHG) emissions scenario considered (shared socioeconomic pathway SSP1-1.9), but between 2.1°C and 3.5°C in the intermediate scenario (SSP2-4.5), and between 3.3°C and 5.7°C in the very high GHG emissions scenario (SSP5-8.5; IPCC 2021). Shared socioeconomic pathways (SSPs) are scenarios of projected socioeconomic global changes up to 2100, and unlike the RCPs describing different levels of greenhouse gases that could occur in the future, but without including any underlying socio-economic connection, SSPs are used to derive greenhouse gas emissions scenarios with different climate policies and socioeconomic trends that could shape future society. The two perspectives (SSP and RCP) were designed to be complementary, since while RCPs set the trajectories of greenhouse gas concentrations and the amount of warming that could occur by the end of the century, SSPs set the socio-economic scenario in which emission reductions will -or will not -be achieved. Thus, SSP1-1.9 represents a scenario with very low GHG emissions in which CO2 emissions are reduced to net zero around 2050, SSP2-4.5 a scenario with intermediate GHG emissions and CO 2 emissions around current levels until 2050, then declining but not reaching net zero by 2100, and SSP5-8.5 represents a scenario with very high GHG emissions, where CO 2 emissions triple by 2075 (IPCC 2021). These present and future scenarios highlight the imminent impact of climate change on biodiversity and the importance of implementing urgent conservation strategies. In addition, current estimates of threat indicate that 39% of global plant species are at risk of extinction . In Argentina, nearly half (47%, 800/1,683 spp.) of endemics are under some category of threat (VU, EN or CR), with the highest number of threatened endemics recorded in the richest families (i.e. Asteraceae, Fabaceae, Poaceae, Cactaceae), but with other families such as Apocynaceae, Bromeliaceae, Iridaceae and Violaceae presenting an elevated proportion on threatened species (Salariato et al. 2021). In this study, hotspots of threatened endemism were detected mainly in valleys and highlands along the Central Andes and the upper portion of the Southern Andes (or Dry Temperate Andes -DTA sensu Ezcurra and Gavini 2020). In particular, transition areas between the biogeographic provinces of the Puna, Prepuna and the Yungas in the Central Andes of Argentina (sensu Cabrera andWillink 1980 andlater Oyarzabal et al. 2018) were found to have the highest number of threatened endemics, and also the most severe declines in future richness associated with climate change. These areas, extended along the Central Andean Puna, Central Andean Dry Puna, Southern Andean Yungas, High Monte and Southern Andean steppe ecoregions (Olson et al. 2001), are characterized by their elevational diversity and the high environmental and climatic heterogeneity (Cabrera and Willink 1980), and have been previously reported as hotspots for endemic vascular plants (Aagesen et al. 2012;Godoy-Bürki et al. 2014;Elías and Aagesen 2019;Salariato and Zuloaga 2020). Furthermore, even though it has been pointed out that global hotspots of species richness may not be congruent with hotspots of endemism and threat (Orme et al. 2005), for the southern portion of the Central Andes areas of high species richness, endemism and threat overlap (Godoy-Bürki et al. 2014), in particular for plant families such as the Poaceae (Aagesen et al. 2009), Cactaceae (Ortega-Baes et al. 2012) and Brassicaceae (Salariato and Zuloaga 2020). The concentration of threatened endemics in these mountainous regions is consistent with previous studies reporting that areas with high topographic and climatic variability contain a disproportionate richness of narrowranging species (Ohlemüller et al. 2008;Luebert  and Weigend 2014; Platts et al. 2014;Perrigo et al. 2020). Over the northern Andes of Argentina, as well as most of the subtropical part of the country, climate warming appears to be accompanied by an increase in precipitation (Barros et al. 2000;Piovano et al. 2002). Conversely, for the Dry Temperate Andes (between 30°S and 40°S) climate change appears to be causing a reduction in precipitation, which together with the increase in temperature leads to glacier retreat, reduced river flows and increased aridification (Barros et al. 2015;Rivera et al. 2017). In particular, for the Andean region located in the 'Cuyan High Andean' biogeographic province (sensu Arana et al. 2017), this precipitation decrease, observed since the end of the 1980s, would represent a risk scenario for biodiversity and water resources of the Andean foothills (Agosta et al. 2014;Cabré and Nuñez 2020).
The decline in the richness of threatened endemics in Argentina was statistically significant also for cells within protected areas. In general, three conservation strategies are used to mitigate the negative effects of climate change on protected areas: to increase the quality of habitats, to increase the quantity of habitats and to increase the connectivity of habitats between the PA network (Mawdsley et al. 2009;Nuñez et al. 2013;Watson et al. 2014). However, it has been reported, from simulation studies, that none of these conservation strategies could fully compensate the negative impact of climate change for vascular plant (Wessely et al. 2017), suggesting that besides the in situ conservation strategies, it is crucially important to halt or at least reduce climate change itself in future years. In addition, global warming is enabling exotic species to invade new regions (Walther et al. 2009), and in many cases favours the extension of crop lands in semi-arid regions (Viglizzo et al. 1995;Magrin et al. 2009;Barros et al. 2015), which increases the levels of threat for threatened endemics by habitat degradation. All these factors urge us to advance in the study of biodiversity in all its dimensions (i.e. taxonomic, phylogenetic and functional) and in the planning of conservation strategies for the vascular flora of Argentina in general, but especially for endemic and threatened species, since this group is the most sensitive to threat and climate change.
The 32 consensus cells identified in this work as the most affected by climate change included nearly half (370 of 785 spp.) of all threatened endemic plant species registered for Argentina. Furthermore, the beta diversity among these cells was high, highlighting the restricted distribution ranges of these threatened species and the importance of each cell in terms of the species that it hosts. The three most vulnerable cells (3, 10, and 12) included mountain environments of the Central Andes along the Southern Andean Yungas, Central Andean Puna, and High Monte ecoregions, as well as their transition zones. This region of the Central Andes harbors the southernmost portion of the Andean Neotropical biota, characterized by a disproportionate species richness with ca 15% of the world's plant species (Myers et al. 2000;Luebert and Weigend 2014). From our occurrence data, these three cells include a total of 118 threatened endemics, about 15% of the total number of threatened endemics registered for Argentina. These results not only highlight the urgency of in situ conservation, but also the critical importance of ex situ (outside natural habitat) conservation programmes for the preservation of these threatened taxa (Havens et al. 2006). Ex situ plant conservation programmes use a diverse set of techniques, such as ex situ cultivation, assisted migration, long-term germplasm storage in seed banks, and plant tissue culture -micropropagation, to conserve threatened species (Li and Pritchard 2009). This is where botanic gardens and arboreta have a major role, being the key to achieving target 8 of the global strategy for plant conservation 2011-2020 (GSPC) (Sharrock 2012), which calls for 'at least 75% of threatened plant species in ex situ collections, preferably in the country of origin, and at least 20% available for recovery and restoration programmes'. In this way, the species list provided here for cells most affected by climate change can serve as a tool for the implementation of ex situ conservation strategies and programmes.
One of the main limitation of the MEM approach, implemented in this work to estimate species richness, remains in that it does not provide information regarding the identity of the species occurring at a given location, so it cannot predict changes in species composition unlike other methods such as S-SDM and SESAM (Guisan and Rahbek 2011;Zurell et al. 2020, Peyre et al. 2020, or J-SDM (Ovaskainen et al. 2017). Additionally, MEMs tend to underpredict the highest values and overpredict the lowest values (Dubuis et al. 2011). However, MEM can accurately predict current and future diversity (Algar et al. 2009;Dubuis et al. 2011;D'Amen et al. 2018;Biber et al.

3
Vol:. (1234567890) 2020) and, in our study, allows us to include into the analyses threatened species with few occurrences. Another important limitation is that the model and results are scale-dependent. Using smaller grid-sizes would allow for finer resolution of species richness patterns, hotspots of threatened species, and most vulnerable cells. Given that environmental variables vary considerably over 0.5 degrees (~ 55 km), especially in mountain areas, analyses at finer scales could reveal new patterns and the influence of additional climatic and other factors such as habitat-related variables (Baudraz et al. 2018). However, to do so, collection biases must be reduced. Using herbaria data, richness may be biased for some localities / cells, both those sampled more intensively (e.g. due to proximity to roads or the vicinity of research centres and herbaria) or, conversely, sparsely collected (due to difficult access; Oliveira et al. 2016). In order to reduce collection biases, poorly explored areas, such as the high peaks of the southernmost Andes or the southern portion of the Patagonian Steppe, need greater sampling effort, as well as incorporation of records from the local herbaria to different biodiversity databases (e.g. GBIF, DFA). After reducing this 'Wallacian shortfall', further analyses might be conducted at finer scales. Richness models might also be adjusted for each particular ecoregion to compare the results of the local models with those of the global models.
The results presented here on the impact of climate change on threatened endemic vascular plants of Argentina, together with the identification of the most affected geographic areas and the species they harbor, are intended to serve as a tool for planning in situ and ex situ conservation strategies in the country, and as a basis for future comprehensive flora conservation assessments. In addition to the analyses presented here, molecular phylogenies that include as many representatives of the Argentine flora as possible are also urgently needed, in order to consider phylogenetic diversity (Swenson 2014) as a complement to the number of species (taxonomic diversity) and to prevent its loss (homogenization of phylogenetic diversity of cells; Saladin et al. 2020). Another aspect worthy of study is functional diversity comprised by threatened endemics with respect to the remaining flora, as rare and / or threatened species may exhibit unique combinations of traits (Mouillot et al. 2013) and have important functional roles for ecosystems (Dıáz and Cabido 2001;Lyons et al. 2005;Dee et al. 2019). Additional steps regarding conservation actions for vascular plants in Argentina should include full assessment of the impact of climate change for the entire vascular flora (all native species, not only endemic ones), SDMs to study the effect of global warming on the species range (size but also fragmentation and elevation shift at the landscape level; Opdam and Wascher 2004;Liang et al. 2018). Complementary to MEM analyses, richness patterns for the native flora should be also estimated by other methods such S-SDM (Ferrier and Guisan 2006), SESAM (Guisan and Rahbek 2011), or J-SDM (Ovaskainen et al. 2017), and the dispersal capacity of species and landscape fragmentation could also be incorporated into hybrid mechanistic models for more accurate estimations of future species distributions (Engler et al. 2009;Zurell et al. 2016;Peyre et al. 2020;Zanatta et al. 2020).