Global drivers of population density in terrestrial vertebrates

Aim: Although the effects of life history traits on population density have been investigated widely, how spatial environmental variation influences population density for a large range of organisms and at a broad spatial scale is poorly known. Filling this knowledge gap is crucial for global species management and conservation planning and to understand the potential impact of environmental changes on multiple species. Location: Global. Time period: Present. Major taxa studied: Terrestrial amphibians, reptiles, birds and mammals. Methods: We collected population density estimates for a range of terrestrial vertebrates, including 364 estimates for amphibians, 850 for reptiles, 5,667 for birds and 7,651 for mammals. We contrasted the importance of life history traits and environmental predictors using mixed models and tested different hypotheses to explain the variation in population density for the four groups. We assessed the predictive accuracy of models through cross-validation and mapped the partial response of vertebrate population density to environmental variables globally. Results: Amphibians were more abundant in wet areas with high productivity levels, whereas reptiles showed relatively higher densities in arid areas with low productivity and stable temperatures. The density of birds and mammals was typically high in temperate wet areas with intermediate levels of productivity. The models showed good predictive abilities, with pseudo- R 2 ranging between 0.68 (birds) and 0.83 (reptiles). Main conclusions: Traits determine most of the variation in population density across species, whereas environmental conditions explain the intraspecific variation across populations. Species traits, resource availability and climatic stability have a different influence on the population density of the four groups. These models can be used to predict the average species population density over large areas and be used to explore macroecological patterns and inform conservation analyses.

Across species, body mass is known to scale negatively with population density, presumably because larger organisms have higher per-capita energy requirements (Blackburn et al., 1993;Currie & Fritz, 1993;Damuth, 1981;Silva & Downing, 1995). Nonetheless, the shape of the relationship with body mass has been extensively debated because a number of confounding factors might play a role, such as sampling area (Blackburn & Gaston, 1996), species co-occurrence and thus resource partitioning, and bias towards high density estimates in the literature (Lawton, 1990;White et al., 2007). Some studies supported a nonlinear relationship with body mass, with density reaching its maximum at intermediate masses and then decreasing asymptotically (Lawton, 1990;Silva & Downing, 1995;White et al., 2007). Trophic levels can alter the body mass-density relationship by influencing the energy conversion efficiency, the availability and accessibility of resources, and the energy investments in foraging, resulting overall in lower densities in higher trophic levels (Carbone & Gittleman, 2002;Silva et al., 1997).
What is less known is the role of environmental conditions and resource availability in shaping population densities (Table 1; Carbone & Gittleman, 2002;Currie & Fritz, 1993;Pettorelli, Bro-Jørgensen, Durant, Blackburn, & Carbone, 2009;Silva, Brimacombe, & Downing, 2001). Resource availability, as expressed by primary productivity for Body mass BM Larger body mass requires larger energy investment; therefore, body mass is expected to decrease with density (Blackburn et al., 1993;Currie & Fritz, 1993;Damuth, 1981;Silva & Downing, 1995). The relationship has been observed to be nonlinear in birds and mammals, with shallower relationships at one or both extremes of the body mass range (Silva et al., 1997) Diet category Diet Diet influences energy availability and accessibility to different trophic levels, and the energy transformation efficiency (Carbone & Gittleman, 2002;Silva et al., 1997). At equal body mass, lower trophic levels are expected to show higher density. Including diet as a categorical variable is thus expected to influence the intercept of the relationship Richness in species Rich Species richness directly affects the energy available for individual species (distant proxy of competition). Thus, it is expected to show a negative relationship with population density. In principle, it also determines the number of prey/predators (prey availability/predatory pressure), but this is likely to be compensated by a lower relative abundance and greater richness of prey Net primary productivity NPP NPP is directly or indirectly related to energy available for individuals (Blackburn & Gaston, 2001;Carbone & Gittleman, 2002;Pettorelli et al., 2009). It is therefore expected to show a positive relationship with population density. The relationship can be nonlinear because resources can be a limiting factor only when they are scarce, whereas population density can be top-down controlled when resources are abundant Annual mean temperature T m Ambient temperatures are usually below the physiological optima of ectotherms; therefore, an increase in temperature is expected to increase their fitness. To a certain extent, the increase in temperature can benefit endotherms by reducing their energy expenditure for thermoregulation, but very high temperatures and heatwaves can increase their mortality (Bowler et al., 2017;Currie & Fritz, 1993). The relationship is thus expected to be positive, and possibly asymptotic or hump shaped in endotherms Temperature or precipitation seasonality T sd or P cv Seasonality in temperature and precipitation can increase mortality either directly or indirectly through resource fluctuations (Adler & Levins, 1994;Williams S. E. & Middleton, 2008). The relationship is expected to be negative Precipitation of the warmest quarter P wq Wet areas are generally suitable to species, whereas reptile species are well adapted to arid conditions (Powney, Grenyer, Orme, Owens, & Meiri, 2010). Precipitation of the warmest quarter better captures water availability to populations, as represented by the seasonal bottleneck in which water can become a limiting factor. The relationship between precipitation and density is expected to be negative, except for reptiles Note. 'Retained' indicates whether the variable was eventually used or excluded because of multicollinearity.
instance, is assumed to allow for populations with higher densities (Currie & Fritz, 1993;Pettorelli et al., 2009). In contrast, at large geographical scales net primary productivity (NPP) is also related to species richness and thus potentially to resource partitioning (Cusens, Wright, McBride, & Gillman, 2012), with unclear effects on species population density (but see Novosolov et al., 2016). Climate can also influence population density. Temperature and precipitation can affect population density by determining resource fluctuations and extreme climatic conditions, which in turn can affect population density by regulating mortality rates (Novosolov et al., 2016;Williams S. E. & Middleton, 2008). As a consequence, we are not yet able to formulate clear hypotheses on many of the factors involved, but we must limit our research questions to broad expectations that deserve further exploration (Table 1).
Most published studies have focused on specific regions and individual taxonomic groups, using different methods or exploring different predictors (Buckley & Jetz, 2007;Carbone & Gittleman, 2002;Cotgreave & Harvey, 1994;Novosolov et al., 2016;Pettorelli et al., 2009;Silva et al., 1997). This has limited the comparability of studies and synthesis. How the ecological drivers of population density vary among groups characterized by different physiology and ecology is still to be explored. The relative influence of traits and environment on population density, for example, is likely to differ between ectotherms and endotherms. For example, energy expenditure increases with temperature in ectotherms but decreases in endotherms (Dillon, Wang, & Huey, 2010). A global and comprehensive comparative study is thus needed to disentangle the biological and environmental drivers of population density.
Indeed, obtaining a better understanding of the drivers of population density for a large range of organisms should lead to a gain in predictive ability and might ultimately allow the mapping of species population density over broad geographical scales. This would be an important step forward to investigate macroecological patterns (McGill & Collins, 2003), such as species abundance distributions (McGill et al., 2007;Xiao, O'Dwyer, & White, 2015), biomass distribution (Hatton et al., 2015) or range size-abundance relationships (Gaston et al., 2000). Understanding the spatial patterns of abundance would allow conservation biologists to identify regions within species distributions where viability is low (Di Marco et al., 2016;Hilbers et al., 2016) and contribute to forecasting the responses of species to climate change, assuming space for time substitution as done in species distribution modelling (Ashcroft et al., 2017). Finally, estimates of species abundance over large geographical scales and for multiple species would permit description of the change in a number of biodiversity measures (e.g., heterogeneity indices, functional or phylogenetic diversity metrics weighted by species relative abundance) over time and space . Although models quantifying the relationship between biological traits and environmental conditions with population densities are urgently needed, their predictive ability needs to be tested carefully before applying them in a conservation context.
Our objective here is 2-fold. First, we aim to understand the relative role of different biological and environmental drivers on the population density in the four terrestrial vertebrates (amphibians, reptiles, birds and mammals) and how they influence geographical population density patterns (Table 1). Second, we aim to evaluate the predictive ability of statistical models at predicting the spatio-temporal distribution of species population density. By quantifying the predictive accuracy of these models, we can shed light on the possibility of predicting abundance for conservation analyses, and assess for which applications density predictions can be considered an option.

| Collection of density estimates
We obtained population density estimates for terrestrial amphibians, reptiles, birds and mammals from the TetraDENSITY database (Santini, Isaac, & Ficetola, 2018). We did not consider density estimates for reptiles living in ponds or lakes (whose densities are strictly dependent on the water surface area and shore length) and for breeding aggregations of amphibians (which may last only a few days or weeks  (Figure 1; Supporting Information Data S1).

| Predictor variables
We considered seven variables as potential predictors of population density, including life history traits, resource availability, local environmental conditions and biotic interactions (Table 1). For life history traits, we focused on body mass and diet. Body mass is expected to show a negative relationship with population density because increasing masses require larger energy investments. We expected more carnivorous diets to be associated with lower population densities because of lower energy availability and accessibility. Body masses for amphibians and reptiles were estimated using length-mass allometric models to have comparable estimates for different morphotypes (Novosolov et al., 2016;Santini, Benítez-L opez, Ficetola, & Huijbregts, 2017 (Wilman et al., 2014). We classified as herbivores species feeding on 80% on plant material, as carnivores species feeding on 80% on animal material, and all the rest as omnivores.
Amphibians were all classified as carnivores.
We considered NPP (Imhoff et al., 2004) as a proxy of resource availability. Local environmental conditions were represented using the following climatic variables: annual mean temperature (T m ), precipitation of the warmest quarter (P wq ), precipitation and temperature seasonality (P cv and T sd , respectively) (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005). Annual mean temperature is expected to show a positive relationship because it is related to energy expenditure for endotherms and to energy availability for ectotherms. Seasonality in temperature and precipitation are expected to decrease population densities because they can increase mortality either directly or indirectly through resource fluctuations (Table 1) Figure S1 in Appendix S1); therefore, to avoid collinearity and maintain comparability among the four groups, we retained only NPP, P wq and P cv and excluded T sd , T m and Rich. However, because these excluded variables underlie different hypotheses, we also present the results for three alternative sets of models in the Supporting Information (Appendix S2): one including P cv , P wq and T m ; one including P cv , P wq and T sd ; and one including Rich, P cv and P wq .
To account for potential location errors in many study sites, all predictor variables were upscaled at 18 resolution. This resolution allowed us to uniform location uncertainty in the original studies, but inevitably masked the effect of species-specific habitat preferences and anthropogenic disturbance. Human pressures, however, are not distributed randomly and are assumed to influence species abundance greatly, so they have the potential to alter the global pattern observed for the environmental variables (Santini, Di Marco, et al., 2017). The scale of our analysis did not permit the testing of the precise effect of human pressures, but we showed that accounting for human influence did not alter the pattern observed for the other variables (see Supporting Information Appendix S3).
Previous studies have pointed out that macroecological investigation of density estimates may be confounded by the effect of sampling area (Blackburn & Gaston, 1996;Gaston, Blackburn, & Gregory, 1999), which scales inversely with population density estimates. This relationship is thought to arise for several reasons, including an increasing amount of unsuitable habitat in large areas, lower sampling efficiency, and a high number of different animals passing through in small areas.
As a consequence, some authors have included sampling area as a predictor (e.g., Novosolov et al., 2016;Pettorelli et al., 2009;Silva et al., 2001). We tested the extent to which our results were sensitive to sampling area and evaluated how the variance associated with sampling area is distributed across taxonomic levels.

| Statistical models
A schematic representation of the analytical framework is presented in Supporting Information Figure S2. Modelling population density for many species across the globe presents several challenges, because density estimates include spatial and temporal pseudo-replicates and can be phylogenetically autocorrelated. To account for all these factors, we used a linear mixed effect model for each vertebrate class. We retained all estimates available for a given species in a given site (replicates were not averaged) to consider the variability in population density within certain species and ecological contexts. As a consequence, one random effect ('Location', i.e., grid cell identity) was used to control for pseudo-replication within grid cells and to control for local effects influencing density other than climate and NPP that we could not otherwise control. One random effect was used to control for different sampling methods ('Method'), which are known to provide different estimates when used to estimate density for the same populations.
Finally, a nested random effect ('Order/Family/Species') was used to account for taxonomy. Genus level was not considered because of its high taxonomic instability. The random effect for the method was classified in broad categories, including 'Incomplete counts' (any incomplete count though plots/transects that is extrapolated to a larger area), for coefficient estimation (Pinheiro & Bates, 2000). As a first step, we compared the relative distribution of variance across the random effect levels (three nested taxonomic levels and the location) to assess whether density varied primarily between high or low taxonomic levels, and the relative effect of geographical space (McGill, 2008). For each taxon, we ran an intercept-only model plus the random effect structure described above (Supporting Information Figure S2a). The variance associated with each random effect was then expressed as a percentage of the total variance. As a second step, to estimate the influence and to tease apart the relative importance of the different retained predictors of population density, we ran a model selection for each taxon by testing all combinations of the variables considered (Supporting Information Figure S2b). After preliminary data exploration that revealed some nonlinear relationships, we also included a cubic term for body mass (also observed in previous studies; Silva & Downing, 1995) and a quadratic term for NPP, P wq and P cv . Models were compared using the Bayesian information criterion (BIC) and model posterior probabilities. The BIC tends to select more parsimonious models compared with the commonly used Akaike information criterion (AIC) and is thus more appropriate when sample size is large and AIC tends to select overfitted models (Raffalovich, Deane, Armstrong, & Tsao, 2008). Model posterior probabilities (equivalent to Akaike weights) indicate the weight of evidence for a given model within a set of competitive models (Burnham & Anderson, 2002). To calculate variable importance, we summed the posterior probabilities of the models in which each variable was included (Supporting Information Figure S2c).
We then obtained an average model by taking the average of models' coefficients weighted by the models' posterior probabilities, and assuming a posterior probability of zero for the models in which a given variable was not included (Supporting Information Figure S2d; Burnham & Anderson, 2002). Individual model performance was evaluated using the R 2 partitioned in its marginal (fixed effects alone) and conditional (fixed and random effects) components.
Here, we did not use a phylogenetic generalized least squares approach because this would not have allowed us to retain multiple density estimates for species while controlling pseudo-replication.
However, because our nested structure assumes that species are equally non-independent within Families, and Families equally nonindependent within Orders, we checked for the presence of phylogenetic autocorrelation in the residuals of the models by calculating Pagel's k (Supporting Information Figure S2e). Pagel's k is a measure of phylogenetic signal in models' residuals (Freckleton, Harvey, & Pagel, 2002) and ranges between zero (no phylogenetic signal) and one ( Previous studies indicated higher densities on islands owing to the so-called 'island compensation effect' (Adler & Levins, 1994;Buckley & Jetz, 2007;MacArthur, Diamond, & Karr, 1972;Novosolov et al., 2016). The definition of island from an ecological point of view is far from trivial (and possibly taxon specific), and the threshold chosen might affect coefficient estimations. Additionally, irrespective of the definition of island, our data were extremely imbalanced between the mainland and islands for birds and mammals (most data are from the mainland), and were insufficient on islands to be modelled for amphibians. To assess whether coefficient estimation was sensitive to islands and their definition, we provide an additional analysis in the Supporting Information ( Figure S3). We tested a model with all variables for reptiles, birds and mammals, including a factor classifying cells as mainland and islands. We repeated the analysis using different area thresholds of landmasses to categorize islands (corresponding to all landmasses where data were collected), and ensured that mainland and islands were balanced by randomly sampling an equal amount of mainland and island records. The whole procedure was replicated 100 times. We then related the coefficient estimates with the threshold area used.

| Predictions of models
To evaluate the predictive accuracy of the models, we ran a 5-fold cross-validation, repeated 10 times for each of the models (Supporting Information Figure S2f). Predictive performances were evaluated using the maximum absolute error (MAE). We also squared the correlation coefficient between the 'Predicted' and the 'Observed' data obtained from all repetitions of the 5-fold cross-validation to calculate a pseudo-R 2 . We compared the MAE with that obtained by randomly sampling species in the datasets 1,000 times and calculating MAE using different records of density estimates for the same species. This procedure allowed us to compare our prediction error with the expected natural intraspecific variability.
We then extracted the model-averaged parameter estimates to map population density globally. Based on the environmental conditions of each grid cell, we calculated the fitted population density for two hypothetical species in each of the four classes: one for a highdensity species and another for a low-density species, with 'high' and 'low' defined by both the fixed effect biological traits and the taxonomic random effects (Figure 4). This is not meant to represent the actual expected density because it does not account for species distribution, but it represents the average species response to environmental conditions.
We provide the R workspace, including the best predictive models, with this manuscript (Supporting Information Data S3). All analyses were conducted in R 3.0.3 (R Core Team, 2016) using several R packages. 'lme4' was used to run the linear mixed models, 'MuMIn' to calculate the marginal and conditional R 2 of the models, 'phytools', 'geiger' and 'ape' working with phylogenetic trees and calculating Pagel's k, 'raster' for all spatially explicit operations.

| RE S U LT S
The residuals of the models did not show any phylogenetic correlation (Supporting Information Table S3), indicating that the random effect specification removed the potential phylogenetic effect that could have influenced parameter estimation and increased type I error. Interceptonly mixed models found marked heterogeneity among the four classes in how variance in population density was distributed among taxonomic levels (Figure 2a). Density varied greatly among Orders in birds (23.6%) and mammals (48.5%) but not within amphibians and reptiles (c. 0%), where density mostly varied across Families (11.8 and 24.5%, respectively). Density varied across species in a similar manner between amphibians (12.1%), reptiles (8.7%), birds (10.8%) and mammals (7.3%). Location explained a large part of the variance in amphibians (56.3%), reptiles (39.4%) and birds (35.6%), which was similar to the total variance attributable to taxonomy (24.9-36.9%). Location was a smaller contributor to the variance in mammals (8.9%). Overall, the variance explained by the sampling method was low, ranging from c. 0% in amphibians to 9.9% in birds. The residual variance was similar between the four classes (12.7-20.4%; Figure 2a).
When we related population densities to life history traits and environmental predictors, most of our expectations were met, but some variables showed unexpected relationships (Figure 3; Supporting Information Table S1). Body mass showed a negative relationship with population density and was an important predictor of density in the four classes (Figures 2b and 3). In amphibians, the relationship showed a positive relationship for small body masses (0.3-1.5 g), whereas in mammals the relationship was curvilinear, with an asymptote in small body masses (c. 10-100 g) and one in large body masses (> 1,000 kg).
Diet was an important predictor of bird density, with carnivorous diets being associated with lower densities. Net primary productivity was important for all groups, especially for birds, and showed a positive linear relationship with amphibian density, negative with reptiles, and hump shaped in birds and mammals (Figures 2b and 3). The 5-fold cross-validation indicated good predictive abilities of the models (Supporting Information Table S2), with pseudo-R 2 between 0.67 (birds) and 0.84 (reptiles), MAE between 0.35 (birds) and 0.48 (reptiles), and without any consistent bias in the estimation (Figure 4).
The forecasting errors were within the range of errors calculated using the conspecific density records within the dataset (Supporting Information Figure S4). All models are available as part of the Supporting Information (Data S2).
Projected densities differed between the four groups of vertebrates, with amphibians displaying higher densities in tropical and highly productive areas, reptiles in tropical arid and non-seasonal environments, birds and mammals in temperate environments, and with the mammals more associated with wet areas when compared with birds ( Figure 5).
In islands, densities were on average higher than those collected on the mainland. However, the difference was sensitive to the area threshold used to define islands (Supporting Information Figure S3). Mammals showed lower densities on very small islands (< 10,000 km 2 ) than on the mainland, but the opposite was true for larger islands. Reptiles showed particularly higher densities on islands, with the highest difference with an island area threshold of c. 100,000 km 2 . The effect of some predictive variables was also sensitive to the area threshold (i.e., body mass in all groups; T sd in reptiles; P wq in mammals), but the sign of the coefficient did not change (Supporting Information Figure S3).
When controlling for variation in sampling area among studies, we found that the parameter estimates were quantitatively different but qualitatively unchanged (Supporting Information Figure S9), confirming that our conclusions are not biased by variation in sampling area among studies (Supporting Information Appendix S4). The results were mostly sensitive to the reduced sample size owing to the consideration of only data associated with sampling area estimates. Furthermore, the variance associated with sampling area is largely explained by the higher taxonomic levels (Orders and Families; Supporting Information Table   S5), suggesting that any artefact attributable to sampling area is limited in a cross-species model.

| Interpretation
Our results confirm that population density is strongly related to, and can be predicted by, a combination of species traits and environmental conditions (Currie & Fritz, 1993;Silva et al., 2001). The major contribution of this study is the finding that population density in the four major groups responds differently with respect to species traits and environmental conditions.
The variance of population densities distributed differently among taxonomic levels and locations in the four groups, with amphibians being particularly influenced by environmental conditions, and mammals being mostly influenced by traits. A number of sampling methods FIG URE 3 Partial response curves of the trait and environmental predictors. BM 5 body mass; NPP 5 net primary productivity; P cv 5 precipitation seasonality; P wq 5 precipitation of the warmest quarter exist to estimate population density in animal populations, and the choice of the method is influenced by the study species and habitat (Williams B. K., Nichols, & Conroy, 2002). Interestingly, although a number of studies demonstrated a relevant difference between the density estimates provided by diverse methods when compared on single species (e.g., Gottschalk & Huettmann, 2011;Seddon, Ismail, Shobrak, Ostrowski, & Magin, 2003), the sampling method explains only a minority of the variance on large-scale analysis across a variety of lifeforms.
Most of our expectations concerning environmental drivers were met. Among the set of traits we considered, body mass was by far the most important predictor for densities of reptiles, birds and mammals (Silva et al., 1997). As observed in previous studies (Silva & Downing, 1995), the relationship was nonlinear (almost asymptotic at extreme body mass values) in mammals, reflecting some constraints at the extremes of the body mass range. The smallest mammals may be limited in population density by the high energetic costs associated with endothermy. In contrast, Allee effects probably prevent large-bodied mammals from persisting at the very low population densities predicted by a linear relationship. Several interpretations have been proposed to explain the high densities of large mammals. These species might have lower rates of predation and might dominate in interspecific aggressions; besides, they might be more efficient in exploiting lower quality food sources, thus extracting more energy from the environment (Silva & Downing, 1995). Diet was a relevant predictor only in birds, whereas surprisingly, it was not important in mammals, for which previous studies identified diet as a major predictor (Silva et al., 1997). It is possible that, when taxonomic levels are treated as random effects, the effect of diet becomes negligible because diet and taxonomy are strictly linked in mammals.
Density showed a hump-shaped relationship with NPP in birds and mammals. This positive relationship is expected as long as resource availability is limiting, but when other factors may come into play, NPP might become unimportant. For example, NPP is positively correlated with species richness (Cusens et al., 2012); thus, although we observe an increase of biomass in high-productive environments, the density of individual species might actually decrease because of higher competition and predation rates, therefore reversing the relationship. The negative relationship of NPP in reptiles can be explained by the lower density of bird and mammal predators in low-productive environments (Janzen, 1976;Kreulen, 1979); it might also be attributable to the higher density associated with low productivity on islands (Novosolov  Figure S3).
Increased precipitation seasonality decreases population density, probably by contributing to high mortality rates in periods of water scarcity.
As expected, the relationship with precipitation of the warmest quarter was negative for reptiles and positive for amphibians and mammals, whereas surprisingly, it had no effect on birds. This might reflect the different adaptations of these taxa to arid conditions; for example, the excretory systems of reptiles and birds (producing uric acid) allow them to retain more water than the excretory system in mammals (producing urea). As for NPP, the negative relationship between reptile densities and precipitation is consistent with previous studies suggesting that the lower density of bird and mammal predators in arid environments may release reptiles from predatory pressures (Janzen, 1976;Kreulen, 1979). On the contrary, this relationship in amphibians was hump shaped, possibly reflecting the increased interspecific competition in high-precipitation areas where amphibian species richness is high (Buckley & Jetz, 2007).

| Predictions and applications
Our global predictions of potential densities depict different spatial patterns in the four taxa ( Figure 5). Overall, tropical and subtropical areas were characterized by a high density of reptiles and amphibians, whereas temperate and northern latitudes were characterized by a high density of birds and mammals. Birds and reptiles, however, showed higher densities in arid and semi-arid environments, whereas amphibians and mammals showed higher densities in high-productivity and humid environments, respectively. However, the interpretation of these predictions is not straightforward. These predictions indicate how population density is expected to respond to environmental gradients but ignore how species traits vary across environmental and geographical space. Traits such as body mass, for example, have been observed to vary along environmental gradients (e.g., Bergmann's rule) in both ectotherms and endotherms (Blackburn, Gaston, & Loder, 1999;Olalla-T arraga & Rodríguez, 2007 analyses where hundreds or thousands of species are jointly considered, such as in conservation planning, but also to derive useful biodiversity metrics that rely on relative densities (Santini, Belmaker et al., 2017) or to forecast population trends under future scenarios of environmental change (Ashcroft et al., 2017). At the same time, the way we built the models makes it possible to map information about the reliability of the predictions. Of course, a multispecies model cannot be particularly informative at the level of a single species, as it assumes the same scaling of specific environmental variables across all species. Nonetheless, data to develop species-specific models of population density are available for only a bunch of species (e.g., Lewis et al., 2017); therefore, at present, multispecies models can be the best surrogates for most species thus far.
As for most macroecological analyses, we were strongly limited by data availability, which led us to a final dataset clearly biased toward certain taxonomic groups and geographical regions (Figure 1). This probably influences our current state of knowledge on large-scale patterns in population density. Nonetheless, the lack of phylogenetic autocorrelation suggests that the control for the taxonomic levels in the random effect part of the models corrects for taxonomic bias. Geographical bias is relevant only when it corresponds to an environmental bias, which is clearly the case for amphibians and reptiles, and to a far lesser extent in birds and mammals ( Figure 5). Therefore, spatial predictions of population density should not be extended beyond interpolation areas.
Our results contribute to the general understanding of population density patterns of the four classes of terrestrial vertebrates and demonstrate a potential for predictions, which can certainly be improved as more data become available.

DATA ACCESSIBILITY
Population density estimates used in the paper are available from the TetraDENSITY database (Santini et al., 2018). The dataset used for the analyses is available as Supporting Information Data S1.