This is the accepted version of the article: Pérez Navarro, María Ángeles; Sapes, Gerard; Batllori, Enric; [et al.]. Climatic suitability derived from species distribution models captures community responses to an extreme drought episode. DOI 10.1007/s10021-018-0254-0 This version is avaible at https://ddd.uab.cat/record/203570 under the terms of the license 1 Climatic suitability derived from species distribution models captures 2 community responses to an extreme drought episode 3 Pérez-Navarro M.A.1*, Sapes G.2, Batllori E.1, Serra-Diaz J.M.3,4, Esteve M.A.5, Lloret F.1 4 1 5 CREAF, Universitat Autònoma Barcelona, Cerdanyola del Vallès, 08193 Barcelona, Spain, 0034 6 935868087, e-mail: m.angeles582@gmail.com, Orcid id: 0000-0001-5553-995X 2 7 Division of Biological Sciences, The University of Montana, Missoula, MT 59812, USA 3 8 Section for Ecoinformatics and Biodiversity, Department of Bioscience, Aarhus University, 9 8000, Aarhus C, Denmark 4 10 UMR Silva, AgroParisTech, Université de Lorraine, INRA, 54000 Nancy, France. 5 11 Department of Ecology and Hydrology, University of Murcia, Campus de Espinardo, 30100, 12 Murcia, Spain * 13 Corresponding author 14 Pérez-Navarro M.A. e-mail: m.angeles582@gmail.com. Ocid id: 0000-0001-5553-995X 15 1 16 Abstract 17 The differential responses of co-occurring species in rich communities to climate change - 18 particularly to drought episodes - have fairly been unexplored. Species Distribution Models 19 (SDMs) are used to assess changes in species suitability under environmental shifts, but whether 20 they can portray population and community responses is largely undetermined, especially in 21 relation to extreme events. Here we studied a shrubland community in SE Spain since this region 22 constitute an ecotone between the Mediterranean biome and subtropical arid areas, and it has 23 recently suffered its driest hydrological year on record. We used four different modelling 24 algorithms (Mahalanobis distance, GAM, BRT and MAXENT) to estimate species’ climatic 25 suitability before (1950-2000) and during the extreme drought. For each SDM, we correlated 26 species’ climatic suitability with their remaining green canopy as a proxy for species resistance 27 to drought. We consistently found a positive correlation between remaining green canopy and 28 species’ climatic suitability before the event. This relationship supports the hypothesis of a 29 higher vulnerability of populations living closer to their species’ limits of aridity tolerance. 30 Contrastingly, climatic suitability during the drought did not correlate with remaining green 31 canopy, likely because of the exceptional episode led to almost zero suitability values. Overall, 32 our approach highlights climatic niche modelling as a robust approach to standardizing and 33 comparing the behavior of different co-occurring species facing strong climatic fluctuations. 34 Although many processes contribute to resistance to climatic extremes, the results confirm the 35 relevance of populations’ position in the species’ climatic niche for explaining sensitivity to 36 climate change. 37 38 Key words: climatic suitability, SDMs, extreme climatic events, drought resistance, niche, 39 dieback. 40 2 41 Introduction 42 The climatic trends observed over the last decades are promoting vegetation shifts (Parmesan 43 and Yohe 2003), phenological changes (Zavaleta and others 2003) and modifications to 44 disturbance regimes (Mouillot and others 2002; Allen and others 2015), as well as altering the 45 interactions between these processes (Franklin and others 2016). However, the adjustment of 46 populations to changing climatic conditions may be more influenced by the extremes of climatic 47 variability than by average climate trends. For instance, mortality and recruitment processes – 48 which shape species' distributions and ranges – may be largely conditioned by pulses of extreme 49 climatic conditions such as extreme drought events (del Cacho and Lloret 2012; Greenwood and 50 others 2017). 51 Vegetation mortality and die-off processes associated with climatic warming have often been 52 observed at ecotones corresponding to the rear edge of species’ distributions (Allen and 53 Breshears 1998; Bigler and others 2006; Jump and others 2006; Lesica and Crone 2016), 54 supporting the assumption that a decline in plant populations may be more significant at their 55 equatorial latitudinal or lowland altitudinal margins (Thomas and others 2004). Translated into 56 the perspective of a plant community, marked by the coexistence of species that have adapted 57 differently to environmental conditions, mortality processes would have a greater influence on 58 the populations of species located close to their tolerance limits, to the benefit of other species 59 that find the new environment more suitable (Martínez-Vilalta and Lloret 2016). This hypothesis 60 implicitly correspond to the biogeographic paradigm that species perform better in their 61 geographical center of distribution than they do in the margins (Centre-Periphery hypothesis, 62 see Sexton and others 2009; but see Dallas and others 2017), with the further assumption that 63 geographical and environmental spaces are mostly concordant (Pironon and others 2015). 64 Species Distribution Models (SDMs) have been used to test the relationship between species’ 65 climatic niche and their physiological or demographic performance (Serra-Diaz and others 2013; 66 Pironon and others 2015; van der Maaten and others 2017). These are statistical models that 3 67 relate the location of species occurrences to the environmental data on these sites (Franklin 68 2010). The SDM approach is based on the assumption that species occurrence portrays the 69 environmental and biotic conditions that are suitable for species to survive and reproduce (i.e. 70 the realized niche; Pulliam 2000; Soberón 2007; Peterson and others 2011), and so model 71 outputs are interpreted as a species-specific index of relative suitability or habitat suitability. 72 Accordingly, SDMs have been widely used under average climatic conditions (climatic norms) to 73 predict past or future distributional changes of species (Elith and Leathwick 2009). However, it 74 is not yet known whether these models are able to capture the impact of extreme climatic 75 events (e.g. droughts), especially in relation to community dynamics. 76 In the Mediterranean basin, vegetation has developed different strategies for dealing with 77 variable rainfall, such as shallow roots and deciduous summer leaves (Valladares and others 78 2004; Zunzunegui and others 2005). Nevertheless, these adaptive syndromes may not be 79 enough under the predicted scenarios of increased climatic extremes (IPCC 2013). In this 80 context, the southeast of the Iberian Peninsula has recently experienced the driest year on 81 record (AEMET 2014), causing an extensive vegetation die-off event in areas dominated by 82 shrubland (Esteve-Selma and others 2015). This Iberian region represents the ecotone between 83 the Mediterranean biome and subtropical shrublands of arid lands (Esteve-Selma and others 84 2010). This recent drought-induced mortality event therefore offers the possiblity to assess 85 community dynamics in relation to biogeographical paradigms at the limits – in this case, the 86 aridity margin - of the biome's distribution, which are areas considered to be very sensitive to 87 climate change (Guiot and Cramer 2016). 88 In this study, we use a shrubland community at the arid southern limit of the Mediterranean 89 biome to assess the differential response of coexisting species to an extreme drought event 90 according to species’ climatic suitability, as determined from SDMs. Specifically, we test whether 91 populations living close to the edge of their species’ climatic niche (i.e. lower suitability values 92 compared to the niche's optimal value) are more vulnerable to such extreme events than 4 93 populations living closer to their niche center. We use the remaining green canopy of species 94 after the drought event to examine the correlation between drought-induced die-off and 95 species’ climatic suitability, considering both the historical suitability, as inferred from historical 96 climatic series, and the drought episode suitability, as reflected by the conditions during the 97 drought event. Given the considerable amount of uncertainty existing with respect to the 98 various SDM techniques, which use model-specific algorithms (Araújo and New 2007), we also 99 test whether the relationship between species' die-off and their climatic suitability depends on 100 the SDM algorithm applied. For this purpose, we applied four SDMs with highly contrasting 101 approaches (Mahalanobis distance, Generalized Additive Models –GAM–, Boosted Regression 102 Trees –BRT–, and Maximum Entropy approaches –MaxEnt–) to determine species suitability, 103 that was later correlated to species' die-off. 104 Material and methods 105 Study area 106 The study was carried out in two semi-arid shrubland areas in the province of Murcia (southeast 107 of the Iberian Peninsula) (Figure 1), Campo de Cagitán (38° 06’ N, 1° 32’ W) and Oro Mountain 108 (38° 11’ N, 01° 30’W), 10 km apart but with similar soil characteristics and climatic conditions. 109 The Campo de Cagitán site was covered by a small expanse of scrubland embedded in an 110 agricultural matrix, and the Oro Mountain site was occupied by shrubland close to an open pine 111 forest on a hill slope. The overall sampled area amounted to 19 hm2. 112 The potential vegetation comprises an open forest of Pinus halepensis L. and a sclerophyllous 113 shrubland (garrigue) dominated by Quercus coccifera L., Pistacia lentiscus L., Olea europea L., 114 Rhamnus lycioides L. and Juniperus oxycedrus L., along with a highly diverse range of small 115 shrubs, such as Thymus hyemalis Lange and Helianthemum spp. (Braun-Blanquet and Bolòs 116 1957). The current landscape in these regions is a highly fragmented cropland that is either in 117 use or recently abandoned, with small patches of forest or shrubland interspersed between the 118 crops. The steep slopes in some areas preclude the presence of agricultural crops but they are 5 119 instead covered by scrublands, often containing the tussock grass Stipa tenacissima L. (which 120 was cultivated for fibers until the 1960s) and occasional open pine forests. 121 The region is included within the Mesomediterranean thermoclimatic belt and the 122 Mediterranean xeric bioclimatic region (Rivas-Martínez and others 2011), characterized by 123 annual mean temperatures of 18.5 °C and an annual rainfall of 200-350 mm. Precipitation in the 124 area is low and mainly concentrated in the fall, with great variability between years. During the 125 hydrological year 2013-2014 the Region of Murcia suffered the worst drought on record since 126 1941, demonstrating the extreme conditions of the event. During the drought event a mean 127 regional rainfall of 146.5 mm was recorded; this corresponds to just 46% of the average value 128 for the period 1971-2000 (Figure 1, AEMET 2014). This episode led to high mortality and 129 defoliation in different forests and shrublands (Esteve-Selma and others 2015). 130 Field sampling 131 In March 2015, a set of ten 50m2 replicate plots were established in the study region, three in 132 the Campo de Cagitán and seven in the Oro Mountain, according to the shrubland surface area 133 available on each site. This sampling design reflects the region's highly fragmented habitat, 134 which prevented us from establishing ten replicates in a single location. 135 Sampling plots were established with the following criteria: 1) shrubland with no signs of recent 136 disturbance, with high species richness, different life forms, and low S. tenacissima density; and 137 2) low pine presence, in order to avoid wetter microenvironments caused by the shade of tree 138 canopies, which could affect the moisture in the air and soil. Each plot consisted of two linear 139 transects of 25 meters long by 1 wide. On each plot, we recorded the total number of individuals 140 per woody species, estimated their size by measuring two perpendicular diameters crossing at 141 the center of each individual and visually estimated the proportion of remaining green canopy 142 (RGC) per individual. A total of 22 species were sampled (Table 2). RGC levels were visually 143 estimated as a proxy for the species response to drought(die-off) as the percentage of green 144 leaves present relative to the amount in healthy individuals found in the study area (Sapes and 6 145 others 2017). To ensure that the green cover loss resulted from the drought of the previous 146 year, we avoided individuals with signs of older decay (e.g., stumps, decomposed stems, 147 branches with no thin tips). To determine the reliability of the visual RGC estimate, we also 148 measured the length (cm) of the segments occupied by green leaves and dry leaves (including 149 segments with no leaves) along a linear path from the tip to the base of two representative 150 branches of ten individuals per species, on every plot where a species was present. Then, a 151 directly measured RGC value was calculated per individual as [branch length with green leaves 152 (cm)/total branch length (cm)]*100. When ten individuals per species were not found within a 153 given plot, we measured the closest individuals to the sampled transects until ten replicates 154 were attained. Pearson’s correlation between direct and visually estimated values of RGC was 155 calculated for plants from each species, always resulting in values higher than 0.7 156 (Supplementary material Appendix 1 Table 1). These analyses support the use of visual estimate 157 of RGC as a proxy for die-off (Sapes and others 2017, see Supplementary material Appendix 3). 158 Considering this high correlation and the limited number of individuals with real measures of 159 defoliation (10 per species and plot), we used the visually estimated RGC (made for every 160 individual) for the statistical analyses. 161 The variables described above were then scaled to the landscape level. First, we calculated the 162 following information for each plot: average percentage of visually estimated RGC for each 163 species; species frequency (the number of plots where each species was found over the total 164 number of plots) and the average size of each species (as the product of the two diameters 165 measured in the plants), since RGC could be affected by species size or relative abundance 166 (Lloret and others 2016; Sapes and others 2017). The values of species' RGC and size were then 167 averaged across plots. Finally, to account for the different species’ strategies in relation to leaf 168 longevity and annual seasonality (Valladares and others 2004), all the species were classified 169 into one of the following foliar strategies: 1-evergreen, 2- semi-deciduous, 3-summer deciduous, 170 and 4-retamoid or leafless species. 7 171 Climatic suitability modelling 172 For the 22 sampled species, we built SDMs using four different algorithms - Mahalanobis 173 distance, GAM, BRT, and MaxEnt - to assess the robustness of the potential relationships 174 between visual estimates of RGC and the climatic suitability output of these models. These 175 models represent four highly differentiated modelling methods: distance-based models 176 (Mahalanobis distance), regression-based models (GAM), decision tree-based methods (BRT) 177 and a machine-learning technique based on the principle of maximum entropy (MaxEnt). 178 Therefore, they represent a gradient of complexity, where some models such as Mahalanobis 179 distance only consider linear relationships between predictors (Franklin 2010) whereas others 180 such as MaxEnt and BRT can include high-order interaction terms (Elith and others 2008; Phillips 181 and Dudík 2008). 182 The geographical occurrence data for each species were obtained from GBIF (Global Biodiversity 183 Information Facility: ). Occurrence data were filtered in order to remove 184 taxonomic and geographic inconsistencies and reduce dense local sampling by randomly 185 thinning species' records to one observation per 1x1 km grid cell. We considered the whole 186 distributional range of species; thus, the available number of occurrences was considerably 187 different from one species to the next, ranging from 200 to 6,000 after filtering. For each set of 188 filtered occurrences, 70% of presences were reserved for fitting the model (training data) and 189 the other 30% for the validation set (testing data), according to the number of environmental 190 predictors selected in our models and following the rule described by Huberty 1994 for 191 determining the optimum partitioning of training and test data. To improve the models' 192 performance accuracy (Barbet-Massin and others 2012), 100,000 random background points 193 were simulated for each species to fit both GAM and MaxEnt algorithms, and a random set of 194 pseudo-absences equivalent to the number of each species occurrences was simulated for BRT. 195 Since Mahalanobis distance works without simulated absences, no background points were 196 used in this case. The background extension was delineated in order to represent the current or 8 197 past available geographical space for the selected species (M dimension sensu Soberón 2007). 198 In our case, the geographical region used to establish species background was the 199 Mediterranean basin. 200 Six climatic variables representative of Mediterranean climate were used as predictors to 201 calibrate the suitability models: isothermality (mean diurnal temperature range/temperature 202 annual range), temperature seasonality, mean temperature of wettest quarter, mean 203 temperature of driest quarter, annual precipitation and precipitation seasonality, all of them 204 with 1 by 1 Km resolution. These variables were selected from the 19 bioclimatic variables 205 available in Worldclim.org (version 1.4) for the period 1950-2000 (Hijmans and others 2005), 206 according to the knowledge of the species’ ecological requirements and in order to reduce 207 variables’ collinearity. Pearson correlation and variation inflation factor (VIF) among variables 208 were always less than 0.75 and 5, respectively. Additionally, we used monthly precipitation and 209 maximum, minimum and mean temperature records over the 2013-2014 period from between 210 68 and 114 weather stations of the Spanish Meteorological Agency (AEMET) to elaborate the 211 climatic layers during the drought event (also in 1 by 1 km resolution), following Ninyerola and 212 others (2000), and using the ‘biovars’ function (dismo package; Hijmans and others 2016). To 213 minimize differences in the climatic interpolation methods between Worldclim and Ninyerola 214 and others (2000), only latitude, longitude and elevation were used as explanatory variables for 215 climatic data. In addition, we applied MESS analysis between these two data set over the Spanish 216 territory (where AEMET data are available) during the 1950-2000 period to assess dataset 217 dissimilarities, showing the high concordance and comparability of both climatic interpolation 218 methods over the entire extension and particularly over the study region (Supplementary 219 material Appendix 1 Figure A1). Finally, species’ historic climatic suitability (HCS) was estimated 220 projecting the models over the climatic layers for the period 1950-2000, whereas species’ 221 climatic suitability during the drought event (episode climatic suitability, ECS) was estimated by 222 projecting the calibrated models over the climate layers of the anomaly period 2013-2014. 9 223 For the four implemented SDM algorithms, model settings were selected following 224 recommendations from the literature, and partial dependence plots and predictive maps were 225 assessed to exclude those settings that produced unreliable response curves or distribution 226 maps. The final selected settings and main literature are showed in Table 1. To make all the 227 model outputs comparable (between 0 and 1 probability values), log raw output transformation 228 was applied for MaxEnt models (Phillips and Dudík 2008) and distance transformation into p- 229 values was applied for Mahalanobis distances (Clark and others 1993). Each model's predictive 230 performance was assessed by comparing model predictions with testing data, using the Area 231 Under Receiver Operating-characteristic Curve (AUC, Fielding and Bell 1997) and the Boyce index 232 (Boyce and others 2002; Hirzel and others 2006). These evaluation methods are considered a 233 reliable approach for our models and allowing comparison among them, since all the models 234 were fitted with the same species data set and environmental extension of layers (Hirzel and 235 others 2006; Franklin 2010). Finally, multivariate environmental similarity surface (MESS, Elith 236 and others 2010) analyses were carried out to measure the similarity between historical climate 237 and the extreme drought period in the occurrence locations; these analyses allowed to identify 238 extrapolation for predictions during the anomalous period as areas with high climatic 239 dissimilarities. 240 Statistical analyses 241 Generalized Linear Models (GLM) with normal distributions were performed to assess the 242 relationship between SDM-inferred HCS and ECS for each species and their die-off recorded in 243 the field. The visually estimated species RGC was used as a response variable whereas HCS, ECS, 244 the interaction between HCS and ECS, species size (logarithmically transformed), species 245 frequency, and species foliar category were introduced as explanatory variables. Difference 246 between HCS and ECS was discarded as an explanatory variable in the models, due to the high 247 correlation with ECS which produced same models results. 10 248 The final models were selected according to stepwise selection based on AIC (Akaike Information 249 Criterion). In addition, Phylogenetic Generalized Linear Models (PGLS) using Phylomatic 250 distances (Webb and Donoghue 2005) were performed in order to assess the potential effect of 251 phylogenetic species relationships in the selected model (Freckleton and others 2002). This PGLS 252 was finally discarded from the final analyses since phylogeny was not significant in any case 253 (lambda =0). Finally, consistency in the climatic suitability estimates (HCS, ECS) obtained from 254 each of the four implemented SDMs was tested by pairwise comparisons, using Pearson 255 correlation tests. All statistical analyses were carried out with R version 3.3.2 (R Core Team 256 2016). 257 Results 258 All the four SDMs algorithms developed showed high performance accuracy values with AUC 259 values higher than 0.75 and Boyce index’ values being always positive and higher than 0.5 (Elith 260 and others 2002; Hirzel and others 2006; Supplementary material Appendix 1 Figure A2). 261 Particularly AUC mean values were 0.96±0.02 and Boyce index mean values were 0.93±0.07 262 MESS analyses showed that precipitation seasonality exhibited high dissimilarity between 263 extreme event climatic data and training predictor data near the coastal region. However, the 264 values corresponding to the drought episode were never outside training boundary values for 265 the study locations (Supplementary material Appendix 1 Figure A3). Climatic suitability dropped 266 dramatically during the drought episode for all species, irrespective of the SDM method, as 267 shown by the comparison between the respective HCS and ECS values (Figure 2, Supplementary 268 material Appendix 1 Table A2). 269 For the majority of SDM algorithms the stepwise GLM model selection determined that the most 270 parsimonious models explaining species’ RGC were those including HCS, foliar category and 271 species size as explanatory variables. Only in the case of Mahalanobis distance did the stepwise 272 GLM model selection fail to remove any explanatory variables from the saturated model. There 273 was some variation in the significant variables associated with RGC in the different SDM 11 274 algorithms. All the selected variables were significant in BRT, while in MaxEnt HCS and foliar 275 category were significant, in GAM foliar category was significant and HCS was only marginally 276 significant, and in Mahalanobis only HCS and foliar category were marginally significant (Table 277 3). 278 In all cases RGC was positively related to HCS (Figure 2 and 3). Foliar categories 3 (summer 279 deciduous species) and 4 (leafless species) presented a significant, negative correlation with RGC 280 in BRT, GAM and MaxEnt models, while in the Mahalanobis distance model, foliar category 3 281 was only marginally significant. Finally, species size was only significantly negatively related with 282 RGC in BRT models (Table 3). Species frequency, ECS or the interaction between HCS and ECS 283 were not significant in any model (Table 3). 284 The values of species’ HCS varied from model to model. BRT predicted the highest suitability 285 values for a given species and MaxEnt the lowest ones. This pattern was consistent for almost 286 all species (Supplementary material Appendix 1 Figure A4). For all species, the majority of 287 pairwise Pearson correlations between the HCS values inferred from the different SDM 288 algorithms were significant, with correlation values ranging between 0.45 and 0.78 289 (Supplementary material Appendix 1 Figure A5). However, the correlation between 290 Mahalanobis distance and BRT was not significant (Supplementary material Appendix 1 Figure 291 A5). In contrast, most pairwise correlations between models were not significant for ECS, likely 292 due to the extremely low ECS values exhibited by most species; in this case, when a significant 293 correlation was found, the relationship was driven by a single outlier value (Supplementary 294 material Appendix 1 Figure A6). 295 Discussion 296 We found a clear relationship between field measurements of species performance under an 297 extreme drought episode and the historical climatic suitability (HCS) of species derived from 298 SDMs. Within the studied community, co-occurring species living closer to their climatic 299 tolerance limit -identified by low HCS values compared to the optimal value of the distribution 12 300 range- proved more vulnerable to the extreme drought episode. This climatic limit corresponds 301 to the aridity margin of species' climatic niche (Supplementary material Appendix 2). These 302 results are consistent with the relationship observed between the decay of shrubland and 303 woodland species and the decrease in climatic suitability in other semi-arid areas in Spain (Sapes 304 and others 2017) and Southwestern North America (Lloret and Kitzberger 2018). This 305 relationship is also consistent with other studies which suggest that species' sensitivity to 306 climate change is related to niche characteristics such as mean niche position and niche breadth 307 (Thuiller and others 2005; Broennimann and others 2006). Species in the climatic niche margins 308 are generally assumed to exhibit lower survivorship and recruitment and higher extinction risk 309 because of the less favorable environmental conditions (Weber and others 2016). Precisely in 310 these situations of the environmental space closest to the physiological tolerance limits of the 311 species, the effect of climate variability is probably more severe (Zimmermann and others 2009), 312 promoting species’ decline or range shifts at the trailing edge of species distribution (Bigler and 313 others 2006; Walther and others 2009). 314 Species’ drought responses and climatic suitability 315 Studies that compare habitat suitability with different species’ performances (population 316 density, growth, recruitment, fecundity, etc.) along the species distribution gradient are scarce 317 and still not fully conclusive with respect to general biogeographic paradigms (Centre-Periphery 318 hypothesis, Wright and others 2006; Sexton and others 2009; Thuiller and others 2010; Abeli 319 and others 2014; Csergő and others 2017). Likely species interaction, local variables or 320 adaptation mechanisms underlie the limited evidence of the relationship between species’ 321 performance and climate suitability (Sexton and others 2009, 2014; Dallas and others 2017; 322 Lloret and Kitzberger 2018). Our results throw some light in this sense as they support the 323 relationship between species’ performance when climatic conditions are extreme and 324 climatically-based descriptions of their suitability (i.e. HCS). 13 325 We also expected that populations experiencing higher displacement of climatic suitability 326 during the extreme event (low ECS) would experience greater leaf loses and higher mortality 327 rates. However, contrary to our expectations, we found that species' suitability during the 328 extreme episode (ECS) did not significantly explain species leaf loses in the studied community. 329 Our extremely low ECS levels observed for all species probably made it impossible to obtain 330 contrasted values of ECS among them. These low values indicate that the climatic episode was 331 extreme enough to displace all the studied populations far from their climatic optimum, even 332 for those species that were closer to this optimum during the historical period (Figure 2). In 333 addition to the exceptionality of the extreme event, the extremely low ECS values may derive 334 from 1) the averaged climatic data used for calibrating the models, which does not reflect the 335 variability or annual extremes during the considered 50-year period and 2) the limited ability of 336 models to predict suitability under climatic scenarios that are highly different from the period 337 used to fit the models (Elith and others 2010), as shown by the low MESS values, particularly for 338 precipitation seasonality –bio15- (although these were not negative in the study site) 339 (Supplementary material Appendix 1 Figure A3). This situation may also amplify the differences 340 between different algorithms and species’ prevalence data in the predictions (Thuiller 2004; 341 Pearson and others 2006), as supported by the low correlations between the ECS predicted by 342 the different models (Supplementary material Appendix 1 Figure A6). 343 In addition to HCS, foliar categories were also significant in explaining the observed species’ RGC, 344 suggesting that leaf strategy and seasonal senescence play a major role in understanding species 345 performance under strong drought conditions, at least in Mediterranean type ecosystems. Our 346 results show that summer-deciduous and leafless species always present significantly lower 347 values of RGC. This result is in part expected due to the general strategy in Mediterranean 348 species of dropping leaves during the dry, hot season to limit evapotranspiration and water loss 349 (de la Riva and others 2016a). In the Mediterranean basin this semideciduous mechanism 350 typically appears in combination with shallow roots and low water potentials as an anisohydric 14 351 syndrome, in contrast to species with hydrostable syndromes, which present sclerophyllous 352 leaves, more sensitive stomatal control and deeper roots (Zunzunegui and others 2005; de la 353 Riva and others 2016b). Thus, estimates of RGC as a proxy of drought resistance can be 354 misleading if these foliar strategies are not considered (Lloret and others 2016). In addition to 355 foliar strategies, other physiological features and local factors may modulate the interspecific 356 variability of responses to a given drought episode. For instance, species-specific resistance to 357 hydraulic failure and carbon economy (McDowell and others 2008; Anderegg and others 2012; 358 Adams and others 2017) and mutualistic and antagonistic biotic interactions (Lloret and others 359 2012; Valladares and others 2014). 360 SDM algorithms and demographic performance 361 Despite all these potential sources of variability, the four different SDM algorithms used in our 362 approach highlighted the positive relationship between climatic suitability (HCS) and resistance 363 to drought (RGC). These results were consistent across species, as shown by the high correlation 364 between the different models' HCS values (Supplementary material Appendix 1 Figure A5). The 365 agreement holds despite the wide variety of the modeling approaches. However, the ‘simplest’ 366 models (Mahalanobis distance and GAM) showed the lowest performance explaining species’ 367 RGC compared to ‘complex’ ones (BRT and MaxEnt). This difference in algorithm performance 368 highlight the importance of the interactions between climatic variables and non-linear 369 relationships when assessing species’ responses to climate, and thus, supporting the use of 370 SDMs versus simpler approaches based on univariate or multivariate correlations of 371 demographic performance with climatic variables. Moreover, algorithms that are generally 372 calibrated to produce smoother response curves, such as GLMs and GAMs, would be more 373 accurate to predict habitat suitability under new conditions (Elith and others 2010; Merow and 374 others 2014), while models based on presence-only data are more appropriate for predicting 375 the lowest suitability values in these scenarios (Pearson and others 2006). There is no general 376 agreement, however, about the most accurate algorithm in relation to situations of range shift 15 377 because even simpler models can lead to erroneous outputs (Elith and others 2010; Merow and 378 others 2014). We therefore urge ecologists to assess the degree of model complexity needed to 379 use SDMs as a proxy of ecological mechanisms, such as defoliation in this case. 380 Caution should also be taken when interpreting SDM predictions, given the assumptions that 381 these kind of models implicitly include (Pearson and Dawson 2003). Among other limitations, 382 these models commonly use only climatic predictors with a broad resolution (~1 km2), 383 disregarding other meaningful abiotic factors, and they are also unable to capture microclimatic 384 effects at small spatial scales (Franklin and others 2013; Lenoir and others 2013; D’Amen and 385 others 2017). Microsite factors could be particularly important for our study, given that soil 386 features and depth, slope and orientation are especially relevant to species survival under 387 extreme drought (Colwell and others 2008; Hamerlynck and McAuliffe 2008). Furthermore, 388 SDMs assume that species respond homogeneously to climate change across their range, not 389 including intra-specific genetic variability and phenotypic plasticity, which may also favor 390 species' local adaptation under unfavorable conditions (Benito Garzón and others 2011; Lloret 391 and García 2016). 392 The drought episode experienced in the Region of Murcia in 2013-2014 was extraordinary in 393 historical terms, but these climatic situations are expected to become more frequent in the 394 future (Sheffield and Wood 2008). The ability of plant communities to withstand these events 395 and subsequently recover their green canopy will depend on both physiological traits related to 396 the adaptive syndromes of Mediterranean species (Peñuelas and others 2001) and the balance 397 between demographic processes such as mortality, growth, and recruitment (Lloret and others 398 2012). Under drier climatic scenarios, leafless and semi-deciduous species with shallow roots 399 (xerophytic malacophyllous) would be expected to be to take more advantage of scarce and 400 irregular rainfalls than sclerophyllous species with deeper roots. This is consistent with the 401 particularly high HCS values obtained for malacophyllous species (Figure 2 and 3). These 402 potential changes in species dominance within the community will likely lead to less productive 16 403 shrublands dominated by smaller species (Valladares and others 2004). Since species' climatic 404 suitability is broadly related to both physiological and demographic species performance 405 (Martinez-Meyer and others 2013), indexes describing climatic suitability can provide rough 406 estimates of species' vulnerability to extreme climatic episodes. While community resistance 407 and resilience could minimize ecological changes, the recurrence of these extreme drought 408 events could lead these Mediterranean communities to cross thresholds beyond which they 409 could collapse (Vicente-Serrano and others 2013; Valladares and others 2014). This depletion of 410 resilience in semi-arid shrubland communities could promote transitions to desert-like 411 ecosystems, as has been predicted by some climate change scenarios for southern areas of 412 Europe (Guiot and Cramer 2016). Accordingly, this study shows the impact of extreme drought 413 events even on communities supposedly well adapted to drought conditions (Lázaro and others 414 2001; Sapes and others 2017). 415 Conclusion 416 This study confirms the role of population position within its species climatic niche in explaining 417 populations' vulnerability to extreme climatic events. In the studied semi-arid shrubland, species 418 closer to their climatic tolerance limit were more vulnerable to extreme drought. Thus, the 419 predicted recurrence of severe drought events could reduce the community resilience, 420 increasing the risk of desertification in these arid lands. Our study empirically concurs with the 421 trends foreseen by theoretical models, based on predicted suitability and correlations with 422 drought response. This concurrence supports the use of SDMs to assess the impact of climate 423 change on plant communities, particularly in extreme climatic conditions. This approach, which 424 links species performance with regional biogeographic patterns, can probably be applied to 425 other processes heavily determined by strong climatic fluctuations. 426 Acknowledgements 427 The research underlying this work has been supported by the Spanish Ministry of Education 428 through a doctoral grant (FPU14/03519) and by the Spanish Ministry of Economy and 17 429 Competitiveness through the BIOCLIM project (CGL2015-67419-R). We thank GBIF and 430 WorldClim for making their data freely available online. We are also grateful to our colleagues 431 at the University of Murcia and several friends who helped during field sampling. Finally, we 432 want to thank the two anonymous reviewers for their helpful comments. 433 18 434 References 435 AEMET, Spanish Weather Agency. 2014. Avance climatológico mensual mes de septiembre. 436 Murcia. 437 Abeli T, Gentili R, Mondoni A, Orsenigo S, Rossi G. 2014. Effects of marginality on plant 438 population performance. J Biogeogr 41:239–49. 439 Allen CD, Breshears DD. 1998. Drought-induced shift of a forest – woodland ecotone : Rapid 440 landscape response to climate variation. Ecology 95:14839–42. 441 Allen CD, Breshears DD, McDowell NG. 2015. On underestimation of global vulnerability to tree 442 mortality and forest die-off from hotter drought in the Anthropocene. Ecosphere 443 6:art129. 444 Anderegg WRL, Berry JA, Smith DD, Sperry JS, Anderegg LDL, Field CB. 2012. The roles of 445 hydraulic and carbon stress in a widespread climate-induced forest die-off. Proc Natl 446 Acad Sci 109:233–7. 447 Araújo MB, New M. 2007. Ensemble forecasting of species distributions. Trends Ecol Evol 448 22:42–7. 449 Barbet-Massin M, Jiguet F, Albert CH, Thuiller W. 2012. Selecting pseudo-absences for species 450 distribution models: How, where and how many? Methods Ecol Evol 3:327–38. 451 Benito Garzón M, Alía R, Robson TM, Zavala MA. 2011. Intra-specific variability and plasticity 452 influence potential tree species distributions under climate change. Glob Ecol Biogeogr 453 Bigler C, Bräker OU, Bugmann H, Dobbertin M, Rigling A. 2006. Drought as an inciting mortality 454 factor in scots pine stands of the Valais, Switzerland. Ecosystems 9:330–43. 455 Boyce MS, Vernier PR, Nielsen SE, Schmiegelow FK. 2002. Evaluating resource selection 456 functions. Ecol Modell 157:281–300. 19 457 Braun-Blanquet J, Bolòs O. 1957. The plant communities of the Central Ebro Basin and their 458 dynamics. An la Estac Exp Aula Dei 5:1–266. 459 Broennimann O, Thuiller W, Hughes G, Midgley GF, Alkemade JMR, Guisan A. 2006. Do 460 geographic distribution, niche property and life form explain plants’ vulnerability to 461 global change? Glob Chang Biol 12:1079–93. 462 del Cacho M, Lloret F. 2012. Resilience of Mediterranean shrubland to a severe drought 463 episode: The role of seed bank and seedling emergence. Plant Biol 14:458–66. 464 Calenge C. 2015. Home Range Estimation in R : the adehabitatHR Package. R vignette:1–60. 465 Clark JD, Dunn JE, Smith KG. 1993. A Multivariate Model of Female Black Bear Habitat Use for a 466 Geographic Information System. J Wildl Manage 57:519. 467 Colwell RK, Brehm G, Cardelus CL, Gilman AC, Longino JT. 2008. Global Warming, Elevational 468 Range Shifts, and Lowland Biotic Attrition in the Wet Tropics. Science (80- ) 322:258–61. 469 Csergő AM, Salguero-Gómez R, Broennimann O, Coutts SR, Guisan A, Angert AL, Welk E, Stott I, 470 Enquist BJ, McGill B, Svenning J-C, Violle C, Buckley YM. 2017. Less favourable climates 471 constrain demographic strategies in plants. Gurevitch J, editor. Ecol Lett 20:969–80. 472 D’Amen M, Rahbek C, Zimmermann NE, Guisan A. 2017. Spatial predictions at the community 473 level: From current approaches to future frameworks. Biol Rev 92:169–87. 474 Dallas T, Decker RR, Hastings A. 2017. Species are not most abundant in the centre of their 475 geographic range or climatic niche. Anderson M, editor. Ecol Lett 20:1526–33. 476 Elith J, Burgman MA, Regan HM. 2002. Mapping epistematic uncertainties and vague conceps 477 in predictions of species distribution. Ecol Modell 157:313–29. 478 Elith J, Kearney M, Phillips S. 2010. The art of modelling range-shifting species. Methods Ecol 20 479 Evol 1:330–42. 480 Elith J, Leathwick JR. 2009. Species Distribution Models: Ecological Explanation and Prediction 481 Across Space and Time. Annu Rev Ecol Evol Syst 40:677–97. 482 Elith J, Leathwick JR, Hastie T. 2008. A working guide to boosted regression trees. J Anim Ecol 483 77:802–13. 484 Esteve-Selma MA, Mart??nez-Fern??ndez J, Hern??ndez I, Mont??vez JP, Lopez JJ, Calvo JF, 485 Robledano F. 2010. Effects of climatic change on the distribution and conservation of 486 Mediterranean forests: The case of Tetraclinis articulata in the Iberian Peninsula. 487 Biodivers Conserv 19:3809–25. 488 Esteve-Selma MA, Martínez-Fernández J, Hernández I, Robledano F, Pérez-Navarro MA, Lloret 489 F. 2015. Cambio climático y biodiversidad en el contexto de la Región de Murcia. In: 490 Consejería de Agua Agricultura y Medio Ambiente, editor. Cambio climático en la Región 491 de Murcia. Evaluación basada en indicadores. Murcia. pp 105–32. 492 Fielding AH, Bell JF. 1997. A review of methods for the assessment of prediction errors in 493 conservation presence/absence models. Environ Conserv 24:38–49. 494 Franklin J. 2010. Mapping species distributions. Spatial inference and prediction. New York: 495 Cambridge University Press. 496 Franklin J, Davis FW, Ikegami M, Syphard AD, Flint LE, Flint AL, Hannah L. 2013. Modeling plant 497 species distributions under future climates: how fine scale do climate projections need to 498 be? Glob Chang Biol 19:473–83. 499 Franklin J, Serra-Diaz JM, Syphard AD, Regan HM. 2016. Global change and terrestrial plant 500 community dynamics. Proc Natl Acad Sci 113:3725–34. 501 Freckleton RP, Harvey PH, Pagel M. 2002. Phylogenetic Analysis and Comparative Data : 21 502 160:712–26. 503 Greenwood S, Ruiz-Benito P, Martínez-Vilalta J, Lloret F, Kitzberger T, Allen CD, Fensham R, 504 Laughlin DC, Kattge J, Bönisch G, Kraft NJB, Jump AS. 2017. Tree mortality across biomes 505 is promoted by drought intensity, lower wood density and higher specific leaf area. Ecol 506 Lett 20:539–53. 507 Guiot J, Cramer W. 2016. Climate change: The 2015 Paris Agreement thresholds and 508 Mediterranean basin ecosystems. Science (80- ) 354:4528–32. 509 Hamerlynck EP, McAuliffe JR. 2008. Soil-dependent canopy die-back and plant mortality in two 510 Mojave Desert shrubs. J Arid Environ 72:1793–802. 511 Hijmans ARJ, Phillips S, Leathwick J, Elith J, Hijmans MRJ. 2016. Package ‘ dismo ’ Species 512 Distribution Modeling. 513 Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. 2005. Very high resolution interpolated 514 climate surfaces for global land areas. Int J Climatol 25:1965–78. 515 Hirzel AH, Le Lay G, Helfer V, Randin C, Guisan A. 2006. Evaluating the ability of habitat 516 suitability models to predict species presences. Ecol Modell 199:142–52. 517 Huberty CJ. 1994. Applied discriminant analysis. New York: Wiley 518 IPCC, 2013. Climate change 2013: the physical science basis. contribution of working group I to 519 the fifth assessment report of the intergovernmental panel on climate change. In: Stocker 520 TF, Qin D, Plattner G-K, Tignor M, Allen SK, Boschung J, Nauels A, Xia Y, Bex V, Midgley 521 PM, Eds. Cambridge University Press, Cambridge, United Kingdom and New York, NY, 522 USA, 1535 pp. 523 Jump AS, Hunt JM, Peñuelas J. 2006. Rapid climate change-related growth decline at the 524 southern range edge of Fagus sylvatica. Glob Chang Biol 12:2163–74. 22 525 de la Riva EG, Olmo M, Poorter H, Ubera JL, Villar R. 2016a. Leaf Mass per Area (LMA) and Its 526 Relationship with Leaf Structure and Anatomy in 34 Mediterranean Woody Species along 527 a Water Availability Gradient. PLoS One 11:e0148788. 528 de la Riva EG, Pérez-Ramos IM, Tosto A, Navarro-Fernández CM, Olmo M, Marañón T, Villar R. 529 2016b. Disentangling the relative importance of species occurrence, abundance and 530 intraspecific variability in community assembly: a trait-based approach at the whole-plant 531 level in Mediterranean forests. Oikos 125:354–63. 532 Lázaro R, Rodrigo FS, Gutiérrez L, Domingo F, Puigdefábregas J. 2001. Analysis of a 30-year 533 rainfall record (1967–1997) in semi–arid SE Spain for implications on vegetation. J Arid 534 Environ 48:373–95. 535 Lenoir J, Graae BJ, Aarrestad PA, Alsos IG, Armbruster WS, Austrheim G, Bergendorff C, Birks 536 HJB, Bråthen KA, Brunet J, Bruun HH, Dahlberg CJ, Decocq G, Diekmann M, Dynesius M, 537 Ejrnaes R, Grytnes J-A, Hylander K, Klanderud K, Luoto M, Milbau A, Moora M, Nygaard B, 538 Odland A, Ravolainen VT, Reinhardt S, Sandvik SM, Schei FH, Speed JDM, Tveraabak LU, 539 Vandvik V, Velle LG, Virtanen R, Zobel M, Svenning J-C. 2013. Local temperatures inferred 540 from plant communities suggest strong spatial buffering of climate warming across 541 Northern Europe. Glob Chang Biol 19:1470–81. 542 Lesica P, Crone EE. 2016. Arctic and boreal plant species decline at their southern range limits 543 in the Rocky Mountains. Ecol Lett:166–74. 544 Lloret F, Escudero A, Iriondo JM, Martínez-Vilalta J, Valladares F. 2012. Extreme climatic events 545 and vegetation: The role of stabilizing processes. Glob Chang Biol 18:797–805. 546 Lloret F, García C. 2016. Inbreeding and neighbouring vegetation drive drought-induced die-off 547 within juniper populations. Field K, editor. Funct Ecol 30:1696–704. 23 548 Lloret F, Kitzberger T. 2018. Historical and event-based bioclimatic suitability predicts regional 549 forest vulnerability to compound effects of severe drought and bark beetle infestation. 550 Glob Chang Biol. 551 Lloret F, de la Riva EG, Pérez-Ramos IM, Marañón T, Saura-Mas S, Díaz-Delgado R, Villar R. 552 2016. Climatic events inducing die-off in Mediterranean shrublands: are species 553 responses related to their functional traits? Oecologia 180:961–73. 554 van der Maaten E, Hamann A, van der Maaten-Theunissen M, Bergsma A, Hengeveld G, van 555 Lammeren R, Mohren F, Nabuurs G-J, Terhürne R, Sterck F. 2017. Species distribution 556 models predict temporal but not spatial variation in forest growth. Ecol Evol 7:2585–94. 557 Martinez-Meyer E, Diaz-Porras D, Peterson AT, Yanez-Arenas C. 2013. Ecological niche 558 structure and rangewide abundance patterns of species. Biol Lett 9:20120637–20120637. 559 Martínez-Vilalta J, Lloret F. 2016. Drought-induced vegetation shifts in terrestrial ecosystems: 560 The key role of regeneration dynamics. Glob Planet Change 144:94–108. 561 Mcdowell N, Pockman WT, Allen CD, Breshears DD, Cobb N, Kolb T, Plaut J, Sperry J, West A, 562 Williams DG, Yepez EA, Mcdowell N, Pockman WT, Allen CD, David D, Mcdowell N, Cobb 563 N, Kolb T, Plaut J, Sperry J. 2008. Mechanisms of Plant Survival and Mortality during 564 Drought : Why Do Some Plants Survive while Others Succumb to Drought ? New Phytol 565 178:719–39. 566 Merow C, Smith MJ, Edwards TC, Guisan A, Mcmahon SM, Normand S, Thuiller W, Wüest RO, 567 Zimmermann NE, Elith J. 2014. What do we gain from simplicity versus complexity in 568 species distribution models? Ecography (Cop) 37:1267–81. 569 Mouillot F, Rambal S, Joffre R. 2002. Simulating climate change impacts on fire frequency and 570 vegetation dynamics in a Mediterranean ecosystem. Glob Chang Biol 8:423–37. 24 571 Ninyerola M, Pons X, Roure JM. 2000. A methodological approach of climatological modelling 572 of air temperature and precipitation through GIS techniques. Int J Climatol 20:1823–41. 573 Parmesan C, Yohe G. 2003. A globally coherent fingerprint of climate change impacts across 574 natural systems. Nature 421:37–42. http://www.ncbi.nlm.nih.gov/pubmed/12511946 575 Pearson RG, Dawson TP. 2003. Predicting the Impacts of Climate Change on the Distribution of 576 Species : Are Bioclimate Envelope Models Useful ? Glob Ecol Biogeogr 12:361–71. 577 Pearson RG, Thuiller W, Araújo MB, Martinez-Meyer E, Brotons L, McClean C, Miles L, 578 Segurado P, Dawson TP, Lees DC. 2006. Model-based uncertainty in species range 579 prediction. J Biogeogr 33:1704–11. 580 Peñuelas J, Lloret F, Montoya R. 2001. Severe drought effects on mediterranean woody flora in 581 Spain. For Sci 47:214–8. 582 Peterson AT (Andrew T, Soberón J, Pearson RG, Anderson RP, Martinez-Meyer E, Nakamura M, 583 Araújo MB. 2011. Ecological Niches and Geographic Distributions. Princeton. 584 Phillips SJ, Dudík M. 2008. Modeling of species distributions with Maxent: new extensions and 585 a comprehensive evaluation. Ecography (Cop) 31:161–75. 586 Pironon S, Villellas J, Morris WF, Doak DF, García MB. 2015. Do geographic, climatic or 587 historical ranges differentiate the performance of central versus peripheral populations? 588 Glob Ecol Biogeogr 24:611–20. 589 Pulliam HR. 2000. On the relationship between niche and distribution. Ecol Lett 3:349–61. 590 Ridgeway G. 2007. Generalized Boosted Models : A guide to the gbm package. Compute 1:1– 591 12. 592 Rivas-Martínez S, Rivas-Sáenz S, Penas-Merino A. 2011. Worldwide bioclimatic classification 25 593 system. Glob Geobot 1:1–638. 594 Sapes G, Serra-Diaz JM, Lloret F. 2017. Species climatic niche explains drought-induced die-off 595 in a Mediterranean woody community. Ecosphere 8:e01833. 596 Serra-Diaz JM, Keenan TF, Ninyerola M, Sabaté S, Gracia C, Lloret F. 2013. Geographical 597 patterns of congruence and incongruence between correlative species distribution 598 models and a process-based ecophysiological growth model. J Biogeogr 40:1928–38. 599 Sexton JP, Hangartner SB, Hoffmann AA. 2014. Genetic isolation by environment or distance: 600 which pattern of gene flow is most common? Evolution (N Y) 68:1–15. 601 Sexton JP, Mcintyre PJ, Angert AL, Rice KJ. 2009. Evolution and Ecology of Species Range Limits. 602 Annu Rev Ecol Evol Syst 40:415–36. 603 Sheffield J, Wood ÆEF. 2008. Projected changes in drought occurrence under future global 604 warming from multi-model , multi-scenario , IPCC AR4 simulations. Clim Dyn 31:79–105. 605 Soberón J. 2007. Grinnellian and Eltonian niches and geographic distributions of species. Ecol 606 Lett 10:1115–23. 607 Thomas CD, Cameron A, Green RE, Bakkenes M, Beaumont LJ, Collingham YC, Erasmus BFN, de 608 Siqueira MF, Grainger A, Hannah L, Hughes L, Huntley B, van Jaarsveld AS, Midgley GF, 609 Miles L, Ortega-Huerta MA, Townsend Peterson A, Phillips OL, Williams SE. 2004. 610 Extinction risk from climate change. Nature 427:145–8. 611 Thuiller W. 2004. Patterns and uncertainties of species’ range shifts under climate change. 612 Glob Chang Biol 10:2020–7. 613 Thuiller W, Cile C, Albert H, Dubuis A, Randin C, Guisan A. 2010. Variation in habitat suitability 614 does not always relate to variation in species’ plant functional traits. Biol Lett 6:120–3. 26 615 Thuiller W, Lavorel S, Araújo MB. 2005. Niche properties and geographical extent as predictors 616 of species sensitivity to climate change. Glob Ecol Biogeogr 14:347–57. 617 Valladares F, Benavides R, Rabasa SG, Diaz M, Pausas JG, Paula S, Simonson WD. 2014. Global 618 change and Mediterranean forests : current impacts and potential responses. In: Forests 619 and Global Change. pp 47–75. 620 Valladares F, Vilagrosa A, Peñuelas J, Ogaya R, Julio J, Corcuera L, Sisó S. 2004. CAPÍTULO 6 621 Estrés hídrico : ecofisiología y escalas de la sequía. In: Valladares F, editor. Ecologia del 622 bosque mediterráneo en un mundo cambiante. Ministerio. Madrid. pp 163–90. 623 Vicente-Serrano SM, Gouveia C, Camarero JJ, Beguería S, Trigo R, López-Moreno JI, Azorín- 624 Molina C, Pasho E, Lorenzo-Lacruz J, Revuelto J, Morán-Tejeda E, Sanchez-Lorenzo A, 625 Dickinson RE. 2013. Response of vegetation to drought time-scales across global land 626 biomes. Proc Natl Acad Sci U S A 110:52–7. 627 Walther G-R, Roques A, Hulme PE, Sykes MT, Pyšek P, Kühn I, Zobel M, Bacher S, Botta-Dukát 628 Z, Bugmann H, Czúcz B, Dauber J, Hickler T, Jarosík V, Kenis M, Klotz S, Minchin D, Moora 629 M, Nentwig W, Ott J, Panov VE, Reineking B, Robinet C, Semenchenko V, Solarz W, 630 Thuiller W, Vilà M, Vohland K, Settele J. 2009. Alien species in a warmer world: risks and 631 opportunities. Trends Ecol Evol 24:686–93. 632 Webb CO, Donoghue MJ. 2005. Phylomatic: Tree assembly for applied phylogenetics. Mol Ecol 633 Notes 5:181–3. 634 Weber MM, Stevens RD, Diniz-Filho JAF, Grelle CE V. 2016. Is there a correlation between 635 abundance and environmental suitability derived from ecological niche modelling? A 636 meta-analysis. Ecography (Cop) 40:817–28. 637 Wood SN. 2011. Fast stable restricted maximum likelihood and marginal likelihood estimation 27 638 of semiparametric generalized linear models. J R Stat Soc 73:3–36. 639 Wright JW, Davies KF, Lau JA, Mccall A, Mckay JK. 2006. Experimental verification of ecological 640 niche modeling in a heterogeneous environment. Ecology 87:2433–9. 641 Zavaleta ES, Shaw MR, Chiariello NR, Thomas BD, Cleland EE, Field CB, Mooney HA. 2003. 642 Grassland Responses To Three Years of Elevated Temperature, Co 2 , Precipitation, and N 643 Deposition. Ecol Monogr 73:585–604. 644 Zimmermann NE, Yoccoz NG, Edwards TC, Meier ES, Thuiller W, Guisan A, Schmatz DR, 645 Pearman PB. 2009. Climatic extremes improve predictions of spatial patterns of tree 646 species. Proc Natl Acad Sci U S A 106 Suppl:19723–8. 647 Zunzunegui M, Díaz Barradas MC, Ain-Lhout F, Clavijo A, García Novo F. 2005. To live or to 648 survive in Doñana dunes: Adaptive responses of woody species under a Mediterranean 649 climate. Plant Soil 273:77–89. 650 28 651 Table captions 652 Table 1 Main settings used in the different SDM modeling approaches. 653 Table 2 Main species information and attributes used in statistical analyses. See methods for 654 details about foliar strategy, RGC (Remaining Green Canopy), Size and Frequency. 655 Table 3 Results of GLMs explaining remaining green canopy as a function of foliar strategy, size, 656 frequency, HCS (Historical Climatic Suitability), ECS (Episode Climatic Suitability) and the 657 interaction between the latter two (HCS:ECS) calculated from four different SDMs (Mahalabonis 658 distance, GAM, BRT, MaxEnt). AIC stepwise selection was applied to obtain the final models. 659 Significant or marginally significant P values are highlighted in bold. 29 660 Method Key reference and main settings Mahal adehabitatHS package (Calenge 2015). Modification of the original function Mahalanobis distance mahasuhab in order to obtain distance and probability values from other layers not used to calculate variables’ mean vector. This function makes it possible to determine probability values, assuming that under multivariate normality, squared Mahalanobis distance is approximately distributed as Chi-square with n- 1 degrees of freedom, which makes it possible to calculate p-value maps (Clark and others 1993). GAM mgcv package v.1.8-16 (Wood 2011). Weighted background number: 100,000 Generalized Additive points. The optimal number of edf for each variable was selected between 1 and Models 4 by cross-validation, using gam function. The number of knots for those species where response curves being biologically counterintuitive was also reduced. BRT gbm package v. 2.1.1 (Ridgeway 2007). Pseudoabsence number equivalent to Boosted Regression Tree presence species data. Tree complexity of 5 for those species with more than 250 occurrences and 3 for those species with less than 250. The learning rate of 0.005 was chosen because it made it possible to achieve at least 1,000 trees in every case, following Elith and others (2008). MaxEnt MaxEnt v. 3.3.1 (Phillips and Dudík 2008) used with default setting with the Maximum entropy exception of: 100,000 background points, 10-fold cross-validation, regularization multiplier of 3, and threshold feature unselect in order to produce smoother response curves. 661 662 30 Family Species Code Life form Foliar strategy RGC (%) Size (cm2) Frequency Anacardiaceae Pistacia lentiscus PLE Microphanerophyte Evergreen 25.33 0.0828 0.5 Asparagaceae Asparagus horridus AHO Chamaephyte Leafless 15.00 0.0291 0.1 Asteraceae Artemisia barrelieri ABA Chamaephyte Semideciduous 44.29 0.0043 0.9 Artemisia campestris ACA Chamaephyte Semideciduous 13.33 0.0115 0.1 Boraginaceae Lithodora fruticosa LFR Nanophanerophyte Semideciduous 47.27 0.0291 0.1 Chenopodiaceae Salsola genistoides SGE Nanophanerophyte Leafless 6.67 0.0795 0.2 Cistaceae Cistus albidus CAL Nanophanerophyte Semideciduous 77.50 0.0452 0.2 Cistus clusii CCL Nanophanerophyte Semideciduous 44.77 0.0456 0.9 Helianthemum syriacum HSY Chamaephyte Semideciduous 61.39 0.0018 0.3 Cupressaceae Juniperus oxycedrus JOX Microphanerophyte Evergreen 56.69 0.5503 0.6 Juniperus phoenicea JPH Microphanerophyte Evergreen 60.00 0.0736 0.1 Fagaceae Quercus coccifera QCO Microphanerophyte Evergreen 33.10 0.0091 0.7 Lamiaceae Rosmarinus officinalis ROF Nanophanerophyte Semideciduous 53.00 0.2353 1.0 Sideritis leucantha SLE Nanophanerophyte Semideciduous 44.77 3.0292 0.2 Teucrium capitatum TCA Chamaephyte Semideciduous 68.00 2.0186 0.6 Thymus hyemalis THY Chamaephyte Semideciduous 45.28 0.3894 1.0 Leguminosae Anthyllis cytisoides ACY Nanophanerophyte Semideciduous 21.46 0.3779 0.6 Dorycnium DPE Nanophanerophyte Semideciduous 21.54 0.5769 0.5 pentaphyllum Ononis fruticosa OFR Nanophanerophyte Semideciduous 10.11 0.4742 0.4 Poaceae Stipa tenacissima STE Hemicryptophyte Evergreen 61.89 0.0088 1.0 Rhamnaceae Rhamnus lycioides RLY Microphanerophyte Semideciduous 8.75 0.0560 0.3 Timeleaceae Daphne gnidium DGN Nanophanerophyte Semideciduous 15.00 0.0045 0.2 663 664 31 665 MAHALANOBIS DISTANCE GAM BRT MAXENT AIC 197.97 AIC 196.10 AIC 187.33 AIC 194.15 R2 adj 0.32 R2 adj 0.33 R2 adj 0.55 R2 adj 0.39 β Std.Error P.value β Std.Error P.value β Std.Error P.value β Std.Error P.value INTERCEPT 0.000 0.000 0.939 0.000 0.000 0.610 0.000 0.000 0.324 0.000 0.000 0.051 (FOLIAR.STR 1) HCS 0.382 0.213 0.096 0.362 0.196 0.083 0.561 0.157 0.003 0.446 0.199 0.039 FOLIAR.STR 2 0.062 0.299 0.838 -0.317 0.241 0.207 -0.317 0.198 0.130 -0.364 0.234 0.139 FOLIAR.STR 3 -0.435 0.230 0.081 -0.563 0.210 0.016 -0.618 0.172 0.002 -0.706 0.223 0.006 FOLIAR.STR 4 -0.254 0.254 0.336 -0.688 0.218 0.006 -0.690 0.173 0.001 -0.628 0.199 0.006 LOG(SIZE) -0.244 0.212 0.271 -0.318 0.201 0.133 -0.373 0.167 0.039 -0.330 0.1934 0.107 FREQUENCY 0.369 0.225 0.125 ECS 27.618 18.640 0.162 HCS:ECS -27.360 18.590 0.165 666 667 32 668 Figure captions 669 Figure 1 Study region within the Iberian Peninsula (Murcia region: black shading in the inset, 670 upper left map) and monthly temperature (lines) and precipitation (bars) data during the 671 historical 1950-2000 period (light grey bars and dots) and the 2013-2014 hydrological year 672 (anomaly period; dark grey bars and triangles). 673 Figure 2 Averaged proportion of Remaining Green Canopy (RGC) in shrub species with foliar 674 strategies 1 and 2 (top graph) for the ten studied plots. Species’ RGC values are sorted in 675 increasing order (x-axis) and error bars are shown. Red and blue dots represent foliar strategies. 676 The lower graph shows the Historical Climatic Suitability (HCS) and the drought Episode Climatic 677 Suitability (ECS) values for each species. Blue and red dots represent median suitability values 678 of HCS and ECS, respectively, and the error bars correspond to the range between maximum 679 and minimum suitability values predicted by the four applied SDM algorithms (see text for 680 details). 681 Figure 3 Partial residual plot of RGC (Remaining Green Canopy) in relation to HCS (Historical 682 Climatic Suitability) obtained for each SDM model. Β (standardized estimate value) and P values 683 for HCS in these models are shown in the left corner of each plot. Species codes are shown in 684 Table 2. 685 33 686 687 34 688 689 35 690 36