Understanding the global geographical distribution of extinction risk is a key challenge in conservation biology. It remains controversial, however, to what extent areas become threat hotspots simply because of high human impacts or due to predisposing ecological conditions. Limits to the taxonomic and geographical extent, resolution and quality of previously available data have precluded a full global assessment of the relative roles of these factors. Here, we use a new global database on the geographical distributions of birds on continents and continental islands to show that, after controlling for species richness, the best predictors of the global pattern of extinction risk are measures of human impact. Ecological gradients are of secondary importance at a global scale. The converse is true for individual biogeographic realms, within which variation in human impact is reduced and its influence on extinction risk globally is therefore underestimated. These results underline the importance of a global perspective on the mechanisms driving spatial patterns of extinction risk, and the key role of anthropogenic factors in driving the current extinction crisis.
Understanding the geographical distribution of extinction risk and its causes are key challenges in conservation biology, and are central to determining spatial priorities for the focus of conservation responses. Major determinants of extinction risk across space include not only anthropogenic environmental impacts but also variation in predisposing ecological conditions (Forester & Machlis 1996). The former include human population density, agricultural and urban land-use, species exploitation, introduced species and disease, and anthropogenic climate change (Soulé 1991; Forester & Machlis 1996). Predisposing ecological conditions include, but are not exclusive to, the availability of ambient environmental energy which is thought to influence speciation rates and thus the occurrence of neoendemics (Rohde 1992), the availability of productive environmental energy which is thought to limit overall and individual species population numbers (Wright 1983), absolute species numbers which may influence food web structure and thus the likelihood of extinction cascades (Gaston 2002), and surface topography which influences the occurrence of narrowly distributed species (Richerson & Lum 1980).
Despite a number of valuable regional studies, restrictions to the taxonomic and geographical extent, resolution and quality of previously available data, have thus far largely precluded a full global assessment for a major taxon of the relative roles of human impact and predisposing ecological factors in determining threatened species richness (Kerr & Currie 1995; Balmford et al. 2001; McKinney 2001; Norris & Pain 2002; Luck et al. 2004; Scharlemann et al. 2005). Given that the sensitivity of individual species to human population density has been shown to vary within and between biogeographical regions (Woodroffe 2000), it remains an open question whether the relative importance of factors indicated by individual regional studies will generalize to other regions or globally. The frequently incomplete representation, within individual continents, of global variation in human impact and ecological gradients, as well as the distinct evolutionary histories of species occurring in different regions, further contribute to this uncertainty. The relevance to conservation policy-making of the answer lies both in its indication of the wider applicability of regional findings and in the confirmation of ultimate causes operating globally. Here, we present an analysis of human and ecological determinants of spatial patterns of extinction risk at continental and global scales. We use a new database on the geographical distribution of the breeding ranges of extant bird species on continents and continental islands on an equal-area grid at a resolution comparable to 1° latitude×longitude (Orme et al. 2005).
We tested equal numbers of predictors of human impact and ecological condition so as to avoid a priori skewing analyses in favour of finding the greater importance of predictors from one category over the other. Building on previous demonstrations of their potential importance in shaping spatial patterns of threatened species richness, for indices of human impact we used human population density, economic activity (purchase power parity gross domestic product, GDP), and extent of agricultural and urban land-area (Kerr & Currie 1995; Balmford et al. 2001; McKinney 2001; Norris & Pain 2002; Luck et al. 2004; Scharlemann et al. 2005). For ecological gradients, we used mean annual temperature as a measure of available ambient energy, while for productive energy we used the Normalized Difference Vegetation Index (NDVI). In addition to topographical variability (elevation range), we used number of land-cover types (habitat diversity) as an alternative measure of habitat heterogeneity. To minimize the risk of including spurious variables in our analyses, we built a multivariate minimum adequate model (MAM) based on regression methods that accounted for spatial autocorrelation in the response variable. To test whether the results of our global model could have been predicted by analyses conducted at a smaller geographical extent, we used the same methodology to construct models separately for six major biogeographic realms (Olson et al. 2001).
2. Material and methods
(a) Species data
The analyses presented here are based on a database of distribution maps for 9626 extant, recognized bird species constructed using a variety of published sources (for details see Orme et al. 2005, in press) and following a standard avian taxonomy (Sibley & Monroe 1990). The polygon breeding range vector maps were converted to an equal-area grid using a Behrmann projection at a cell resolution of 96 486.2 m. This grid cell size is equivalent to 1° longitude and 1° latitude at 30° latitude N/S (1/360th of the width of the globe under a Behrmann projection using the WGS84 datum). The global grid therefore contains 360 by 152 cells, omitting the partial cells at latitudes higher than 87.13°. Species were scored as present in a grid cell if any of the available sources indicated that the breeding range fell within the cell boundaries. Threatened species were those classified as Critical, Endangered and Vulnerable, but not those in lower risk categories (Near Threatened, Least Concern) or other categories (Data Deficient and Not Evaluated; BirdLife International 2000). Where necessary, we converted the taxonomy used in BirdLife International (2000) back to the standard avian taxonomy (Sibley & Monroe 1990), and calculated the number of threatened species in each grid cell. Biogeographic realms were delimited using the World Wildlife Fund ecoregions map (Nearctic, Palaearctic, Neotropical, Afrotropical, Indo-Malayan, and Australasian; Olson et al. 2001). The final dataset used for analyses omitted grid cells falling within Oceania or Antarctica, since environmental data were not available for these realms. Remaining true oceanic islands, defined as any land area located further than 200 km from the edge of continental shelf, were also omitted since these are known to differ markedly in the kinds and intensities of evolutionary and threatening processes affecting constituent avifauna, in comparison with continental locations (Manne et al. 1999; Blackburn et al. 2004; Duncan & Blackburn 2004). Finally, so as to avoid bias in terms of the contribution of coastal land-area to the regression models, grid cells with less than 50% land-cover were omitted from the final dataset.
(b) Environmental data
Data for the eight selected environmental and human impact predictors (see above) were each re-projected and re-sampled to the same equal-area grid as the species richness data. Human population density, GDP, NDVI, agricultural and urban land-area, and elevation range were all log10-transformed for the analysis.
Sources and raw resolutions of the eight selected environmental variables are as follows: (a) human population density (people km−2) for 1995 at 2.5 arc-min resolution (CIESIN 2003) derived from human population census data for 127 105 administrative units, and based on national population estimates that have been adjusted to match the UN national estimated population for each country; (b) purchase power parity GDP data (US $) for 1990 at 0.25° resolution (CIESIN 2005); (c) agricultural land-area (km2); (d) urban and built-up land-area (km2); and (e) number of land-cover types (habitat diversity) occurring in a grid cell, all computed using remotely sensed data for the 12-month period between April 1992 and March 1993 at 30 arc-s resolution classified to the US Geological Survey (USGS) 25-category land-cover classification (USGS 2003a); (f) elevation range (m), maximum minus minimum elevation within each grid cell, from 30 arc-s resolution data (USGS 2003b); (g) mean annual temperature data (°C) for the period 1961–1990 at 10 min resolution interpolated from station means (New et al. 2002); (h) mean annual remotely sensed NDVI for the period 1982–1996 at 0.25° resolution (ISLSCP Initiative II 2005). Agricultural land-area (above) was computed as the sum of all agricultural land-use classes from the USGS data (USGS 2003a) (2, dryland cropland and pasture; 3, irrigated cropland and pasture; 4, mixed dryland/irrigated cropland and pasture; 5, cropland/grassland mosaic; 6, cropland/woodland mosaic). In order to standardize the definition of terrestrial land-area across raw environmental datasets, each was overlaid with a high-resolution terrestrial areas map (ESRI 1993) prior to re-sampling to the 1° Behrmann gird. Raw-data cells, or portions of cells, falling outside this definition of land-area were excluded from re-sampling calculations, and the latter were weighted by the land-area associated with each remaining raw-data cell.
(c) Statistical analyses
To deal simultaneously with spatial autocorrelation and a response variable that was not normally distributed, analyses were based on a Poisson errors generalized linear mixed modelling (GLMM) method (SAS, Littell et al. 1996) in which an exponential spatial covariance structure is fitted with longitudinal and latitudinal cell centroid values as spatial variables using Proc Glimmix v. 1.0 add-in in SAS v. 9.1.3. The choice of the exponential, over other spatial covariance structures, was based on inspection of semi-variograms of non-spatial Poisson error model residuals. Spatial GLMMs took account of the differences among major biogeographical realms in spatial autocorrelation by estimating the maximum geographic distance or range parameter (ρ), measured in degrees, over which spatial autocorrelation in equivalent independent errors model residuals was observed to occur. This involved estimating ρ from the semi-variogram of residuals of non-spatial Poisson errors models that included the relevant combination of predictors, separately for each realm. All six estimates of ρ were then entered as spatial covariance parameters in the global model, with spatial autocorrelation assumed for observations within the same realm. GLMMs used the pseudo-likelihood (PL) procedure (Wolfinger & O'Connell 1993) that obtains a maximum-likelihood-like estimate of the scale parameter (φ) (Littell et al. 1996). PL does not compute a true log-likelihood, precluding use of model selection procedures based on Akaike's Information Criterion, and forward stepwise model-building procedures were employed to determine MAMs. The fit of quadratic as well as linear terms for predictors was tested in order to allow for nonlinear relationships. Estimates of variance explained cannot be derived from spatial models that use PL, so we used percentage of total deviance explained from equivalent non-spatial models as an indication. For all geographical areas modelled, we explored collinearity among predictor variables using tolerance levels (Quinn & Keough 2002). Tolerance levels were sufficiently high (i.e. greater than 0.1, following Quinn & Keough 2002) in all cases except for the Nearctic and Neotropical realms where some redundancy was observed between human population density, GDP and urban land-area (for tolerance values see table 1 in electronic supplementary material). Hence, this had the relatively minor consequence that significance of human population density in the Neotropical MAM could not be separated from the effect of the other two socio-economic predictors.
3. Results and discussion
The global distribution of threatened avian species richness exhibits marked large-scale spatial heterogeneity (figure 1a), being highest across much of the Indo-Malayan realm and parts of the Neotropics, including areas of the Andes, Amazonia and the Atlantic coastal forests. Previous work has shown that the geographical distribution of threatened species richness is, to some extent, dependent on that of overall species richness (Kerr & Currie 1995; McKee et al. 2003). However, the proportion of species threatened (figure 1b) is far from constant and does not simply mirror the patterns for absolute numbers of either threatened (figure 1a) or non-threatened species (figure 1c). Instead, the richest continental or larger-island areas with respect to proportional threat are the southern Palaearctic, Madagascar, and New Zealand.
Our global model revealed that, after controlling for the effect of spatial variation in non-threatened species richness, human population density (figure 1d) was the primary global driver of geographical patterns of numbers of threatened bird species, followed by extent of agricultural activity as a secondary human influence (table 1). Environmental factors played a more minor role, with elevation range and NDVI entering as subsidiary factors in the global MAM. Accounting for the effects of human factors and non-threatened avian richness, we might expect gradients such as elevation range and primary productivity to be inversely related to the remaining variation in numbers of threatened species, since we could expect fewer threatened species in pristine areas that coincide with areas of lower human impact. However, the importance of elevation range as a positive predictor of threatened species richness globally may result from the influence of topography on the occurrence of restricted-range species (Jetz et al. 2004) and their inherent associated vulnerability to population decline (Stattersfield et al. 1998; Manne et al. 1999). Alternatively, in some cases this occurrence may result from range contractions from human-impacted lower elevation areas leaving remnant populations in more mountainous regions. Similarly, the positive influence of NDVI may be linked to contraction of the geographic distribution of threatened species to remaining areas of high plant productivity. The global MAM also included the significant negative influence of GDP indicating that, having accounted for other factors, areas of high economic development are coincident with lower numbers of threatened species. This is more likely to result from local extinction resulting in threatened species range contractions from areas of highest economic activity (including urban areas), rather than from any positive influence of economic development on conservation investment and effectiveness.
In contrast to the global model, our analyses of threat within biogeographic realms generally suggested that human impact variables ranked relatively lower in importance compared with ecological predictors (table 2). The only exception to this was Australasia, which showed human population density to be the primary predictor of threatened species richness, ahead even of non-threatened species richness. However, the latter was the primary predictor of numbers of threatened species for all remaining realms, and was the only predictor to enter MAMs for the Nearctic and Palaearctic, the two higher latitude realms with the lowest absolute numbers of threatened species. For the major tropical realms, subsequent predictors entered MAMs in different rank orders, with NDVI, elevation range, and temperature being the second strongest predictors, respectively, for Afrotropical, Neotropical, and Indo-Malayan realms. Likely reasons for the positive slopes for NDVI and elevation range are the same as those proposed for the global MAM, whereas the importance of temperature over NDVI or elevation range in Indo-Malaya is indicative of the impact of widespread lowland deforestation on the large-scale distribution of plant productivity.
These results have important consequences for understanding the mechanisms that underlie contemporary extinction processes. For instance, the differences between the results of our global analyses and those of the realm-specific analyses indicate that the combination of predictors of threatened species richness observed for one realm cannot be assumed to apply either to other realms or globally. One reason for this appears to be that individual realms often contain relatively limited geographical variation in human impacts (figure 2a) and/or the numbers of threatened species (figure 2b). For example, the Indo-Malayan realm shows the maximum global average levels both of threatened species richness (figure 2b) and of human population density, these being an order of magnitude higher than in any other realm. In spite of this, human population density does not even enter the MAM for Indo-Malaya since variation in human population density within the realm is relatively limited at this spatial resolution (figure 2a). Hence, the weak predictive strength of population density within Indo-Malaya belies the fact that this variable is largely driving the peaks in global avian threatened species richness. The same pattern is repeated across several of the other realms, with intra-realm variation typically being relatively minor in comparison to the global patterns, resulting in the regional models having considerably less power than the full global analysis.
Our analyses also indicate likely directions for future work in this area. In particular, the estimates of the percentage of deviance explained by models range from 21 to 66%, indicating that considerable variation in the number of threatened species remains to be explained (tables 1 and 2). This is expected, since time lags in species' responses to threatening processes on the one hand, and removal by local extinction of threatened species from high-threat areas (extinction filtration) on the other, will tend to reduce variation in numbers of threatened species present along spatial gradients in the duration and/or intensity of a given driver of threat (Pimm et al. 1995; Balmford 1996; Forester & Machlis 1996). Even globally, therefore, the full impact of human activities will be largely underestimated in the continuing absence of comprehensive historical species distribution data (Cardillo et al. 2004). Our findings thus lend support to conservation prioritization initiatives focused on areas of high human population density and impact, as well as those with the best remaining natural habitat. Moreover, lack of marked congruence in the intra-regional spatial distribution of numbers of threatened species across major taxa (Kerr & Currie 1995; Dobson et al. 1997) may belie a stronger emergent signal of taxonomic congruence globally. Hence, global studies of other tractable taxa are urgently needed to test the wider validity of these models. Nevertheless, being the first global models of extinction risk for the best-known vertebrate taxon, our findings represent one of the highest-resolution analyses of human impacts on global biodiversity currently available to inform conservation policy.
We thank T. Allnutt, B. Beehler, T. Brooks, B. Coates, J. Cromie, H. Fry, P. Higgins, D. McNicol, D. Mehlman, C. Perrins, R. Porter, H. Pratt, N. Redman, R. S. Ridgely, C. Robertson, A. Silcocks, M. Strange, M. Unwin, M. Weston, M. Whitby, P. Williams, D. Wynn, B. Young, J. Zook, A. & C. Black, Academic Press, BirdGuides Ltd, Birds Australia, Christopher Helm, Conservation International, NatureServe, Oxford University Press, the Ornithological Society of New Zealand and Princeton University Press for access to data; L. Birch, R. Prys-Jones, B. Sheldon, the Alexander Library (Oxford University), and the Natural History Museum (Tring) for access to libraries; M. Burgess, F. Eigenbrod and N. Pickup for technical assistance; three anonymous reviewers for comments; and O. Schabenberger for analytical advice. This work was funded by The Natural Environment Research Council (grant nos NER/O/S/2001/01258, NER/O/S/2001/01257, NER/O/S/2001/01230, and NER/O/S/2001/01259).