Mammal population densities at a global scale are higher in human-modified areas

We fitted a mixed to the in based on a as a of the index (HFI), and indirect while accounting for trophic level and primary productivity (normalized index; We found a significant positive relationship between population density and HFI, where population densities were higher in areas with a higher HFI (e.g. agricultural or suburban areas – no populations were located in very high HFI urban areas) compared to areas with a low HFI (e.g. wilderness areas). We also tested the effect of the individual components of the HFI and still found a consistent positive effect. The relationships remained positive even across populations of the same species, although variability among species was high. Our results indicate shifts in mammal population densities in human modified landscapes, which is due to the combined effect of species filtering, increased resources and a possible reduction in competition and predation. Our study provides further evidence that macroecological patterns are being altered by human activities, where some species will benefit from these activities, while others will be negatively impacted or even extirpated.


Introduction
Population density of terrestrial mammals span from < 0.0001 to > 25 000 individuals km −2 , and can vary over several orders of magnitude across conspecific populations (Santini et al. 2018a, Stephens et al. 2019) Understanding the drivers of this variation is a key objective of macroecology (Brown andMaurer 1989, Stephens et al. 2019). There has been a wide range of research examining the non-anthropogenic drivers of terrestrial mammal population density including the effects of species traits and environmental conditions (Table 1). Overall, there is a negative relationship between body mass and population density (i.e. larger species tend to have lower population densities) and this is a reflection of the variation in space use patterns and the relative energy use of species across the body size spectrum (Damuth 1981, Fa and Purvis 1997. In addition, population densities tend to be higher in temperate areas with intermediate Table 1. Hypotheses and predictions of anthropogenic drivers of population density. 'Human pressure' represents any of the anthropogenic variables.

Predictor
Expected overall relationship Rationale Reference Human footprint index (HFI) − Assuming a negative relationship between mammal population density and species richness, we expect that areas that are more impacted by humans will have fewer species, but higher average population densities. Berger 2007, Šálek et al. 2015, Venter et al. 2016, Santini et al. 2018b Night-time lights − Night-time lights represent areas with human infrastructure of varying density (low to high). We assume that higher values of night-time lights reflect areas with higher human activities. We expect a negative relationship between mammal population density and night-time lights. Davies et al. 2013, Gaston et al. 2013 Human population density +/− Human population density reflects human pressures within close proximity to human settlements (e.g. disturbance and hunting), but also infrastructure and anthropogenic resources. We might expect a negative relationship between human population density and mammal population densities due to the negative impacts of human disturbance. A positive relationship is also possible due to anthropogenic resources supporting higher mammal densities. Saari et al. 2016, Venter et al. 2016 Cropland +/− Croplands reflect land transformation, where cropland structure can vary depending on the management strategy (e.g. large-scale monoculture versus small-scale mosaic agriculture). We might expect a positive relationship between cropland and mammal population density due to croplands providing additional food resources. A negative relationship is also possible due to increasing croplands leading to fragmentation and habitat degradation or loss and therefore lower mammal densities. Schley et al. 2008, Bleier et al. 2012, Venter et al. 2016, Lewis et al. 2017 Pasture + Pastures also reflect land transformation, but pastures are more widely distributed than croplands. We expect a positive relationship between pastures and mammal population density due to pastures providing an additional food source for grazing herbivores (and higher densities of prey species for carnivores). Ramankutty et al. 2008, Venter et al. 2016 Accessibility − Accessibility is the ease with which humans can access cities and reflects infrastructure (e.g. roads and cities). We would expect a negative relationship between accessibility and mammal population density due to regions with low accessibility values being closer to cities and higher human pressures (e.g. increased habitat fragmentation and increasing access to nature). Saari et al. 2016, Venter et al. 2016, Tucker et al. 2018 Human pressure:mass − There may be an interaction between the anthropogenic variables and mass, where larger species are more susceptible to human impacts than smaller species and may have lower densities when human impacts are high. Cardillo et al. 2005, Benítez-López et al. 2017, Crees et al. 2019 Human pressure:diet +/− There might also be an interaction with diet, as studies suggest that densities of some carnivores increase along urbanisation gradients, while others are sensitive. Herbivore densities may increase with human impacts as some species may use humans as protection against predators. Omnivore densities may also increase due to anthropogenic food resources. Fedriani et al. 2001, Berger 2007, Pita et al. 2009, Craigie et al. 2010, Šálek et al. 2015, Saari et al. 2016, Venter et al. 2016, Santini et al. 2018b, Tucker et al. 2018 productivity levels, mediated by food availability and competition (Santini et al. 2018b). Carnivores, herbivores and omnivores forage on different food resources resulting in variable home range sizes and spacing patterns of individuals (Tucker et al. 2014, Tamburello et al. 2015, Stephens et al. 2019. As a consequence, herbivores generally occur in higher densities followed by omnivores and carnivores (Fa andPurvis 1997, Silva et al. 2001). These general density patterns across body mass, resources and diet were based on global datasets, however different trends (e.g. non-linear or reverse) have been found at finer spatial or taxonomic scales (e.g. regions or ecosystems; White et al. 2007, Pettorelli et al. 2009, Hempson et al. 2017, Pedersen et al. 2017. Therefore, terrestrial mammal population density patterns are not only highly variable across species, but also across spatial scales ranging from regional to global (Santini et al. 2018b). The expansion of human activities such as habitat fragmentation, land-use change and hunting have had important consequences for mammal abundance and distribution, and have altered global macroecological patterns. Previous studies have reported human impacts on species traits, ecosystem-level processes and ecological patterns across largespatial scales including global reductions in average body mass (Diniz-Filho et al. 2009, Rapacciuolo et al. 2017, Santini et al. 2017, declining species richness (Torres-Romero and Olalla-Tárraga 2015) and homogenisation of the zoogeographic regions (Bernardo-Madrid et al. 2019). In addition, human activities have reduced mammalian movements globally (Tucker et al. 2018). It has been proposed that this reduction in movement may not only be a consequence of movement barriers and habitat fragmentation, but also of higher resource availability in human-modified areas (Tucker et al. 2018). Higher resource availability in human-modified areas may come from various sources such as the presence of crops, supplemental feeding and artificial water sources (Prange et al. 2004, Jones et al. 2014. Species ranging patterns and population density are often linked (Makarieva et al. 2005), and we may expect that the proposed mechanisms of higher resource availability may also apply to population density. For example, the presence of abundant resources may not only reduce the need to move, but may also support higher population densities. In contrast, we would expect that movement barriers, such as roads or fences, would limit access to resources and lead to a decline in population density.
However, so far evidence that human activities have had consequences for mammal population densities consists of a mixture of studies that focus on a single species or single system, or has focussed on species abundances across large scales. The results of these studies are mixed, with localised studies identifying both negative and positive effects on population density (Šálek et al. 2015(Šálek et al. , Said et al. 2016, while large-scale studies exploring abundance trends has identified general patterns of declines in response to humans (Estes et al. 2011, Newbold et al. 2014, Purvis et al. 2018. Negative effects leading to declining densities include habitat degradation, direct disturbance and exploitation (Benítez-López et al. 2010, Newbold et al. 2015. Positive effects leading to increasing densities may be due to various processes including increased resource availability (e.g. food/water supplementation), humans acting as predator deterrents (i.e. known as 'human shield' effect: Berger 2007, Wang et al. 2015, Steyaert et al. 2016, or extinction of predators and/or competitors (Crooks and Soulé 1999). For example, small carnivores such as the red fox Vulpes vulpes and raccoon Procyon lotor can have higher population densities along urbanisation gradients (Šálek et al. 2015), which may result from the absence of large carnivores (mesopredator release hypothesis; Crooks and Soulé 1999). Conversely, some species such as wildebeest Connochaetes taurinus and impala Aepyceros melampus tend to have lower population densities in human-modified landscapes (Said et al. 2016). This differential response of species to humans is most likely due to differences in behaviour or diet plasticity, and life history traits that enable some species to adapt to human-modified landscapes (Šálek et al. 2015(Šálek et al. , Santini et al. 2019). There may also be an effect of species richness, where a negative relationship between population density and species richness could result in human-modified areas having fewer species, but higher average population densities (Santini et al. 2018b).
Previous studies have shown that human impacts tend to have a greater impact on species with larger body size (Cardillo et al. 2005, Diniz-Filho et al. 2009). Larger species tend to have slower rates of population growth (Fenchel 1974), have lower population densities (Damuth 1981) and are often persecuted (e.g. hunting; Benítez-López et al. 2017), which increases their susceptibility to human impacts. However, studies have also shown that small and large mammals may be vulnerable to different threats, with small mammals being more sensitive to habitat degradation and loss, while large mammals are more sensitive to over-exploitation and persecution (González-Suárez et al. 2013. There is evidence that carnivores, omnivores and herbivores have mixed responses to humans. Studies suggest that some carnivore densities increase along urbanisation gradients (Šálek et al. 2015), although other carnivores are more sensitive (Pita et al. 2009). Similarly, some herbivore densities may increase with human impacts by living close to humans as protection against predators (Berger 2007). There is also evidence that some omnivorous species are attracted to areas where humans are present due to the accessibility of anthropogenic food resources that may support higher densities (Fedriani et al. 2001). However, declines in densities across the diet classes have also been described (Craigie et al. 2010), suggesting that species responses are dependent upon whether they are able to adapt to human modified areas (Santini et al. 2019).
With these conflicting results, it remains unclear how the results of these small-scale phenomena translate into global patterns of population density in response to human impacts. Here we examined whether there is -on a global scale -an overall positive or negative relationship between mammal population density and human impacts. We related 6729 population density estimates for 468 species with a number of intrinsic and extrinsic drivers explored in previous studies (Table 1; Santini et al. 2018b) and with the Human footprint index (HFI) (Venter et al. 2016), a composite measure of human impacts on ecosystems. To disentangle the various sources of human impact on mammal population density, we repeated the analysis on the specific components of the HFI including human population density, pastures, croplands, night-time lights and accessibility (Table 1). For species with sufficient sample sizes we also examined the intra-specific relationships and summarized their effects using a metaanalytical technique. If movement barriers are the key factor influencing population density in high HFI areas, we would expect a decline in population density overall (Tucker et al. 2018). However, if food resources are a key factor shaping population density, then we would expect an increase in population density (Santini et al. 2018b, Tucker et al. 2018). However, we expect different species to have mixed responses to some of the anthropogenic pressure variables due to differences in diet and body size (Table 1).

Population density data
We used terrestrial mammal population density estimates from an unpublished extended version of the TetraDENSITY database (Santini et al. 2018a). TetraDENSITY includes spatially explicit terrestrial vertebrate population density estimates (number of individuals km −2 ) from around the globe. The density estimates were extracted from the literature, therefore, the sampling methods and timing of data collection varied according to the study. Population density was estimated using six broad sampling methods: incomplete counts (any incomplete count that is extrapolated to a larger area), censuses ('complete' counts, which assume full detection of individuals), distance sampling (including different algorithms and sampling designs), home range extrapolation (derived from home range area estimation), mark-recapture (including different algorithms and capture approaches) and trapping (removal methods, indicate the minimum number known to be alive) (Santini et al. 2018a). From this database, we extracted population density, sampling method, locality (study site), continent and year that the estimate was collected in. Due to the limited temporal span of the environmental and anthropogenic correlates, we only used population density estimates collected between 1990 and 2019. We also excluded estimates that referred to non-native species, and for which spatial coordinates were highly uncertain (> 1 degree). After filtering, the extended dataset included 6729 estimates for 468 species (17 orders and 68 families) globally, excluding Antarctica ( Fig. 1).

Human impact, primary productivity and life history data annotation
As a measure of human impact, we used the Human footprint index (HFI), which includes a range of human disturbances including roads, croplands, night-time lights and built environments (Venter et al. 2016). HFI has been calculated for two time periods (1993 and 2009) at a 1 km resolution and values range from 0 (low footprint; e.g. Greenland) to 50 (high footprint; e.g. London) (Venter et al. 2016). The HFI values have previously been assigned to five human pressure categories including 'no pressure' (HFI = 0), 'low pressure' (HFI = 1-2), 'moderate pressure' (HFI = 3-5), 'high pressure' (HFI = 6-11) and 'very high pressure' (HFI = 12-50) (Venter et al. 2016). Some of the longitude and latitude positions had unknown spatial uncertainty as the positions were estimated based on names of locations provided in the Figure 1. Spatial distribution of the density estimates in our sample. Point transparency is used to aid with the visualisation of overlapping study locations (i.e. darker points indicate a higher number of estimates for the same location). study (Santini et al. 2018a). This might result in a spatial mismatch, where the longitude and latitude assigned to a population might place the population in a low HFI cell, but the population were actually located in a neighbouring cell with a different HFI value. Furthermore, population density estimates in large mammals are often performed over large areas (tens of km 2 ; e.g. spatially explicit capture recapture methods or aerial censuses). As a consequence, we resampled the 1-km resolution HFI data to a 1-degree resolution by averaging the values across the cells, and this was performed prior the HFI data annotation process. We then annotated the data with HFI using the date of collection, longitude and latitude of each estimate. If a population density estimate was collected between 1990 and 2000, we assigned a HFI value from 1993, and density estimates from 2001 to 2019 were assigned a HFI value from 2009.
As the HFI is a composite measure of different human pressures that may have conflicting impacts of mammals, we collected additional anthropogenic variables including human population density, accessibility, night-time lights, croplands and pastures. With these additional data, we tested specific hypotheses (Table 1) about how humans impact mammal population densities. All anthropogenic data were first resampled the data to a 1-degree resolution by averaging the values, prior to the annotation process. Human population density, croplands and pastures data were extracted from the History database of the global environment (HYDE, ver. The cropland and pasture measures represent the percentage of the cell covered by croplands and pastures, respectively. We calculated the average human population density, percentage of croplands and percentage of pasture using the data from all years (i.e. [2000][2001][2002][2003][2004][2005][2006][2007][2008][2009][2010][2011][2012][2013][2014][2015] and annotated each density estimate with this averaged value based on the longitude and latitude of the population. Accessibility data for 2000 were extracted from Nelson (2008) and for 2015 from Weiss et al. (2018). Population density estimates collected between 1990 and 2007 were assigned an accessibility value from 2000, and density estimates collected from 2008 and 2017 were assigned an accessibility value from 2015. Nighttime light data (digital number; DN) spanning 1992-2013 were obtained from the National Geophysical Data Center of National Oceanic and Atmospheric Administration (NOAA/ NGDC) (<https://ngdc.noaa.gov/eog/dmsp/downloadV-4composites.html>, accessed 15/10/2019). We averaged the night-time light 1-degree data across all years and annotated the density data based on longitude and latitude position.
We used an average of the normalized vegetation index (NDVI) from 2000 to 2019 to account for the effects of mean productivity at each sampling location (<https:// urs.earthdata.nasa.gov>, accessed 12/09/2019). Similar to the HFI, the NDVI data were first resampled to a 1-degree resolution by averaging the values to account for spatial uncertainty of the population locations. We also annotated the density estimates with mammal species richness data (1-degree resolution; IUCN 2017) to control for the negative relationship between population density and species richness (Santini et al. 2018b). We annotated each population density estimate with the respective species-level body mass and diet category from the PanTHERIA (Jones et al. 2009) and EltonTraits (Wilman et al. 2014) databases. We used specieslevel body mass values because the population-level values were not available.

Human footprint index
Using all 6729 density estimates, we ran a linear mixed effects model with log 10 population density as the response variable and body mass, HFI, NDVI (quadratic term), trophic level and two interaction terms (HFI:mass and HFI:diet) as predictor variables. We included an interaction term between HFI and body mass to test for any changes in the response of different sized species to human-modified landscapes. The interaction between HFI and diet enabled us to examine whether carnivores respond differently to human impacts compared with omnivores and herbivores ( Table 1). The model also included random effects for taxonomy (order/family/species), sampling location, sampling method, sampling year and continent to account for various inherent methodological structures and spatial bias in the data. It was important to include sampling method as a random effect because the different methods tend to have a different species coverage due to sampling bias, and the random effect accounts for this heterogeneity by adjusting the intercept. We did not include genus in our models due to its high taxonomic instability (Santini et al. 2018b). In addition to the nested taxonomy random effect included in our models, we checked for the model residuals for the presence of phylogenetic autocorrelation using Pagel's λ, where 0 indicates no phylogenetic signal and 1 indicates a strong phylogenetic signal (i.e. distribution expected under the Brownian motion model of evolution) (Freckleton et al. 2002) using the mammalian supertree (Fritz et al. 2009).
We checked for multicollinearity among the variables included in the models using correlation coefficients and variance inflation factors (VIFs). The correlation coefficients among the predictor variables within each model were |r| ≤ 0.30, which is below the common cut-off value of 0.7 (Dormann et al. 2013), and the VIFs were below the commonly accepted cut-off value of 4.0 (Zuur et al. 2010). Mammal population density and body mass values were log 10 -transformed to normalise the model's residual distribution and reduce heteroscedasticity.
The full dataset was skewed towards populations that occur in HFI areas between 0 and 15 (~84%), so we also resampled the data so that the same number of data points were selected across the range of HFI. HFI was split into five human pressure categories: 'no pressure', 'low pressure', 'moderate pressure', 'high pressure' and 'very high pressure', with each category representing ~20% of the Earth's surface (Venter et al. 2016). The minimum sample size across 6 these categories was 34 in the 'very high pressure' category (HFI = 12-50). We repeated this sampling 1000 times and then examined the proportion of the time that a positive relationship between density and HFI occurred. We also ran the HFI model using the same data restricted to HFI values between 0 and 15 to ensure that our results were not driven by those points > HFI 15.

Other anthropogenic variables
We ran additional models for the other anthropogenic variables (i.e. human population density, accessibility, night-time lights, etc.). Each anthropogenic variable was examined separately (except for croplands and pastures) due to significant correlations (Supplementary material Appendix 1 Table A1, Fig. A1). In each of these models, we included the same predictor variables as the human footprint model: body mass, HFI, NDVI, trophic level and two interaction terms (human pressure:mass and human pressure:diet). Mammal population density, body mass, night-time lights, accessibility and human population density values were log 10 -transformed to normalise the model's residual distribution and reduce heteroscedasticity. We note that the quadratic NDVI term was not significant in any of the models and the models including the quadratic term had VIF values larger than 4. To reduce multicollinearity in the models, we excluded the quadratic NDVI term from all of the final models.

Sensitivity to spatial resolution
As the 1-degree resolution data has a large grid size (~110 × 110 km) thus increasing the risk of averaging the environmental values experienced by the populations with others that are actively avoided, we ran sensitivity analyses to ensure that the data resolution did not influence our results. We averaged the HFI, NDVI and anthropogenic variables to 10 and 50 km. We then annotated the mammal density data with the 1, 10 and 50 km data based on the longitude and latitude positions. In addition, we also annotated the 10 km data using a buffer around the coordinates, where the buffer size was equal to the area needed to host 100 individuals according to the local density estimate (i.e. density × 100). This enabled us to account for range of environmental conditions experienced by a population proportional to its density, and the varying size of sampling areas that are inversely proportional to the population density. We then ran the same linear models with mammal density as the response variable and body mass, HFI, NDVI, trophic level and two interaction terms (human pressure:mass and human pressure:diet) as the predictor variables. We also included random effects for taxonomy (order/family/species), sampling location, sampling method, sampling year and continent.

Cross-validation
We assessed the robustness of our inferences using spatial block cross-validation (Roberts et al. 2017, Currie 2019. This enabled us to test a) if our models could predict population density patterns using independent data from different spatial locations across the globe, and b) if the direction of the relationship between mammal population density and human impacts was consistent across space. We ran cross-validation with 1-degree spatial blocks and a systematic sample fold for each of the five anthropogenic variables (i.e. HFI, accessibility, cropland, pasture, night-time lights and human population density). For each of the anthropogenic variables, we ran the complete model including both interactions, a model with only the mass interaction and a model with only the diet interaction. We also ran a null model without any anthropogenic variables to test whether the models with the anthropogenic variables had better predictive capabilities than the null model. We evaluated the performance of the cross-validation using the mean absolute error (MAE), selecting the best predictive model based on the lowest MAE (i.e. averaged across the spatial blocks per model). We also calculated a pseudo-R 2 by squaring the correlation coefficient between the 'Predicted' and the 'Observed' data obtained from the cross-validation (including all spatial blocks).

Taxonomic diversity across human impact gradients
To assess potential biases in the data and/or filtering effects on species composition along the human impact gradients, we examined the taxonomic distinctness of density estimates (Clarke and Warwick 1998) along gradients of each of the anthropogenic variables. Each anthropogenic variable was split into 10 quantile bins (i.e. 0-1, by 0.1 increments) and the taxonomic distinctness was calculated for each quantile bin. If species filtering is occurring, then we would expect a negative relationship between taxonomic distinctness and the anthropogenic variables.

Intraspecific patterns
We also examined intraspecific patterns using a meta-analysis approach to quantify the direction, strength and consistency of the association between density and the anthropogenic variables across conspecific populations. First, we estimated the correlation coefficient between density and each anthropogenic variable for all species with a n ≥ 5 using a Spearman rank correlation. We limited the correlations per species to estimates sampled using the same methodology. A minimum sample size of 5 is required to estimate the variance of the correlation. This resulted in a sample size of 85 correlation coefficients for 80 species spanning 12 orders (47.0% Cetartiodactyla, 22.35% Primates, 15.3% Carnivore, 4.7% Rodentia and 10.6% others). Then we transformed the correlation coefficients to Fisher's z scores using the correlation sample size to obtain the effect size and variance of each correlation. To estimate a summary effect size across species, we ran a mixed-effect meta-analysis on the z scores and their variance (Borenstein et al. 2009). We included a nested random effect for taxonomic orders, families and species. We tested the residual heterogeneity using the Q-statistic. A significant Q test indicates that a significant amount of variability exists between the effect sizes. Then, we ran the same meta-analyses using diet and body mass as moderators (meta-regression).
This allowed us to test whether the variance across effect sizes could be explained by these two traits. We then ran the Omnibus test to assess if the residual heterogeneity of the meta-analyses is significantly explained the by the body mass and diet of the species.
All analyses were carried out in R ver. 3.5.1 (R Core Team) and details on the R packages used in the analyses can be found in the Supplementary material Appendix 1 Table A2.

Human footprint index
The final dataset included body mass values from 7 g to 3940 kg, HFI values from 0 to 33, species richness values from 1 to 193 species and NDVI values from 0 to 0.83. Areas with a HFI value of ~33 represent regions of very high human pressure that include a mix of agricultural landscapes and cities such as the area surrounding Cambridge (UK) with an average human population density of 776 individuals km -2 .
Based on this dataset, we found a significant positive relationship between terrestrial mammal population density and the HFI (Fig. 2). On average, population densities in areas of high HFI were up to 25 times higher than areas of low HFI (HFI = 0: 1.02 individuals km -2 , HFI = 15: 4.53 individuals km -2 , HFI = 33: 26.29 individuals km -2 , based on a 15 kg herbivorous mammal). We did not detect any significant effects of primary productivity (NDVI; p = 0.938) or the interaction between HFI and diet for carnivores (HFI:Diet (Carnivore); p = 0.543), but there was a significant interaction between HFI and diet for omnivores (HFI:Diet (Omnivore); p ≤ 0.001). There was a significant effect of diet, where omnivores have significantly lower population densities than herbivores (p = 0.011). The interaction term between HFI and body mass was also significant (p = 0.01), where there was a positive relationship between population density and HFI across different body masses, but the intercept is lower for large species (Supplementary material Appendix 1 Fig. A2). As predicted, we found a significant negative relationship between mammal population density and body mass ( Table 2).
The models using resampled data across the human pressure bins so that the same number of data points were selected across the range of HFI showed a positive relationship for 93.3% of all models, further suggesting the positive relationship between population density and HFI is not an artefact of uneven sampling of density estimates across the distribution of the HFI (Supplementary material Appendix 1 Fig. A3). After running the HFI model focusing on the data between HFI 0-15, we still found a significant relationship between density and HFI, providing further support for the increasing population density relationship (Supplementary material Appendix 1 Fig. A4, Table A3).

Other anthropogenic variables
Consistent with the HFI results, we found a significant positive relationship between terrestrial mammal population density and night-time lights (Supplementary material Appendix 1 Fig. A5 Figure 2. Partial response plot for the relationship between mammal population density (individuals km -2 ) and the Human footprint index (n = 6729) with 95% confidence intervals (shaded area). This is a visualisation of the model results presented in Table 2. significant negative relationship between population density and body mass across the models (Supplementary material  Appendix 1 Table A4-A7). We did not detect any significant effects of primary productivity suggesting that anthropogenic resources might play a role mammal density patterns, and the results were mixed for both diet, and the interaction term between the human pressure and diet (Supplementary material Appendix 1 Table A4-A7). We consistently found a significant interaction between human pressure and omnivory, suggesting that the slope of the relationship between population density and human pressure was shallower than the slope for herbivores and carnivores. The interaction term between the anthropogenic variables and body mass were all significant, with a similar effect as HFI (i.e. positive relationships for all body masses, but a higher intercept for smaller species (Supplementary material Appendix 1 Table A4-A7)). For all of the models, we did not find a phylogenetic signal in the residuals, with λ < 0.001 (Supplementary material Appendix  1 Table A8).

Sensitivity to spatial resolution
The results from models including the 1, 10 and 50 km data were qualitatively similar to the 1-degree results, with a consistently positive relationship between terrestrial mammal population density and HFI (Supplementary material Appendix 1 Fig. A6a) The relationship between population density and accessibility was also consistently negative (Supplementary material Appendix 1 Fig. A6f ).

Cross-validation
The cross-validation results showed that the pseudo-R 2 of the models ranged between 0.506 and 0.564 and MAE between 0.663 and 0.787 (Supplementary material Appendix 1 Table  A9). The best predictive model included mass, HFI, NDVI, species richness, diet and an interaction term between diet and HFI, with a MAE value of 0.663 and a pseudo-R 2 of 0.559 (Supplementary material Appendix 1 Table A9). This suggests that the inclusion of an anthropogenic variable, in this case HFI, improved the model predictability compared to the null model. We did not detect any consistent bias in the population density estimation (Supplementary material Appendix 1 Fig. A7).

Intraspecific patterns
The mixed-effects meta-regression results and Spearman's correlation coefficients suggested qualitatively similar trends to the mixed effects models, with a general trend of increasing densities with increasing human pressure (Supplementary  material Appendix 1 Table A10, Fig. A8). However, there was a significant amount of variation across the populations based on the significant Q-test results (Supplementary material Appendix 1 Table A10). The Spearman's correlation coefficients further supported this variation across populations, where 60% of species had a correlation coefficient > 0 for HFI and 54-64% > 0 for the remaining anthropogenic variables (Supplementary material Appendix 1 Fig. A7). When including body mass and diet as moderators in the metaregression, the results again suggested similar patterns to the main analyses, however the significance of these relationships and the amount of variance explained by diet and mass were mixed (Supplementary material Appendix 1 Table A11). As predicted, there were differences in the response of different diet categories to human impacts, with carnivores being the only group that had significant results for the meta-regression analyses, and also had consistently higher intercepts (Supplementary material Appendix 1 Table A11). These results suggest that carnivore population densities consistently increase with human impacts, excluding accessibility, which was predicted to have a negative effect on density.

Taxonomic diversity across human impact gradients
The taxonomic distinctness results were mixed, with a general trend toward high distinctness at intermediate levels of human pressure, however distinctness decreased with increasing cropland and night-time lights (Supplementary material Appendix 1 Fig. A9). These results suggest that species composition in our sample changes along human pressure gradients, with a general trend towards hump-shaped relationships, however there was some variation in the diversity  A9). There was a negative relationship between diversity and night-time lights, and a positive relationship with percentage of pasture (Supplementary material Appendix 1 Fig. A9).

Discussion
Overall, we found a significant positive relationship between terrestrial mammal population densities and the human footprint based on the density estimates of 468 species spanning HFI values 0-33 (no human pressure to very high human pressure). When examining the individual components of the HFI, we were not able to completely tease apart the mechanisms underlying the density-HFI relationship. Our results consistently show increasing mammal densities with each of the anthropogenic variables, however, we did not find a significant relationship with productivity, indicating that anthropogenic resources may be a contributing factor. Combined, our results suggest that mammal population densities are on average higher in areas with higher human disturbance.
Our results may seem contrary to what one might expect, particularly given that human-modified landscapes are seen as undesirable and can have negative consequences for many species (Seto et al. 2011). However, with landscapes changing globally and the resulting decreases in global biodiversity (e.g. the Living planet index (McRae et al. 2017)), the underlying trends shaping these patterns are often mixed depending on the scale at which the trends are examined (i.e. populationlevel versus assemblage-level metrics) (Dornelas et al. 2019). In many cases, there are both species that can adapt and thrive in modified landscapes ('winners') and species that are unable to adapt and disappear ('losers'). This process is also known as species filtering, where a relatively small number of species can adapt and thrive in human-modified landscapes while most cannot (Brashares 2010, Riggio et al. 2018, Dornelas et al. 2019, Santini et al. 2019. It may also be that the overall positive relationship is due to the decline or loss of rare sensitive species from human-modified landscapes, while common species increase in density (Davies et al. 2004). This differential response between common and rare species has been demonstrated for birds and mammals (Lennon et al. 2004, Vázquez andGaston 2004).
Previous studies have shown that in human-modified landscapes those species able to adapt and thrive show signs of increasing population densities due to a reduction in competition for resources (Peres and Dolman 2000, Ruscoe et al. 2011, Wolf and Ripple 2017. This leads to the increasing density of successful human-adapted species, but an overall decrease in species diversity as non-adapted species disappear (Parsons et al. 2018). In the case of our results, it is possible that while we see an overall positive trend for mammal population density, there are changes in the species composition underlying this trend, with some species increasing and others decreasing (Dornelas et al. 2019).
There are potential benefits associated with these humanmodified landscapes, which span a range of habitats such as urban and agricultural areas. For example, human-modified landscapes may provide supplemental food and water resources (e.g. anthropogenic food such as garbage or food supplementation), shelter and altered climate (e.g. light and temperature) (Bateman andFleming 2012, Newsome et al. 2015) that may support higher population densities. Studies have shown that mammals may use human-modified landscapes as a way of minimising predation risk, known as the 'human shield effect' (Berger 2007). Additionally, agricultural regions may include benefits such as consistent and easily digestible food resources (McLennan and Hockings 2014).
Our taxonomic distinctness analysis results suggest that changes in species composition or species filtering may be occurring, with high species richness generally occurring at intermediate levels of human pressure. This could be a reflection of the intermediate disturbance hypothesis, where the diversity of competing species is expected to be high at intermediate intensities of disturbance (Connell 1978). Interestingly, there was a negative relationship between taxonomic distinctness and cropland, suggesting that croplands may support fewer species. Changes in taxonomic distinctness could also be due to sampling bias, where population density studies are predominantly conducted in areas that are at intermediate distances from human settlements. The remaining population density studies are conducted on human-adapted species near human settlements, or on rare and sensitive species in more remote areas, but these are not as common. As we still see increasing densities with HFI within species' (i.e. the meta-analysis models), it is possible that the species for which we have density estimates across the HFI continuum are already filtered (i.e. species that are already lost may not be included in our analysis). However, as our cross-validation results confirm that regardless of the mechanism (e.g. competition or predation release, resource availability or human-shield effect), areas with higher HFI have higher densities of terrestrial mammals compared to areas with low HFI.
Increasing mammal population density may also be a response to habitat or landscape simplification, where the reduction in complexity provides a reduced number of niches for a few species, which then become highly abundant (Coleman andBarclay 2012, Hohnen et al. 2016). For example, deer mice Peromyscus maniculatus were found to be more abundant and had higher rates of reproduction in burnt (i.e. more simplistic) habitats compared to unburnt (i.e. more complex) habitats (Zwolak et al. 2012). Our results provide some support that the simplification of landscapes due to human modification is altering the abundance and density of mammals, where some species are becoming more abundant and others are disappearing.
Our results provide strong support for one of the hypotheses raised by Tucker et al. (2018), who suggested that increased resource availability in human-modified landscapes as a potential mechanisms underlying mammal movement reductions. Species that are able to adapt and thrive in human-modified landscapes may have higher population densities due to the increased availability of resources and reduced predation or competition, and in turn the reduced need to move extensively (see also, Stephens et al. 2019). Barriers (e.g. roads) and habitat fragmentation on the other hand, may impact human-sensitive species, leading to decreasing densities and local extinctions (Buchmann et al. 2013, Chase et al. 2020), our results show higher mammal densities in areas with human pressure. Taken together, higher resource availability in high HFI areas, competition/ predation release mechanisms or land simplification, may underlie our finding that animal population densities are higher in human-modified areas.
Our results also provide additional evidence that carnivores may benefit from human-modified landscapes due to their higher population densities compared with herbivores and omnivores. Šálek et al. (2015) also found that carnivorous mammal population densities were higher and home range sizes were smaller in human-modified landscapes. They attributed these patterns to the behavioural flexibility and life history adaptions of carnivores, in addition to anthropogenic resources, favourable abiotic conditions and reduced predation that enable carnivores to adapt to human-modified environments.
Our results differ to Saari et al. (2016) who examined the abundance of species based on a meta-analysis and found that the overall abundance of terrestrial animals (i.e. birds, mammals and arthropods) was lower in urban areas. We note that our results do not reflect those populations living in highdensity urban areas (i.e. HFI > 40). The populations with the highest HFI value in our study were in areas of high human footprint with various configurations such as patches of natural vegetation (e.g. forests) interspersed with urban areas (e.g. roads and housing). Other areas were on the outskirts of cities interspersed with intensive agricultural landscapes. This could be one reason that our results differ from the findings of Saari et al. (2016) and there is a clear need for further work focusing on compiling population density data for species in urban areas (i.e. > 40 HFI) to examine effects of humans across the entire spectrum of HFI.
There may be various consequences associated with the trend of increasing population densities with higher human footprint. For example, higher population density can lead to changes in disease dynamics (Pongsiri et al. 2009), increased human-wildlife interactions (Soulsbury and White 2016, König et al. 2020, Pozo et al. 2020, habitat changes (Wallach et al. 2010) and declines of other species (e.g. native species; Hollings et al. 2016). Future work should explore the relationship between modified landscapes and population density to tease apart what these consequences are for populations and ecosystems.
The positive relationship between population density and HFI we reveal here was robust against various potential biases in the data that we accounted for (e.g. sampling). However, there are other potential correlates of HFI that may positively influence mammal population densities that we could not account for here, such as difficult to measure covariates like freshwater availability. It is also possible that the effect of humans on mammal population density may relate to a potential overlap between the environmental preferences of humans and mammals, where humans and mammals prefer to colonise areas with favourable characteristics, such as adequate resource availability and mild climatic conditions, could provide additional explanation for the positive relationship between density and the HFI. With the inclusion of NDVI in our models, we should have accounted for at least some of the impacts of resource availability.
We note that the validity of the patterns described here is scale-dependent. Our results reflect the relationship between human impact and mammal population density globally. Our main results and sensitivity analysis provide evidence that population densities tend to be higher in human-dominated regions, however, the relationship between human impact and mammal population density certainly could differ at a local scale.
Our findings -based on a rather coarse 1-degree resolution and not including population estimates in very high HFI areas such as cities -suggest that terrestrial mammal populations are higher in areas modified by humans and future research should disentangle the mechanisms shaping mammal population density and examine what the consequences of these patterns are. For example, it would be important to perform longitudinal studies that track abundances of mammals under varying human pressures (also see Barnes et al. 2016, Kiffner et al. 2020. It would be also important to determine what the effects of increasing abundance for community composition (e.g. species richness) and ecological functions (e.g. predation) are. Furthermore, if there is a species filtering effect due to humans, then it would be useful to understand which traits enable species to tolerate humans. As human populations continue to expand, understanding biodiversity in human-modified landscapes becomes more important. Management strategies that foster human-wildlife coexistence will be critical for species survival and human wellbeing.
Author contributions -MAT, LS, CC and TM conceived the manuscript. MAT and LS conducted the analyses. MAT and TM wrote the first manuscript draft. All authors contributed to interpretation of results and writing the final version of the manuscript.