Scholarly article on topic 'Wildlife population trends in protected areas predicted by national socio-economic metrics and body size'

Wildlife population trends in protected areas predicted by national socio-economic metrics and body size Academic research paper on "Biological sciences"

Share paper
Academic journal
Nature Communications
OECD Field of science

Academic research paper on topic "Wildlife population trends in protected areas predicted by national socio-economic metrics and body size"



Received 28 Aug 2015 | Accepted 29 Jul 2016 | Published 1 Sep 2016^^HDOI:iai038/ncomms12747l OPEN

Wildlife population trends in protected areas predicted by national socio-economic metrics and body size

Megan D. Barnes1,2'*, Ian D. Craigie3'*, Luke B. Harrison4, Jonas Geldmann5,6, Ben Collen7, Sarah Whitmee8, Andrew Balmford6, Neil D. Burgess5,9, Thomas Brooks10,11,12, Marc Hockings1,9,10 & Stephen Woodley13

Ensuring that protected areas (PAs) maintain the biodiversity within their boundaries is fundamental in achieving global conservation goals. Despite this objective, wildlife abundance changes in PAs are patchily documented and poorly understood. Here, we use linear mixed effect models to explore correlates of population change in 1,902 populations of birds and mammals from 447 PAs globally. On an average, we find PAs are maintaining populations of monitored birds and mammals within their boundaries. Wildlife population trends are more positive in PAs located in countries with higher development scores, and for larger-bodied species. These results suggest that active management can consistently overcome disadvantages of lower reproductive rates and more severe threats experienced by larger species of birds and mammals. The link between wildlife trends and national development shows that the social and economic conditions supporting PAs are critical for the successful maintenance of their wildlife populations.

1 School of Geography Planning and Environmental Management, the University of Queensland, St Lucia, Queensland 4067, Australia. 2 Australian Research Council Centre of Excellence for Environmental Decisions, the University of Queensland, St Lucia, Queensland 4072, Australia. 3 Australian Research Council Centre of Excellence for Coral Reef Studies, James Cook University, Townsville, Queensland 4811, Australia. 4Redpath Museum, McGill University, 859 Sherbrooke Street West, Montreal, Quebec H3A 0C4, Canada. 5 Center for Macroecology, Evolution and Climate, Natural History Museum of Denmark, University of Copenhagen, Universitetsparken 15, Copenhagen E 2100, Denmark. 6 Conservation Science Group, Department of Zoology, University of Cambridge, Cambridge CB2 3EJ, UK. 7Centre for Biodiversity and Environment Research, University College London, Gower Street, London WC1E 6BT, UK. 8 Indicators and Assessment Unit, Institute of Zoology, Zoological Society of London, Regent's Park, London NW1 4RY, UK. 9 United Nations Environment Programme World Conservation Monitoring Centre (UNEP-WCMC), 219 Huntington Road, Cambridge CB3 0DL, UK.10 International Union for Conservation of Nature, 28 rue Mauverney, Gland 1196, Switzerland. 11 World Agroforestry Center (ICRAF), University of the Philippines Los Banos, Laguna 4031, Philippines. 12 School of Geography and Environmental Studies, University of Tasmania, Hobart, Tasmania TAS 7001, Australia. 13 Woodley and Associates, Chelsea, Quebec J9B 1T3, Canada. * These authors contributed equally to this work. Correspondence and requests for materials should be addressed to M.D.B. (email: or to I.D.C. (email:

Biodiversity is in crisis1,2. A key response to global biodiversity declines3 and the associated threatening processes4, has been the establishment of protected areas (PAs). PAs underpin most global and national conservation strategies5, covering at least 15.4% of global land surface area6. The importance of PAs is set to increase further given the latest Convention on Biological Diversity targets to increase global land coverage of PAs to 17% by 2020 (ref. 7). Ensuring that biodiversity is maintained within-PA boundaries is consequently fundamental to achieving global conservation goals.

Despite a central objective of PAs being to conserve wildlife populations within their boundaries, many PAs are experiencing undesirable wildlife population declines8,9. Worldwide, wildlife population changes in PAs are patchily documented10, unquantified and poorly understood. While a commonly held conception is that PAs are effective at maintaining wildlife populations within their borders, this assumption has not been widely tested. Conversely, the perception that some PAs are failing, or at least performing inadequately, has precipitated calls to radically change both conservation decision-making and PA management11, and emphasized the need to ensure that PAs are effectively managed in the long term12.

It is thus vital to quantify how well PAs are conserving wildlife13, and to identify enabling conditions and barriers to effective conservation. By identifying properties of those PAs more likely to maintain wildlife populations, it will be possible to

ensure new PAs are established in a more spatially and financially efficient configuration and maximize biodiversity outcomes within resource constraints. Without a better understanding of those factors that contribute to wildlife outcomes in PAs, then their selection, design and management are likely to remain sub-optimal.

Wildlife population change is an important and useful metric for evaluating wildlife conservation outcomes in PAs. It is sensitive to long-term environmental change14, often directly linked to PA objectives, and valuable in diagnosing extinction risk15. Critically, population trends can also quantify biodiversity change in a variety of habitat types, including savannah and other non-forested habitats, complementing studies of PA impacts on maintaining forest cover16. We compiled an extensive data set of 1,902 vertebrate population abundance time series from 447 terrestrial PAs, and calculated bird and mammal population trends for 556 species as a metric of PA effectiveness in meeting conservation goals. Direct counterfactual data17 from similar but unprotected populations were not available due to insufficient monitoring effort outside PAs. We conducted a broad-scale evaluation of within-PA wildlife abundance trends (change over time) to identify properties of PAs that contribute to variation in trends among PAs. Learning from conservation outcomes for wildlife among PAs has the capacity to dramatically improve policy and management for native wildlife in PAs by promoting the propagation of enabling practices globally.

Figure 1 | Locations of our PA population time series. Countries in grey are those included in the analysis. (a) Proportionally sized pie charts indicate the number of bird (red) and mammal (blue) time series in each PA. (b) Mean per cent annual change in population abundances for each PA represented. Lighter (more yellow) dots represent greatest declines, and darker (more blue) dots greatest increases.

° 100 -

Decreasing abundance

Increasing abundance

-50 -40 -30 -20 -10 0 10 20 30 40 50

Decreasing abundance

Increasing abundance

-50 -40 -30 -20 -10 0 10 20 30 40 50

Annual change in population size (%)

Figure 2 | Frequency distribution of wildlife abundance changes in protected areas. (a) Changes by species type. (b) Changes by location. In a, grey shows all species, green shows all mammals and blue shows all birds. In b, grey shows all sample populations, green shows sampled populations in Africa and blue shows sampled populations in Europe.

We find that on average PAs are maintaining the abundance of populations of monitored birds and mammals within their boundaries, and that wildlife population trends are more positive in countries with higher development scores, as well as for larger bodied species. This body mass finding suggests active management can overcome disadvantages of more severe threats experienced by larger species18,19. It also suggests there is a need to manage smaller species directly, rather than assuming that management actions targeted at conservation of iconic taxa will lead to effective conservation of all species. Our results also underscore the need to address social and economic conditions that support PA management—as these appear to be critical for maintaining wildlife populations within PA borders.


Dataset. The full data set contains 1,902 population time series for 556 species in 447 PAs (Supplementary Table 1, Supplementary Data 1 and 2) across 72 countries (Fig. 1 and Supplementary Fig.1), from time periods between 1970 and 2010. The species in our data set are dominated by large mammalian herbivores and waterfowl (body mass distributions— Supplementary Fig. 2), reflecting taxonomically uneven global monitoring efforts. Note our data end before the recent poaching crisis that has significantly reduced populations of African elephants (Loxodonta africana) and also affected rhinoceros species (Diceros bicornis and Cerotherium simium)20.

Overall Trends. Overall, the mean percentage annual change in population size within PAs was near zero (slightly positive: mean 0.52%, median 0.81%, s.d. 12.7, Fig. 2). Overall, bird trends were marginally positive (mean annual change 1.72%, median 1.71%, s.d. 12.45, Table 1), whereas mammal trends were slightly negative (mean — 1.00%, median — 0.62%, s.d. 12.45, Fig. 2a). Trends in Europe were more positive than those in Africa (Fig. 2b and Table 1).

Population trends showed substantial variation across species and PAs (Fig. 1b). To explore this variation, we tested factors previously identified as likely to influence biodiversity outcomes in PAs10, and collated data sets addressing those factors (Supplementary Tables 2-4). We examined six groups of possible influences: (1) PA design (for example, size, shape, IUCN management category), (2) socio-economic context of the region and country in which the PA is located (for example, wealth, corruption), (3) species' traits which might determine

Table 1 | Annual percentage population change in each

data set.

Data set Mean Median s.d.

Global 0.52 0.81 12.72

Mammal -1.00 - 0.62 12.45

Bird 1.72 1.71 12.45

Africa -1.79 - 1.67 13.68

Europe 2.05 2.15 12.05

response (for example, body mass), (4) local human impacts (for example, road density, land-use change), (5) biophysical context (for example, PA elevation) and (6) time series characteristics (for example, length). We used linear mixed effect models to account for the hierarchically nested data structure, and in addition to a global model, produced separate models for mammals and birds, and for the two most data-rich regions (Europe and Africa).

Model Results. Correlates of species population trends in PAs were identified across four of the six groups of factors (Figs 3-5, Supplementary Fig. 3 and Supplementary Tables 5 and 6). In order of importance: first, in most models, population trends were more positive in areas with higher national Human Development Index scores (global, mammal, bird and Europe models), and greater Gini indices (that is, level of income inequality, in global and bird models; see also Supplementary Note 1). Second, population trends increased with body mass (in all models) and differed among taxonomic classes (all models where tested, except Europe). Third, among our measures of anthropogenic impact, population trends were positively correlated with local road density (mammal and Africa models) and local Human Influence Index (Europe model). Fourth, population trends in African PAs were more positive in later years. Neither of the other two groups of factors, PA design or biophysical context, emerged as significant fixed effects in any models: in all models species and socio-economic factors were more important. Models were found to be robust to the effects of phylogenetic influence as reflected by taxonomy, and exhibited low sensitivity to statistical outliers such as the positive outliers of African elephants and rhinoceroses. Note our abundance data pre-dates the post-2008 surge in illegal hunting of elephants and rhinoceroses20,21.

Birds Mammals

-1-1-1- U_1-1-1-1-r

40 50 60 0.3 0.5 0.7 0.9

Gini index HDI


Figure 3 | Partial-effects plots for variables in the most parsimonious global model. Partial-effects plots showing fitted relationships between change in population size and (a) body mass, (b) taxonomic class, (c) Gini index and (d) Human Development Index (HDI). In a, c, and d, dashed lines are 95% credibility intervals. In b, the circles indicate the estimated partial effect size for each factor level with credible intervals displayed as error bars.


Given the central role of PAs as global conservation tools, it is reassuring that, on average, monitored populations within PAs are stable or marginally positive (mean 0.52% annual increase). The differences between birds and mammals (birds more positive), and between Europe and Africa (Europe more positive) are likely an effect of differing regional histories, as well as current pressures. Large-scale land conversion outside the tropics has led to broad-scale historical extirpations22. More recent policy changes in Europe have led to widespread improvements in biodiversity management23, and the recovery and reintroduction of many wildlife populations, including birds in PAs24. African wildlife populations are typically more intact in absolute terms, but are under increasing anthropogenic pressure; causing abundance declines8.

Strikingly, larger-bodied species had more positive population trends in all models except Europe, indicating that PAs are more likely to maintain populations of larger-bodied wildlife than smaller bodied species. This finding is consistent across geographic realms, and taxonomic class, so it was not driven by the difference in body mass between birds and mammals. In the mammal and Africa models the relationship of body mass to trend was u-shaped (Fig. 5), suggesting perhaps that the smallest species are more resilient because of their high reproductive rates25; while intermediate sized species, lacking active management and having slower reproductive rates, are experiencing greatest decline. These findings have not previously been quantitatively demonstrated, but are supported by previous anecdotal evidence from Kenyan PAs26,

which posited a shift from elephant and rhino to poaching smaller species in response to increased penalties, before the recent ivory and horn poaching crisis.

One explanation for these findings is that threat processes of high severity (for example, hunting) are impacting intermediate bodied species particularly badly—a pattern previously detected among threat processes for mammals and birds27-29. Management effort and external project funding is commonly prioritized towards large-bodied flagship and charismatic species30, and larger species are key to tourist revenues and public/political priorities31. As a result, monitoring and management actions focus on the needs of these species and it is likely that smaller species do not receive the same benefits from PAs. We detail interactions between body mass, threatening processes and influential factors in a conceptual diagram (Supplementary Fig. 4). Larger species also tend to be preferred for ecological study and monitoring. Consequently, population declines are likely to be noticed sooner and their causes better understood, leading to more effective management responses. This result has substantial implications for future allocation of conservation effort among species within PAs.

Indicators of greater human wealth (gross domestic product) and development (Human Development Index) were associated with more positive population trends. Wealth and development have been shown to have a complex relationship with conservation. Poor outcomes have been associated with both historical and accelerating threatening processes, while increasing wildlife trends can be associated with greater awareness and

Body mass HDI Gini index

Body mass HDI

Road density in 25 k Body mass2 HDI2 HDI3

I +++ +++

Road density in 25 k Mid year Body mass2 Mid year2 Body mass

+++ +++

Body mass HDI

Land use change (Hll in 25 k)

HDI ++ c

Gini index ++

Body mass 1 1 ++ 1

1 1 ++ i

-1.0 -0.5 0 0.5 1.0

Parameter estimates for the most parsimonious (preferred) model in each dataset

Figure 4 | Informative continuous fixed effects in preferred models for each subset. Data are shown for (a) global, (b) mammals, (c) birds, (d) Africa, (e) Europe. Bar lengths indicate the size of the parameter estimates and their explanatory power in the model. Significance levels derived from highest posterior credibility intervals from Markov-Chain Monte Carlo methods: NS P>0.05; + P<0.05; + + P<0.01; + + + P<0.001. Colours show groups of potential explanatory factors: dark blue, species' traits; light blue, socio-economic context; green, local human impacts; and yellow, time series characteristics. Superscripts indicate power (that is, squared, cubed) of the variable for those variables which exhibit higher order relationships.

management capacity, that is, wealthier countries may have more resources available for PA management, and spend more on conservation32-34. Moreover, human populations in wealthier areas have less need for resources directly extracted from PAs to support livelihoods (for example, bushmeat, firewood). A final contributing factor to this pattern could be historical species loss (an extinction filter effect)—with the biota of wealthier countries having been more extensively purged of species which are more sensitive to anthropogenic threats35.

This finding of an association with development is promising, as it shows that additional capacity or decreased dependence on natural resources, can help ameliorate wildlife declines in PAs. Thus, economic development and the associated improvement in food security and governance may lead to effective conservation. In reality, outcomes are likely related to a balance of multiple socio-economic factors. Regardless of which of these explanations is dominant, extra effort will be required to retain species in developing regions. However, to avoid the extinction filter in developing regions resulting in wildlife loss catching up with historic losses in developed countries, extra effort will be required to retain species in these countries. The marginal finding that higher national Gini indices predict more positive population trends is unexpected. The result can be explained by a small number of countries (for example, South Africa) that have both very high Gini indices and relatively positive wildlife trends.

Counter-intuitively, more positive trends in PAs were correlated with increased anthropogenic landscape changes,

indicated by locally denser road networks in the mammals and Africa models, and land-use change in the Europe model. Development of new roads may open up areas and cause wildlife and habitat declines36, but in areas of historical road construction there may be a filter effect with PAs currently experiencing a wildlife recovery, especially for certain robust species37. Wildlife populations surrounded by extensive historical land clearing and roads perhaps experienced their declines before our data set but are now stable, whereas places with fewer roads have recently experienced development and clearing, driving declines38. Conversely, higher road density may correlate with increased levels of resources for management, improving PA outcomes: road density was found to be positively correlated with greater PA management effectiveness by Geldmann et al.37

Several well-studied ecological factors that conservation theory predicts should be important determinants of wildlife trends, including PA size39 and PA shape40, have no explanatory power in our models. We do not interpret the lack of significance of these variables to mean they are not important. We suggest that over the timescales addressed by our data their influence is overwhelmed by either the priorities of managers, the more general landscape scale recovery of large species23, or the socio-economic context of PAs, but we cannot not easily differentiate between these drivers. Over multi-decadal timescales the ecological drivers are still likely to be influential, but our results show the importance of managing PAs for more immediate threats. A substantial quantity of research effort examines optimal design of PAs and PA networks with respect to

■ ■ « « ■ ..... < -0.05 0.00 0.05 0.10 0.70 0.75 0.80 0.85 0.90 0.95 Birds Mammals -0.04 0.00 0.04

Body mass (grams) (log10) HDI C|ass Land use change

(Hll in 25 km)

Figure 5 | Partial-effects plots for the most parsimonious model of each subset modelled. Data are shown for (a) mammals, (b) birds, (c) Africa, (d) Europe. Dashed lines are 95% credibility intervals based on MCMC sampling with 10,000 samples. They show the relationship between each fixed effect: (a) body mass, HDI, road density; (b) body mass, Gini index, HDI; (c) body mass, class, road density, mid year; and (d) body mass, HDI, class, land-use change (Human Impact Index (HII) in 25 km), and population trends when all others are held constant. For categorical variables (class) dots indicate the estimated partial effect sizes of each factor level, and error bars show 95% credible intervals. HDI, Human Development Index; MCMC, Markov-Chain Monte Carlo.

ecological processes, but our findings suggest that greater conservation benefit would result from a focus on optimizing network design considering management and human influences on PAs.

The ability of PAs to maintain species populations is critical to global conservation. Our finding that abundances of monitored birds and mammals are maintained within PA boundaries is encouraging. However, our results show that PAs do not work

equally well for all species or in all circumstances. Moreover the recent poaching crisis in Africa shows that population gains can be rapidly reversed, if the threatening processes grow too large to mitigate20, emphasizing the need to scale management effort with threat intensity41,42. In addition, the time series data for this analysis are from PAs that are older, larger, and further removed from humans than most PAs (see Supplementary Fig. 5). Hence, we should not be complacent. These findings are likely to represent a best-case scenario for protected populations as they are sourced from PAs experiencing lower than average anthropogenic threat. If we expect PAs to act as refuges for all species in perpetuity, then a wider range of species must be targeted for management, and a particular focus on conserving medium-sized species may be required.

Further, it is clear that PAs do not exist in a vacuum. Our results show that the social and economic conditions that support PA management are critical for the maintenance of wildlife populations within PA boundaries. Anthropogenic drivers of wildlife abundance change appear to have more influence than those drivers affecting ecological processes such as PA size, at least over the period of a few decades. Much of the research effort into PAs targets ecological processes but it seems greater return would come from focusing efforts more on anthropogenic drivers of PA performance.

Managing PAs' socio-political, rather than simply ecological, dimensions is pivotal to wildlife conservation in PAs. Human dependence on PA resources in poorer countries must be addressed if existing PAs are to retain their contents as these nations continue to develop. Finally, to understand the return on our investment in PAs in the long-term, best practice adaptive management and systematic monitoring of biological outcomes, including appropriate counterfactual monitoring, is essential43. The tools to understand impact and improve outcomes exist—but we must strengthen the will and capacity to implement them.


Wildlife population time series data. A global database of population abundance time series of all available time series for native birds and mammals in terrestrial PAs worldwide was compiled from sources including the Living Planet Database44, PA agencies, published literature, grey literature and non-governmental organizations. Time series in the data set represent the majority of the data available globally to address this topic. Time series consisted of population abundance count estimates, or proxies of abundance such as nest density, mark-recapture or density estimates. Marine species, other than those with at least one critical life history stage on land (for example, breeding colonies), were excluded. We used population time series that met a number of criteria based on Collen et al.44: (1) the technique used to measure abundance was comparable over the length of the time series; (2) the geographic location of the population was provided; (3) the majority (> 50%) of the measured population was within a PA; and (4) population time series were a minimum of 5 years in length between 1970 and 2010, with at least three measures of abundance within that time period (that is, an estimate was not required for every year within a time series). Where data were available for multiple time points in a single year (for example, wet season and dry season), data were standardised to obtain a single per annum abundance estimate for each population. Standardization was carried out using the most appropriate and comparable method given data type and species ecology to obtain an estimate of species population changes most likely to remain consistent through time, and accurately represent population change, in accordance with established LPI database practice44,45. For instance, monthly abundance counts were averaged (using mean estimates) for resident non-irruptive species.

Estimation of wildlife population trends. Population trends were estimated by fitting a generalized linear regression model on time with a log-link function to each population time series46,47. The trend was taken to be the slope value from each regression. Fitting a log-link function assumes the response variable (population abundance) has a Poisson error distribution, and that the logarithm of its expected value can be modelled by a linear combination of unknown parameters, and is most appropriate for zero-inflated data, such as count and catch per-unit effort data46. Leading and trailing zero values (that is, zeros which occurred at the beginning or end of time series) were excluded from each time series before calculating trends; such zero values generally occur when populations are present at a level below that at which the sampling method is able detect

individuals, rather than when a population has been extirpated or introduced. These zeros are therefore inaccurate estimates of abundance, yet they exert undue influence on the estimated slope values (population trends). Slope values exceeding ± 0.5 on the natural log scale were excluded following inspection of the data distribution. Such values are equivalent to annual rates of population growth and decline of more than 50% per annum and biologically implausible over periods greater than 5 years (equivalent of a population of 1,000 individuals increasing to 7,594 in 5 years or 57,665 in 10 years). Such values are likely the result of errors in surveys or data entry48-50.

Explanatory variable selection and preparation. Variables were selected to represent characteristics of PAs considered most likely to be important determinants in maintaining wildlife populations in PAs for vertebrate species, based on both theoretical supposition and empirical observation as identified via a comprehensive literature review10 and can be regarded as belonging to six groups of possible influences: PA design (for example, size, shape, IUCN management category), socio-economic context of the region and country in which the PA is located (for example, wealth, corruption), species' traits which might determine response (for example, body mass), local human impacts (for example, road density, land-use change), biophysical context (for example, PA elevation), and time series characteristics (for example, length). Several variables suggested by the literature to be important for determining PA outcomes at the site scale (for example, PA-specific management budgets, threat intensity) were unavailable for most sites, and therefore could not be tested in this analysis. The key socio-economic variables used in the analysis were only consistently available at the national scale, not at finer spatial scales relative to the PAs. A table summarising the justification for each variable, a priori hypotheses, and references underpinning the selection of each factor are provided in Supplementary Table 2. For each explanatory variable, information was collected at a scale based on a combination of availability, standard practice and relevance for population dynamics. Generally, this meant that the finest resolution data available at a global scale was used. Descriptions of the explanatory variables and data sources and resolution are given in Supplementary Table 3.

Spatial data were analysed using ArcGIS 10.0 and R 2.15.0 (ref. 51) using the packages raster52, geosphere53, maptools54, rgeos55, rgdal56 and sp57. PA boundaries were calculated using spatial information from the World Database of Protected Areas58. Spatial data relating to PA context (for example, human population density) were calculated in buffers of three different sizes (5, 10 and 25 km) around each PA polygon. Multiple buffer sizes were used because it was not known a priori over what distance potential correlates would be most likely to be acting. PAs represented only by point information in the WDPA were included in the analysis by creating appropriately sized circular buffers using the ArcGIS buffer tool under an equal-area (Mollewide) projection59. The size of the buffer was given by the area (km2) recorded for the individual PAs in the WDPA. In cases where the projection of the PA data set differed from the explanatory data layer, the PAs were reprojected using ArcGIS. Details of the preparation of the individual explanatory variables are provided in Supplementary Tables 3 and 4. Several explanatory data sets were available as raster layers. The R raster package52 was used to overlay the PA polygon onto those cells and calculate metrics based on the underlying raster cell values. However, using the raster package, overlayed polygons must cover the centre of a raster cell to be considered as inside the polygon. Small, spatially complex PA polygons may only encompass a few raster cells, particularly when overlayed on raster layers with coarse spatial resolution. The selected cells may therefore poorly represent the overall cell area truly overlapped by the polygon. For this reason, we disaggregated each raster layer by a factor of 10 (30 for the spatially coarse agricultural suitability layers) using the R raster package for small PAs (<2,000 cells overlapped by the PA polygon). We also enforced a minimum overlapping area for very small PAs before performing all calculations (noted in Supplementary Table 3 in terms of the equivalent number of raster cells before disaggregation). For some variables (noted in Supplementary Table 3), PA polygons were clipped of adjoining marine areas in ArcGIS using the World Vector Shoreline Plus layer (

Explanatory variables were transformed to normalize distributions, where necessary, and standardised using the R function 'scale' so all variables in the data set had equal means and standard deviations but different ranges. Normalized and standardised variables were evaluated to determine collinearity by visual inspection of the data and by calculating Pearson's correlation coefficient (see Supplementary Fig. 6).

Population trend modelling relative to explanatory variables. The slope of the generalized linear regressions for all populations was used as the response variable to address the question: what factors predict trend in abundance for bird and mammal species in terrestrial PAs? We applied a linear mixed-effects modelling approach to explore the key correlates of wildlife population of birds and mammals through time in PAs. Explanatory variables were hierarchically spatially structured; at the species, PA and national levels. We used linear mixed-effects modelling to account for the data structure and investigate the relationship between population trends and the suite of potential explanatory variables. Mixed-effects models allow partitioning of the variance in population trends in a nested hierarchy60,61. In this analysis, the data were structured such that each population trend referred to a

particular species within a given PA. Most sites contain several species, and most species occur at several sites so there are multiple observations of individual species across a suite of site subsets, such that individual population trends are not mutually independent. Further, populations were distributed non-randomly across continents and countries. All models were implemented in R 2.15.0 (ref. 51) using the packages lme4 (ref. 62) and MuMin63. Random effects were kept consistent in all models with species, site (PA), and country fitted as random effects in every model based on the a priori understanding of the data structure. In models where both birds and mammals were present, taxonomic Class was included as a fixed effect in the models. Models were generated for the entire global data set and subsets for taxonomic class and geographic realm. Subset models were generated for: birds, mammals, Africa and Europe. Although some population data was available outside Europe and Africa, other geographic realms had insufficient sample sizes to construct subset models.

Modelling procedure. Model selection was made following an information-theoretic approach, using the corrected Akaike Information Criterion AICc46,64. Forward and backward stepwise model selections were conducted to narrow the set of candidate models. Thereafter all possible subsets of candidate models were compared using the function dredge63. This included testing all plausible interactions, and polynomials (orthogonal squares and cubes). Some variables could not be fitted simultaneously as they were highly collinear (Pearson's correlation coefficient >0.5). In these cases the variable with the greatest explanatory power, as defined by the best AICc value, was found by substitution. Substitution was conducted by exchanging variables in the model that were collinear and expected to explain the same component of the variation to assess which of them provided the best fit to the data. Fit was assessed by change in AICc, and during substitution the remainder of the model specification was held constant. Substitution was carried out both during initial data exploration and before final model selection for each data set.

For each model with a DAICc of less than four the fitted residuals were examined using qq-plots, Cook's distance leverage plots and histograms. Models selection was tested for sensitivity to outliers through systematic removal and replacement of outlying data. Outliers were identified using qqplots and Cooks distance leverage plots of the model residuals. When outliers were removed the preferred models and their parameter estimates were remarkably consistent, and additionally exhibited low sensitivity the removal of extremely large species such as Africa Elephant (L. africana) and rhinoceroses (D. bicornis and C. simium). The elephant and rhino outliers were highly positive; if data including the recent years' poaching-related population declines of these species had been available these species may not have been outliers. Fitting genus and order as random effects, to test for phylogenetic influence as reflected by taxonomy, did not improve model fit significantly for any models. When they were fitted, body mass parameter estimates remained stable and significant.

Effect sizes. Markov-Chain Monte Carlo Highest Posterior Density estimates with 10,000 samples were calculated to estimate effect sizes and 95% credibility intervals46 for the most parsimonious models. Partial-effects sizes were calculated and plots produced for parameters of the best-fit models (Fig. 5). Variable relative importance was calculated following Zuur et al.65 using Akaike weights and relative frequency standardised across all models with a DAICc <4 (Supplementary Fig. 3).

Data availability. All relevant data are available from the authors upon request. The population time series are available from the Living Planet Database to registered users: Links and citations for all the freely available predictor data sets are available in Supplementary Table 3.


1. Butchart, S. H. M. et al. Global biodiversity: indicators of recent declines. Science 328, 1164-1168 (2010).

2. Tittensor, D. P. et al. A mid-term analysis of progress toward international biodiversity targets. Science 346, 241-244 (2014).

3. Barnosky, A. D. et al. Has the Earth's sixth mass extinction already arrived? Nature 471, 51-57 (2011).

4. Geldmann, J., Joppa, L. N. & Burgess, N. D. Mapping change in human pressure globally on land and within protected areas. Conserv. Biol. 28, 1604-1616 (2014).

5. Secretariat of the Convention on Biological Diversity. Global Biodiversity Outlook 4: A Mid-Term Assessment of Progress Towards the Implementation of the Strategic Plan for Biodiversity 2011-2020 (Secretariat of the Convention on Biological Diversity, Montreal, 2014).

6. Juffe-Bignoli, D. et al. Protected Planet Report 2014 (UNEP-WCMC, 2014).

7. Secretariat of the Convention on Biological Diversity. Strategic Plan for Biodiversity 2011-2020, Including Aichi Biodiversity Targets (Secretariat to the Convention on Biological Diversity, 2010).

8. Craigie, I. D. et al. Large mammal population declines in Africa's protected areas. Biol. Conserv. 143, 2221-2228 (2010).

9. Laurance, W. F. et al. Averting biodiversity collapse in tropical forest protected areas. Nature 489, 290-294 (2012).

10. Geldmann, J. et al. Effectiveness of terrestrial protected areas in reducing habitat loss and population declines. Biol. Conserv. 161, 230-238 (2013).

11. Fuller, R. A. et al. Replacing underperforming protected areas achieves better conservation outcomes. Nature 466, 365-367 (2010).

12. Watson, J. E. M., Dudley, N., Segan, D. B. & Hockings, M. The performance and potential of protected areas. Nature 515, 67-73 (2014).

13. McDonald-Madden, E. et al. ENVIRONMENT: 'True' Conservation Progress. Science 323, 43-44 (2009).

14. Buckland, S. T., Magurran, A. E., Green, R. E. & Fewster, R. M. Monitoring change in biodiversity through composite indices. Philos. Trans. R. Soc. Lond. B Biol. Sci. 360, 243-254 (2005).

15. Hoffmann, M. et al. The difference conservation makes to extinction risk of the world's ungulates. Conserv. Biol. 29, 1303-1313 (2015).

16. Andam, K. S., Ferraro, P. J., Pfaff, A., Sanchez-Azofeifa, G. A. & Robalino, J. A. Measuring the effectiveness of protected area networks in reducing deforestation. Proc. Natl Acad. Sci. USA 105, 16089-16094 (2008).

17. Ferraro, P. J. & Pattanayak, S. K. Money for nothing? A call for empirical evaluation of biodiversity conservation investments. PLoS Biol. 4, e105 (2006).

18. Cardillo, M. et al. Multiple causes of high extinction risk in large mammal species. Science 309, 1239-1241 (2005).

19. Bennett, P. M. & Owens, I. P. F. Evolutionary Ecology of Birds - Life histories, Mating Systems and Extinction (Oxford University Press, 2002).

20. Wittemyer, G. et al. Illegal killing for ivory drives global decline in African elephants. Proc. Natl Acad. Sci. USA 111, 13117-13121 (2014).

21. Biggs, D., Courchamp, F., Martin, R. & Possingham, H. P. Legal trade of Africa's Rhino horns. Science 339, 1038-1039 (2013).

22. Fritz, S. A., Bininda-Emonds, O. R. P. & Purvis, A. Geographical variation in predictors of mammalian extinction risk: big is bad, but only in the tropics. Ecol Lett. 12, 538-549 (2009).

23. Chapron, G. et al. Recovery of large carnivores in Europe's modern human-dominated landscapes. Science 346, 1517-1519 (2014).

24. Eaton, M. A. et al The State of the UK's Birds 2012 (RSPB, BTO, WWT, CCW, NE, NIEA, SNH and JNCC, Sandy, Bedfordshire, 2012).

25. Brook, B. & Bowman, D. J. S. One equation fits overkill: why allometry underpins both prehistoric and modern body size-biased extinctions. Popul. Ecol. 47, 137-141 (2005).

26. Gibson, C. C. Politicians and Poachers: The Political Economy of Wildlife Policy in Africa (Cambridge University Press, 1999).

27. Owens, I. P. F. & Bennett, P. M. Quantifying biodiversity: a phenotypic perspective. Conserv. Biol. 14, 1014-1022 (2000).

28. Isaac, N. J. B. & Cowlishaw, G. How species respond to multiple extinction threats. Proc. R. Soc. Biol. Sci. Ser. B 271, 1471-2954 (2004).

29. Brashares, J. S. et al. Bushmeat hunting, wildlife declines, and fish supply in west Africa. Science 306, 1180-1183 (2004).

30. Smith, R., Verissimo, D., Isaac, N. & Jones, K. Identifying Cinderella species: uncovering mammals with conservation flagship appeal. Conserv. Lett. 5, 205-212 (2012).

31. Hilborn, R. A. Y. et al. Effective enforcement in a conservation area. Science 314, 1266 (2006).

32. McCreless, E., Visconti, P., Carwardine, J., Wilcox, C. & Smith, R. J. Cheap and nasty? The potential perils of using management costs to identify global conservation priorities. PLoS ONE 8, e80893 (2013).

33. Balmford, A., Gaston, K. J., Blyth, S., James, A. & Kapos, V. Global variation in terrestrial conservation costs, conservation benefits, and unmet conservation needs. Proc. Natl Acad. Sci. USA 100, 1046-1050 (2003).

34. Waldron, A. et al. Targeting global conservation funding to limit immediate biodiversity declines. Proc. Natl Acad. Sci. USA 110, 12144-12148


35. Balmford, A. Extinction filters and current resilience: the significance of past selection pressures for conservation biology. Trends Ecol. Evol. 11, 193-196 (1996).

36. Laurance, W. F. et al. A global strategy for road building. Nature 513, 229-232


37. Geldmann, J. et al. Changes in protected area management effectiveness over time: a global analysis. Biol. Conserv. 191, 692-699 (2015).

38. Laurance, W. F. & Balmford, A. Land use: a global map for road building. Nature 495, 308-309 (2013).

39. Maiorano, L., Falcucci, A. & Boitani, L. Size-dependent resistance of protected areas to land-use change. Proc. R. Soc. Biol. Sci. Ser. B 275, 1297-1304 (2008).

40. Ewers, R. M. & Didham, R. K. The effect of fragment shape and species sensitivity to habitat edges on animal population size. Conserv. Biol. 21, 926-936 (2007).

41. Ihwagi, F. W. et al. Using poaching levels and elephant distribution to assess the conservation efficacy of private, communal and government land in Northern Kenya. PLoS ONE 10, e0139079 (2015).

42. Wasser, S. K. et al. Genetic assignment of large seizures of elephant ivory reveals Africa's major poaching hotspots. Science 349, 84-87 (2015).

43. Craigie, I. D., Barnes, M. D., Geldmann, J. & Woodley, S. International funding agencies: potential leaders of impact evaluation in protected areas? Phil. Trans. R. Soc. B 370, 20140283 (2015).

44. Collen, B. et al. Monitoring change in vertebrate abundance: the Living Planet Index. Conserv. Biol. 23, 317-327 (2009).

45. McLellan, R., Iyengar, L., Jeffries, B. & Oerlemans, N. Living planet report 2014: species and spaces, people and places, all_publications/living_planet_report/ (WWF International, Gland, 2014).

46. Zuur, A. F., Ieno, E. N., Walker, N. J., Saveliev, A. A. & Smith, G. M. Mixed Effects Models and Extentions in Ecology with R (Springer Science, 2009).

47. Zuur, A., Ieno, E. N. & Smith, G. M. Analysing ecological data (Springer Science & Business Media, 2007).

48. Norton-Griffiths, M. Counting Animals. 2nd Edn (African Wildlife Leadership Foundation, 1978).

49. Thomas, L. et al. Distance software: design and analysis of distance sampling surveys for estimating population size. J. Appl. Ecol. 47, 5-14 (2010).

50. McCarthy, M. A. et al. The influence of abundance on detectability. Oikos 122, 717-726 (2013).

51. R Core Team. R: A Language and Environment for Statistical Computing, (Foundation for Statistical Computing, Vienna, Austria, 2012).

52. Hijmans, R. J. et al. Raster: Geographic analysis and modeling with raster data. R package version 2.2-31 (2013).

53. Hijmans, R. J., Williams, E. & Vennes, C. Geosphere: Spherical Trigonometry. R package version 1.2-28 (2012).

54. Lewin-Koh, N. J. et al. Maptools: Tools for reading and handling spatial objects. R package version 0.8-16 (2012).

55. Bivand, R. & Rundel, C. Rgeos: Interface to Geometry Engine-Open Source (GEOS). R package version 0.2-7 (2012).

56. Keitt, T. H., Bivand, R., Pebesma, E. & Rowlingson, B. Rgdal: Bindings for the Geospatial Data Abstraction Library. R package version 0.7-20 (2012).

57. Pebesma, E. J. & Bivand, R. S. Classes and methods for spatial data in R. R News 5, (2005).

58. WDPA. World Database of protected areas (WDPA) Release June 2012 (web download version). The WDPA is a joint product of UNEP and IUCN, prepared by UNEP-WCMC, supported by IUCN WCPA and working with Governments, the Secretariats of MEAs and collaborating NGOs (2012).

59. UNEP-WCMC. World Database on Protected Areas User Manual 1.1., http:// (UNEP-WCMC, 2015).

60. Pinheiro, J. C. & Bates, D. M. Mixed-Effects Models in S and S-PLUS (US Government Printing Office, 2000).

61. McCulloch, C., Searle, S. & Neuhaus, J. Generalized, Linear, and Mixed Models, (Wiley, 2001).

62. Bates, D. & Maechler, M. R package 'lme4' (2010).

63. MuMIn: Multi-model inference (CRAN, 2012).

64. Burnham, K. P. & Anderson, D. R. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach (Springer-Verlag, 2002).

65. Zuur, A., Ieno, E. N., Walker, N., Saveliev, A. A. & Smith, G. M. Mixed Effects Models and Extensions in Ecology with R (Wiley International, 2009).


We are members of the IUCN WCPA-SSC Joint Taskforce on Biodiversity and Protected Areas, which supported this project. We thank W.N. Venables for input on the design of the statistical analysis, We thank other members of the taskforce who contributed initial thoughts and feedback for this project: S. Andelman, N. Dudley, E. Enkerlin, N. Lopoukhine, K. MacKinnon, S. Stuart, K. Redford and S. Stolton. We thank L. McRae for her assistance with LPI data and J. Ringma for graphics assistance. We thank WWF and ZSL who funded and compiled the living planet database and everyone who contributed to the database. We thank H. Possingham and C. Carbone for manuscript comments and suggestions, and E. McDonald Madden for advice. Funding was provided by IUCN's Global Protected Areas Program and Parks Canada (S. Woodley), European Union's Development Fund through the Biodiversity and Protected Areas Management (BIOPAMA) Programme, School of Geography Planning and Environmental Management (M.D.B., M.H.) and the Centre of Excellence for Environmental Decisions at the University of Queensland (M.D.B.) and Natural Environmenlt Research Council studentship NER/S/A/2006/14094 with CASE support from UNEP-WCMC (I.D.C.).

Author contributions

M.D.B., I.D.C., A.B., S. Wood. and M.H conceived the project. I.D.C., M.D.B. designed analysis. M.D.B., I.D.C., S. Wood. M.H., L.H. J.G. prepared the data. M.D.B., I.D.C., L.H. analysed the data M.D.B., I.D.C., L.H. and S. Whit. displayed the items. M.D.B. and I.D.C. wrote the manuscript, with input from all authors. All authors contributed to the paper.

Additional information

Supplementary Information accompanies this paper at naturecommunications

Competing financial interests: The authors declare no competing financial interests.

Reprints and permission information is available online at reprintsandpermissions/

How to cite this article: Barnes, M. D., Craigie, I. D. et al. Wildlife population trends in protected areas predicted by national socio-economic metrics and body size. Nat. Commun. 7:12747 doi: 10.1038/ncomms12747 (2016).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit

© The Author(s) 2016