A comparison of macroecological and stacked species distribution models to predict future global terrestrial vertebrate richness

Predicting future changes in species richness in response to climate change is one of the key challenges in biogeography and conservation ecology. Stacked species distribution models (S‐SDMs) are a commonly used tool to predict current and future species richness. Macroecological models (MEMs), regression models with species richness as response variable, are a less computationally intensive alternative to S‐SDMs. Here, we aim to compare the results of two model types (S‐SDMS and MEMs), for the first time for more than 14,000 species across multiple taxa globally, and to trace the uncertainty in future predictions back to the input data and modelling approach used.


| INTRODUC TI ON
One of the current major challenges in biogeography is to understand and predict the potential impacts of global change on the distribution of biological diversity. In addition to land-use change and its consequences for natural habitats, climate change has been identified as one of the most prominent drivers of biodiversity change (IPBES, 2019;Sala et al., 2000). Changes in biological systems in response to climate change are frequently documented, including shifts in species distribution (Chen, Hill, Ohlemüller, Roy, & Thomas, 2011). Such changes in species distributions ultimately result in changes of biodiversity patterns, such as the geographical variation in species richness.
Species richness is a simple albeit important measure for biodiversity and has been identified as one of the Essential biodiversity variables (EBVs; Pereira et al., 2013). EBVs function as an interface between raw data and indicators and are meant to provide robust and coordinated data about biodiversity change on a global scale in order to inform policy makers (Brummitt et al., 2017;Geijzendorffer et al., 2016). EBVs require representative sampling across taxonomic groups especially for assessing changes in ecosystem services. Field studies and monitoring schemes provide data for a wide range of EBVs, but are often limited in spatial or temporal coverage, whereas large-scale data sources, such as the Global Biodiversity Information Facility (GBIF), are inherently biased (Meyer, Kreft, Guralnick, & Jetz, 2015) and not representative on a global scale (Proença et al., 2017).
Global climate models predict an increase in global mean surface temperature of up to 4.5°C by the end of this century compared to today (IPCC, 2013). This would be comparable to the difference between the last glacial maximum and the pre-industrial climate (Otto-Bliesner et al., 2006;Shakun & Carlson 2010). Given that species distributions and thus species richness have changed over the last century in response to an increase in the global average temperature of about 1°C, these changes are likely to continue in response to future climate change as well, even under optimistic (but currently unlikely) scenarios such as the ones in line with the 2°C or 1.5°C targets of the Paris Agreement (Hof et al., 2018). Thus, reliable predictions of future species richness under different climate change scenarios are of great need.
Species range shifts in response to climate change have been found across a wide array of taxa, with the majority of species shifting their distribution towards higher latitudes and altitudes (Chen et al., 2011); but also idiosyncratic responses have been observed (Dobrowski et al., 2013). To estimate potential shifts in species distributions as a response to climate change, current species distributions and climatic data can be fed into statistical models to infer the climatic niche of a species (species distribution models (SDMs)) or the relationship between climate and species richness (macroecological models (MEMs)).
Modelling distributions for all taxa present in a region of interest and aggregating them to a single species richness layer, an approach termed 'stacked SDMs' (S-SDMs), allows predicting current and future species richness (Ferrier & Guisan, 2006). However, besides the various assumptions inherent to SDMs, which have been discussed elsewhere in more detail (Elith & Leathwick, 2009); we see two major drawbacks of S-SDMs as a tool for species richness projections. First, the potential errors of SDMs add up when stacking the output of multiple SDMs, which may increase the potential error of species richness estimates (Pineda & Lobo, 2009). Second, SDMs are, by definition, species-specific and thus require occurrence information for each of the investigated species. Furthermore, they require a minimum number of occurrence points, making them unsuitable for modelling small-ranging species Wisz et al., 2008), which constitute significant amounts of the species numbers of many taxonomic groups (Platts et al., 2014). Thus, sufficient species-specific information is lacking to reliably project species richness based on species richness. Model type by far contributes to most of the variation in the different future species richness predictions, indicating that the two model types should not be used interchangeably. Nevertheless, both model types have their justification, as MEMs can also include species with a restricted range, whereas S-SDMs are useful for looking at potential species-specific responses.

K E Y W O R D S
biodiversity, climate change, cluster analysis, macroecological model, richness model, species distribution model, species richness, terrestrial vertebrates, variance partitioning S-SDMs for the vast majority of the world's taxonomic groups, with only a few exceptions (such as terrestrial birds or mammals; Beck et al., 2012).
They further only require total richness values and no species-specific information. This is particularly useful in cases where morphospecies or inaccurate occurrence data preclude the application of species-specific distribution models, but allow rough estimates of the spatial variation in species richness. MEMs have a long history in macroecological research (Currie, 1991) and are often associated with a discussion on the factors that drive large-scale species richness patterns (Currie, 1991;Currie et al.,. 2004;Francis & Currie, 2003;Hawkins et al., 2003;Jetz & Rahbek, 2002;Rahbek & Graves, 2001;Rahbek et al., 2007;Rangel et al., 2018;Thuiller, Midgley, Rouget, & Cowling, 2006;Wright, 1983).
In this study, we aim to assess whether MEMs may be a reliable tool for predicting current as well as future species richness patterns in comparison to S-SDMs. Several studies have addressed this question under current conditions, based on different approaches. For example Guisan and Rahbek (2011) found that S-SDMs consistently over-predict species richness and suggested to combine S-SDMs with MEMs to derive more reliable species richness estimates.
Inspired by the findings of Rahbek (2011), Calabrese, Certain, Kraan, andDormann (2014) showed that the over-prediction of S-SDMs is due to the thresholding of individual occurrence probabilities and can be avoided by stacking the non-thresholded probability values. More recently, Harris et al. (2018) showed that MEMs perform better than S-SDMs when predicting future breeding bird richness of North America.
While a comparison of MEM and S-SDM predictions is rather simple for current conditions, as they can be compared with observed richness values, quality assessments of future predictions remain challenging, as there are no future reference values to compare with (but see Harris et al., 2018). However, an exploration of the commonalities and differences of predictions rendered by dif-

| MATERIAL S AND ME THODS
The methods for species data, climate data and S-SDMs were performed following the same procedure as Hof et al. (2018) and thus are only described briefly here.

| Species data
Species presence data for three taxonomic groups (amphibians, birds and mammals) were derived from expert extent-of-occurrence (EOO) range maps provided by BirdLife International and NatureServe (2015) and the International Union for Conservation of Nature (2016). Range maps were gridded to a spatial resolution of 0.5° using the raster package (Hijmans 2017) in R (R Core Team, 2018). Polygons were first rasterized as lines and then as polygons in order to include also cells that are touched by the polygon boundaries, which is particularly important when considering coastal areas and islands. Only polygons where a species was extant or probably extant, occurred natively and was resident or occurred regularly during the breeding season were considered (this information is part of the expert range map data).
Pseudo-absence data for each species was generated by randomly selecting absences, equal to the number of presences or 1,000 absences for all species with less than 1,000 presences, using a distance-weighted approach, where the probability of randomly selecting a point decreases by 1/(De^2) where De is the distance from the range edge. This way the risk of only sampling absences close to the range of a species was reduced (Barbet-Massin, Jiguet, Albert, & Thuiller, 2012;Thuiller, 2004), whereas at the same time avoiding to extensively sample too far beyond the range of a species, where absences are likely to occur due to non-bioclimatic reasons (Anderson & Raza, 2010).

| Bioclimatic variables
Current bioclimatic variables were derived from 30-year (1980-2009) monthly means of minimum temperature, maximum temperature and precipitation extracted from the meteorological forcing dataset 'EartH2Observe, WFDEI and ERA-Interim data Merged and Bias-corrected for ISIMIP' (EWEMBI; Lange, 2016) using the 'dismo' package (Hijmans, Phillips, Leathwick, Elith & Hijmans, 2017) in R (R Core Team, 2018). The EWEMBI dataset provides bias-corrected global daily climate data at a spatial resolution of 0.5° and was specifically compiled for impact assessments of a 1.5°C global warming above pre-industrial levels (Lange, 2016). For a representative subset of 10% of the species of each taxon we built Generalized Additive Models (GAMs) for each of these variable combinations and for each taxon selected the variable combination, resulting in the largest number of species with SDM models of high accuracy (highest area under the curve (AUC; Fielding & Bell, 1997)). For birds and mammals, the final variables were temperature seasonality (bio4), maximum temperature of the warmest month (bio5), annual precipitation (bio12) and precipitation seasonality (bio15). For amphibians, the variables were temperature seasonality (bio4), maximum temperature of the warmest month (bio5), precipitation of the warmest quarter (bio18) and precipitation of the coldest quarter (bio19; Figure S1).
Species distribution data often exhibit spatial autocorrelation, which can bias parameter estimates and error probabilities (Kühn, 2007). We reduced the effect of spatial autocorrelation in the SDMs by applying two different methods. For species with more than 50 presences, we divided the world into 10 blocks, based on a representative subset of the climatic space of each of the world's ecoregions, as defined by Olson et al. (2001) and built 10 models leaving one block out at a time, using the left out block for model evaluation (Bagchi et al., 2013). For range-restricted species (≤50 presences), we split the data into 10 datasets by repeatedly randomly selecting 70% of the data, using the left-out 30% for model evaluation.
Species occurring in less than 10 grid cells were not considered in this analysis (3,318 amphibians, 896 birds, 968 mammals).
The performance of the fitted SDMs was evaluated by calculating the overall AUC for each species (the average AUC across the 10 blocks and the 10 sets of pseudo-absences) and models with an overall AUC smaller than 0.7 were dropped (64 amphibians, 149 birds, 332 mammals). This left us with SDMs for 2,964 amphibian, 8,493 terrestrial bird and 4,039 terrestrial mammal species, which represents more than 85% of the entire species set across more than 80% of the considered area (amphibians = 96.8%, birds = 83.3%, mammals 90%, Figure S5). The use of AUC as a correct metric for model evaluation is debatable (Hirzel, Le Lay, Helfer, Randin, & Leroy et al., 2018;Lobo, Jiménez-Valverde, & Real, 2008), as is the use of other metrics (Fourcade, Besnard, & Secondi, 2018). However, given that we have performed a rigorous variable selection procedure and a spatial blocking cross-validation, as suggested by Fourcade et al. (2018) to avoid the inflation of SDM performance scores, we are confident that the selected models are robust.
Future species distributions were derived by predicting the models using future bioclimatic variables. They were limited by the extent of neighbouring zoogeographic realms, as defined by Holt et al. (2013), to avoid predictions for areas that mirror analogue climatic conditions.
To account for the unlikely assumption of unlimited dispersal of each species within realms, we also applied a species-specific dispersal buffer. Species-specific dispersal distances are still unavailable for most species considered here (Nathan, Klein, Robledo-Arnuncio, & Revilla, 2012), thus we limited the dispersal of each species by applying a buffer of d/4 to the range polygon, where d is the diameter of the largest range polygon of a species, and by clipping the current and future predictions of each species by the buffered range polygon. This approach builds upon previous studies (i.e. Barbet-Massin & Jetz, 2015), but rather than taking one or multiple buffer distances which are kept constant for all species, a species-specific buffer distance according to the species' largest range extent was used (also see Hof et al., 2018). This is in line with a study by Whitmee and Orme (2013), who found that dispersal distance is species-specific and together with home range area and body mass is mostly explained by geographical range size. To verify this approach, we explored the impact of choosing a specific dispersal correction factor on the results by comparing the low dispersal scenario (d/4) to several larger dispersal scenarios (d/2, d and 2*d) as well as to a no dispersal and a full dispersal scenario ( Figures S2, S3, S7 and S8).
The current and future raw probabilities of occurrence of the individual SDMs were then stacked without thresholding, following the procedure suggested by Calabrese et al. (2014), to derive current and future predictions of species richness for each of the three taxa.

| MEMs
For the MEMs, we used the same gridded range data, but only of the species that we used for the S-SDMs, and combined these to create gridded species richness data for each taxon. The resulting species richness values per grid cell were then used as response variable for the MEMs using the same explanatory variables as for the SDMs to allow for a direct comparison. In order to avoid the violation of key statistical assumptions (Dormann, 2007), spatial autocorrelation was again reduced by applying the ecoregion-blocking approach described above (Bagchi et al., 2013;Hof et al., 2018).
MEMs were fitted using the same modelling approaches as for the SDMs. GAMs were this time fitted using a Poisson response with a logit link using thin-plate regression splines (Wood, 2003(Wood, , 2006(Wood, , 2011. GBMs were again conducted using the 'gbm' package in R (R Core Team, 2018; Ridgeway, 2017) following the same procedure as for the SDMs, but using a Poisson distribution. The model parameters were again optimized using cross-validation, but this time using three different learning rates (0.1, 0.01, 0.001), as otherwise optimal models had too few trees (number of trees <1,000).
MEMs were fitted and predicted using the same current and future climatic data as for the SDMs. Dispersal scenarios do not apply to MEMs, as they are inherently included in the models.
To evaluate the model fit of the MEMs we calculated the Root Mean Square Error (RMSE; Wilmott, 1981) between EOO-based and modelled species richness (see Table S1).
In addition to the MEMs using species richness based on the species sets modelled in the S-SDMs, we also ran MEMs using species richness of each taxon with the respective entire species set (6,381 amphibians, 9,885 birds, 5,276 mammals) as explanatory variable to see how incorporating all species would influence the output of our MEMs (see Figures S6, S15 and S16).

| Model comparison
While predictions of current species richness can be evaluated by the observed data used for fitting the models, it is impossible to compare the performance of future predictions of S-SDMs and MEMs, due to the lack of future data for validation. To further identify the sources that contribute to the variance in the patterns of predicted future species richness across each grid cell, we performed a three-way Analysis of Variance (ANOVA) without replication (Legendre & Legendre, 1998;Sokal & Rohlf, 1995)  separately. The proportion of the sum of squares from each of these sources (and their interactions) of the total sum of squares is an estimate of the variance that can be attributed to one of the sources (see Table 1). Note that the variance determined by the full interaction cannot be differentiated from the residual variance (the part of variance that is not explained by any of the factors or their interactions).
Doing this for every grid cell individually allows us to display spatial differences and so identify regions of low and high variance among the different sources.
TA B L E 1 Relative contributions to the overall variation in predicted future species richness from different sources of variance Note: Values represent median proportions (%) of the total sum of squares from the three-way ANOVA performed for each grid cell evaluating the relative contributions of model type, model algorithm, general circulation model (GCM), as well as their interactions to the variance of predicted future species richness for the two different RCPs, separately for each taxon. Note that the variance determined by the full interaction cannot be differentiated from the residual (unexplained) variance.
One major benefit of MEMs is that they can be applied for all species (see Figures S6, S15 and S16), including rare and smallranging ones, which cannot be modelled using SDMs. Therefore, we performed a sensitivity analysis on how the species coverage (see Figure S5) affects the similarity and variance contribution of the different future predictions (see Figures S21 and S22, Table S2).

| Current patterns of species richness
Current ensemble predictions of species richness based on S-SDMs and MEMs both had a high correlation with EOO-based species richness across all three taxa (R 2 = 0.77-0.92; Figure 1a-f and Table S1).
F I G U R E 1 EOO-based versus current predicted species richness for each grid cell per taxon for ensemble models (average among GCMs and model algorithms) of (a), (b), (c) MEMs and (d), (e), (f) S-SDMs and (g), (h), (i) the correlation between S-SDMs and MEMs. S-SDMs were performed using a low dispersal scenario (d/4). MEMs were performed using the same species as for the S-SDMs. Black lines represent the fit of the respective regression, with the model equation and the R 2 value given in black as well. Perfect fits would have intercepts of 0, slopes of 1 and a R 2 value of 1 (dashed grey line) However, both model types over-predicted richness in species-poor sites and under-predicted richness in species-rich sites (Figure 1).
The relationship between EOO-based and predicted current richness varied among the different model algorithms used. For the S-SDMs, the GAM algorithm showed a better fit (intercept closer to 0 and slope closer to 1) than the GBM algorithm, whereas for the MEMs, it was vice versa (Table S1).
For S-SDMs the relationship between EOO-based and predicted richness under current conditions varied with the dispersal assumption used. Under a no-dispersal scenario S-SDMs always under-predicted richness, whereas dispersal scenarios with a large dispersal buffer (d/2, d, 2*d and full dispersal) always over-predicted richness ( Figures S2 and S3). Under a low dispersal scenario (d/4), current predictions derived from S-SDMs provided a better fit with EOO-based species richness than MEMs. This was consistent across all three taxa ( Figure 1) and across the different model algorithms (Table S1).
Nevertheless, current predictions of MEMs and S-SDMs showed a strong correlation (R 2 > 0.8; Figure 1g-i).
When we compared the spatial patterns of EOO-based and predicted current richness, we found that both MEMs and S-SDMs largely reproduced the global variation in vertebrate richness ( Figure 2, Figure S4). S-SDMs provided a better fit with EOO-based richness patterns (53.3%-77.6%) compared to MEMs (22.4%-46.7%). S-SDMs were in particular much better in areas of high species richness, that is across most parts of South America ( Figure 2).
Running MEMs for all vertebrate species overall resulted, as expected, in on average about 8.55% (2.27% amphibians, 12.4% birds, 9.64% mammals) higher richness estimates, although in some areas predicted current species richness got smaller ( Figure S6).

| Future patterns of species richness
The correlations between future richness predictions from S-SDMs and MEMs were weaker (R 2 = 0.78-0.81) than those between current richness predictions. This was also reflected by the spatial pattern of future predictions, which showed large deviations in particular in areas of high species richness, i.e. central South America and central Africa, for both RCPs (Figure 3, Figures S7 and S8). Future predictions derived from S-SDMs varied with the dispersal assumption used and under a low dispersal scenario (d/4) provided the best fit with predictions derived from MEMs ( Figures S9 and S10). The correlation and spatial consistency of future richness predictions from MEMs and S-SDMS was even weaker when looking at the absolute (R 2 < 0.55) and relative change (R 2 < 0.5) in species richness between current and future predictions (see Figures S11-14). Future predictions from MEMs for the entire species set differed on average by around 8% (2% amphibians, 12% birds, 9% mammals) from MEMs that used only the species that could also be modelled with SDMs ( Figures S15 and S16). MEMs showed a stronger similarity in predictions across model algorithms and GCMs than S-SDMs. All these patterns were consistent across the two RCPs (Figure 4, Figure S17). GCM and the full interaction was negligible. All patterns were consistent across the two RCPs considered (Table 1).
Looking at the sources of variance in species richness from a spatial perspective, we found that the used model type mostly explained the variance in the majority of areas across the globe.
However, there were also some areas, that is along the Andes, east  Table S2).

| Current patterns of species richness
MEMs and S-SDMs both provided reliable predictions of EOO-based global vertebrate richness at a 0.5° resolution. This is in line with previous studies (Calabrese et al., 2014;Distler et al., 2015;Dubuis et al., 2011), which compared the fit of MEMs and S-SDMs across various scales and taxonomic groups.
MEMs and S-SDMs both under-predicted areas of low species richness and over-predicted areas of high species richness ( Figure 1a-f), which has also been found by Harris et al. (2018). Guisan and Rahbek (2011) argued that MEMs could be used to improve S-SDM predictions, as the latter consistently over-pre- MEMs performed worse for taxa where species richness was low (in our case amphibians) than for taxa with higher species richness values (birds and mammals). This appears to be in line with a study by Rahbek et al. (2007), which showed that models for the highest richness-quartile of birds performed better than models for the lower richness-quartiles. This might also be the reason why Da It has to be noted that for our analyses we had to exclude a large number of species which could not be modelled by SDMs due to their small geographical ranges, that is their low number of occurrence records. We also excluded these species from the MEM analyses to allow for consistency in the comparison between MEMs and S-SDMs, but additionally also provided the MEM results for the entire species set (Figures S6, S15 and S16). Even though S-SDMs outperformed MEMs in the prediction of species richness especially to include these range-restricted species is an advantage, as they can capture the impacts of climate change on entire taxa. However, we stress that range-restricted species are more affected by habitat type, whereas large-ranging species are more limited by climatic zones and biome types (Brown & Maurer, 1989); thus applying climate-only models to small-ranging species might not reflect biological patterns.
Spatial variation in species richness is not only driven by climate, which we considered here, but also by other factors such as topography (Davies et al., 2007;Rahbek & Graves, 2001), productivity (Coops, Kearney, Bolton, & Radeloff, 2018) and land use (Kehoe et al., 2017). The factors which explain the variation in species richness depend strongly on the spatial scale considered (i.e. see Chase et al., 2018). Both, MEMs and S-SDMs could potentially be improved using high-resolution climatic and non-climatic factors. Baudraz et al. (2018) showed that land use and very highresolution topo-climatic factors can improve MEMs for predicting mountain grassland species richness in the Swiss Alps, whereas Seo, Thorne, Hannah and Thuiller (2008) highlighted the influence of spatial scale on the accuracy of SDM predictions. However, for SDMs the appropriate spatial resolution depends on the species considered, as large and mobile organisms might be well-represented by large-scale climatic conditions, whereas small and less mobile species might not (Nadeau, Urban, & Bridle, 2017). This also strongly depends on the data source considered, as range maps at a high resolution typically result in incorrect spatial patterns of species richness (Hurlbert & Jetz, 2007). Thus, models based on expert range maps should ideally be performed on a rather coarse resolution, and high-resolution global occurrence data for a representative number of vertebrate species is currently still unavailable.

| FUTURE PAT TERN S OF S PECIE S RICHNE SS
Future predictions of MEMs and S-SDMs also showed a high correlation (R 2 > 0.75, Figure 3c,f,i), but were quite different when looking at the predicted future change in species richness ( Figures   S11-14). This is in contrast to the results from Distler et al. (2015), who found that predicted changes in summer and winter bird species richness of North America are consistent across S-SDMs and MEMs. The difference between their and our study might be due to the fact that Distler et al. (2015) did not use range maps, but point occurrence data from the Audubon Christmas Bird Count and North American Breeding Bird Survey. This leads to the question whether range maps should be used in SDMs, which has been extensively discussed elsewhere (e.g. Ficetola et al., 2014;Fourcade, 2016;Gaston & Fuller, 2009;Herkt, Skidmore, & Fahr, 2017;Pineda & Lobo, 2012). Despite their limitations, expert range maps provide one of the most comprehensive global biodiversity datasets with little spatial and taxonomic bias, unlike global data collection initiatives, like GBIF (Meyer et al., 2015).
There is a growing body of literature comparing range maps with different sources of occurrence records for different regions and taxa (i.e. see Barbosa, Estrada, Márquez, Purvis, & Orme, 2012;Fourcade, 2016;Herkt et al., 2017;Meyer et al., 2015), coming to mixed conclusions; however, EOO range maps are still widely used in macroecological research (Belmaker & Jetz, 2015;Hof et al., 2018;Slavenko & Meiri, 2015;Thuiller et al., 2019;Torres-Romero & Olalla-Tárraga, 2015;Zurell et al., 2018). Being aware of these recent discussions and acknowledging the limitations and caveats when using EOO range maps, we think that they remain valid for our purpose of comparing different modelling types for coarse species richness predictions.
Evaluating whether S-SDMs or MEMs provide better predictions of future species richness is a challenge that will remain unresolved, due to the lack of future data for validation. However, Harris et al. (2018) used observed data from two time spans (1982-2003 vs. 2004-2013)  often termed niche-conservatism, but also that the distribution of a species is well-known and near equilibrium (Holt, 2009). However, niche-conservatism is strongly debated (Pearman et al., 2008;Wiens & Graham, 2005), especially the question whether niches will remain stable under novel future climates (Veloz et al., 2012).In contrast to S-SDMs, MEMs try to establish a direct relationship between environmental variables and species richness and assume that the environment constrains the number of species that can co-exist in an area (Guisan & Rahbek, 2011). This relationship is mainly driven by climate, productivity, environmental heterogeneity, disturbance and history (Currie et al., 2004;Field et al., 2009) and strongly depends on the spatial scale considered (Field et al., 2009;Whittaker, Willis, & Field, 2001). The direct underlying mechanisms behind this are still discussed, but are obviously driven by evolutionary and biogeographic processes, that is speciation, extinction and dispersal, which again assume niche conservatism (Wiens & Donoghue, 2004).
These inherent differences between S-SDMs and MEMs are also reflected in the variable contribution to the model fit, which is quite different among the two model types (Distler et al., 2015). While both approaches do not consider population dynamics, S-SDMs further have the advantage that they consider species-specific responses to climatic change rather than one general response as in MEMs.
When partitioning the variance of the future predictions of species richness into different sources, we found that model type was the largest variance source, followed by GCM and model algorithm (

| CON CLUS IONS
With this study, we aimed to compare patterns of species richness predicted by macroecological models (MEMs) and stacked species distribution models (S-SDMs). Our findings show that under current conditions the patterns rendered by the two model types do closely resemble the global variation in EOO species richness. However, the two model types produce less collinear predictions of future species richness and model type by far contributes to most of the variation among future predictions. This suggests that both approaches are valid for coarse estimates of species richness at large geographic scales. However, the failure of both approaches to capture the full spectrum of global species richness variation calls for caution when trying to predict future richness, especially for the regions harbouring the highest number of species.
MEMs are a potentially attractive alternative to S-SDMs because of the lower number of model-inherent assumptions, the drastic reduction in computational time compared to S-SDMs, and their ability to be applied to taxa and regions for which species-richness estimates, but no species-specific occurrence data are available. As MEMs do not perform much worse than S-SDMs under current conditions, this suggests that they may be useful as a first quick explorative analysis or for outlining future species richness patterns, where running multiple S-SDM is unfeasible.
Since, as for any modelling exercise, it is impossible to assess the quality of future predictions directly, and predictions between model types vary considerably when looking at predicted changes in species richness, specific findings need to be interpreted with great care.

ACK N OWLED G EM ENTS
We thank the Inter-Sectoral Impact Model Intercomparison Project