This is the accepted version of the journal article: Lerner, David; Fernández-Martínez, Marcos; Livne- Luzon, Stav; [et al.]. «A biome-dependent distribution gradient of tree species range edges is strongly dictated by climate spatial heterogeneity». Nature Plants, Vol. 9, issue 4 (April 2023), p. 544–553. DOI 10.1038/s41477-023-01369-1 This version is available at https://ddd.uab.cat/record/284235 under the terms of the license 10 11 A biome-dependent distribution gradient of tree species range edges is 12 strongly dictated by climate spatial heterogeneity 13 David Lerner1, Marcos Fernández Martínez2, Stav Livne-Luzon1, Jonathan Belmaker3,4, 14 Josep Peñuelas2,5, Tamir Klein1 15 1 Department of Plant & Environmental Sciences, Weizmann Institute of Science. Rehovot 16 76100, Israel. 17 2 CREAF, Cerdanyola de Valles, Catalonia E08193, Spain. 18 3 School of Zoology, Faculty of Life Sciences, Tel Aviv University, Tel Aviv, Israel. 19 4 Steinhardt Museum of Natural History, Tel Aviv University, Tel Aviv, Israel. 20 5 CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Catalonia E08193, Spain 21 22 Abstract 23 Understanding the causes of the arrest of species distributions has been a fundamental 24 question in ecology and evolution. These questions are of particular interest for trees due to 25 their long lifespan and sessile nature. A surge in data-availability evokes for a macro- 26 ecological analysis to determine the underlying forces limiting distributions. Here we analyse 27 the spatial distribution of >3600 major tree species to determine geographical areas of range- 28 edge hotspots and find drivers for their arrest. We confirmed biome edges to be strong 29 delineators of distributions. Importantly, we identified a stronger contribution of temperate 30 than tropical biomes to range edges, adding strength to the notion that tropical areas are 31 centers of radiation. We subsequently identified a strong association of range-edge hotspots 32 with steep spatial climatic gradients. We linked spatial and temporal homogeneity and high 33 potential evapotranspiration in the tropics as the strongest predictors for this phenomenon. 34 We propose that the poleward migration of species in light of climate change might be 35 hindered due to steep climatic gradients. 36 37 Main Text 38 Introduction 39 The geographical distributions of species are marked by their range limits. Understanding the 40 causes of distribution arrest has been a fundamental question in ecology and evolution1–4. 41 Given the strong interplay between biotic, abiotic, demographic, physical and historical 42 forces in predicting range-edges, it has been challenging to find underpinnings for their 5 43 formation. Two main environmental forces seem to play a major role in the formation of 44 range edges; spatial environmental heterogeneity and habitat quality4,54,5. Most models on the 45 formation of species-range edges rely on the interplay between either one of these two forces 46 with non-climatic pressures to explain their formation. For example, steep climate gradients 47 combined with high dispersal and gene flow reduces species’ fitness and genetically 48 constrains their evolution into novel environments6,7. Likewise, low habitat quality reduces 49 population size2,8, increasing drift and migration load9–11. Nevertheless, the significance of 50 climate in the interplay between these two environmental components in defining range edges 51 for an array of species or on a wide biogeographical scale remains elusive12,13. 52 53 The field of biogeography has long sought understanding of species ranges despite having 54 limited tools (e.g. 3,14). Given the surge of large-scale datasets, it is now possible to better 55 identify the underpinnings of species distributions by studying the macro eco-evolutionary 56 processes involved in their formation11,12,15–17. Although methods and results are disparate 57 between studies, there is almost a consensus that the presence of large-scale biogeographical 58 units confine species with climate as their primary predictor. For example, Bontager et al.11 59 suggested distinct characteristics for range edge populations dependent on their latitude. 60 Likewise, niche conservatism and strong beta-diversity patterns seem to withhold at large 61 macro-ecological scales16–18. This is generally true for plant species, with biomes being the 62 most consistent classifier based on structural and functional similarity19,20. Recent efforts 63 have been made to understand how accurate and substantial biome entities are at defining 64 species distributions 20,21. Nevertheless, it remains an open question if, and to what extent, the 65 intersection between biomes is a source of species range edge hotspots. Deciphering such 66 patterns will enable the proper understanding of how communities of species redistribute and 67 are structured geographically, and if similar biomes in distinct geographic areas have similar 68 effects on the distribution of species. 69 70 However, although there has been an increased interest in defining the biogeographical 71 underpinnings of species distributions, most techniques have used species relatedness and 72 diversity metrics to test for the existence of shared niche space between species and 73 communities. The direct analysis of the ecological and climatic limitations to geographic 74 space, although trivial, remain elusive. Here, we look at the universal set of climatic factors 75 and geographical patterns of species distributions by focusing directly on species’ range edge 76 distributions. 6 77 78 We present the first-to-date global study of tree species range edges, applying a novel, 79 simple, yet effective method of delineating range edges (REs) in order to: (1) identify 80 deterministic patterns of REs, as seen by RE-dense areas (RE hotspots: REH); (2) determine 81 whether the classification of biomes as distinct community-level patterns of biodiversity 82 properly delineate the niche of species; (3) identify global-scale REH patterns, and; (4) 83 discern the underlying niche factors responsible for RE formation. In particular, we 84 investigate whether spatial heterogeneity of abiotic factors or a universal predictor for habitat 85 quality are determinants of RE formation. We focus specifically on tree species, as they are 86 an exemplary group of species in the study of the ecological changes predicted to occur at the 87 peripheries of distributions, given their long-lived characteristics and fundamental role in 88 many ecosystems25,26, specifically biomes 18,27. 89 90 Although we do not study the interplay between climate and other ecological and 91 evolutionary limiting factors to species distributions (e.g., seed dispersal, plasticity, and 92 adaptation), discerning these patterns and the climatic components leading to such 93 distributions will enable a better understanding of interplays between biotic and abiotic 94 factors in future studies. A better understanding of the climatic factors affecting dispersal 95 enable for better predictions of the success of species to track changing climates, and in turn, 96 if they will be subject of migration lags22–24. 97 98 Results 99 Global data set and REH distribution 100 We present the first report of the global distribution of tree range edge hotspots (REH) (Fig. 101 1), marked by the hexagons with significant clustering of REs (Supplementary Fig. 1). We 102 did not identify any significant coldspots, given the baseline presence of REs around the 103 globe. The visual patterns emerging from these distributions indicated that distributions 104 stopped disproportionally more at the edges of biomes than within them. For example, 105 northern REH occurred mostly at the intersection between a montane (Himalayas) and a 106 desert biome (Gobi Desert) or at the edge of the tundra in North America. Southern REH in 107 Africa and southern Eurasia tended to stop at the edges of desert biomes (e.g. Sahara Desert) 108 but stopped mainly at the edges of temperate and montane grasslands in the Southern 109 Hemisphere (pampas and Andes, respectively). Eastern and western REH were notably 7 110 mostly in similar geographical locations, e.g. at the intersections between the Himalayas and 111 central Asian deserts or at the edge of the Atacama Desert in South America. 112 113 In order to identify the underlying niche factors that define REs, we focused only on inland 114 REs, since REs at the edge of a water source is probably due to an obvious geographical 115 barrier rather than an ecological effect. The fraction of inland REs was globally similar 116 between the continents (between 65 and 75%), except for Australia (~45%, Supplementary 117 Fig. 2a). The large fraction of inland REs (Supplementary Fig. 2b) was mainly associated 118 with the edges of biomes (e.g. northern REs in Africa, Asia, and North America were 119 associated with the Sahara Desert, montane grasslands of the Himalayas, and the tundra 120 biome, respectively; Fig. 1). A significant fraction of eastern inland REs were in Europe, but 121 almost no significant REH were identified, implying the sparse distribution of REs 122 throughout the continent, or as an effect from a smaller area for distribution compared to 123 other continents. 124 125 Contribution of biome-biome intersections to REHs 126 We quantitatively identified the global patterns of arrest by analysing the fraction of REH 127 that stopped at biome-biome intersections (14 central global ecological regions best 128 distinguished by their climate, fauna, and flora obtained from World Wild Life (WWF, 129 http://www.worldwildlife.org/; Fig. 1) and by delineating a buffer zone at points of 130 intersection between two or more biomes. Although not significant, the number of REH was 131 strongly associated with the intersection between biomes (Fig. 2a, first panel), indicating that 132 climatic conditions were a probable cause for the REs of tree species at biome edges. Our 133 results, however, identified an unequal contribution of the different biome edges to the 134 fraction of REH (Fig 2a). We then analysed the individual biome-biome intersections 135 normalised by a global permutation (i.e., from permutation of the global distribution of REH, 136 see Materials and Methods). Here we identified (1) a strong contribution of REH at the 137 intersections between temperate and desert biomes (Fig. 2b, left panel) in comparison to (2) a 138 weaker contribution at the intersections of tropical and subtropical biomes (between 139 themselves and with temperate biomes). The significance of the contribution of REH 140 (p<0.05) (see Materials and Methods) was attributed almost exclusively to the intersections 141 within temperate biomes (asterisks in Fig. 2b). This unequal contribution between the 142 temperate vs tropical and subtropical biomes was globally consistent. Nevertheless, under a 8 143 per-biome permutation (i.e., normalising each intersect by a selective permutation from the 144 respective biome combinations separately; Fig 2b right panel), we observed a much weaker 145 contribution of range edges to the formation of biomes. Only a selected number of biome- 146 biome intersections had a significant contribution to the formation of REH, yet no specific 147 pattern between temperate and tropical biomes was observed. Indeed, the strong positive 148 correlation between the number of REs between biomes with the number of REs within 149 biomes (Fig. 2c; R = 0.9, linear regression, p < 0.001) reflects such discrepancy between 150 panels in Fig 2b, as the number of REH at the edge of biomes had a strong linear association 151 with the number of REH within that biome. We identified five biomes, however, outside the 152 95% confidence interval (CI) of the regression (Fig. 2c). The two biomes above the 153 regression CI (desert; and boreal) were biomes where the number of REH at the edges was 154 larger than predicted (Fig 2c). In contrast, temperate coniferous forest, mangrove and flooded 155 grasslands fell below the CI of the regression, indicating a much larger number of REH 156 within the biome compared to the edge. 157 158 A parallel analysis using the distribution of REs (rather than REH) was conducted to identify 159 similarities and differences between the distributions. The results indicated a similar pattern 160 of distribution (Supplementary Fig. 4a-b compared to Fig. 2b-c), suggesting global forces 161 associated with REH. 162 163 Climatic predictors of RE formation 164 Intersections between biomes were a significant cause of RE formation, so we also 165 investigated the dependence of RE formation on climate. We tested both the ‘absolute 166 climate’ (i.e. annual and seasonal average temperature and precipitation) and the ‘spatial 167 heterogeneity of climate’ (SH, i.e. spatial variability of absolute climate) at each of the global 168 hexagons using 19 bioclimatic variables obtained from WorldClim (see Materials and 169 Methods). Elevation (absolute and SH) and latitude were also accounted for at each hexagon 170 resolution. Generalised linear models of regression between each climatic variable as an 171 independent predictor indicated that all 40 climatic variables were significant predictors of 172 RE formation (p<0.05) (Supplementary Fig5a). Interestingly, most absolute climatic 173 variables were negatively associated with the global RE distribution (Fig. 3a and b), implying 174 a general prediction of REs occurring in climates with low temperatures and low 9 175 precipitation. A positive association was attributed to the four climatic variables that defined 176 the temporal heterogeneity of temperature and precipitation (mean diurnal range, temperature 177 seasonality temperature annual range and precipitation seasonality; BioClim variables 2, 4, 7, 178 and 15, respectively; Fig. 3a, Supplementary Fig. 5) and all of the SH climatic variables. A 179 mixed model with biomes and continents as random factors gave a reduced number of 180 variables that significantly predicted the formation of RE (Fig. 3a). A crossed model with 181 both continents and biomes consistently gave a stronger fit (AIC values) than only 182 considering either of these categorical random factors independently (Supplementary Fig 6b). 183 In this more stringent global analysis, latitude, absolute temperature and SH temperature, 184 precipitation and elevation were significant predictors of REH. The goodness of fit (R2) of 185 each model, as a measure of predictive strength, indicated that spatial heterogeneity 186 accounted for RE formation better than did their absolute equivalents. A model selection was 187 carried out to identify the most important factors associated with RE formation, followed by a 188 model averaging of the models with a ∆AIC < 2, (Figs. 3c). The SH climatic variables again 189 defined RE better than the absolute climatic variables. Although absolute climatic variables 190 such as isothermality (BioClim3), temperature of the wettest quarter (BioClim8) (both 191 characteristic of tropical and subtropical climates), and annual precipitation (BioClim12) 192 were strongly associated with REs in a generalised linear model, temperature of the warmest 193 quarter (BioClim10) was the only predictor strongly associated with RE under the mixed 194 model (Figs. 3a and S6). This difference in results can be visually seen when comparing 195 between continents (Fig. 3b) e.g. in panel 3; with temperature of wettest quarter (BioClim8) 196 having partial dissociations in Africa and South America, i.e. the two continents with the 197 most tropical biomes. Furthermore, SH isothermality (BioClim3), temperature of coldest 198 month (SH BioClim 6), precipitation (SH BioClim13 and 16) and elevation change (SH 199 elevation) were all predictors of REs. Interestingly, the absolute climatic predictors 200 (BioClim8 and BioClim10) were strongly negatively correlated with spatial heterogeneity at 201 the minimum temperature of the coldest month (SH BioClim6) (Supplementary Fig. 7b), 202 indicating its representation of a tropical biome climate. The weaker relative importance of 203 SH precipitation (Fig. 3c, left panel) is due to its dissociation with the continents that most 204 strongly represent temperate regions; Europe and North America (Fig 3b, panel 4). 205 206 In parallel, we ran models with the ENVIREM dataset, a dataset of environmental variables 207 complementary to WorldClim that are more ecophysiologically meaningful for plant 208 species28. Most ENVIREM variables associated with REH are spatially heterogeneous 10 209 variables (Fig. 3c, Supplementary Fig. 7). Absolute potential evapotranspiration (PET) of the 210 coldest quarter was strongly correlated with both absolute temperature of the warmest quarter 211 (BioClim10) and SH temperature of the coldest month (BioClim6); the two predictors 212 indicative of a transition between tropical and temperate biomes (Supplementary Fig. 7b). SH 213 EmbergerQ and SH PET of warmest quarter were strongly associated with SH precipitation 214 (BioClim13 and 16) and SH temperature of the warmest quarter (BioClim10), respectively. 215 Albeit the strong correlation of all absolute ENVIREM variables with absolute BioClim 216 variables28, we identified several SH ENVIREM variables to be weakly correlated to BioClim 217 (e.g., SH moisture index and PET of driest month). All results in Fig. 3 were robust against 218 spatial autocorrelation. (Supplementary Fig. 5b, see Materials and Methods for further 219 information). 220 221 The REH distribution from our generated polygons was compared to a REH distribution from 222 expert polygons (see Materials and Methods), in order to test for the accuracy of our 223 generated dataset. Indeed, results of GLMs converge, indicating the robustness of our 224 generated dataset to describe REs (Supplementary Fig. S8). 225 226 Discussion 227 The results support our hypothesis of a nonrandom distribution of REs. We were able to 228 confirm that many tree species range edges were clustered rather than sparsely (stochastic) 229 distributed, by obtaining a larger number of significant REH. This finding suggests the 230 underlying presence of ecological and evolutionary forces governing REH formation. 231 Similarly, the matching results obtained from the biome analyses when accounting for REs 232 (Supplementary Fig. 4) or REH (Fig. 2) was also indicative of the deterministic clustering of 233 species REs at these specific ecological barriers (biome edges). Nevertheless, we identified 234 case-specific exceptions such as the scattered distribution of eastern REs throughout Europe, 235 indicated by the strong identification of internal REs (Supplementary Fig. 2b) but not REH 236 (Fig. 1c). This could potentially have occurred as a consequence of either the intense 237 anthropogenic activity throughout29,30 that prevented the distribution of tree species to reach 238 their natural REs, or the effect of smaller land area compared to other continents, altering the 239 effect of how species interact with abiotic factors, thus affecting their distribution and 240 adaptation. 241 11 242 Our results also indicated a strong dependence of REH on the edges of biomes, strongly 243 supporting the many efforts to determine whether the division of the planet into discrete 244 geographical units has been appropriately delineated21,31–34. Our findings, however, were 245 unexpectedly biome-specific, identifying the biome borders that most defined tree 246 distributions (Fig. 2). There is no seemingly obvious pattern of differential contribution to 247 REH when analysing each biome independently (per-biome bootstrap; Fig 2a). Yet, a clear 248 distinction between a tropical/temperate distribution of REH is strongly observed under a 249 biome-biome pairs analysis (global bootstrap; Fig 2b). Specifically, REs were strongly 250 dependent on desert, temperate, and montane edges compared to the weak dependence on 251 tropical and subtropical biome edges. Despite the marginally significant association of range 252 edges to biomes globally (Fig 2a), a tropical-temperate REH distribution is not specifically 253 dependent on biome edges, but rather, it represents a global biome trend (Fig 2b, left panel; 254 Fig. 2c). In addition, the importance of biomes in describing climate’s association to REs, as 255 seen from a mixed effects model (Fig. 3a, S6c), as well as the absence of latitude in 256 predicting REs in a global model (Fig 3c) is further indicative of the importance of biomes 257 (rather than a latitudinal effect) in defining REs. 258 259 The large differences in the relative number of REs between the tropical and temperate 260 biomes can indicate adaptive mechanisms between species residing in either of these two 261 types of biomes. These nontrivial results may have been related to the Latitude Diversity 262 Gradient35–37, a well-established pattern in which biodiversity is higher in the tropics than in 263 temperate regions. Firstly, the notion that these differences surged through the effective 264 evolutionary time hypothesis36,38 suggests that genetic diversity is higher in the tropics due to 265 the possible longer times needed for adaptation and expansion. Similarly, long-term climatic 266 oscillations have been suggested to reduce cladogenesis at higher latitudes, consistent with 267 the observation that tropical areas are centers of evolutionary novelty36,39,40 and evolve faster 268 than temperate regions35, leading to the notion of tropics as centers for radial expansion of 269 clades and species. Recent studies have identified mechanisms for this ‘out of the tropics 270 expansion’ model41,42 and have reported a higher fraction of bridge species (species violating 271 niche conservatism) from the tropics compared to temperate regions. We conclude that the 272 lower contribution to REs in tropical vs temperate regions is consistent with stronger radial 273 expansions from the tropics than temperate regions. 274 12 275 However, other mechanisms could explain the distinctive REH patterns between the tropic 276 and temperate regions. For example, based on several studies demonstrating the large 277 physiological and evolutionary effects of forest fragmentation29,43,44, the increased long-term 278 anthropogenic activity, and consequently excessive fragmentation, in forests in temperate 279 regions compared to tropical and subtropical regions may have hindered the adaptation of 280 species to novel climates. The lack of association of REH specifically to temperate 281 grasslands, temperate conifers and Mediterranean biomes specifically (Fig 2a) could be 282 indicative of such an effect, as these biomes have been historically subject to strong 283 anthropogenic activity. 284 285 Our analyses confirmed the strong dependence of RE on climate. Lower temperature, higher 286 climatic heterogeneity (both temporal and spatial), and elevation changes were the strongest 287 climatic predictors of REs. The analysis of the best predictive variables indicated a noticeable 288 weak correlation between most of the leading factors determining REs (Supplementary Fig. 289 7b), suggesting that these factors were site-specific predictors of REs (e.g. desert biome with 290 low annual precipitation, montane grassland with high elevation, and tropical and subtropical 291 biomes with temperature homogeneity (isothermality, temperature of the wettest quarter and 292 spatial heterogeneity of low temperature – SH6)). We suggest that the strong negative 293 correlation between these absolute climatic variables and the spatial heterogeneity of low 294 temperature is an indication of different temperature patterns between tropical and temperate 295 regions, i.e. a buffered heterogeneous temperature in the tropics in contrast to the latitudinal 296 effect of decreasing temperatures in the temperate biomes. There was an overall positive 297 trend for the effect of precipitation (BC12-BC19) when controlling for biome since the signal 298 appeared only when biome is included as a random factor (Supplementary Fig. 6). This 299 means there must be other non-accounted confounding factors altering the relationship when 300 analysing all biomes together, suggesting that the effect is not biome dependent. Although 301 this study does not make a note of the complex interplay between biotic and abiotic forces or 302 the ecological traits of tree species (e.g. seed dispersal, phenology) in the formation of RE, 303 we find a consistent global effect of temperature and spatial heterogeneity of temperature and 304 precipitation to predict their formation throughout the different models. This is indicative of 305 their principal role as universal predictors of species distributions. 306 307 The strong prediction of REH formation from PET of the coldest quarter and its covariance to 308 the tropical-temperate transition variables (SH BioClim6 and absolute BioClim10), could 13 309 reflect a possible mechanistic evolutionary constrain for the distribution of woody species’ at 310 temperate biomes. At higher latitudinal temperate biomes, where evapotranspiration is 311 strongly reduced, especially during the coldest quarter, there is a strong limitation to 312 photosynthesis and growth by the significant time reduction of stomatal conductance45,46. 313 Similarly, given the strength of its prediction of REH, SH embergerQ (pluviothermic 314 quotient) could indicate a more refined mechanism for RE formation than its covariates of 315 precipitation (SH BioClim 13 and 16). This index describes mean annual precipitation in 316 relation to annual changes in temperature. EmbergerQ thus increases the predictability of 317 how precipitation also dictates the formation of REH in the more temperate biomes (Fig 3b, 318 panel 4). The consistency of our results indicates spatial and temporal heterogeneity of 319 climate and topography as overwhelmingly stronger predictors of RE formation than their 320 absolute climatic counterparts (Fig. 3) and, in turn, frail evidence for a universal poor habitat 321 quality. Even in the cases where mean climatic variables strongly predict REH, these were 322 strongly associated with this transition between the spatial and temporal climatically 323 homogenous tropical and subtropical biomes to the more heterogeneous temperate biomes. 324 Nevertheless, we note the importance of a lack of evapotranspiration, particularly in cold 325 climates, as a main predictor of RE formation. 326 327 These observations have substantial implications for the effects of climate change on tree 328 distributions and its effects on tree migration. Although predictions for future steeper 329 temperature gradient as a result of greenhouse gas emission and climate change has not been 330 trivial47–49, such an increase in temperature gradients could vastly affect the distribution of 331 species. In particular, our results strengthen the growing understanding that the predicted 332 poleward migration of tree species might not be as successful as previously predicted23,24. The 333 increase of stronger spatial gradients (especially in the lower latitudes)24,48,50 or extreme and 334 spontaneous events might all be causes of migration lags, despite the suitable temperatures at 335 higher latitudes and altitudes. Likewise, the importance of PET from temperate biomes on the 336 formation of REH presented here could also suggest a possible migration lag or loss of 337 adaptation due to the predicted reduction in PET at higher latitudes51. Our results thus 338 highlight the importance of accounting for more precise spatial heterogeneity of climate as a 339 critical feature in future models of species distribution and the development of more precise 340 conservation efforts such as assisted migration. 341 342 Materials and Methods 14 343 Data acquisition and polygon formation 344 Supplementary Fig. 9 visually summarises the methodologies used to obtain global range 345 edge hotspot (REH) distributions. We downloaded a data set of tree species from the open- 346 source data set using R packages rgbif and taxize. Global Biodiversity Information Facility 347 (GBIF; 05 July 2021, https://doi.org/10.15468/dl.ajen6k) using the Botanic Gardens 348 Conservation list of 60000 tree species. We downloaded occurrences with entries from 1980 349 onwards, removing any occurrence reported with a geospatial issue, species not belonging to 350 the kingdom Plantae (in case of mismatched species names), and any occurrence marked as 351 unlikely, mismatched, or invalid. We removed occurrences that had reported uncertainties of 352 >100 km and records based on fossils and unknown sources. We then used the 353 CoordinateCleaner R package52 to remove any occurrences with zero coordinates, equal x 354 and y coordinates, duplicates, occurrences at sea, coordinates at capitals, and centroids. To 355 finalise, we again removed species with <300 occurrences. The data we used undoubtedly 356 contained sampling bias53, probably overrepresenting the number of REs in some regions 357 with a reduced or negligible sampling effort. We tried to overcome this issue by basing our 358 filtering steps on several previous studies52–54. The strength of the critical filtering steps 359 applied in our analysis resembled those previously presented53. 360 We converted the georeferenced species occurrences (x and y coordinates) into distributional 361 polygons in parallel using two independent techniques; through concave-hull 362 (Supplementary Methods) and multivariate kernel-density estimation (described below). 363 Given the strong similarity between the two methods (Supplementary Fig. 8), we discuss the 364 methods and results in detail only for the kernel-density estimated polygons. We created 365 polygons using two-dimensional kernel density estimations. We first divided the extent of all 366 the coordinates into 800 grid points in each dimension (longitude and latitude) in order to 367 produce a matrix of 640,000 grid cells for each of the species. Subsequently, we selected for 368 the grid cells with the highest 99% estimation of the species’ occurrence and subsequently 369 rasterised these. Polygons were then delineated around the contour of the rasters. 370 371 Polygon groupings 372 All polygons belonging to the same species were grouped based on absolute distance from 373 one another. Polygons separated by ≤500 km were grouped together, with the assumption that 374 fragmentation, gene flow, and unreported data could all warrant two nearby populations to be 375 considered as one. We used the hclust function (package stats, agglomeration method: 15 376 complete) to hierarchically cluster populations from a sequence of three or more populations 377 by their distances, also using a cutoff of 500 km (cutree function, package stats) for 378 determining the clusters. The final data set comprised >3600 tree species, ranging from one to 379 nine populations (polygons) per species, for a total of 8500 populations. All spatial data was 380 analysed using R packages sf and raster. 381 382 RE determination 383 RE dense areas were determined by (1) defining distinct global units, (2) identifying the RE 384 of each species and (3) map species’ RE to the global units to calculate the density of 385 REs/unit. In detail: (1) We rasterised the world map to spatially bin the density of REs. To 386 overcome the problem of spatial distortion, we used hexagonal bins with the dggridR R 387 package, developed using the ISEA Discrete Global Grids system, a repetition of polygons on 388 the surface of an icosahedron, allowing for the projection of equal sized bins onto a 2D plane. 389 We defined the size of each hexagon as ~23 000 km2 (with an average spacing between center 390 nodes of 165 km). (2) We used coordinates of the cardinal directions (north, south, east, and 391 west) to represent species REs by subdividing each polygon cluster into four quartiles in the 392 four cardinal directions (NE, NW, SE, and SW). The REs for each quartile were determined 393 as the two most-outward coordinates of the corresponding cardinal directions (e.g. north and 394 east cardinal coordinates for the NE quartile). Eight REs were thus determined for each 395 population (two for each cardinal direction). As a filtering step, we accounted for both REs 396 from the same cardinal direction if they were >20 arc-degrees apart, otherwise we only 397 accounted for the farthest point from the centroid. (3) The total number of REs obtained using 398 this method was normalised by the total number of species intersecting its respective 399 hexagon. In parallel, we also defined REs by accounting for the perimeter of the polygon for 400 each species (Supplementary Fig. 10). The ‘perimeter’ system may be a more realistic and 401 complete system for identifying REs, but the ‘cardinal coordinate’ system, although more 402 simplistic in nature, (1) provides a clearer visual representation of the distribution of REs, and 403 (2) allows for the directionality of REs to be compared, essential farther along the pipeline by 404 distinguishing between coastline and inland REs and identifying hotspots in the four cardinal 405 directions. 406 407 We then classified coastline and inland REs, assuming that the arrest of species distribution at 408 the edge of a water source was probably due to an obvious geographical barrier rather than an 409 ecological effect. Coastline REs were determined by creating a semicircle (buffer of 3 arc- 16 410 degrees) around each RE in the direction of its cardinal coordinate and measured the 411 percentage overlap with water. Cardinal coordinates with >50% overlap were considered a 412 “coastal RE”. 413 414 In order to find the probability distribution of REs at the edge of water source 415 (Supplementary Fig. 2), we permutated the global population of REs (except for the 416 Australian population) and used the mean of this permutation to compare to the number of 417 REs in each of the continents. 418 419 Hotspot analysis 420 Hotspots were identified using the Getis-Ord Gi* hotspot analysis55,56 to find spatial 421 correlations between hexagon (inland and normalised) RE densities. We initially compiled a 422 list of neighbors between all hexagons using the poly2nb function and then obtained the local 423 G statistic using the weighted density (normalised number of REs) of the global hexagons 424 and their relative distance from each other. The G statistic calculates a Z-score (measure of 425 standard deviation) for each hexagon. P values were then determined using the critical Z- 426 scores at 95% confidence levels followed by a Bonferroni correction using the p.adjustSP 427 function (using the number of neighbors between hexagons rather than the total number of 428 hexagons). All analyses were carried out using the R spdep package57. 429 430 We compared the analysis from the linear models with expert based polygons from three 431 different sources – IUCN, BIEN and EUFORGEN. A randomised weighted sample of all of 432 this dataset was used to generate a global distribution of REH by running this sample through 433 our pipeline. GLMs were run on the global distribution of expert-based REH in the same way 434 as with our generated polygons. Given the uneven distribution of expert-based polygons 435 globally (Supplementary Fig. 8a), we ran models excluding Asia and Africa, in order to 436 account for this bias. As seen by the strong similarity between the GLMs of expert polygons 437 in a global and filtered model (Supplementary Fig 8b and Fig 8c, panel 3), we observed an 438 overrepresentation of the expert polygons for these continents. Likewise, the GLMs from our 439 generated polygons are much similar to those obtained from the filtered model. In this case, 440 practically all variables showed the same relationship with range edges (either positive or 441 negative β values) as well as similar magnitudes. 442 443 Statistical analyses 17 444 Contribution of biome edge to RE. We used the 14 biomes defined by the World Wildlife 445 Fund (WWF) for our analyses. The distributions were downloaded from the WWF webpage 446 (http://www.worldwildlife.org/). To identify the intersection between biomes, we reduced the 447 complexity of the polygon edge using the rmapshaper package, which can perform 448 topologically aware polygon simplifications, thus maintaining the intersection between 449 biomes upon reduction of “edginess” of the polygons. The intersections were delineated and 450 subsequently enlarged (with a buffer distance of 0.1 arc-minute, ~185 m at the equator). A: 451 452 i. Global permutation assay was carried out by randomising (1000 iterations) the 453 global distribution of hexagons with REH (absolute Z-score) 454 ii. Per-biome stratified permutation was carried out by randomising the distribution 455 of hexagons with REH within each biome independently. 456 457 The averaged global bootstrap shown in panel 1 of Fig 2a, was calculated using a per-biome 458 bootstrap to obtain the probability distribution of REH at biome edges (Supplementary Fig. 459 3a). Distributions that were not normally distributed as a result of their small size (flooded 460 grassland, mangrove and tropical and subtropical coniferous forest; Supplementary Fig 3b), 461 were removed from the analysis. A general trend for the probability of range edges falling at 462 the intersection of biomes was therefore measured as a unified standardised z-distribution, 463 and compared to the median z-score from the actual percent overlap for each biome. 464 465 The density of hotspots at the intersection between biomes was calculated using the sum of 466 Z-scores of the hotspots at that intersection, and the percentage contribution was then 467 calculated using this value over the total Z-score at all biome intersections. The global or per- 468 biome 1000 permutation means were used as a normalising denominator for the values 469 obtained from our data set (Fig 2b, Supplementary Fig 4a). The denominator could be either 470 larger or smaller than the numerator, so we log-transformed the outcome to obtain a linear- 471 like relationship. Contribution within a biome was calculated the same way as for the 472 contribution at the edge, using the mean from a permutated assay to normalise for the 473 absolute value. 474 475 Climatic dependency of RE. We tested the relationship between RE density and climatic 476 features by assigning a set of environmental variables to each hexagon. We used the 477 bioclimatic attributes downloaded from WorldClim Global Climate Data58 at a resolution of 5 18 478 arc-minutes. The 19 BIOCLIM variables and elevations for each hexagon were extracted 479 using the R raster package. We also used the 16 ENVIREM variables described by 28, 480 downloaded from their website, at a resolution of 2.5 arc-minutes. Absolute climate for each 481 hexagon was obtained using the mean over all pixels. Climatic spatial heterogeneity (SH) 482 was calculated using the proportional variability index59,60 (PV) over all pixels in each 483 hexagon. 484 485 We used an array of linear mixed models (Fig 3, Supplementary Fig 5 and 6) to test for the 486 dependence of REs to climate and the robustness of the results. Linear mixed models were 487 carried out to account for biomes and continents. Both of these variables were introduced as 488 random effects in random intercept models. A model selection analysis was used to 489 determine the models that best predicted the formation of RE. Random intercept models 490 using both continent and biome as random variables were run, and models with ∆AIC <2 491 were selected for. Given the strong correlation between different predictor variables, we ran 492 models only with variable combinations that had a Pearson’s correlation value r < 0.7. The 493 relative contribution of the variable included in the model were calculated from the selected 494 models. Analyses were run using R package MuMIn61. 495 496 All statistical analyses (individual GLMs, and multiple-predictor GLMs) were tested for their 497 robustness to spatial autocorrelation by creating a spatial autocovariate (autocov_dist 498 function, spdep package), calculated as the distance-weighted average of neighboring 499 dependent variables62, so hexagons in proximity were averaged and those farther away 500 received a lower weighting average. We set the predetermined distance to 200 km based on 501 the average distance between cells. The spatial autocovariate was then included in the 502 regression model as a dependent variable. 503 504 Data Availability 505 The occurrence points used from GBIF can be found in the GBIF webpage 506 (https://doi.org/10.15468/dl.ajen6k). Polygons generated from occurrence points are provided 507 in the public Zenodo repository 10.5281/zenodo.7613535. Biome polygons were obtained 508 from the WWF webpage (http://www.worldwildlife.org/). Bioclimatic attributes were 509 downloaded from WorldClim Global Climate Data58. ENVIREM variables were downloaded 510 from their webpage (https://envirem.github.io/). 511 19 512 Code Availability 513 Custom codes related to this paper can be found in a GitHub repository at 514 https://github.com/dlernerg/Global-Range-edges 515 516 Acknowledgments 517 TK wishes to thank Edith and Nathan Goldenberg Career Development Chair; Mary and Tom 518 Beck-Canadian Center for Alternative Energy Research; Larson Charitable Foundation New 519 Scientist Fund; Yotam Project; Dana and Yossie Hollander; Estate of Emile Mimran; and the 520 Estate of Helen Nichunsky. DL was supported by the Sustainability and Energy Research 521 Initiative PhD grant. MF-M and JP were supported by the PID2019-110521GB-I00 and 522 TED2021-132627B-I00 grants funded by MCIN, AEI/10.13039/501100011033 and the 523 NextGeneration EU/PRTR. MF-M. was supported by a postdoctoral fellowship from ”la 524 Caixa” Foundation (ID 100010434), code: LCF/BQ/PI21/11830010. 525 Author Contributions: D.L., J.B. and T.K. designed the research; D.L. performed the 526 research and analysed the data; M.F.M, J.P., T.K., J.B., and S.L.L. provided scientific advice; 527 M.F.M. and J.B. advised on statistical aspects; D.L. wrote the paper with special contribution 528 from J.P., J.B., M.F.M. and T.K. 529 Competing Interest Statement: The authors declare no conflict of interest. 530 531 Figure 1. Range-edge hotspots. Hotspots were identified using Getis-Ord Gi* analysis, 532 which returns a Z-score for each hexagon in the world. Only hexagons with p<0.05 are 533 considered hotspots and subsequently shown here. The Z-score of each hexagon is 534 represented by the color gradient. Biomes (as defined by the World Wildlife Fund) are 535 marked by colors. ‘T&sT’ and ‘Temp’ stand for ‘Tropic and Subtropic’ and ‘Temperate’, 536 respectively. 537 538 Figure 2. Range-edge hotspots at intersections between biomes. (a) Modeled distribution 539 of the percentage of hotspots at the edges of biomes from a permutated (randomised) per- 540 biome distribution of hotspots. The first panel represented the median value of all the other 541 biomes in the figure over a standardised z-distribution of biomes. The arrow marks the 542 percentage of RE hotspots at biome intersections in the data set (one-sided p-values) (b) Heat 543 maps of the percentage of RE hotspots (relative to the total number of hotspots) at the 20 544 intersection between two biomes. A biome-pair intersection with a significant number of 545 hotspots (p-value < 0.08) is marked with an asterisk. Biome-pairs that have no intersections 546 are gray. Panel 1 is normalised over a global bootstrap and panel 2 over a per-biome 547 bootstrap. (c) Correlations (and regression lines) of the relationships between the number of 548 REs at the edges of biomes and the number of REs within the biomes. The shaded area 549 represents the 95% confidence interval around the regression line. P-values are calculated 550 using a two-sided Student’s T-test (degrees of freedom = 12). See Methods and Materials for 551 further information on the methodology for obtaining significance levels (for (b)) and hotspot 552 permutations. ‘T&sT’ and ‘Temp’ stand for ‘Tropic and Subtropic’ and ‘Temperate’, 553 respectively. 554 555 556 Figure 3. Climatic predictors of range-edge formation. Models of RE formation using the 557 absolute climate (mean) and the spatial heterogeneity (PV index) of the 19 WorldClim 558 variables and elevation. (a) Individual (binomial) mixed regression models between each of 559 the predictor variables and number of REs (1-19 are BioClim variables) and accounting for 560 continents and biomes as random effects. The estimated coefficients of the explanatory 561 variables (β) are represented by the color gradient. ‘*’ and ‘**’ represent p-values <0.05 and 562 <0.01, respectively (two-sided Student T-test). 563 BioClim1 – Annual Mean Temperature, BioClim2 - Mean Diurnal Range (Mean of monthly 564 (max temp - min temp), BioClim3 – Isothermality, BioClim4 - Temperature Seasonality, 565 BioClim5 - Max Temperature of Warmest Month, BioClim6 - Min Temperature of Coldest 566 Month, BioClim7 - Temperature Annual Range, BioClim8 - Mean Temperature of Wettest 567 Quarter, BioClim9 - Mean Temperature of Driest Quarter, BioClim10 - Mean Temperature of 568 Warmest Quarter, BioClim11 - Mean Temperature of Coldest Quarter, BioClim12 - Annual 569 Precipitation, BioClim13 – Precipitation of Wettest Month, BioClim14 – Precipitation of 570 Driest Month, BioClim15 – Precipitation Seasonality (PV), BioClim16 – Precipitation of 571 Wettest Quarter, BioClim17 – Precipitation of Driest Quarter, BioClim18 – Precipitation of 572 Driest Quarter, BioClim19 – Precipitation of Coldest Quarter. 573 (b) Violin plots depicting the results from (a) for four predictor variables. The distribution of 574 continental climates is shown in gray, contrasted with the climatic distribution specific to the 575 RE hotspots (scaled to the intensity of the hotspot, i.e. the Z- score). (c) Forest plot of a 576 model average from the highest predicting LMM with both BioClim and ENVRIEM 577 variables (identified with a model selection). Beta values (log-odds) are shown for each 21 578 predictor. Absolute and SH climate was obtained from all of the inland global hexagonal 579 units (n = 5851). Error bars represent 95% confidence interval around the average effect. 580 581 References 582 1. Soule, M. The epistasis cycle: A theory of marginal populations. Annu Rev Ecol Syst 4, 165– 583 187 (1973). 584 2. Brown, J. H. On the Relationship between Abundance and Distribution of Species. Am Nat 585 124, 255–279 (1984). 586 3. Gaston, K. J. The Structure and Dynamics of Geographic Ranges. Oxford University Press 587 (2003). 588 4. Sexton, J. P., McIntyre, P. J., Angert, A. L. & Rice, K. J. Evolution and ecology of species range 589 limits. Annu Rev Ecol Evol Syst 40, 415–436 (2009). 590 5. Gaston, K. J. Geographic range limits: Achieving synthesis. Proceedings of the Royal Society B: 591 Biological Sciences 276, 1395–1406 (2009). 592 6. Goldberg, E. E. & Lande, R. Notes and Comments Species’ Borders and Dispersal Barriers. Am. 593 Nat vol. 170 (2007). 594 7. Bachmann, J. C., Rensburg, A. J. van, Cortazar-Chinarro, M., Laurila, A. & Buskirk, J. van. Gene 595 flow limits adaptation along steep environmental gradients. American Naturalist 195, E67– 596 E86 (2020). 597 8. Hargreaves, A. L., Samis, K. E. & Eckert, C. G. Are species’ range limits simply niche limits writ 598 large? A review of transplant experiments beyond the range. American Naturalist 183, 157– 599 173 (2014). 600 9. Henry, R. C., Bartoń, K. A. & Travis, J. M. J. Mutation accumulation and the formation of range 601 limits. Biol Lett 11, 11DUUMY (2015). 602 10. Perrier, A., Sánchez-Castro, D. & Willi, Y. Environment dependence of the expression of 603 mutational load and species’ range limits. J Evol Biol 35, 731–741 (2022). 604 11. Bontrager, M. et al. Adaptation across geographic ranges is consistent with strong selection in 605 marginal climates and legacies of range expansion. Evolution (N Y) 1–18 (2021) 606 doi:10.1111/evo.14231. 607 12. Santini, L., Pironon, S., Maiorano, L. & Thuiller, W. Addressing common pitfalls does not 608 provide more support to geographical and ecological abundant-centre hypotheses. 609 Ecography 42, 696–705 (2019). 610 13. Oldfather, M. F., Kling, M. M., Sheth, S. N., Emery, N. C. & Ackerly, D. D. Range edges in 611 heterogeneous landscapes: Integrating geographic scale and climate complexity into range 612 dynamics. Glob Chang Biol 26, 1055–1067 (2020). 613 14. Janzen, D. H. Why Mountain Passes are Higher in the Tropics WHY MOUNTAIN PASSES ARE 614 HIGHER IN THE TROPICS*. Source: The American Naturalist vol. 101. 22 615 15. Maxwell, M. F., Leprieur, F., Quimbayo, J. P., Floeter, S. R. & Bender, M. G. Global patterns 616 and drivers of beta diversity facets of reef fish faunas. J Biogeogr 49, 954–967 (2022). 617 16. Roy, K., Hunt, G., Jablonski, D., Krug, A. Z. & Valentine, J. W. A macroevolutionary perspective 618 on species range limits. Proceedings of the Royal Society B: Biological Sciences 276, 1485– 619 1493 (2009). 620 17. Loiseau, N. et al. Global distribution and conservation status of ecologically rare mammal and 621 bird species. Nat Commun 11, (2020). 622 18. Kerkhoff, A. J., Moriarty, P. E. & Weiser, M. D. The latitudinal species richness gradient in New 623 World woody angiosperms is consistent with the tropical conservatism hypothesis. Proc Natl 624 Acad Sci U S A 111, 8125–8130 (2014). 625 19. Donoghue, M. J. & Edwards, E. J. Biome shifts and niche evolution in plants. Annu Rev Ecol 626 Evol Syst 45, 547–572 (2014). 627 20. Ringelberg, J. J., Zimmermann, N. E., Weeks, A., Lavin, M. & Hughes, C. E. Biomes as 628 evolutionary arenas: Convergence and conservatism in the trans-continental succulent 629 biome. Global Ecology and Biogeography 29, 1100–1113 (2020). 630 21. Smith, J. R. et al. A global test of ecoregions. Nat Ecol Evol 2, 1889–1896 (2018). 631 22. Zhu, K., Woodall, C. W. & Clark, J. S. Failure to migrate: Lack of tree range expansion in 632 response to climate change. Glob Chang Biol 18, 1042–1052 (2012). 633 23. Corlett, R. T. & Westcott, D. A. Will plant movements keep up with climate change? Trends in 634 Ecology and Evolution vol. 28 482–488 Preprint at https://doi.org/10.1016/j.tree.2013.04.003 635 (2013). 636 24. la Sorte, F. A., Butchart, S. H. M., Jetz, W. & Böhning-Gaese, K. Range-Wide Latitudinal and 637 Elevational Temperature Gradients for the World’s Terrestrial Birds: Implications under 638 Global Climate Change. PLoS One 9, e98361 (2014). 639 25. Paquette, A. & Messier, C. The effect of biodiversity on tree productivity: From temperate to 640 boreal forests. Global Ecology and Biogeography 20, 170–180 (2011). 641 26. Pichancourt, J. B., Firn, J., Chadès, I. & Martin, T. G. Growing biodiverse carbon-rich forests. 642 Glob Chang Biol 20, 382–393 (2014). 643 27. Pennington, R. T., Lavin, M. & Oliveira-Filho, A. Woody plant diversity, evolution, and ecology 644 in the tropics: Perspectives from seasonally dry tropical forests. Annu Rev Ecol Evol Syst 40, 645 437–457 (2009). 646 28. Title, P. O. & Bemmels, J. B. ENVIREM: an expanded set of bioclimatic and topographic 647 variables increases flexibility and improves performance of ecological niche modeling. 648 Ecography 41, 291–307 (2018). 649 29. Veresoglou, S. D. & Peñuelas, J. Variance in biomass-allocation fractions is explained by 650 distribution in European trees. New Phytologist 222, 1352–1363 (2019). 651 30. Grantham, H. S. et al. Anthropogenic modification of forests means only 40% of remaining 652 forests have high ecosystem integrity. Nat Commun 11, (2020). 23 653 31. Holdridge, L. R. Determination of world plant formations from simple climatic data. Science 654 (1979) 105, 367–368 (1947). 655 32. Whittaker, R. H. Classification of Natural Communities. Botanical Review 28, 1–239 (1962). 656 33. McDonald, R. et al. Species compositional similarity and ecoregions: Do ecoregion boundaries 657 represent zones of high species turnover? Biol Conserv 126, 24–40 (2005). 658 34. von Humboldt, A. & Bonpland, A. Essay on the Geography of Plants. (2013). 659 35. Cardillo, M. Latitude and rates of diversifcation in birds and butterfies. Proc. R. Soc. Lond. B 660 266, 1221–1225 (1999). 661 36. Hillebrand, H. On the Generality of the Latitudinal Diversity Gradient. Am. Nat vol. 163 (2004). 662 37. Mittelbach, G. G. et al. Evolution and the latitudinal diversity gradient: Speciation, extinction 663 and biogeography. Ecology Letters vol. 10 315–331 Preprint at 664 https://doi.org/10.1111/j.1461-0248.2007.01020.x (2007). 665 38. Hewitt, G. M. Genetic consequences of climatic oscillations in the Quaternary. in 666 Philosophical Transactions of the Royal Society B: Biological Sciences vol. 359 183–195 (2004). 667 39. Crane, P. & Scott, L. Angiosperm Diversification and Paleolatitudinal Gradients in Cretaceous 668 Floristic Diversity. Science (1979) 246, 675–678 (1989). 669 40. Jablonski, D. The tropics as a source of evolutionary novelty through geological time. P. J. IEEE 670 Trans. Geosci. Remote Sensing 361, 180–191 (1993). 671 41. Jablonski, D. et al. Out of the tropics, but how? Fossils, bridge species, and thermal ranges in 672 the dynamics of the marine latitudinal diversity gradient. Proc Natl Acad Sci U S A 110, 673 10487–10494 (2013). 674 42. Antonelli, A. et al. An engine for global plant diversity: Highest evolutionary turnover and 675 emigration in the American tropics. Front Genet 6, (2015). 676 43. Jump, A. S. & Peñuelas, J. Running to stand still: Adaptation and the response of plants to 677 rapid climate change. Ecol Lett 8, 1010–1020 (2005). 678 44. Morreale, L. L., Thompson, J. R., Tang, X., Reinmann, A. B. & Hutyra, L. R. Elevated growth and 679 biomass along temperate forest edges. Nat Commun 12, 7181 (2021). 680 45. Wilkinson, S., Clephan, A. L. & Davies, W. J. Rapid Low Temperature-Induced Stomatal Closure 681 Occurs in Cold-Tolerant Commelina communis Leaves But Not in Cold-Sensitive Tobacco 682 Leaves, via a Mechanism That Involves Apoplastic Calcium But Not Abscisic Acid. Plant Physiol 683 126, 1566–1578 (2001). 684 46. Brodribb, T. J. & Holbrook, N. M. Stomatal protection against hydraulic failure: a comparison 685 of coexisting ferns and angiosperms. New Phytologist 162, 663–670 (2004). 686 47. Davis, B. A. S. & Brewer, S. Orbital forcing and role of the latitudinal insolation/temperature 687 gradient. Clim Dyn 32, 143–165 (2009). 688 48. Seager, R. et al. Strengthening tropical Pacific zonal sea surface temperature gradient 689 consistent with rising greenhouse gases. Nat Clim Chang 9, 517–522 (2019). 24 690 49. Xu, Y. & Ramanathan, V. Latitudinally asymmetric response of global surface temperature: 691 Implications for regional climate change. Geophys Res Lett 39, n/a-n/a (2012). 692 50. Colwell, R. K., Brehm, G., Cardelús, C. L., Gilman, A. C. & Longino, J. T. Global Warming, 693 Elevational Range Shifts, and Lowland Biotic Attrition in the Wet Tropics. vol. 322 694 www.sciencemag.org (2008). 695 51. Basso, B., Martinez-Feria, R. A., Rill, L. & Ritchie, J. T. Contrasting long-term temperature 696 trends reveal minor changes in projected potential evapotranspiration in the US Midwest. 697 Nat Commun 12, 1476 (2021). 698 52. Zizka, A. et al. CoordinateCleaner: Standardized cleaning of occurrence records from 699 biological collection databases. Methods Ecol Evol 10, 744–751 (2019). 700 53. Serra-diaz, J. M., Enquist, B. J., Maitner, B., Merow, C. & Svenning, J. Big data of tree species 701 distributions : how big and how good ? (2018) doi:10.1186/s40663-017-0120-0. 702 54. Zizka, A. et al. No one-size-fits-all solution to clean GBIF. PeerJ 8, (2020). 703 55. Getis, A. & Ord, J. K. The Analysis of Spatial Association by Use of Distance Statistics. Geogr 704 Anal (1992). 705 56. Mendez, C. Spatial autocorrelation analysis in R. R Studio/RPubs. https://rpubs.com/quarcs- 706 lab/spatial-autocorrelation (2020). 707 57. Bivand, R. S., Pebesma, E. & Gómez-Rubio, V. Applied Spatial Data Analysis with R. (2013). 708 58. Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution 709 interpolated climate surfaces for global land areas. International Journal of Climatology 25, 710 1965–1978 (2005). 711 59. Heath, J. P. Quantifying temporal variability in population abundances. Oikos vol. 115 573– 712 581 Preprint at https://doi.org/10.1111/j.2006.0030-1299.15067.x (2006). 713 60. Fernández-Martínez, M. et al. The consecutive disparity index, D: a measure of temporal 714 variability in ecological studies. Ecosphere 9, (2018). 715 61. Bartoń, K. MuMIn: Multi-model inference. R package version 1.10.0.. 1. . (2013). 716 62. F. Dormann, C. et al. Methods to account for spatial autocorrelation in the analysis of species 717 distributional data: A review. Ecography vol. 30 609–628 Preprint at 718 https://doi.org/10.1111/j.2007.0906-7590.05171.x (2007). 719 720 25 North Hotspots South Hotspots East Hotspots West Hotspots Hotspots Z-score 30 20 T&sT Moist broadleaf Mediterranean T&sT Dry broadleaf Temp conif Temp grasslands Montane grasslands Desert T&sT coniferous Boreal Flooded grasslands 10 Tundra Mangrove Temp broadleaf T&sT grasslands (a) 0.4 T&sT Moist broadleaf T&sT grasslands Temp grasslands % RE at 0.4 biome p= 0 p = 0.067 0.4 p = 0.902 0.3 0.3 Density intersect 0.3 0.3 0.2 p = 0.06 0.2 0.2 0.2 0.1 0.1 0.1 0.1 0 0 0 0 -4 -2 0 2 -2 0 2 4 -2 0 2 -2 0 2 Montane grasslands Tundra Mediterranean Desert 0.4 0.4 p = 0.023 p = 0.927 0.4 p = 0.946 p= 0 Density 0.3 0.3 0.3 0.3 0.2 0.2 0.2 0.2 0.1 0.1 0.1 0.1 0 0 0 0 -2 0 2 -3 -2 -1 0 1 2 -4 -2 0 2 -3 0 3 6 T&sT Dry broadleaf Temp broadleaf Temp conif Boreal 0.4 0.4 0.4 p = 0.113 p = 0.006 0.4 p= 1 p= 0 0.3 0.3 Density 0.3 0.3 0.2 0.2 0.2 0.2 0.1 0.1 0.1 0.1 0 0 0 0 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 3 -5 -2.5 0 -2 0 2 4 (b) Montane grasslands * Temp grasslands * * * Temp conif * * Temp broadleaf * * Desert * Mediterranean * * Boreal * Tundra * * Mangrove * Flooded grasslands T&sT grasslands * 1 T&sT coniferous * 0 0 T&sT Dry broadleaf * -1 -1 T&sT Moist broadleaf -2 ds ds n sl s sl s s de ras ve rr ra f ds Tu s n de Ma nds sl s sl s l f as ve p rt ne ass if ra af ni af rr al M mp mp af ea br rea oa rt ne ass if d d T& M rou ea ea sT nif af d d u ea n Te Te ese n ite nd an an gr an gr lan br se g gro nd le ite re co dle an an gr lan gr ero Te Te dle gr gro gr co T& co dle gr co dl dl an a an o fe ad ed Bo ed Tu e sl oa sl oa B D oa an D sT oa n as as as ro as m br br r sT t b T& y b st p ry sT p T& ois d Fl sT M mp m d oi m r M D M D T& Te M ta Te M ta sT oo oo on on sT sT T& Fl T& T& (c) 6.0 R = 0.89 p = 1.6e-05 Tundra , Boreal log(mean hotspots at biome edge) 5.0 Temp grasslands Temp broadleaf 4.0 Temp conif Desert Mediterranean T&sT Dry broadleaf Montane grasslands 3.0 T&sT grasslands T&sT Moist broadleaf T&sT coniferous Mangrove Flooded grasslands 0.5 1.0 1.5 2.0 log(hotspot within biome) (a) ** ** ** ** 0.3 ** ** * ** * ** * 0.2 * β R2 * ** ** 0.4 ** ** ** 0.2 ** * ** 0.1 * 0.0 -0.2 -0.4 0.0 -0.6 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 El El La . ev ev (b) t. Mean Climate Spatial Hererogeneity Climate (SH climate) . SH BioClim16: Precipitation of Wettest Q (mm) mean BioClim10: Temperature Warmest Q (°C) mean BioClim8: Temperature Wettest Q (°C) mean BioClim1: Mean Annual Temperature (°C) 40 30 30 20 20 10 0 10 -10 0 -20 -10 0.6 20 0.4 0 0.2 -20 0.0 a a lia a a e lia ric ric a ia op e ric ric rica ia op tra me me Af ric As r tra me me Af As Eu r s A A Eu Au s A A Au rth uth rth uth No So No So (c) BioClim Variables SH climate T. annual range Min T. coldest M Isothermality Elevation PPT. wettest M PPT. wettest Q T. warmest Q T. warmest Q -0.002 -0.001 0 0.001 Mean 0.00 0.25 0.50 0.75 1.00 ENVIREM variables SH climate embergQ Continentally Moisture Index PET wettest Q PET warmest Q PET driest Q PET coldest Q Monthly county by T. PET coldest Q -0.004 -0.003 -0.002 -0.001 0 0.001 Mean 0.00 0.25 0.50 0.75 1.00 1 A biome-dependent distribution gradient of tree species range edges is 2 strongly dictated by climate spatial heterogeneity 3 David Lerner1, Marcos Fernández Martínez2, Stav Livne-Luzon1, Jonathan Belmaker3,4, 4 Josep Peñuelas2,5, Tamir Klein1 5 1 Department of Plant & Environmental Sciences, Weizmann Institute of Science. Rehovot 6 76100, Israel. 7 2 CREAF, Cerdanyola de Vallès, Catalonia E08193, Spain. 8 3 School of Zoology, Faculty of Life Sciences, Tel Aviv University, Tel Aviv, Israel. 9 4 Steinhardt Museum of Natural History, Tel Aviv University, Tel Aviv, Israel. 10 5 CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Catalonia E08193, Spain 11 12 Abstract 13 Understanding the causes of the arrest of species distributions has been a fundamental 14 question in ecology and evolution. These questions are of particular interest for trees due to 15 their long lifespan and sessile nature. A surge in data-availability evokes for a macro- 16 ecological analysis to determine the underlying forces limiting distributions. Here we analyse 17 the spatial distribution of >3600 major tree species to determine geographical areas of range- 18 edge hotspots and find drivers for their arrest. We confirmed biome edges to be strong 19 delineators of distributions. Importantly, we identified a stronger contribution of temperate 20 than tropical biomes to range edges, adding strength to the notion that tropical areas are 21 centers of radiation. We subsequently identified a strong association of range-edge hotspots 22 with steep spatial climatic gradients. We linked spatial and temporal homogeneity and high 23 potential evapotranspiration in the tropics as the strongest predictors for this phenomenon. 24 We propose that the poleward migration of species in light of climate change might be 25 hindered due to steep climatic gradients. 26 27 Main Text 28 Introduction 29 The geographical distributions of species are marked by their range limits. Understanding the 30 causes of distribution arrest has been a fundamental question in ecology and evolution1–4. 31 Given the strong interplay between biotic, abiotic, demographic, physical and historical 32 forces in predicting range-edges, it has been challenging to find underpinnings for their 33 formation. Two main environmental forces seem to play a major role in the formation of 1 34 range edges; spatial environmental heterogeneity and habitat quality4,54,5. Most models on the 35 formation of species-range edges rely on the interplay between either one of these two forces 36 with non-climatic pressures to explain their formation. For example, steep climate gradients 37 combined with high dispersal and gene flow reduces species’ fitness and genetically 38 constrains their evolution into novel environments6,7. Likewise, low habitat quality reduces 39 population size2,8, increasing drift and migration load9–11. Nevertheless, the significance of 40 climate in the interplay between these two environmental components in defining range edges 41 for an array of species or on a wide biogeographical scale remains elusive12,13. 42 43 The field of biogeography has long sought understanding of species ranges despite having 44 limited tools (e.g. 3,14). Given the surge of large-scale datasets, it is now possible to better 45 identify the underpinnings of species distributions by studying the macro eco-evolutionary 46 processes involved in their formation11,12,15–17. Although methods and results are disparate 47 between studies, there is almost a consensus that the presence of large-scale biogeographical 48 units confine species with climate as their primary predictor. For example, Bontager et al.11 49 suggested distinct characteristics for range edge populations dependent on their latitude. 50 Likewise, niche conservatism and strong beta-diversity patterns seem to withhold at large 51 macro-ecological scales16–18. This is generally true for plant species, with biomes being the 52 most consistent classifier based on structural and functional similarity19,20. Recent efforts 53 have been made to understand how accurate and substantial biome entities are at defining 54 species distributions 20,21. Nevertheless, it remains an open question if, and to what extent, the 55 intersection between biomes is a source of species range edge hotspots. Deciphering such 56 patterns will enable the proper understanding of how communities of species redistribute and 57 are structured geographically, and if similar biomes in distinct geographic areas have similar 58 effects on the distribution of species. 59 60 However, although there has been an increased interest in defining the biogeographical 61 underpinnings of species distributions, most techniques have used species relatedness and 62 diversity metrics to test for the existence of shared niche space between species and 63 communities. The direct analysis of the ecological and climatic limitations to geographic 64 space, although trivial, remain elusive. Here, we look at the universal set of climatic factors 65 and geographical patterns of species distributions by focusing directly on species’ range edge 66 distributions. 67 2 68 We present the first-to-date global study of tree species range edges, applying a novel, 69 simple, yet effective method of delineating range edges (REs) in order to: (1) identify 70 deterministic patterns of REs, as seen by RE-dense areas (RE hotspots: REH); (2) determine 71 whether the classification of biomes as distinct community-level patterns of biodiversity 72 properly delineate the niche of species; (3) identify global-scale REH patterns, and; (4) 73 discern the underlying niche factors responsible for RE formation. In particular, we 74 investigate whether spatial heterogeneity of abiotic factors or a universal predictor for habitat 75 quality are determinants of RE formation. We focus specifically on tree species, as they are 76 an exemplary group of species in the study of the ecological changes predicted to occur at the 77 peripheries of distributions, given their long-lived characteristics and fundamental role in 78 many ecosystems25,26, specifically biomes 18,27. 79 80 Although we do not study the interplay between climate and other ecological and 81 evolutionary limiting factors to species distributions (e.g., seed dispersal, plasticity, and 82 adaptation), discerning these patterns and the climatic components leading to such 83 distributions will enable a better understanding of interplays between biotic and abiotic 84 factors in future studies. A better understanding of the climatic factors affecting dispersal 85 enable for better predictions of the success of species to track changing climates, and in turn, 86 if they will be subject of migration lags22–24. 87 88 Results 89 Global data set and REH distribution 90 We present the first report of the global distribution of tree range edge hotspots (REH) (Fig. 91 1), marked by the hexagons with significant clustering of REs (Supplementary Fig. 1). We 92 did not identify any significant coldspots, given the baseline presence of REs around the 93 globe. The visual patterns emerging from these distributions indicated that distributions 94 stopped disproportionally more at the edges of biomes than within them. For example, 95 northern REH occurred mostly at the intersection between a montane (Himalayas) and a 96 desert biome (Gobi Desert) or at the edge of the tundra in North America. Southern REH in 97 Africa and southern Eurasia tended to stop at the edges of desert biomes (e.g. Sahara Desert) 98 but stopped mainly at the edges of temperate and montane grasslands in the Southern 99 Hemisphere (pampas and Andes, respectively). Eastern and western REH were notably 100 mostly in similar geographical locations, e.g. at the intersections between the Himalayas and 101 central Asian deserts or at the edge of the Atacama Desert in South America. 3 102 103 In order to identify the underlying niche factors that define REs, we focused only on inland 104 REs, since REs at the edge of a water source is probably due to an obvious geographical 105 barrier rather than an ecological effect. The fraction of inland REs was globally similar 106 between the continents (between 65 and 75%), except for Australia (~45%, Supplementary 107 Fig. 2a). The large fraction of inland REs (Supplementary Fig. 2b) was mainly associated 108 with the edges of biomes (e.g. northern REs in Africa, Asia, and North America were 109 associated with the Sahara Desert, montane grasslands of the Himalayas, and the tundra 110 biome, respectively; Fig. 1). A significant fraction of eastern inland REs were in Europe, but 111 almost no significant REH were identified, implying the sparse distribution of REs 112 throughout the continent, or as an effect from a smaller area for distribution compared to 113 other continents. 114 115 Contribution of biome-biome intersections to REHs 116 We quantitatively identified the global patterns of arrest by analysing the fraction of REH 117 that stopped at biome-biome intersections (14 central global ecological regions best 118 distinguished by their climate, fauna, and flora obtained from World Wild Life (WWF, 119 http://www.worldwildlife.org/; Fig. 1) and by delineating a buffer zone at points of 120 intersection between two or more biomes. Although not significant, the number of REH was 121 strongly associated with the intersection between biomes (Fig. 2a, first panel), indicating that 122 climatic conditions were a probable cause for the REs of tree species at biome edges. Our 123 results, however, identified an unequal contribution of the different biome edges to the 124 fraction of REH (Fig 2a). We then analysed the individual biome-biome intersections 125 normalised by a global permutation (i.e., from permutation of the global distribution of REH, 126 see Materials and Methods). Here we identified (1) a strong contribution of REH at the 127 intersections between temperate and desert biomes (Fig. 2b, left panel) in comparison to (2) a 128 weaker contribution at the intersections of tropical and subtropical biomes (between 129 themselves and with temperate biomes). The significance of the contribution of REH 130 (p<0.05) (see Materials and Methods) was attributed almost exclusively to the intersections 131 within temperate biomes (asterisks in Fig. 2b). This unequal contribution between the 132 temperate vs tropical and subtropical biomes was globally consistent. Nevertheless, under a 133 per-biome permutation (i.e., normalising each intersect by a selective permutation from the 134 respective biome combinations separately; Fig 2b right panel), we observed a much weaker 4 135 contribution of range edges to the formation of biomes. Only a selected number of biome- 136 biome intersections had a significant contribution to the formation of REH, yet no specific 137 pattern between temperate and tropical biomes was observed. Indeed, the strong positive 138 correlation between the number of REs between biomes with the number of REs within 139 biomes (Fig. 2c; R = 0.9, linear regression, p < 0.001) reflects such discrepancy between 140 panels in Fig 2b, as the number of REH at the edge of biomes had a strong linear association 141 with the number of REH within that biome. We identified five biomes, however, outside the 142 95% confidence interval (CI) of the regression (Fig. 2c). The two biomes above the 143 regression CI (desert; and boreal) were biomes where the number of REH at the edges was 144 larger than predicted (Fig 2c). In contrast, temperate coniferous forest, mangrove and flooded 145 grasslands fell below the CI of the regression, indicating a much larger number of REH 146 within the biome compared to the edge. 147 148 A parallel analysis using the distribution of REs (rather than REH) was conducted to identify 149 similarities and differences between the distributions. The results indicated a similar pattern 150 of distribution (Supplementary Fig. 4a-b compared to Fig. 2b-c), suggesting global forces 151 associated with REH. 152 153 Climatic predictors of RE formation 154 Intersections between biomes were a significant cause of RE formation, so we also 155 investigated the dependence of RE formation on climate. We tested both the ‘absolute 156 climate’ (i.e. annual and seasonal average temperature and precipitation) and the ‘spatial 157 heterogeneity of climate’ (SH, i.e. spatial variability of absolute climate) at each of the global 158 hexagons using 19 bioclimatic variables obtained from WorldClim (see Materials and 159 Methods). Elevation (absolute and SH) and latitude were also accounted for at each hexagon 160 resolution. Generalised linear models of regression between each climatic variable as an 161 independent predictor indicated that all 40 climatic variables were significant predictors of 162 RE formation (p<0.05) (Supplementary Fig5a). Interestingly, most absolute climatic 163 variables were negatively associated with the global RE distribution (Fig. 3a and b), implying 164 a general prediction of REs occurring in climates with low temperatures and low 165 precipitation. A positive association was attributed to the four climatic variables that defined 166 the temporal heterogeneity of temperature and precipitation (mean diurnal range, temperature 5 167 seasonality temperature annual range and precipitation seasonality; BioClim variables 2, 4, 7, 168 and 15, respectively; Fig. 3a, Supplementary Fig. 5) and all of the SH climatic variables. A 169 mixed model with biomes and continents as random factors gave a reduced number of 170 variables that significantly predicted the formation of RE (Fig. 3a). A crossed model with 171 both continents and biomes consistently gave a stronger fit (AIC values) than only 172 considering either of these categorical random factors independently (Supplementary Fig 6b). 173 In this more stringent global analysis, latitude, absolute temperature and SH temperature, 174 precipitation and elevation were significant predictors of REH. The goodness of fit (R2) of 175 each model, as a measure of predictive strength, indicated that spatial heterogeneity 176 accounted for RE formation better than did their absolute equivalents. A model selection was 177 carried out to identify the most important factors associated with RE formation, followed by a 178 model averaging of the models with a ∆AIC < 2, (Figs. 3c). The SH climatic variables again 179 defined RE better than the absolute climatic variables. Although absolute climatic variables 180 such as isothermality (BioClim3), temperature of the wettest quarter (BioClim8) (both 181 characteristic of tropical and subtropical climates), and annual precipitation (BioClim12) 182 were strongly associated with REs in a generalised linear model, temperature of the warmest 183 quarter (BioClim10) was the only predictor strongly associated with RE under the mixed 184 model (Figs. 3a and S6). This difference in results can be visually seen when comparing 185 between continents (Fig. 3b) e.g. in panel 3; with temperature of wettest quarter (BioClim8) 186 having partial dissociations in Africa and South America, i.e. the two continents with the 187 most tropical biomes. Furthermore, SH isothermality (BioClim3), temperature of coldest 188 month (SH BioClim 6), precipitation (SH BioClim13 and 16) and elevation change (SH 189 elevation) were all predictors of REs. Interestingly, the absolute climatic predictors 190 (BioClim8 and BioClim10) were strongly negatively correlated with spatial heterogeneity at 191 the minimum temperature of the coldest month (SH BioClim6) (Supplementary Fig. 7b), 192 indicating its representation of a tropical biome climate. The weaker relative importance of 193 SH precipitation (Fig. 3c, left panel) is due to its dissociation with the continents that most 194 strongly represent temperate regions; Europe and North America (Fig 3b, panel 4). 195 196 In parallel, we ran models with the ENVIREM dataset, a dataset of environmental variables 197 complementary to WorldClim that are more ecophysiologically meaningful for plant 198 species28. Most ENVIREM variables associated with REH are spatially heterogeneous 199 variables (Fig. 3c, Supplementary Fig. 7). Absolute potential evapotranspiration (PET) of the 200 coldest quarter was strongly correlated with both absolute temperature of the warmest quarter 6 201 (BioClim10) and SH temperature of the coldest month (BioClim6); the two predictors 202 indicative of a transition between tropical and temperate biomes (Supplementary Fig. 7b). SH 203 EmbergerQ and SH PET of warmest quarter were strongly associated with SH precipitation 204 (BioClim13 and 16) and SH temperature of the warmest quarter (BioClim10), respectively. 205 Albeit the strong correlation of all absolute ENVIREM variables with absolute BioClim 206 variables28, we identified several SH ENVIREM variables to be weakly correlated to BioClim 207 (e.g., SH moisture index and PET of driest month). All results in Fig. 3 were robust against 208 spatial autocorrelation. (Supplementary Fig. 5b, see Materials and Methods for further 209 information). 210 211 The REH distribution from our generated polygons was compared to a REH distribution from 212 expert polygons (see Materials and Methods), in order to test for the accuracy of our 213 generated dataset. Indeed, results of GLMs converge, indicating the robustness of our 214 generated dataset to describe REs (Supplementary Fig. S8). 215 216 Discussion 217 The results support our hypothesis of a nonrandom distribution of REs. We were able to 218 confirm that many tree species range edges were clustered rather than sparsely (stochastic) 219 distributed, by obtaining a larger number of significant REH. This finding suggests the 220 underlying presence of ecological and evolutionary forces governing REH formation. 221 Similarly, the matching results obtained from the biome analyses when accounting for REs 222 (Supplementary Fig. 4) or REH (Fig. 2) was also indicative of the deterministic clustering of 223 species REs at these specific ecological barriers (biome edges). Nevertheless, we identified 224 case-specific exceptions such as the scattered distribution of eastern REs throughout Europe, 225 indicated by the strong identification of internal REs (Supplementary Fig. 2b) but not REH 226 (Fig. 1c). This could potentially have occurred as a consequence of either the intense 227 anthropogenic activity throughout29,30 that prevented the distribution of tree species to reach 228 their natural REs, or the effect of smaller land area compared to other continents, altering the 229 effect of how species interact with abiotic factors, thus affecting their distribution and 230 adaptation. 231 232 Our results also indicated a strong dependence of REH on the edges of biomes, strongly 233 supporting the many efforts to determine whether the division of the planet into discrete 234 geographical units has been appropriately delineated21,31–34. Our findings, however, were 7 235 unexpectedly biome-specific, identifying the biome borders that most defined tree 236 distributions (Fig. 2). There is no seemingly obvious pattern of differential contribution to 237 REH when analysing each biome independently (per-biome bootstrap; Fig 2a). Yet, a clear 238 distinction between a tropical/temperate distribution of REH is strongly observed under a 239 biome-biome pairs analysis (global bootstrap; Fig 2b). Specifically, REs were strongly 240 dependent on desert, temperate, and montane edges compared to the weak dependence on 241 tropical and subtropical biome edges. Despite the marginally significant association of range 242 edges to biomes globally (Fig 2a), a tropical-temperate REH distribution is not specifically 243 dependent on biome edges, but rather, it represents a global biome trend (Fig 2b, left panel; 244 Fig. 2c). In addition, the importance of biomes in describing climate’s association to REs, as 245 seen from a mixed effects model (Fig. 3a, S6c), as well as the absence of latitude in 246 predicting REs in a global model (Fig 3c) is further indicative of the importance of biomes 247 (rather than a latitudinal effect) in defining REs. 248 249 The large differences in the relative number of REs between the tropical and temperate 250 biomes can indicate adaptive mechanisms between species residing in either of these two 251 types of biomes. These nontrivial results may have been related to the Latitude Diversity 252 Gradient35–37, a well-established pattern in which biodiversity is higher in the tropics than in 253 temperate regions. Firstly, the notion that these differences surged through the effective 254 evolutionary time hypothesis36,38 suggests that genetic diversity is higher in the tropics due to 255 the possible longer times needed for adaptation and expansion. Similarly, long-term climatic 256 oscillations have been suggested to reduce cladogenesis at higher latitudes, consistent with 257 the observation that tropical areas are centers of evolutionary novelty36,39,40 and evolve faster 258 than temperate regions35, leading to the notion of tropics as centers for radial expansion of 259 clades and species. Recent studies have identified mechanisms for this ‘out of the tropics 260 expansion’ model41,42 and have reported a higher fraction of bridge species (species violating 261 niche conservatism) from the tropics compared to temperate regions. We conclude that the 262 lower contribution to REs in tropical vs temperate regions is consistent with stronger radial 263 expansions from the tropics than temperate regions. 264 265 However, other mechanisms could explain the distinctive REH patterns between the tropic 266 and temperate regions. For example, based on several studies demonstrating the large 267 physiological and evolutionary effects of forest fragmentation29,43,44, the increased long-term 268 anthropogenic activity, and consequently excessive fragmentation, in forests in temperate 8 269 regions compared to tropical and subtropical regions may have hindered the adaptation of 270 species to novel climates. The lack of association of REH specifically to temperate 271 grasslands, temperate conifers and Mediterranean biomes specifically (Fig 2a) could be 272 indicative of such an effect, as these biomes have been historically subject to strong 273 anthropogenic activity. 274 275 Our analyses confirmed the strong dependence of RE on climate. Lower temperature, higher 276 climatic heterogeneity (both temporal and spatial), and elevation changes were the strongest 277 climatic predictors of REs. The analysis of the best predictive variables indicated a noticeable 278 weak correlation between most of the leading factors determining REs (Supplementary Fig. 279 7b), suggesting that these factors were site-specific predictors of REs (e.g. desert biome with 280 low annual precipitation, montane grassland with high elevation, and tropical and subtropical 281 biomes with temperature homogeneity (isothermality, temperature of the wettest quarter and 282 spatial heterogeneity of low temperature – SH6)). We suggest that the strong negative 283 correlation between these absolute climatic variables and the spatial heterogeneity of low 284 temperature is an indication of different temperature patterns between tropical and temperate 285 regions, i.e. a buffered heterogeneous temperature in the tropics in contrast to the latitudinal 286 effect of decreasing temperatures in the temperate biomes. There was an overall positive 287 trend for the effect of precipitation (BC12-BC19) when controlling for biome since the signal 288 appeared only when biome is included as a random factor (Supplementary Fig. 6). This 289 means there must be other non-accounted confounding factors altering the relationship when 290 analysing all biomes together, suggesting that the effect is not biome dependent. Although 291 this study does not make a note of the complex interplay between biotic and abiotic forces or 292 the ecological traits of tree species (e.g. seed dispersal, phenology) in the formation of RE, 293 we find a consistent global effect of temperature and spatial heterogeneity of temperature and 294 precipitation to predict their formation throughout the different models. This is indicative of 295 their principal role as universal predictors of species distributions. 296 297 The strong prediction of REH formation from PET of the coldest quarter and its covariance to 298 the tropical-temperate transition variables (SH BioClim6 and absolute BioClim10), could 299 reflect a possible mechanistic evolutionary constrain for the distribution of woody species’ at 300 temperate biomes. At higher latitudinal temperate biomes, where evapotranspiration is 301 strongly reduced, especially during the coldest quarter, there is a strong limitation to 302 photosynthesis and growth by the significant time reduction of stomatal conductance45,46. 9 303 Similarly, given the strength of its prediction of REH, SH embergerQ (pluviothermic 304 quotient) could indicate a more refined mechanism for RE formation than its covariates of 305 precipitation (SH BioClim 13 and 16). This index describes mean annual precipitation in 306 relation to annual changes in temperature. EmbergerQ thus increases the predictability of 307 how precipitation also dictates the formation of REH in the more temperate biomes (Fig 3b, 308 panel 4). The consistency of our results indicates spatial and temporal heterogeneity of 309 climate and topography as overwhelmingly stronger predictors of RE formation than their 310 absolute climatic counterparts (Fig. 3) and, in turn, frail evidence for a universal poor habitat 311 quality. Even in the cases where mean climatic variables strongly predict REH, these were 312 strongly associated with this transition between the spatial and temporal climatically 313 homogenous tropical and subtropical biomes to the more heterogeneous temperate biomes. 314 Nevertheless, we note the importance of a lack of evapotranspiration, particularly in cold 315 climates, as a main predictor of RE formation. 316 317 These observations have substantial implications for the effects of climate change on tree 318 distributions and its effects on tree migration. Although predictions for future steeper 319 temperature gradient as a result of greenhouse gas emission and climate change has not been 320 trivial47–49, such an increase in temperature gradients could vastly affect the distribution of 321 species. In particular, our results strengthen the growing understanding that the predicted 322 poleward migration of tree species might not be as successful as previously predicted23,24. The 323 increase of stronger spatial gradients (especially in the lower latitudes)24,48,50 or extreme and 324 spontaneous events might all be causes of migration lags, despite the suitable temperatures at 325 higher latitudes and altitudes. Likewise, the importance of PET from temperate biomes on the 326 formation of REH presented here could also suggest a possible migration lag or loss of 327 adaptation due to the predicted reduction in PET at higher latitudes51. Our results thus 328 highlight the importance of accounting for more precise spatial heterogeneity of climate as a 329 critical feature in future models of species distribution and the development of more precise 330 conservation efforts such as assisted migration. 331 332 Materials and Methods 333 Data acquisition and polygon formation 334 Supplementary Fig. 9 visually summarises the methodologies used to obtain global range 335 edge hotspot (REH) distributions. We downloaded a data set of tree species from the open- 336 source data set using R packages rgbif and taxize. Global Biodiversity Information Facility 10 337 (GBIF; 05 July 2021, https://doi.org/10.15468/dl.ajen6k) using the Botanic Gardens 338 Conservation list of 60000 tree species. We downloaded occurrences with entries from 1980 339 onwards, removing any occurrence reported with a geospatial issue, species not belonging to 340 the kingdom Plantae (in case of mismatched species names), and any occurrence marked as 341 unlikely, mismatched, or invalid. We removed occurrences that had reported uncertainties of 342 >100 km and records based on fossils and unknown sources. We then used the 343 CoordinateCleaner R package52 to remove any occurrences with zero coordinates, equal x 344 and y coordinates, duplicates, occurrences at sea, coordinates at capitals, and centroids. To 345 finalise, we again removed species with <300 occurrences. The data we used undoubtedly 346 contained sampling bias53, probably overrepresenting the number of REs in some regions 347 with a reduced or negligible sampling effort. We tried to overcome this issue by basing our 348 filtering steps on several previous studies52–54. The strength of the critical filtering steps 349 applied in our analysis resembled those previously presented53. 350 We converted the georeferenced species occurrences (x and y coordinates) into distributional 351 polygons in parallel using two independent techniques; through concave-hull 352 (Supplementary Methods) and multivariate kernel-density estimation (described below). 353 Given the strong similarity between the two methods (Supplementary Fig. 8), we discuss the 354 methods and results in detail only for the kernel-density estimated polygons. We created 355 polygons using two-dimensional kernel density estimations. We first divided the extent of all 356 the coordinates into 800 grid points in each dimension (longitude and latitude) in order to 357 produce a matrix of 640,000 grid cells for each of the species. Subsequently, we selected for 358 the grid cells with the highest 99% estimation of the species’ occurrence and subsequently 359 rasterised these. Polygons were then delineated around the contour of the rasters. 360 361 Polygon groupings 362 All polygons belonging to the same species were grouped based on absolute distance from 363 one another. Polygons separated by ≤500 km were grouped together, with the assumption that 364 fragmentation, gene flow, and unreported data could all warrant two nearby populations to be 365 considered as one. We used the hclust function (package stats, agglomeration method: 366 complete) to hierarchically cluster populations from a sequence of three or more populations 367 by their distances, also using a cutoff of 500 km (cutree function, package stats) for 368 determining the clusters. The final data set comprised >3600 tree species, ranging from one to 11 369 nine populations (polygons) per species, for a total of 8500 populations. All spatial data was 370 analysed using R packages sf and raster. 371 372 RE determination 373 RE dense areas were determined by (1) defining distinct global units, (2) identifying the RE 374 of each species and (3) map species’ RE to the global units to calculate the density of 375 REs/unit. In detail: (1) We rasterised the world map to spatially bin the density of REs. To 376 overcome the problem of spatial distortion, we used hexagonal bins with the dggridR R 377 package, developed using the ISEA Discrete Global Grids system, a repetition of polygons on 378 the surface of an icosahedron, allowing for the projection of equal sized bins onto a 2D plane. 379 We defined the size of each hexagon as ~23 000 km2 (with an average spacing between center 380 nodes of 165 km). (2) We used coordinates of the cardinal directions (north, south, east, and 381 west) to represent species REs by subdividing each polygon cluster into four quartiles in the 382 four cardinal directions (NE, NW, SE, and SW). The REs for each quartile were determined 383 as the two most-outward coordinates of the corresponding cardinal directions (e.g. north and 384 east cardinal coordinates for the NE quartile). Eight REs were thus determined for each 385 population (two for each cardinal direction). As a filtering step, we accounted for both REs 386 from the same cardinal direction if they were >20 arc-degrees apart, otherwise we only 387 accounted for the farthest point from the centroid. (3) The total number of REs obtained using 388 this method was normalised by the total number of species intersecting its respective 389 hexagon. In parallel, we also defined REs by accounting for the perimeter of the polygon for 390 each species (Supplementary Fig. 10). The ‘perimeter’ system may be a more realistic and 391 complete system for identifying REs, but the ‘cardinal coordinate’ system, although more 392 simplistic in nature, (1) provides a clearer visual representation of the distribution of REs, and 393 (2) allows for the directionality of REs to be compared, essential farther along the pipeline by 394 distinguishing between coastline and inland REs and identifying hotspots in the four cardinal 395 directions. 396 397 We then classified coastline and inland REs, assuming that the arrest of species distribution at 398 the edge of a water source was probably due to an obvious geographical barrier rather than an 399 ecological effect. Coastline REs were determined by creating a semicircle (buffer of 3 arc- 400 degrees) around each RE in the direction of its cardinal coordinate and measured the 401 percentage overlap with water. Cardinal coordinates with >50% overlap were considered a 402 “coastal RE”. 12 403 404 In order to find the probability distribution of REs at the edge of water source 405 (Supplementary Fig. 2), we permutated the global population of REs (except for the 406 Australian population) and used the mean of this permutation to compare to the number of 407 REs in each of the continents. 408 409 Hotspot analysis 410 Hotspots were identified using the Getis-Ord Gi* hotspot analysis55,56 to find spatial 411 correlations between hexagon (inland and normalised) RE densities. We initially compiled a 412 list of neighbors between all hexagons using the poly2nb function and then obtained the local 413 G statistic using the weighted density (normalised number of REs) of the global hexagons 414 and their relative distance from each other. The G statistic calculates a Z-score (measure of 415 standard deviation) for each hexagon. P values were then determined using the critical Z- 416 scores at 95% confidence levels followed by a Bonferroni correction using the p.adjustSP 417 function (using the number of neighbors between hexagons rather than the total number of 418 hexagons). All analyses were carried out using the R spdep package57. 419 420 We compared the analysis from the linear models with expert based polygons from three 421 different sources – IUCN, BIEN and EUFORGEN. A randomised weighted sample of all of 422 this dataset was used to generate a global distribution of REH by running this sample through 423 our pipeline. GLMs were run on the global distribution of expert-based REH in the same way 424 as with our generated polygons. Given the uneven distribution of expert-based polygons 425 globally (Supplementary Fig. 8a), we ran models excluding Asia and Africa, in order to 426 account for this bias. As seen by the strong similarity between the GLMs of expert polygons 427 in a global and filtered model (Supplementary Fig 8b and Fig 8c, panel 3), we observed an 428 overrepresentation of the expert polygons for these continents. Likewise, the GLMs from our 429 generated polygons are much similar to those obtained from the filtered model. In this case, 430 practically all variables showed the same relationship with range edges (either positive or 431 negative β values) as well as similar magnitudes. 432 433 Statistical analyses 434 Contribution of biome edge to RE. We used the 14 biomes defined by the World Wildlife 435 Fund (WWF) for our analyses. The distributions were downloaded from the WWF webpage 436 (http://www.worldwildlife.org/). To identify the intersection between biomes, we reduced the 13 437 complexity of the polygon edge using the rmapshaper package, which can perform 438 topologically aware polygon simplifications, thus maintaining the intersection between 439 biomes upon reduction of “edginess” of the polygons. The intersections were delineated and 440 subsequently enlarged (with a buffer distance of 0.1 arc-minute, ~185 m at the equator). A: 441 442 i. Global permutation assay was carried out by randomising (1000 iterations) the 443 global distribution of hexagons with REH (absolute Z-score) 444 ii. Per-biome stratified permutation was carried out by randomising the distribution 445 of hexagons with REH within each biome independently. 446 447 The averaged global bootstrap shown in panel 1 of Fig 2a, was calculated using a per-biome 448 bootstrap to obtain the probability distribution of REH at biome edges (Supplementary Fig. 449 3a). Distributions that were not normally distributed as a result of their small size (flooded 450 grassland, mangrove and tropical and subtropical coniferous forest; Supplementary Fig 3b), 451 were removed from the analysis. A general trend for the probability of range edges falling at 452 the intersection of biomes was therefore measured as a unified standardised z-distribution, 453 and compared to the median z-score from the actual percent overlap for each biome. 454 455 The density of hotspots at the intersection between biomes was calculated using the sum of 456 Z-scores of the hotspots at that intersection, and the percentage contribution was then 457 calculated using this value over the total Z-score at all biome intersections. The global or per- 458 biome 1000 permutation means were used as a normalising denominator for the values 459 obtained from our data set (Fig 2b, Supplementary Fig 4a). The denominator could be either 460 larger or smaller than the numerator, so we log-transformed the outcome to obtain a linear- 461 like relationship. Contribution within a biome was calculated the same way as for the 462 contribution at the edge, using the mean from a permutated assay to normalise for the 463 absolute value. 464 465 Climatic dependency of RE. We tested the relationship between RE density and climatic 466 features by assigning a set of environmental variables to each hexagon. We used the 467 bioclimatic attributes downloaded from WorldClim Global Climate Data58 at a resolution of 5 468 arc-minutes. The 19 BIOCLIM variables and elevations for each hexagon were extracted 469 using the R raster package. We also used the 16 ENVIREM variables described by 28, 470 downloaded from their website, at a resolution of 2.5 arc-minutes. Absolute climate for each 14 471 hexagon was obtained using the mean over all pixels. Climatic spatial heterogeneity (SH) 472 was calculated using the proportional variability index59,60 (PV) over all pixels in each 473 hexagon. 474 475 We used an array of linear mixed models (Fig 3, Supplementary Fig 5 and 6) to test for the 476 dependence of REs to climate and the robustness of the results. Linear mixed models were 477 carried out to account for biomes and continents. Both of these variables were introduced as 478 random effects in random intercept models. A model selection analysis was used to 479 determine the models that best predicted the formation of RE. Random intercept models 480 using both continent and biome as random variables were run, and models with ∆AIC <2 481 were selected for. Given the strong correlation between different predictor variables, we ran 482 models only with variable combinations that had a Pearson’s correlation value r < 0.7. The 483 relative contribution of the variable included in the model were calculated from the selected 484 models. Analyses were run using R package MuMIn61. 485 486 All statistical analyses (individual GLMs, and multiple-predictor GLMs) were tested for their 487 robustness to spatial autocorrelation by creating a spatial autocovariate (autocov_dist 488 function, spdep package), calculated as the distance-weighted average of neighboring 489 dependent variables62, so hexagons in proximity were averaged and those farther away 490 received a lower weighting average. We set the predetermined distance to 200 km based on 491 the average distance between cells. The spatial autocovariate was then included in the 492 regression model as a dependent variable. 493 494 Data Availability 495 The occurrence points used from GBIF can be found in the GBIF webpage 496 (https://doi.org/10.15468/dl.ajen6k). Polygons generated from occurrence points are provided 497 in the public GitHub repository https://github.com/dlernerg/Global-Range-edges. Biome 498 polygons were obtained from the WWF webpage (http://www.worldwildlife.org/). 499 Bioclimatic attributes were downloaded from WorldClim Global Climate Data58. ENVIREM 500 variables were downloaded from their webpage (https://envirem.github.io/). 501 502 Code Availability 503 Custom codes related to this paper can be found in a GitHub repository at 504 https://github.com/dlernerg/Global-Range-edges 15 505 506 Acknowledgments 507 TK wishes to thank Edith and Nathan Goldenberg Career Development Chair; Mary and Tom 508 Beck-Canadian Center for Alternative Energy Research; Larson Charitable Foundation New 509 Scientist Fund; Yotam Project; Dana and Yossie Hollander; Estate of Emile Mimran; and the 510 Estate of Helen Nichunsky. DL was supported by the Sustainability and Energy Research 511 Initiative PhD grant. MF-M and JP were supported by the PID2019-110521GB-I00 and 512 TED2021-132627B-I00 grants funded by MCIN, AEI/10.13039/501100011033 and the 513 NextGeneration EU/PRTR. MF-M. was supported by a postdoctoral fellowship from ”la 514 Caixa” Foundation (ID 100010434), code: LCF/BQ/PI21/11830010. 515 Author Contributions: D.L., J.B. and T.K. designed the research; D.L. performed the 516 research and analysed the data; M.F.M, J.P., T.K., J.B., and S.L.L. provided scientific advice; 517 M.F.M. and J.B. advised on statistical aspects; D.L. wrote the paper with special contribution 518 from J.P., J.B., M.F.M. and T.K. 519 Competing Interest Statement: The authors declare no conflict of interest. 520 521 Figure 1. Range-edge hotspots. Hotspots were identified using Getis-Ord Gi* analysis, 522 which returns a Z-score for each hexagon in the world. Only hexagons with p<0.05 are 523 considered hotspots and subsequently shown here. The Z-score of each hexagon is 524 represented by the color gradient. Biomes (as defined by the World Wildlife Fund) are 525 marked by colors. ‘T&sT’ and ‘Temp’ stand for ‘Tropic and Subtropic’ and ‘Temperate’, 526 respectively. 527 528 Figure 2. Range-edge hotspots at intersections between biomes. (a) Modeled distribution 529 of the percentage of hotspots at the edges of biomes from a permutated (randomised) per- 530 biome distribution of hotspots. The first panel represented the median value of all the other 531 biomes in the figure over a standardised z-distribution of biomes. The arrow marks the 532 percentage of RE hotspots at biome intersections in the data set (one-sided p-values) (b) Heat 533 maps of the percentage of RE hotspots (relative to the total number of hotspots) at the 534 intersection between two biomes. A biome-pair intersection with a significant number of 535 hotspots (p-value < 0.08) is marked with an asterisk. Biome-pairs that have no intersections 536 are gray. Panel 1 is normalised over a global bootstrap and panel 2 over a per-biome 16 537 bootstrap. (c) Correlations (and regression lines) of the relationships between the number of 538 REs at the edges of biomes and the number of REs within the biomes. The shaded area 539 represents the 95% confidence interval around the regression line. P-values are calculated 540 using a two-sided Student’s T-test (degrees of freedom = 12). See Methods and Materials for 541 further information on the methodology for obtaining significance levels (for (b)) and hotspot 542 permutations. ‘T&sT’ and ‘Temp’ stand for ‘Tropic and Subtropic’ and ‘Temperate’, 543 respectively. 544 545 546 Figure 3. Climatic predictors of range-edge formation. Models of RE formation using the 547 absolute climate (mean) and the spatial heterogeneity (PV index) of the 19 WorldClim 548 variables and elevation. (a) Individual (binomial) mixed regression models between each of 549 the predictor variables and number of REs (1-19 are BioClim variables) and accounting for 550 continents and biomes as random effects. The estimated coefficients of the explanatory 551 variables (β) are represented by the color gradient. ‘*’ and ‘**’ represent p-values <0.05 and 552 <0.01, respectively (two-sided Student T-test). 553 BioClim1 – Annual Mean Temperature, BioClim2 - Mean Diurnal Range (Mean of monthly 554 (max temp - min temp), BioClim3 – Isothermality, BioClim4 - Temperature Seasonality, 555 BioClim5 - Max Temperature of Warmest Month, BioClim6 - Min Temperature of Coldest 556 Month, BioClim7 - Temperature Annual Range, BioClim8 - Mean Temperature of Wettest 557 Quarter, BioClim9 - Mean Temperature of Driest Quarter, BioClim10 - Mean Temperature of 558 Warmest Quarter, BioClim11 - Mean Temperature of Coldest Quarter, BioClim12 - Annual 559 Precipitation, BioClim13 – Precipitation of Wettest Month, BioClim14 – Precipitation of 560 Driest Month, BioClim15 – Precipitation Seasonality (PV), BioClim16 – Precipitation of 561 Wettest Quarter, BioClim17 – Precipitation of Driest Quarter, BioClim18 – Precipitation of 562 Driest Quarter, BioClim19 – Precipitation of Coldest Quarter. 563 (b) Violin plots depicting the results from (a) for four predictor variables. The distribution of 564 continental climates is shown in gray, contrasted with the climatic distribution specific to the 565 RE hotspots (scaled to the intensity of the hotspot, i.e. the Z- score). (c) Forest plot of a 566 model average from the highest predicting LMM with both BioClim and ENVRIEM 567 variables (identified with a model selection). Beta values (log-odds) are shown for each 568 predictor. Absolute and SH climate was obtained from all of the inland global hexagonal 569 units (n = 5851). Error bars represent 95% confidence interval around the average effect. 570 17 571 References 572 1. Soule, M. The epistasis cycle: A theory of marginal populations. Annu Rev Ecol Syst 4, 165– 573 187 (1973). 574 2. Brown, J. H. On the Relationship between Abundance and Distribution of Species. Am Nat 575 124, 255–279 (1984). 576 3. Gaston, K. J. The Structure and Dynamics of Geographic Ranges. Oxford University Press 577 (2003). 578 4. Sexton, J. P., McIntyre, P. J., Angert, A. L. & Rice, K. J. Evolution and ecology of species range 579 limits. Annu Rev Ecol Evol Syst 40, 415–436 (2009). 580 5. Gaston, K. J. Geographic range limits: Achieving synthesis. Proceedings of the Royal Society B: 581 Biological Sciences 276, 1395–1406 (2009). 582 6. Goldberg, E. E. & Lande, R. Notes and Comments Species’ Borders and Dispersal Barriers. Am. 583 Nat vol. 170 (2007). 584 7. Bachmann, J. C., Rensburg, A. J. van, Cortazar-Chinarro, M., Laurila, A. & Buskirk, J. van. Gene 585 flow limits adaptation along steep environmental gradients. American Naturalist 195, E67– 586 E86 (2020). 587 8. Hargreaves, A. L., Samis, K. E. & Eckert, C. G. Are species’ range limits simply niche limits writ 588 large? A review of transplant experiments beyond the range. American Naturalist 183, 157– 589 173 (2014). 590 9. Henry, R. C., Bartoń, K. A. & Travis, J. M. J. Mutation accumulation and the formation of range 591 limits. Biol Lett 11, 11DUUMY (2015). 592 10. Perrier, A., Sánchez-Castro, D. & Willi, Y. Environment dependence of the expression of 593 mutational load and species’ range limits. J Evol Biol 35, 731–741 (2022). 594 11. Bontrager, M. et al. Adaptation across geographic ranges is consistent with strong selection in 595 marginal climates and legacies of range expansion. Evolution (N Y) 1–18 (2021) 596 doi:10.1111/evo.14231. 597 12. Santini, L., Pironon, S., Maiorano, L. & Thuiller, W. Addressing common pitfalls does not 598 provide more support to geographical and ecological abundant-centre hypotheses. 599 Ecography 42, 696–705 (2019). 600 13. Oldfather, M. F., Kling, M. M., Sheth, S. N., Emery, N. C. & Ackerly, D. D. Range edges in 601 heterogeneous landscapes: Integrating geographic scale and climate complexity into range 602 dynamics. Glob Chang Biol 26, 1055–1067 (2020). 603 14. Janzen, D. H. Why Mountain Passes are Higher in the Tropics WHY MOUNTAIN PASSES ARE 604 HIGHER IN THE TROPICS*. Source: The American Naturalist vol. 101. 605 15. Maxwell, M. F., Leprieur, F., Quimbayo, J. P., Floeter, S. R. & Bender, M. G. Global patterns 606 and drivers of beta diversity facets of reef fish faunas. J Biogeogr 49, 954–967 (2022). 607 16. Roy, K., Hunt, G., Jablonski, D., Krug, A. Z. & Valentine, J. W. A macroevolutionary perspective 608 on species range limits. Proceedings of the Royal Society B: Biological Sciences 276, 1485– 609 1493 (2009). 18 610 17. Loiseau, N. et al. Global distribution and conservation status of ecologically rare mammal and 611 bird species. Nat Commun 11, (2020). 612 18. Kerkhoff, A. J., Moriarty, P. E. & Weiser, M. D. The latitudinal species richness gradient in New 613 World woody angiosperms is consistent with the tropical conservatism hypothesis. Proc Natl 614 Acad Sci U S A 111, 8125–8130 (2014). 615 19. Donoghue, M. J. & Edwards, E. J. Biome shifts and niche evolution in plants. Annu Rev Ecol 616 Evol Syst 45, 547–572 (2014). 617 20. Ringelberg, J. J., Zimmermann, N. E., Weeks, A., Lavin, M. & Hughes, C. E. Biomes as 618 evolutionary arenas: Convergence and conservatism in the trans-continental succulent 619 biome. Global Ecology and Biogeography 29, 1100–1113 (2020). 620 21. Smith, J. R. et al. A global test of ecoregions. Nat Ecol Evol 2, 1889–1896 (2018). 621 22. Zhu, K., Woodall, C. W. & Clark, J. S. Failure to migrate: Lack of tree range expansion in 622 response to climate change. Glob Chang Biol 18, 1042–1052 (2012). 623 23. Corlett, R. T. & Westcott, D. A. Will plant movements keep up with climate change? Trends in 624 Ecology and Evolution vol. 28 482–488 Preprint at https://doi.org/10.1016/j.tree.2013.04.003 625 (2013). 626 24. la Sorte, F. A., Butchart, S. H. M., Jetz, W. & Böhning-Gaese, K. Range-Wide Latitudinal and 627 Elevational Temperature Gradients for the World’s Terrestrial Birds: Implications under 628 Global Climate Change. PLoS One 9, e98361 (2014). 629 25. Paquette, A. & Messier, C. The effect of biodiversity on tree productivity: From temperate to 630 boreal forests. Global Ecology and Biogeography 20, 170–180 (2011). 631 26. Pichancourt, J. B., Firn, J., Chadès, I. & Martin, T. G. Growing biodiverse carbon-rich forests. 632 Glob Chang Biol 20, 382–393 (2014). 633 27. Pennington, R. T., Lavin, M. & Oliveira-Filho, A. Woody plant diversity, evolution, and ecology 634 in the tropics: Perspectives from seasonally dry tropical forests. Annu Rev Ecol Evol Syst 40, 635 437–457 (2009). 636 28. Title, P. O. & Bemmels, J. B. ENVIREM: an expanded set of bioclimatic and topographic 637 variables increases flexibility and improves performance of ecological niche modeling. 638 Ecography 41, 291–307 (2018). 639 29. Veresoglou, S. D. & Peñuelas, J. Variance in biomass-allocation fractions is explained by 640 distribution in European trees. New Phytologist 222, 1352–1363 (2019). 641 30. Grantham, H. S. et al. Anthropogenic modification of forests means only 40% of remaining 642 forests have high ecosystem integrity. Nat Commun 11, (2020). 643 31. Holdridge, L. R. Determination of world plant formations from simple climatic data. Science 644 (1979) 105, 367–368 (1947). 645 32. Whittaker, R. H. Classification of Natural Communities. Botanical Review 28, 1–239 (1962). 646 33. McDonald, R. et al. Species compositional similarity and ecoregions: Do ecoregion boundaries 647 represent zones of high species turnover? Biol Conserv 126, 24–40 (2005). 19 648 34. von Humboldt, A. & Bonpland, A. Essay on the Geography of Plants. (2013). 649 35. Cardillo, M. Latitude and rates of diversifcation in birds and butterfies. Proc. R. Soc. Lond. B 650 266, 1221–1225 (1999). 651 36. Hillebrand, H. On the Generality of the Latitudinal Diversity Gradient. Am. Nat vol. 163 (2004). 652 37. Mittelbach, G. G. et al. Evolution and the latitudinal diversity gradient: Speciation, extinction 653 and biogeography. Ecology Letters vol. 10 315–331 Preprint at 654 https://doi.org/10.1111/j.1461-0248.2007.01020.x (2007). 655 38. Hewitt, G. M. Genetic consequences of climatic oscillations in the Quaternary. in 656 Philosophical Transactions of the Royal Society B: Biological Sciences vol. 359 183–195 (2004). 657 39. Crane, P. & Scott, L. Angiosperm Diversification and Paleolatitudinal Gradients in Cretaceous 658 Floristic Diversity. Science (1979) 246, 675–678 (1989). 659 40. Jablonski, D. The tropics as a source of evolutionary novelty through geological time. P. J. IEEE 660 Trans. Geosci. Remote Sensing 361, 180–191 (1993). 661 41. Jablonski, D. et al. Out of the tropics, but how? Fossils, bridge species, and thermal ranges in 662 the dynamics of the marine latitudinal diversity gradient. Proc Natl Acad Sci U S A 110, 663 10487–10494 (2013). 664 42. Antonelli, A. et al. An engine for global plant diversity: Highest evolutionary turnover and 665 emigration in the American tropics. Front Genet 6, (2015). 666 43. Jump, A. S. & Peñuelas, J. Running to stand still: Adaptation and the response of plants to 667 rapid climate change. Ecol Lett 8, 1010–1020 (2005). 668 44. Morreale, L. L., Thompson, J. R., Tang, X., Reinmann, A. B. & Hutyra, L. R. Elevated growth and 669 biomass along temperate forest edges. Nat Commun 12, 7181 (2021). 670 45. Wilkinson, S., Clephan, A. L. & Davies, W. J. Rapid Low Temperature-Induced Stomatal Closure 671 Occurs in Cold-Tolerant Commelina communis Leaves But Not in Cold-Sensitive Tobacco 672 Leaves, via a Mechanism That Involves Apoplastic Calcium But Not Abscisic Acid. Plant Physiol 673 126, 1566–1578 (2001). 674 46. Brodribb, T. J. & Holbrook, N. M. Stomatal protection against hydraulic failure: a comparison 675 of coexisting ferns and angiosperms. New Phytologist 162, 663–670 (2004). 676 47. Davis, B. A. S. & Brewer, S. Orbital forcing and role of the latitudinal insolation/temperature 677 gradient. Clim Dyn 32, 143–165 (2009). 678 48. Seager, R. et al. Strengthening tropical Pacific zonal sea surface temperature gradient 679 consistent with rising greenhouse gases. Nat Clim Chang 9, 517–522 (2019). 680 49. Xu, Y. & Ramanathan, V. Latitudinally asymmetric response of global surface temperature: 681 Implications for regional climate change. Geophys Res Lett 39, n/a-n/a (2012). 682 50. Colwell, R. K., Brehm, G., Cardelús, C. L., Gilman, A. C. & Longino, J. T. Global Warming, 683 Elevational Range Shifts, and Lowland Biotic Attrition in the Wet Tropics. vol. 322 684 www.sciencemag.org (2008). 20 685 51. Basso, B., Martinez-Feria, R. A., Rill, L. & Ritchie, J. T. Contrasting long-term temperature 686 trends reveal minor changes in projected potential evapotranspiration in the US Midwest. 687 Nat Commun 12, 1476 (2021). 688 52. Zizka, A. et al. CoordinateCleaner: Standardized cleaning of occurrence records from 689 biological collection databases. Methods Ecol Evol 10, 744–751 (2019). 690 53. Serra-diaz, J. M., Enquist, B. J., Maitner, B., Merow, C. & Svenning, J. Big data of tree species 691 distributions : how big and how good ? (2018) doi:10.1186/s40663-017-0120-0. 692 54. Zizka, A. et al. No one-size-fits-all solution to clean GBIF. PeerJ 8, (2020). 693 55. Getis, A. & Ord, J. K. The Analysis of Spatial Association by Use of Distance Statistics. Geogr 694 Anal (1992). 695 56. Mendez, C. Spatial autocorrelation analysis in R. R Studio/RPubs. https://rpubs.com/quarcs- 696 lab/spatial-autocorrelation (2020). 697 57. Bivand, R. S., Pebesma, E. & Gómez-Rubio, V. Applied Spatial Data Analysis with R. (2013). 698 58. Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution 699 interpolated climate surfaces for global land areas. International Journal of Climatology 25, 700 1965–1978 (2005). 701 59. Heath, J. P. Quantifying temporal variability in population abundances. Oikos vol. 115 573– 702 581 Preprint at https://doi.org/10.1111/j.2006.0030-1299.15067.x (2006). 703 60. Fernández-Martínez, M. et al. The consecutive disparity index, D: a measure of temporal 704 variability in ecological studies. Ecosphere 9, (2018). 705 61. Bartoń, K. MuMIn: Multi-model inference. R package version 1.10.0.. 1. . (2013). 706 62. F. Dormann, C. et al. Methods to account for spatial autocorrelation in the analysis of species 707 distributional data: A review. Ecography vol. 30 609–628 Preprint at 708 https://doi.org/10.1111/j.2007.0906-7590.05171.x (2007). 709 710 21