Synergistic effects of habitat fragmentation and hunting on the extinction risk of neotropical primates

Habitat fragmentation and overexploitation of natural resources are the most prevalent and severe threats to biodiversity in tropical forests. Several studies have estimated the effect of these threats on species extinction risk, however the effect resulting from their interaction remains poorly understood. Here, we assess whether and how habitat area, fragmentation, and hunting can synergistically affect the extinction risk of neotropical primates (Platyrrhine). We use a Random Forest model to estimate the Red List extinction risk category of 147 primate species based on their biological traits and the environmental predictors they are exposed to. We find that environmental variables are better predictors of extinction risk than biological traits, and that hunting and fragmentation interact creating synergistic feedback that lead to higher extinction risk than when considered in isolation. We also show that the effect of environmental predictors is mediated by biological traits, with large species being sensitive to habitat area and fragmentation, and frugivorous species more threatened by hunting. Our results increase the understanding of potentially interactive effects between different threats, habitat area and species traits, supporting the idea that multiple threats can reinforce each other and should be thus addressed simultaneously in conservation agendas.


Introduction
Habitat loss and overexploitation of natural resources are the most prevalent biodiversity threats globally (Maxwell et al. 2016), and are particularly severe in tropical forests (Lewis et al. 2015). In fact, tropical ecosystems are increasingly encroached by crops or livestock (Potapov et al. 2017), timber, energy production (Duden et al. 2020) and infrastructure (Laurance et al. 2009), resulting in intense deforestation. Furthermore, seemingly undisturbed forests may be under hunting pressure, with mammal populations reduced by Communicated by David Hawksworth.
Extended author information available on the last page of the article 1 3 more than 80% due to hunting (Benítez-López et al. 2017). These threats are intertwined and may act in synergy (Brook et al. 2008;Romero-Muñoz et al. 2020), as road building to access natural resources, deforestation and settlement expansion facilitate the access of hunters to forest fragments (Benítez-López et al. 2017, 2019. Most studies have quantified the effects of each threat in isolation, showing that their relative impact depends on species ecology. In fact, while larger species are generally more targeted by hunters (Benítez-López et al. 2017), forest specialist species are more vulnerable to land-use change due to deforestation (Galán-Acedo et al. 2019a, b) and generalist species are usually able to cope with disturbed habitat (Galán-Acedo et al. 2019b). In turn, species with low population density would be highly affected by hunting, in terms of high probability of extirpation due to low number of individuals. Recent studies have evaluated the effect of habitat loss and hunting in combination, either at regional (Symes et al. 2018;Romero-Muñoz et al. 2020) or pantropical scale (Gallego-Zamorano et al. 2020), estimating that tropical mammals have lost about 40% of their distribution due to the combined effects of habitat loss and hunting (Gallego-Zamorano et al. 2020). However, habitat fragmentation has not been formally considered in these studies (but see Peres 2001).
Habitat loss often results in habitat fragmentation, which can cause further decline in biodiversity (Crooks et al. 2017) reducing the abundance of populations making them sensitive to demographic and genetic stochasticity and local extinction, while hampering recolonization due to increased isolation (Haddad et al. 2015).Disentangling the effect of fragmentation from that of habitat loss presents some methodological challenges as these two processes are intimately linked. Studies estimating habitat fragmentation separately from habitat loss showed weak or positive effects on species abundance and richness (Fahrig 2019;Fahrig et al. 2019;Watling et al. 2020). However, these interpretations are currently debated on both empirical and theoretical grounds (Fletcher et al. 2018;Betts et al. 2019). For example, Saura (2020) argued that even assuming the habitat amount hypothesis holds (e.g. habitat fragmentation per se does not have negative effects on biodiversity, but species richness is positively related to the amount of surrounding habitat), changes in habitat configuration are expected to alter species richness (but see Fahrig 2021 andSaura 2021). Recent studies instead highlighted the long-term effect of fragmentation (i.e. extinction debt), therefore questioning the reliability of studies contrasting habitat loss and fragmentation effects within short timeframes (Semper-Pascual et al. 2021;Broekman et al. 2022).
Primates are among the mammal groups most affected by hunting (Ripple et al. 2016), and are particularly sensitive to deforestation which usually results in the fragmentation of continuous habitat in smaller patches (Estrada et al. 2017;Eppley et al. 2020). The loss of primate species is detrimental to the functioning of tropical forest ecosystems, due to their key roles as seed dispersers, folivores, and, in cases, as seed predators influencing the forest regeneration and diversity (Nuñez-Iturri & Howe 2007;Barnett et al. 2012;Bello et al. 2015;Estrada et al. 2017). Approximately 40% of neotropical primates (Platyrrhine) assessed by the International Union for Conservation of Nature (IUCN) are considered at risk of extinction (IUCN 2020), in particular, Lagothrix is the most threatened Platyrrhine genus due to hunting (Stafford et al. 2017), whereas Ateles is threatened by both hunting and deforestation (Stafford et al. 2017;Aquino et al. 2018).
Here, we present a quantitative analysis of the relative effects of area of habitat, habitat fragmentation, hunting pressure, species biological traits (body mass, population density, diet, gestation length and litter size), and their interactions, on the extinction risk of Neotropical primates. We expect extinction risk to increase with decreasing habitat area (Betts et al. 2017), increasing fragmentation (Crooks et al. 2017) and hunting pressure (Estrada et al. 2017), and that the effects of these factors are mediated by species traits. Specifically, we expect large and slow-reproducing species to be more sensitive to hunting (Ripple et al. 2016;Benítez-López et al. 2017, 2019, and species with low population density to be more sensitive to amount of habitat area (Sykes et al. 2020), fragmentation (Eppley et al. 2020) and hunting. We also expect a positive interaction between environmental predictors yielding a multiplicative effect on extinction risk (Romero-Muñoz et al. 2020).

Species and biological traits selection
We extracted the geographic ranges and conservation status of all neotropical primate species from the IUCN Red List database in 2019. We excluded all data deficient (DD) neotropical primates for which threat status is unknown. We only selected species that were not assessed under criterion B of the IUCN Red List, which concerns the extent of occurrence (B1) and the area of occupancy (B2) which may be related with the area of habitat used as model predictor, therefore introducing circularity. All the species we selected were assessed under the Red List criteria A and C, which are based on species population trends and are not directly calculated using species spatial information. None of the species was assessed under criteria D and E.
We used biological traits to account for characteristics that can make species vulnerable to threats (Purvis et al. 2000;Cardillo et al. 2008;Davidson et al. 2009). We excluded biological variables with more than 50% of missing data (see below), therefore selected the following biological traits: body mass, a proxy of vulnerability to hunting, as larger tropical mammals were predicted to be more hunted (Benítez-López et al. 2017, 2019; species' average population density, a measure of population level space use and vulnerability to fragmentation, since species with large spatial requirements at the population level are more sensitive to fragmentation (Santini et al. 2018a, b;Eppley et al. 2020); and percentage of frugivory in the diet (hereafter, frugivory), under the assumption that species with more specialized diet generally need more space to find appropriate resources, being thus more vulnerable to fragmentation (Eppley et al. 2020). We also selected two reproductive traits as measures of population recovery potential: gestation length and litter size, which represent reproductive timing and output, respectively (Bielby et al. 2007). Species traits were extracted from different databases: PanTHERIA (Jones et al. 2009 The biological traits that we selected for our analysis had different proportions of missing values-frugivory diet 4%, body mass 12%, population density 32%, litter size 44%, gestation length 45%. Therefore, we imputed the data following Penone et al. (2014) using the "mice" package (Van Buuren & Groothuis-Oudshoorn 2011) in R and the phylogenetic eigenvectors. We used phylogenetic eigenvectors to account for latent traits and phylogenetic relatedness (Diniz-Filho et al. 1998). We obtained the phylogenetic tree from the PHYLACINE database (Faurby et al. 2018). The source 1 3 phylogeny was derived using a hierarchical Bayesian approach with a posterior distribution of 1,000 trees. We extracted 10 random phylogenetic trees from the phylogeny and extracted 20 eigenvectors from each tree, which we used to test the sensitivity of our imputation to phylogenetic uncertainty. We repeated the imputation using 5, 10 and 20 phylogenetic eigenvectors. We produced 10 different imputed datasets each iteration with the predictive mean matching method. To assess the imputation performance, we calculated the Normalized Root Mean Squared Error (NRMSE), which ranges from 0 to 1, for each variable. Lower values of NRMSE indicate better estimates of the variables.

Area of Habitat and Fragmentation
We computed the amount of area of habitat within the IUCN Red List species ranges (Brooks et al. 2019), to represent the amount of habitat available and potentially occupied by species within their ranges. We estimated the area of habitat by combining the raster layers of forest cover and loss from Hansen et al. (2013), the digital elevation model from Jarvis et al. (2008), and the polygons of the species geographic ranges available from the IUCN database (IUCN 2020). First, as forest layers are expressed as percentage of canopy cover, we binarized these layers into forest and non-forest layers. Although some neotropical primates can perform specific activities at the ground level (Mourthé et al. 2007;Souza-Alves et al. 2019;Eppley et al. 2022), all Plathrrine monkeys are strictly arboreal, thus we followed previous studies and considered only areas where the tree cover was > 75% (Aleman et al. 2018;Vieilledent et al. 2018;Eppley et al. 2020). Then, we overlaid the geographic range areas with the binary forest maps to calculate the amount of forest habitat for each species. Subsequently, we excluded the portions of forest habitat outside the species altitudinal range of presence (Tracewski et al. 2016). The area of habitat was measured at 90 m resolution, matching the coarsest resolution of the two raster layers (i.e. 30 m forest coverage, 90 m elevation model). The values of the area of habitat were reported in km 2 .
We calculated two common fragmentation metrics from the area of habitat maps to estimate the degree of the forest fragmentation: the mean patch area in km 2 (Innes & Koch 1998) and the distance from forest edge in km (Crooks et al. 2017), measured as the average Euclidean distance of all cells within a species area of suitable habitat from the nearest non-forest edge. The former takes into account the area of continuous habitat and thus is a proxy of the potential size of isolated or semi-isolated population, while the latter is a measure of habitat degradation in terms of habitat alteration due to edge effects (Pfeifer et al. 2017). Low values of mean patch area and low values of distance from forest edge indicate a more fragmented habitat.

Defaunation index
As a proxy of hunting pressure, we used the hunting-induced defaunation index from Benítez-López et al. (2019). This index is derived from a model predicting the hunting pressure on tropical mammal species. The defaunation index ranges from 0 (no defaunation) to 1 (local extirpation). We extracted the average defaunation index within the species available habitat in the distribution range, which represents how much a species population size is reduced by hunting pressure.

Modeling
We used the IUCN Red List category of species as the response variable representing extinction risk in our model. Following previous studies (Purvis et al. 2000;Cardillo et al. 2004;Polaina et al. 2016), we converted this categorical variable to numeric assigning a value to each extinction risk category: Least Concern = 1, Near Threatened = 2, Vulnerable = 3, Endangered = 4, Critically Endangered = 5. We also transformed the area of habitat, mean patch area, distance from edge and population density using natural logarithm for graphical purposes.
To quantify the relationship between the predictors and the species extinction risk we used a Random Forest model. Random Forest is a machine learning approach that has been successfully used in other ecological and conservation analyses (Di Marco et al. 2015;Pacifici et al. 2020). These models are more flexible than statistical linear models, are robust to collinearity and structured data, and allow to estimate complex non-linear relationships and interactions among predictors (Cutler et al. 2007). Also, the model does not make any assumption on the distribution of the response variable, which is generally problematic in linear models such as phylogenetic least square models (Lucas et al. 2019;Cazalis et al. 2022), i.e. ordinal distribution with non-even increase in extinction risk between Red List categories. We estimated the relative importance of each variable in predicting extinction risk category by measuring the relative increase in the mean square error of the Random Forest model when the values of the variables are randomly permuted (Cutler et al. 2007).
To estimate the accuracy of the model, we performed a 5-repeated tenfold cross validation with an 80-20% training-testing sets and assessing each model using the Root Mean Square Error (RMSE) at each iteration. We selected the model with the lowest RMSE, which was the most accurate under cross-validation, as our final model. Then, we assessed the performance of the final model using the percentage of variance explained. We repeated the Random Forest model across all imputed datasets in order to assess how uncertainty in the imputation procedure could influence our conclusions (Conenna et al. 2021). To quantify the interaction between environmental predictors and traits, we created partial response plots between the predicted extinction risk as a function of an environmental predictor (area of habitat, mean patch area, distance from edge or defaunation index), while modifying the value of the other interacting environmental predictor or biological trait. The rest of predictors were maintained at average values. For example, to assess the interaction between the defaunation index and the area of habitat, we plotted extinction risk as a function of increasing defaunation, and for area of habitat at high (95 th percentile) and low (5 th percentile) values. To better show the interaction effects, we generated a plot displaying the distribution of the distances between two response curves at the two extremes of the range across the ten imputed datasets (Fig. A.1). The median of the distribution indicates the effect size of the interactions, and distributions not overlapping with zero are an indication of consistent directional interaction effects.

Results
The best imputation procedure included 10 phylogenetic eigenvectors in addition to life history trait variables, resulting in an average NRMSE across all imputed trait variables ranging from 0.08 to 0.24, which indicates good imputation performance (Table A.1). The Random Forest model was able to explain 44.5% (± 2.49%, 95% confidence interval) of the variance. Overall, the environmental predictors were more important than biological predictors in our model, with all environmental variables among the top four important variables. The most important variables in our model were the area of habitat and the defaunation index, followed by the fragmentation variables and body mass (Fig. 1).
The area of habitat and the defaunation index yielded, respectively, a strong negative and strong positive effect on extinction risk (Fig. 2a,b). Thus, extinction risk increased as the amount of habitat decreased, and with increasing hunting-induced defaunation. The two variables of fragmentation also showed a negative influence on extinction risk, with low values of mean patch area and low values of distance from the edge associated with higher extinction risk (Fig. 2c,d). The average extinction risk was higher for primates with large body mass, high gestation length and low population density, while species with higher litter size showed a slight increase in extinction risk (Fig. 2e-h). The percentage of frugivory showed a slightly non-linear trend, with the highest risk of extinction for species consuming very low, or very high percentage of fruit in their diets (Fig. 2i).
The interactions between the amount of area of habitat, fragmentation and hunting variables showed synergistic effects on species extinction risk (Fig. 3). The effect of hunting was stronger with decreasing average distance from edge (Figs. 3a,4) and with smaller values of mean patch areas (Figs. 3b,4), indicating that species living in fragmented areas are also more vulnerable to hunting. The effect of mean patch area was stronger when the area of habitat was low, indicating a stronger effect of fragmentation in species distributed in small areas (Figs. 3d,e,4). Hunting pressure and area of habitat did not show a clear interaction (Figs. 3c,4).
Area of habitat, fragmentation and hunting variables also exhibited variable interaction effects with biological traits (Fig. 4). For example, the extinction risk of large-sized species Fig. 1 Variable importance in the Random Forest model predicting species extinction risk. Variable importance, with 95% confidence intervals, is represented as the increase in Mean Square Error (MSE) associated with random permutation of its values. Red bars indicate environmental predictors, blue bars indicate species biological traits decreased as the mean patch area and distance to edge increased (indicating less fragmentation), but at a lower rate than in the case of small-sized species (Fig. A.2), suggesting higher sensitivity to fragmentation in large species. Similarly, large species showed a stronger positive effect of area of habitat. The effect of area of habitat was more important Fig. 2 Partial responses plots displaying the relationship between extinction risk and each predictor, while all the other variables are kept to their average value. Shaded colours represent 95% confidence intervals. Environmental predictors are represented in red and biological predictors in blue Fig. 3 Interaction between area of habitat and threat variables. The plots show the response of one variable at two fixed values of the interacting variable, corresponding to its 5th percentile (green) and 95th percentile (purple). The shaded areas represent the 95% confidence intervals of the responses calculated across the responses of the ten imputed datasets. Different slopes in the two curves represent the interaction between the two variables 1 3 for species with long gestation length (Fig. A.3) and living at low density (Fig. A.5). Further, the defaunation index had a stronger negative effect on more frugivorous species (Fig.  A.6). Finally, frugivory and species density were also associated with higher sensitivity to distance from edge (Fig. A.5, A.6). Other traits did not show clear interactions with any of the environmental predictors (Fig. A.4, A.5, A.6).

Discussion
In this study we estimated the effect of the area of habitat, fragmentation, and hunting on the extinction risk of Neotropical primates, while considering differences in their biological traits. We found that reduced area of habitat, increased hunting pressure and increased habitat fragmentation were associated with higher extinction risk. Additionally, larger, slow-reproducing, and low-density species were more threatened on average. Furthermore, larger species appeared to respond more negatively to hunting and fragmentation, and be more sensitive to a lower area of habitat than small-sized species. Finally, we report clear Fig. 4 Size of the interaction effects between pairs of variables on the extinction risk of neotropical primates. The effect size of the interactions is calculated by measuring the delta in extinction risk prediction distances between two response curves at the two extremes of their distribution (see Fig. S1). This procedure is repeated across the 10 imputed datasets to estimate a distribution of effect sizes. Distributions whose at least 90% of the values were above or below 0 were represented with labels in bold. Interactions between environmental predictors are represented in red and interactions between an environmental predictor and a biological predictor in blue 1 3 interactions between fragmentation and hunting pressure, and between area of habitat and fragmentation, suggesting their synergistic effect in determining species extinction risk.
Overall, threats and the extent of the area of habitat and body mass were the most important drivers of extinction risk, with other biological predictors showing a minor contribution to extinction risk. This result supports the idea that extrinsic variables (environmental or threat variables) should be included in comparative risk analyses in addition to biological predictors (Murray et al., 2014;Di Marco et al. 2018). This, however, ultimately depends on the group under investigation, the variability in biological traits and knowledge of significant external pressures on the group of species under study. Nonetheless, our study supports previous comparative extinction risk analyses in concluding that larger, slow reproducing, and low-density species are more threatened on average (Cardillo et al. 2005;Fritz et al. 2009;Hilbers et al. 2016). It also supports previous studies on the role of traits in mediating the effects of fragmentation and hunting, showing that larger species are more sensitive to habitat area and fragmentation (Crooks et al. 2017;Ripple et al. 2017), and hunting (Redford 1992;Ripple et al. 2016;Benítez-López et al. 2019). Matching the results of previous studies (Eppley et al. 2020;Sykes et al. 2020), we found that species living at low population density showed higher sensitivity to habitat area and fragmentation. Frugivory showed a weak non-linear effect, with species at a slightly higher risk of extinction being at the extremes of the frugivory continuum. The species considered in this study were mostly folivores or frugivores, so species at the extreme of the frugivory continuum showed a more specialized diet, whereas species at the center had a more diverse diet. Frugivores also exhibited a slightly higher sensitivity to distance from edge. Indeed, frugivorous species require larger areas of continuous habitat to forage as their trophic resources are more sparse and clumped in space compared to those of folivorous and omnivorous species (Milton & May 1976), resulting in a reduced tolerance to fragmented habitats. Our result may thus suggest that species with more diverse diets can better cope with degraded habitat. Finally, we found that frugivorous species were more threatened by hunting than omnivores and folivores, probably because many highly frugivorous primates in the neotropics are large-bodied species that are heavily hunted (e.g. Ateles, Lagothrix, Alouatta).
Several recent studies have questioned the individual role of habitat fragmentation in determining species risk of extinction. Some studies have argued that habitat fragmentation generally has no, or even positive, effects on species richness and abundance (Fahrig 2019;Fahrig et al. 2019), while others support the traditional view that habitat configuration, irrespective of habitat amount, yields negative effects on species persistence (Haddad et al. 2015;Fletcher et al. 2018;Saura 2020). In this work we show that fragmentation, either expressed as mean patch area or distance from edge, is positively related to the extinction risk of neotropical primates, and these measures interact negatively with the amount of habitat area, reinforcing the idea that habitat configuration, not just amount, is a key parameter for species persistence (Ramírez-Delgado et al. 2022) particularly in the case of habitat specialists with limited mobility across non-habitat areas such as arboreal primates.
Our results should be interpreted considering some limitations. First, in this study we measured fragmentation using a binary forest layer, thresholding the forest layer at 75% of canopy cover following previous papers on tropical forests (Aleman et al. 2018;Vieilledent et al. 2018;Eppley et al. 2020). This process resulted in a simplification of the landscape, as species can experience a wide variety of canopy covers within their range and the lower boundary likely differs among species. Second, the response variable in our model, the IUCN Red List extinction risk category, reflects a coarse categorization of species extinction risk (but see Mooers et al. 2008). This likely reduces our ability to capture nuances in the data and is expected to flatten the estimated effects of drivers. Finally, we worked under the assumption that these categories were correctly assigned. It is possible, however, that Red List categories are incorrectly assigned because of uncertainty of data at the time of the assessment ). This may have introduced noise in the model, however it is unlikely to have biased the results as we only focused on one mammalian taxon in one region, which is assessed consistently by the same group of experts during collective workshops (IUCN Primate Specialist Group, Neotropic section).
Overall our study supports the role of habitat area, fragmentation and hunting as important drivers of primate extinction risk (Estrada et al. 2017). However, it also shows that these factors do not act consistently across all species in our sample, but rather their impact may be exacerbated or mitigated by species traits, such as body mass and trophic ecology (Ripple et al. 2015Benítez-López et al. 2017, 2019Eppley et al. 2020;Sykes et al. 2020). Finally, it supports the idea that the amount of area of habitat and threats like fragmentation and hunting are self-reinforcing and may lead to more severe reductions in mammal distributions and their extinction risk than anticipated (Brook et al. 2008;Gallego-Zamorano et al. 2020;Romero-Muñoz et al. 2020). These results suggest that conservation efforts aimed at restoring habitat and connectivity for large species can also indirectly reduce the impact of hunting, provided that human accessibility to restored areas is restricted to local communities, and that these should be actively involved in the design, implementation, monitoring, and evaluation of their hunting activities to ensure sustainability. However, considering the reported interactions among threats, addressing threats in isolation would unlikely result in effective conservation outcomes. Effective conservation of threatened primates will require adopting a more holistic approach.