Half of all artiodactyls (even-toed hoofed mammals) are threatened with extinction, around double the mammalian average. Here, using a complete species-level phylogeny, we construct a multivariate model to assess for the first time which intrinsic (biological) and extrinsic (anthropogenic and environmental) factors influence variation in extinction risk in artiodactyls. Globally artiodactyls at greatest risk live in economically less developed areas, have older weaning ages and smaller geographical ranges. Our findings suggest that identifying predictors of threat is complicated by interactions between both biological and anthropogenic factors, resulting in differential responses to threatening processes. Artiodactyl species that experience unregulated hunting live in significantly less economically developed areas than those that are not hunted; however, hunted species are more susceptible to extinction if they have slower reproductive rates (older weaning ages). In contrast, risk in non-hunted artiodactyls is unrelated to reproductive rate and more closely associated with the economic development of the region in which they live.
Understanding the biological processes underlying species extinction is one of the most important goals for conservation biology, particularly for effective management of populations and predicting species that require immediate conservation measures (Mace & Balmford 2000). It has been shown in a variety of taxa, at both local and global scales, that biological characteristics are important determinants of variation in extinction risk (see review by Reynolds 2002). To date, studies have largely focused on overall correlates of extinction risk of mammals (e.g. Purvis et al. 2000; Cardillo & Bromham 2001; Jones et al. 2003), while not assessing the complexity of how different traits predispose species to extinction via different threatening processes (Jennings et al. 1998; Owens & Bennett 2000; Fisher & Owens 2004; Isaac & Cowlishaw 2004; Keane et al. 2005).
Anthropogenic factors such as human population density and economic development are expected to correlate with extinction risk (e.g. Davies et al. 2006) regardless of threatening process. Higher human densities lead to higher levels of human impact on species and the environment (e.g. McKee et al. 2004; Keane et al. 2005; Davies et al. 2006) and have been found to be important in explaining local and global extinctions (Brashares et al. 2001; Cardillo et al. 2004, respectively). Economic status is also important; it affects how much governments can spend on conservation efforts and how its people exploit natural resources. The situation in economically less developed countries often necessitates opportunistic and short-term exploitation of the local flora and fauna by its citizens, increasing the risk of local extinction (Rodriguez 2000; Matos & Bovi 2002). For example, gross domestic product (GDP) is negatively related to the number of threatened birds globally (Davies et al. 2006).
Habitat loss (including fragmentation and degradation) and exploitation are the two major processes threatening mammals today; they account for 46.6 and 33.9% of all recorded threats, respectively (Mace & Balmford 2000). Each threat process exerts fundamentally different selective pressures: exploitation acts directly upon the species by increasing mortality while habitat loss acts indirectly by reducing the carrying capacity of the environment. The Artiodactyla (even-toed hoofed mammals) provide a unique test case for evaluating the importance of hunting because exploitation is the primary threatening process responsible for 40% of recorded threats, while habitat loss accounts for 36% (www.redlist.org).
Predictions concerning how species traits increase susceptibility to exploitation can be divided into two main categories: those pertaining to reproductive rates and those related to hunter behaviour (see review by Fitzgibbon 1998). As an example, the Quaternary mammalian megafauna extinctions have often been associated with over-hunting and the spread of modern humans (see review by Barnosky et al. 2004), and both hunter selection of large-bodied prey (blitzkrieg hypothesis, Martin 1984) and slow reproductive rates (Johnson 2002) have been proposed as determinants of the loss of mammal species.
Hunters have been shown to prefer to hunt larger-bodied species (Mittermeier 1987; Jerozolimski & Peres 2003), and the large body size has been found to correlate with threat status in exploited birds (Owens & Bennett 2000; Keane et al. 2005) and hunted primates (Isaac & Cowlishaw 2004). However, the harvest rate may reflect the encounter rate, which is controlled by hunter numbers, hunter behaviour and prey biology (Fitzgibbon 1998). Increasing size is also associated with slow reproductive rates, which increase vulnerability to extinction via hunting. The principle of sustainable harvesting (Hartig 1796; Clark 1990) suggests that a population will remain stable even when individuals are harvested as long as off-take does not exceed rN(1−N/K), where N is the population size, K is the carrying capacity and r is the intrinsic rate of population growth. Therefore, species with fast rates of increase (r) will have a larger surplus available for harvesting and should be more resilient to increased mortality (Bodmer 1995; Bodmer et al. 1997; Jennings et al. 1998). Thus, traits that reflect reproductive speed such as weaning age, age at sexual maturity and body mass will potentially show an association with extinction risk via exploitation (Owens & Bennett 2000; Purvis 2001; also see reviews by Cowlishaw & Dunbar 2000; Kokko et al. 2001).
Here, we have constructed a multivariate phylogenetic regression model to identify the biological, ecological and anthropogenic variables that correlate with levels of threat across the Artiodactyla as given by the IUCN Red List (sensu Purvis et al. 2000). We then investigated whether different traits elevate the risk of extinction depending on whether a species is hunted or not. The majority of artiodactyls are hunted for food, often termed bushmeat hunting, while others are hunted as trophies (e.g. Urial sheep, Ovis vignei; sable antelope, Hippotragus niger) or for their soft hair (e.g. chiru, Pantholops hodgsoni). The analyses were repeated with the dataset partitioned into two: the hunted species forming one partition and the species that we could not validate as being hunted forming the other ‘non-hunted’ partition.
2. Material and methods
The IUCN (2006) Red List (IUCN 2006) was used as the measure of current extinction risk. The Red List categories were divided into six levels (following Purvis et al. 2000): least concern, 0; near threatened, 1; near threatened conservation dependent and vulnerable, 2; endangered, 3; critically endangered, 4; extinct in wild and extinct, 5, although no extinct species were added to the analysis. To ameliorate the effect of autocorrelation between species traits and the criteria for IUCN classification, all analyses were performed on species whose Red List classification was based on a decline in population density or geographical range size (criterion A), rather than the absolute measures of those variables. After the data deficient species and those not classified under criterion A were removed, 144 artiodactyl species categorized from least concern to critically endangered were available for analysis (electronic supplementary material).
Eleven predictor variables were chosen to represent the anthropogenic, hunter behaviour and reproductive rates hypotheses (summarized in table 1). Five additional variables were also included to represent other commonly cited hypotheses concerning elevated extinction risk: small population size as indicated by small geographical range and low population density (e.g. Gaston 1994); poor ecological flexibility (e.g. Brown 1971; Laurance 1991) as represented by narrow habitat breadth and dietary specialization; and, finally, large home ranges which make species especially vulnerable to habitat loss (e.g. Woodroffe & Ginsberg 1998). Thirteen quantitative traits were obtained from the PanTheria trait database (K. E. Jones, J. Bielby, A. Purvis, D. Orme, A. Teacher, J. L. Gittleman, R. Grenyer, et al., unpublished manuscript): weaning age; maximum lifespan; gestation length; age at first birth; inter-birth interval; age at sexual maturity; home range; body mass; group size; population density; geographical range; mean human population density; and actual evapotranspiration (AET). AET was included as a measure of primary productivity, which is thought to be a confounding factor when using mean human population density (Balmford et al. 2001). The database is the product of a collaborative effort to construct a comprehensive database of major biological traits for all mammal species; information is drawn from both primary and secondary sources. All quantitative data used in this analysis represent central tendency measures for each species and were +1 natural log transformed.
Gross national income (GNI; US $ millions) taken from UNEP (2006) was added as an indicator of the socio-economic condition experienced by a species across its range. The 2003 estimates of GNI were used unless unavailable; the most recent estimate was then used. Using ArcGIS, a weighted average of GNI across the species range was calculated from country estimates of GNI and the area of the species range (km2) in each country (ranges taken from Sechrest 2003). Three categorical traits were collected from the literature: mating strategy (polygamous, monogamous); diet (specialist (grazer or browser) or generalist (mixed grazer/browser or omnivore)); and habitat breadth. Habitat breadth can vary between the levels 1 and 6 and was calculated from the dataset presented in Caro et al. (2004), which lists a species as living in as many as six different habitats (grassland/scrubland, dense forest, desert, rocky, tundra and swamp). All data were converted to the taxonomy of Grubb (1993) and are available in the electronic supplementary material.
The artiodactyl dataset of 144 species was used to identify traits that correlate with extinction risk across artiodactyls. The complete dataset was then partitioned on the basis of threatening process and reanalysed to look for the differences between traits that predispose hunted and non-hunted species to extinction. Literature searches (Web of Science, Biological Abstracts and Zoological Record) using the terms Artiodactyl and hunting or bushmeat were used to identify species that are known to be hunted, and we also followed up citations in the papers found during the literature search. If the only citations we could find referred to regulated hunting then that species was not added to the hunted partition, e.g. Ovis dalli. It is important to exclude species that only experience regulated hunting from the hunted partition as, theoretically at least, hunting should not be a threatening process if appropriate quotas are set and enforced. We identified 111 artiodactyl species that are cited within the primary and secondary literature as being hunted without regulation; 94 of these were classified under criterion A in the IUCN Red List (IUCN 2006). The list of hunted species is a conservative estimate of the species hunted worldwide. Species that have experienced over-hunting in the past, but are now protected, are not included, and species with small geographical ranges are likely to be under-represented in our sample, as we are restricted to the species that live within the areas where studies of hunting have been undertaken.
All statistical analyses were conducted in R (http://www.r-project.org/) on phylogenetically independent contrasts (PICs; Felsenstein 1985), as the majority of traits showed phylogenetic patterning using Pagel's λ-statistic (results not shown; Pagel 1999). The topology presented in Price et al. (2005) was used to generate the PICs. Independent contrasts, including the modified version for categorical traits, were calculated in R (http://www.r-project.org/) using the Ape package (Paradis) and code provided by Andy Purvis, David Orme and Rich Grenyer (Pendek package in development, available upon request from ). The independent contrasts were standardized using branch lengths (taken from Bininda-Emonds et al. 2007), which were then transformed on a trait-by-trait basis following Garland et al. (1992). Polytomies in the phylogeny were treated as soft; they were resolved arbitrarily and the resulting n contrasts from each down-weighted by 1/n, with each node thereby contributing one degree of freedom (following Purvis & Garland 1993).
We used a three-stage process to identify correlates of extinction risk in all three data partitions (complete artiodactyl dataset, hunted partition and non-hunted partition). The first stage involved regressing IUCN threat rating against each continuous predictor variable; although threat is a discrete variable, a continuous distribution underlies the categories (Purvis et al. 2000). To test the effect of treating threat as a continuous variable, we ran analyses with IUCN threat category converted to a binary response variable (non-threatened, least concern and near threatened; threatened, near threatened conservation dependent through to critically endangered) with each life-history trait in turn as the Y variable. Trait values were not available for every species (electronic supplementary material), which meant that sometimes the degrees of freedom were quite low (10–20 degrees of freedom); however, no analysis was run without at least 10 degrees of freedom. Body mass was added as a covariate to all analyses that included weaning age, home range, population density, age at first birth, inter-birth interval and gestation length, as preliminary analysis confirmed that they covaried with body mass. Collinearity, however, should not be a problem as all correlations between traits and body mass indicated R-squared values of under 0.40. AET, as a measure of primary productivity, was added to all analyses including mean human population density, as preliminary analysis confirmed that it had the potential to be a confounding factor due to its significant association with mean human population density. Categorical data were analysed using the Wilcoxon signed-rank test on the contrasts. The second stage involved repeating the regression analyses including geographical range to allow comparisons to previous mammalian extinction risk analyses (Purvis et al. 2000; Jones et al. 2003). During each of these stages, the regressions were plotted, and contrasts that had undue influence over the regression line were deleted. Finally, the third stage involved building a multivariate regression model for each data partition that included all significant (p<0.05) and marginally significant (p<0.1) continuous traits from the first two stages. This regression model formed the starting set for model simplification to find the minimum adequate model (MAM), following the procedure outlined in Purvis et al. (2000).
Of the 144 artiodactyls that are categorized under criterion A in the IUCN Red List, 67 are currently threatened with extinction. Biological, ecological and anthropogenic correlates of extinction risk across all artiodactyls are presented in table 2. When single traits were regressed against threat status, small geographical range, older weaning age and low GNI were all significantly associated with a higher risk of extinction. When geographical range was added to the regression model, weaning age lost some significance and population density became marginally significant. The MAM (table 3) revealed weaning age, population density and GNI as the important predictors of extinction risk in the Artiodactyla, explaining 36.2% of the variance. Geographical range explains 6.9% but loses significance when placed in a multivariate model with weaning age (d.f.=60, t=−0.885, p=0.38). None of the three categorical traits showed a significant relationship with threat category.
When the effect of threat distribution was removed using the dichotomous threat variable, older weaning was confirmed as the most significant predictor of elevated threat (V=168, p=0.002). GNI was confirmed as an important predictor of threat (V=117, p=0.051); species that live in economically less developed countries are more threatened. Larger body mass and longer gestation length showed associations with higher threat (body mass V=299, p=0.079; gestation V=232, p=0.063), while geographical range (V=190.5, p=1) and population density (V=101, p=0.8) showed no significant predictive ability.
(b) Hunted and non-hunted species
Approximately 50% of the species in each partition are threatened with extinction. Traits that are associated with elevated threat levels in hunted artiodactyls are different from those that increase vulnerability to extinction in artiodactyls that are not known to experience uncontrolled hunting (table 4). In the first two stages of the analysis, weaning age and geographical range were the only significant predictors of threat status in hunted species, while age at first birth and population density were nearly significant (p<0.1). The MAM for the hunted partition contains only weaning age, which explains 19.7% of the variance in threat rating. By contrast, in the analysis of the non-hunted partition, geographical range, population density and GNI were significantly correlated with threat status, although population density lost significance when geographical range was added to the model. The MAM for the non-hunted partition contains both geographical range and GNI and explains 26.8% of the variance in threat rating. None of the categorical traits were significantly related to threat status in either partition (results not shown); the degrees of freedom are very low (d.f.<12).
Our analyses confirm the prediction that the influence of particular threatening processes on extinction risk depends on species-specific biological and ecological traits. Traits that were significant predictors of extinction risk across artiodactyls segregated among the hunted and non-hunted processes; geographical range was the only trait to remain significant regardless of threatening process. This result contrasts the comparative analyses of exploited primates and carnivores, which show the same set of extinction risk correlates across all species regardless of threat process (Purvis 2001).
Although geographical range was the only trait significantly associated with extinction risk in all dataset partitions, weaning age appears to be the key determinant of artiodactyl threat rating. Weaning age explains the greatest amount of variance in artiodactyl threat rating (R2=11%) and, unlike geographical range, weaning age is retained in the artiodactyl MAM and remains significant when threat is treated as a binary variable (threatened or non-threatened). This finding is not consistent with other studies in which geographical range is the single most important predictor of global mammalian extinction risk (carnivores and primates, Purvis et al. 2000; marsupials, Cardillo & Bromham 2001; bats, Jones et al. 2003). The greater importance placed on a reproductive trait in artiodactyls may relate to the fact that artiodactyls are primarily threatened by hunting (Mace & Balmford 2000), and weaning age is only associated with extinction in hunted artiodactyls. The mechanism for this result may be that hunted artiodactyls are more prone to extinction if they wean at older ages, a conclusion which is consistent with the prediction based on sustainable yield theory (Hartig 1796; Clark 1990). Species with older weaning ages have slower reproductive rates and, consequently, have lower sustainable off-takes, making them more prone to extinction via unsustainable harvesting. Our conclusion that the global extinction of hunted artiodactyls is related to slow reproductive rate is congruent with evidence from local extinctions of hunted mammals (Bodmer et al. 1997) and studies of bushmeat hunting, which show that slow-growing species are being hunted unsustainably (e.g. Barnes 2002; Bennett et al. 2002; Fa et al. 2003).
There is no evidence that hunter behaviour is playing a role in determining global artiodactyl extinction risk; hunter preference for large body size and hence slow life-history traits (e.g. Isaac & Cowlishaw 2004) cannot explain our results as body mass is not a significant predictor of threat. Our conclusion that present vulnerability to hunting-induced extinction is related to slow reproductive rates, not hunter preference for body size, is consistent with recent conclusions regarding the determinants of past exploitation-related extinctions (Johnson 2002). Extrapolation from extant relatives of species that went extinct during the Late Quaternary mammalian ‘megafaunal’ extinctions has shown that species with slow reproductive rates were more likely to become extinct regardless of body size (Johnson 2002).
GNI was a significant predictor of threat across all artiodactyls: species that live in areas of low economic development are more threatened. These findings agree with the results from a recent study (Davies et al. 2006) on the global distribution of extinction risk in birds, wherein areas of high economic development (as measured by GDP) are coincident with lower numbers of threatened species worldwide. When the dataset was partitioned, GNI was only associated with threat status in non-hunted artiodactyls, which contradicts our prediction that a weak economy will elevate the threat of extinction for all species regardless of threatening process. Although the threat status of hunted species is not associated with national economic status, risk of being hunted is associated with GNI: artiodactyls that experience uncontrolled hunting live in areas with significantly lower GNI than non-hunted species (Mann–Whitney test, W=5889, p=3.047×10−6). It is perhaps not surprising that species that experience unregulated hunting live in countries with low GNI because regulation (setting and enforcing quotas, etc.) of hunting is costly. Thus, the exclusion of artiodactyls that experience only regulated hunting may actually prevent us from seeing a relationship between threat status and GNI in the hunted partition as we have excluded the species that live in areas with higher GNI.
Species that are in the ‘non-hunted’ partition are likely to be primarily threatened by habitat loss, as it is the second most important threatening process in artiodactyls, accounting for 36% of all threats (www.redlist.org). Several hypotheses are commonly associated with extinction due to habitat loss: those reflecting ecological flexibility (Brown 1971; Laurance 1991; Norris & Harper 2004; but see Vazquez & Simberloff 2002) and those associated with small population size (Terborgh 1974; Simberloff 1986; Owens & Bennett 2000; Koh et al. 2004). Low population densities and small geographical ranges are associated with extinction risk in non-hunted artiodactyls, which is consistent with the small population hypothesis, but neither of the traits that reflect ecological flexibility, dietary specialization and habitat breadth, is associated with threat.
We conclude that different biological traits elevate vulnerability to extinction in artiodactyl species depending on whether a species is hunted. Correlates of extinction risk across all artiodactyls are a composite of the traits that increase vulnerability to different threats. Hunted artiodactyls with slower reproductive rates are more at risk of extinction, even though artiodactyls per se are less vulnerable than primates to extinction via hunting due to their relatively fast rates of reproduction (Bodmer et al. 1997). It is therefore important to know what type of threat a species is facing, particularly whether it is hunted or not, before identifying correlates of extinction risk. This study is an initial step in understanding how artiodactyls respond to anthropogenic extinction processes; the effects of habitat loss and hunting are often synergistic (Peres 2001). To improve the predictive ability of our extinction risk models, future studies must quantify how species respond to multiple extinction threats (Isaac & Cowlishaw 2004) and determine the spatial and/or temporal variation in the different threats (Fisher & Owens 2004).
We thank Andy Purvis, Rich Grenyer and David Orme for supplying the ‘R’ code for the phylogenetic independent contrast analysis, Hilmar Lapp for help with ‘R’, Georgina Mace, Rich Grenyer, Jonathan Davies and Jessica Partain for helpful discussion, Andy Purvis and two anonymous reviewers for comments and suggestions on earlier drafts of the manuscript, and Kate Jones for help with PanTheria and comments on a more recent draft of the manuscript. Financial support was received from the National Science Foundation (DEB/0129009) and latterly a NESCent post-doctoral fellowship to S.A.P. (NESCent NSF no. EF-0423641).