Increasing synergistic effects of habitat destruction and hunting on mammals over three decades in the Gran Chaco

1 –––––––––––––––––––––––––––––––––––––––– © 2020 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Mikko Mönkkönen Editor-in-Chief: Miguel Araújo Accepted 17 March 2020 43: 1–13, 2020 doi: 10.1111/ecog.05053 43 1–13


Introduction
Habitat destruction and overexploitation are the two main drivers of the unfolding sixth mass extinction, and both threats continue to expand (IPBES 2019). On one hand, growing demands for food, livestock feed and biofuels trigger widespread land-use changes, including agricultural expansion into remaining natural habitats in the Global South (Kehoe et al. 2017). On the other hand, overexploitation (i.e. the unsustainable hunting, collection of animals and plants, logging or fishing) (IPBES 2019), expands rapidly as the global human population grows, affluence increases and demand for wild animals and plants (e.g. meat, live specimens) increases (Benítez-López et al. 2017. Therefore, understanding the extent of these threats and how they change over time is critically important to inform conservation actions (Wilson et al. 2005, Pressey et al. 2007.
Habitat destruction and overexploitation may synergise where they act simultaneously, exacerbating their individual impacts on biodiversity (Brook et al. 2008). For instance, deforestation increases hunter access to shrinking habitat and formerly remote areas (Peres 2001) as does road infrastructure development related to expanding agriculture (Laurance et al. 2014). Yet, despite these synergistic effects, the interactions among habitat destruction and overexploitation remain weakly understood, and most studies in conservation and ecology continue to study threats in isolation (Brook et al. 2008, Dirzo et al. 2014, because approaches and datasets to jointly study multiple threats are lacking (Wilson et al. 2005, Pressey et al. 2007, Joppa et al. 2016.
Assessing the spatial footprint of threats to biodiversity, how these footprints overlap, where they remain absent, and how they change over time -which we here collectively refer to as the 'geographies of threat' -can help understand the individual and combined effects of those threats. Understanding geographies of threat is also imperative for guiding conservation planning by identifying where threat-specific conservation actions should take place (Wilson et al. 2005, Pressey et al. 2007). However, mapping the geographies of threat is challenging. Few studies have mapped multiple threats at broad scales; typically within the scope of single-species studies (Bleyhl et al. 2015, Romero-Muñoz et al. 2019b), which has limited value for conservation planning that targets wider biodiversity facets (Nicholson and Possingham 2006). Studies assessing broader groups of species on the other hand, usually rely on IUCN's expert-based threat categorizations and range maps, thereby assuming that threats impact multiple species uniformly (Allan et al. 2019, Gallego-Zamorano et al. 2020, which is too simplistic. In addition, expert-based range maps contain false presences and vary tremendously in quality, depending on regions and taxa, and their use is therefore limited to very coarse resolutions (Ficetola et al. 2014, Di Marco et al. 2017. Consequently, these approaches are insufficient to inform threat-specific management actions on the ground (Wilson et al. 2005, Tulloch et al. 2015. New approaches to map the species-specific responses to threats for multiple species simultaneously and at resolutions useful for practitioners are urgently needed (Wilson et al. 2005, Pressey et al. 2007, Tulloch et al. 2015. Recent advances in remote sensing now allow the reconstruction of detailed land-change histories across several decades and large areas (Hansen et al. 2013, Baumann et al. 2017, Song et al. 2018). This provides opportunities for assessing habitat change dynamically, but few studies to date have made use of these opportunities (Maguire et al. 2015, Oeser et al. 2019, Romero-Muñoz et al. 2019b). Likewise, new approaches for assessing the impact of hunting in space are developed (Benítez-López et al. 2017). Such hunting pressure models synthesise knowledge across local studies, to describe how species-specific responses to hunting vary across landscapes (Benítez-López et al. 2019). Here, we propose to combine habitat suitability and hunting pressure models for characterising the footprints of habitat destruction and hunting, and how they overlap.
Understanding of the interacting footprints of habitat destruction and hunting is particularly poor in tropical deforestation frontiers, where rapid habitat destruction often couples with high hunting pressure (Peres 2001, Benítez-López et al. 2019). This situation is particularly dire in the world's tropical dry forests, which are vanishing quickly across the globe (Hoekstra et al. 2005, Curtis et al. 2018). However, these systems remain weakly protected (Miles et al. 2006) and under-researched (Blackie et al. 2014). The individual and combined impacts of habitat destruction and hunting on biodiversity in these forests are highly unclear, translating into a real barrier towards implementing conservation planning and action.
At 1.1 million km 2 , the Gran Chaco region (hereafter Chaco) in South America, extending into parts of Argentina, Paraguay and Bolivia, is the largest continuous tropical and subtropical dry forest globally, but it has recently turned into a global deforestation hotspot due to rapid agricultural expansion across the several deforestation frontiers that it encompasses (Baumann et al. 2017, Curtis et al. 2018, Le Polain de Waroux et al. 2018. At the same time, hunting is very widespread there, causing massive defaunation (Noss et al. 2005, Altrichter 2006, Periago et al. 2014). Together, these threats render the Chaco a global conservation priority (WWF 2015, Kuemmerle et al. 2017. Increasing evidence suggests important interactions between habitat destruction and hunting in this region. For example, large mammals disappear from remaining forest patches soon after the surrounding areas are deforested because they are easily hunted out (Núñez-Regueiro et al. 2015, Semper-Pascual et al. 2019. Likewise, cattle ranchers in areas where pastures expand often persecute large predators over fears of attacks on cattle (Quiroga et al. 2016, Romero-Muñoz et al. 2019b). Yet, our understanding of how these threats play out and interact in space is very limited.
Here, we reconstruct the individual and combined spatial footprints of habitat destruction and hunting pressure for larger mammals (> 1 kg body weight) across the entire Chaco between 1985 and 2015. We combine satellite-based landuse reconstructions with species-specific, time-calibrated habitat suitability models and hunting pressure models. This allows to assess the footprints of habitat destruction and hunting and to identify threat hotspots as well as how they change over time. Specifically, we aimed to answer the following questions: 1) how have the footprints of habitat destruction and hunting pressure on larger mammals changed in the Chaco since 1985? 2) What is the relative importance, in terms of the share of species' ranges affected and their overall footprints, of these two threats acting alone versus together, and how this has changed over time? 3) How has the distribution of core areas, where threats are absent, changed since 1985, and where are current hotspots of threats and priority areas for conservation?

Study region
The Chaco region is a highly biodiverse region comprising parts of Bolivia, Paraguay and Argentina (Olson et al. 2001, TNC et al. 2005. The climate ranges from tropical (north) to subtropical (south). Precipitation is seasonal and ranges from > 1200 mm yr −1 (east) to < 400 mm yr −1 (west and south). Xerophilous forests are the dominant vegetation, interspersed with gallery forests and savannas (Prado 1993). The Chaco has a long land-use history, with Indigenous Peoples using the area for millennia, and criollo people practicing subsistence ranching for up to 200 years (Camino et al. 2018). Recent expansion of intensified agriculture, mainly driven by largescale, market-oriented agribusiness, has converted more than 142 000 km 2 of forests (> 20% of the Chaco's forests) to pastures and croplands between 1985 and 2015 (Baumann et al. 2017). Hunting is also widespread (see Extended methods in Supplementary material Appendix 1), with many actors hunting for subsistence, commercial, cultural and retaliatory reasons, together producing widespread defaunation across the Chaco (Periago et al. 2014, Torres et al. 2014, Semper-Pascual et al. 2018. Only 9.1% of the Chaco is currently protected (Nori et al. 2016).

Data preparation
We gathered 27 408 presence locations from local surveys and opportunistic observations for 56 larger terrestrial mammals. These records were collected from 1978 to 2018, partly by the authors, and from public (e.g. GBIF), and governmental and non-governmental organisations' databases (Supplementary material Appendix 1 Table A1 for details). To reduce sampling bias, we spatially filtered presence locations by enforcing a minimum distance of 10 km between presence locations (Kramer-Schadt et al. 2013). We only included species with more than 10 points after applying the spatial filtering, resulting in a final list of 48 species, for which we retained a total of 4611 presence locations.
As potential predictors for our habitat suitability and hunting pressure models, we generated 11 variables at a 1-km 2 resolution (Supplementary material Appendix 1 Table A2). All variables covered the entire Chaco plus a 30-km buffer to account for potential border effects (Piquer-Rodríguez et al. 2015). For the habitat model, we included four variables characterizing land cover (%Forest, %Cropland, %Grassland, %Pastures), three variables describing habitat structure (%Forest edge, Distance to water) and two climate variables (mean annual temperature, mean annual precipitation; Supplementary material Appendix 1 Table A2). To assess collinearity among predictors, we calculated Pearson's correlation coefficients for each variable pair and kept the variable with the higher explanatory power for pairs with r > 0.7 (Dormann et al. 2013) (Supplementary material Appendix 1 Fig. A1).
For the hunting-pressure model, we followed Benítez-López et al. (2019) and used three predictors: Distance to hunter access points, Human population density (both indicators of hunting risk), and Species body mass (an indicator of a species' intrinsic vulnerability to population decline as a result of hunting) (Supplementary material Appendix 1 Table A2). We defined spatial features representing hunter access points for each species separately, based on the regional expertise of the authors (Supplementary material Appendix 1 Table A3). Assessment of subsistence ranches involved screen-digitizing > 27 000 individual ranches spread across the Chaco and assessing their persistence over time using high-resolution imagery in Google Earth. Likewise, we reconstructed the evolution of the road network since 1985 based on historic satellite imagery (see Extended methods in the Supplementary material Appendix 1 for details).

Modelling habitat destruction and hunting pressure over time
We parameterized 1) a habitat suitability model, characterizing resource availability, and 2) a hunting pressure model, characterizing species-specific population declines due to hunting. By overlaying the two resulting maps, we then 4 identified four habitat categories for each species individually, according to the prevailing threats: core area (good habitat suitability and low hunting pressure), poor habitat-only (poor habitat suitability, but low hunting pressure), hunting pressure-only (high hunting pressure, but good habitat suitability), and both threats together (poor habitat suitability and high hunting pressure). We tracked these habitat categories across time using time-calibrated models for each species, resulting in time series of the individual and combined threat footprints (Fig. 1).
To represent habitat suitability, we used maximum entropy modelling (Phillips et al. 2017). This is a presence-only, nonparametric species distribution modelling algorithm that performs well in predicting habitat suitability, even for small samples (Elith and Leathwick 2009) and for time-calibrated habitat models (Kuemmerle et al. 2012, Sieber et al. 2015, Romero-Muñoz et al. 2019b). Time-calibrated models have two key advantages: 1) they make use of all available data, across the entire time period studied, and 2) they ensure that observed changes in habitat suitability are solely due Figure 1. Framework for reconstructing 'geographies of threat' due to habitat destruction and hunting pressure for 48 larger mammals in the Chaco from 1985 to 2015. We first modelled the spatial footprint of each threat per species, then stacked these footprints across the community, and then used this information to assess how spatial footprints of threats changed over time (including threat overlaps).
to changes in predictor variables, and not due to uneven distribution of points over time or varying sampling bias (Nogués-Bravo 2009, Sieber et al. 2015. We fitted maximum entropy models for each species using Maxent (v3.4.1) (Phillips et al. 2017) using only hinge features to avoid overfitting (Elith et al. 2010). We tested a range of parameterizations and selected a regularisation multiplier of 1 and a prevalence value of 0.5 (Elith et al. 2010).
As background points, we created sets of points for each species individually to account for differences in species' distribution as well as sampling effort in space and time, which helps to avoid issues arising from sampling biases (Elith et al. 2010, Merow et al. 2013. We used 10 000 background points that we distributed proportionally in time according to the presence points. We then extracted predictor values for each presence and background point from the year each point was sampled (Sieber et al. 2015). This yielded a single, timecalibrated Maxent model per species, which we then projected onto the sets of predictor variables from 1985, 2000 and 2015 (see Extended methods in the Supplementary material Appendix 1 for further details). To assess the robustness of our models, we ran 10-fold cross-validation. We assessed the models' predictive performance with the average area under the curve (AUC) values across the 10 replicates. We defined species' ranges as those areas with habitat suitability values above the 5% quantile in 1985 (Pearson et al. 2004).
To model hunting pressure, we relied on a recently-developed approach to capture hunting-induced defaunation for tropical mammals (Benítez-López et al. 2019). This approach uses a two-stage mixed model that describes a species' population responses to hunting pressure. First, a binomial model was fitted to discriminate extant and locally extinct species. Second, a Gaussian model was fitted to the non-zero response ratios in abundance change due to hunting based on 3281 abundance estimates in hunted and non-hunted areas studies across the tropics (Benítez-López et al. 2019) (see Extended methods). This results in a hunting pressure index ranging from 0 (no decline in abundance) to 1 (total local extirpation). We re-fitted the original global model to Neotropical mammals only (n = 1945 abundance ratios) and then evaluated the predictive accuracy with 5-fold cross-validation with an 80%/20% training/testing set. We split our predictions into two categories of high (> 0.3), and low (≤ 0.3) hunting pressure, based on the International Union for Conservation of Nature (IUCN) criterion of 30% population decline due to threats that have not ceased that renders a species threatened (criterion A4, IUCN 2012). We assessed the accuracy of our model for predicting these hunting pressure categories using sensitivity and specificity.

Mapping the footprints of habitat destruction and hunting pressure
To map the spatial footprints of threats, we first applied thresholds to the habitat suitability maps and the hunting pressure maps to classify good and poor suitability, and high and low hunting pressure (Fig. 1), respectively. For the habitat suitability maps, we used the 'maximum sensitivity plus specificity' threshold (Liu et al. 2013). For the hunting pressure maps, we used the threshold of 0.3 to separate high and low hunting pressure. We then overlaid the two binary maps to identify the four habitat categories according to threat levels (Fig. 1).
We stacked the raster maps across all 48 species to obtain per-pixel species counts for each category for the years 1985, 2000 and 2015 (Fig. 1). We also calculated for each year the overall area affected by poor habitat and hunting pressure, and the share of each species' ranges affected by either threat alone or by both together. In the habitat model, we kept climate conditions constant for the entire study period (by using 30-yr climate averages) but allowed land cover and land use to vary over time. Therefore, expansion of poor habitat over time can only be attributed to impacts of land cover/use change and we refer to this as habitat destruction (Fig. 1). We refer to the increases of hunting pressure over time as 'increasing hunting pressure'.
To identify priority areas (i.e. the most important areas with high-quality habitat and low threat levels) and hotspots of threats (i.e. areas where threats have disproportionally high impacts), we adopted a rarity-weighted richness measure (Kier and Barthlott 2001), which considers both richness (i.e. how many species have their core area in a given cell) and range size (i.e. whether a species has a large or small core area). This approach compares favourably to other prioritisation algorithms (Albuquerque and Beier 2015). Priority areas can guide efforts to expand habitat protection (e.g. via additional protected areas), while threat hotspots can help to spatially target threat-specific conservation action (see Extended methods in the Supplementary material Appendix 1 for further details).

Results
Both our habitat suitability models, and hunting-pressure models performed well. Our habitat suitability models had overall high to very high model fit and discrimination values for all 48 modelled species (AUC consistently > 0.7; Supplementary material Appendix 1 Fig. A2). For the hunting pressure model, overall sensitivity and specificity were 0.9 and 0.5, respectively, indicating good predictive performance.
In terms of the predicted threat footprints, our habitat suitability models showed that by 2015, poor habitat covered on average 49% (± 20% SD) of the ranges of the species we investigated (Fig. 2A). Similarly, hunting pressure was on average high across 45% (± 30% SD) of species' ranges in 2015 (Fig. 2C). Between 1985 and 2015, large areas of the Chaco became affected by habitat destruction and hunting pressure (38 and 41% of the region, respectively; Fig. 2B, D). For some species, hunting pressure expanded over wide areas and even inside protected areas (Fig. 2D).
At the species level, the footprint of habitat destruction showed an average expansion of 9.6% (± 22.7% SD) or 22 000 km 2 (± 51 000 km 2 SD; Fig. 3A). This threat increased for 34 mammals (71%), while it either remained constant or decreased for the remaining 14 species (Fig. 3A). For example, since 1985 land-use change affected over 25% of the high-quality habitat of the jaguar Panthera onca, puma Puma concolor, the white-lipped peccary Tayassu pecari, and the collared peccary Pecari tajacu. In contrast, species such as the maned wolf Chrysocyon brachyurus or the crab-eating fox Cerdocyon thous experienced declining pressure from habitat destruction over time (Fig. 3A). Among countries, the footprint of habitat destruction expanded faster in Paraguay than in Bolivia and Argentina (Fig. 3B).
The footprint of hunting pressure expanded on average by 8.4% (± 6.7% SD) or 23 000 km 2 (± 34 000 km 2 SD; Fig. 3A). Generally, this footprint changed more evenly than the footprint of habitat destruction, with increasing hunting pressure for almost all species (i.e. 44 species = 92%). For instance, hunting pressure on the puma, the jaguar, the giant armadillo Priodontes maximus, and the grey brocket deer Mazama gouazoubira each increased by more than 20%. For some frequently-hunted species, such as the white-lipped peccary and tapir Tapirus terrestris, the footprint of hunting pressure increased only slightly, as this footprint was already large in 1985. Only very few species, such as Geoffroy's cat Leopardus geoffroyi) experienced slightly shrinking hunting pressure (Fig. 3A). The footprint of hunting pressure expanded faster in Paraguay and Bolivia than in Argentina (Fig. 3B).
In addition to the individual expansion of threat footprints, we found a strong increase between 1985 and 2015 in the area where habitat destruction and hunting pressure overlap (Fig. 4). The cumulative area for all mammals affected by both threats expanded by 465 000 km 2 (or 43% of the Chaco) between 1985 and 2015 (Fig. 4). In comparison, habitat destruction-only and hunting pressure-only cumulatively expanded by 300 000 km 2 and 363 000 km 2 (34% and 28% of the Chaco), respectively (Fig. 4). At the species level, the area of both threats acting together increased by 17% (± 20.2%) on average between 1985 to 2015. In contrast, the area where only one threat impacts species decreased (by 39.5% and 6.1%, for habitat destruction and hunting pressure, respectively; Fig. 5).
Regarding core areas (i.e. good habitat suitability and low hunting pressure), 36 species (75%) experienced a contraction (on average 38% ± 62.2% SD) between 1985 and 2015 (Supplementary material Appendix 1 Fig. A3). Contractions were particularly common in northern Paraguay and the northernmost Chaco in Bolivia, where up to 34 species lost core areas in some locations (Supplementary material Appendix 1 Fig. A3). By 2015, remaining core areas were mainly concentrated in southern Bolivia, north-eastern Paraguay and some smaller areas in northern Argentina (Supplementary material Appendix 1 Fig. A3). The cumulative core area lost for all species between 1985 and 2015 was 407 000 km 2 .
Our rarity-weighted richness analyses revealed that priority areas for the community of larger mammals as a whole covered large areas of the northern Chaco in 2015, mainly in Bolivia and northern Paraguay, as well as the eastern-most Chaco in Argentina (Fig. 6A). In contrast, hotspots where habitat destruction and hunting pressure acted simultaneously covered broad areas in north-western Paraguay, northeastern Argentina and south-western Bolivia (Fig. 6B). Hotspots of habitat destruction-only were spread across central and northern Paraguay, southern Bolivia and the central Chaco in Argentina; whereas hotspots of hunting pressure-only were most common in northern Paraguay, south-western Bolivia and northern Argentina (Fig. 6B). For comparison, threat hotspots based on species' global ranges were similar to those based on the Chaco ranges. This was different for priority areas, where calculations based on global ranges revealed priority areas in the Bolivian Chaco (as in the analyses using Chaco ranges), but to a lesser extent in Paraguay and northern Argentina, and not at all in eastern Argentina (Supplementary material Appendix 1 Fig. A6).

Discussion
Understanding the individual and combined effects of different threats to biodiversity is critical for identifying effective conservation interventions to halt biodiversity loss. Yet, we currently lack approaches to map the spatial footprints of threats at resolutions fine enough to be useful for conservation planning. By combining land-cover time series mapped from satellite imagery, habitat suitability models and hunting 8 pressure models, we reconstructed the footprints of habitat destruction and hunting for the entire community of larger mammals of South America's Gran Chaco, a 1.1 million km 2 deforestation hotspot. We found that the footprints of both threats expanded considerably -and much more than deforestation alone -across the Chaco over three decades, producing a widespread loss of core areas. In addition, habitat destruction and hunting pressure acted simultaneously over increasing portions of the Chaco over time, suggesting that synergistic effects are becoming the norm. The priority areas and hotspots of threat that we identified point to key areas for larger mammals, where habitat protection and threatspecific management actions should swiftly be implemented to avoid further biodiversity loss. Overall, our findings suggest increasing synergistic effects between habitat destruction and hunting pressure in the Chaco, a situation likely common in many tropical deforestation frontiers around the world. Our work therefore highlights the urgent need to better understanding how these threats act on species in space and time, in other words, the geographies of threat to biodiversity. We here develop an effective and easily transferable approach to do so.  The footprints of habitat destruction and hunting pressure expanded hugely across the Chaco between 1985 and 2015 for almost all mammals we assessed. This is exemplified by the cumulative footprints of threats expanded over more than double the area of forest and natural grassland loss in that period (142 600 and 31 700 km 2 , respectively) (Baumann et al. 2017). The footprint of hunting pressure penetrated even further into remote areas, including protected areas, than habitat destruction. Hunting is the main cryptic disturbance for mammals, and often extends into otherwise 'intact' forests (Peres et al. 2016, Benítez-López et al. 2019. Similarly, the footprint of habitat destruction was also broader than that of deforestation, likely because small fragments are unsuitable for wide-ranging species, and because edge effects decrease resource availability close to deforested areas (Barlow et al. 2016). Only a few forest patches remain in the Chaco that are large enough to be effectively remote from hunter access points and agricultural lands. Other studies in deforestation frontiers have also reported that anthropogenic disturbance can extend over much larger areas than the area undergoing deforestation alone (Peres et al. 2006, Barlow et al. 2016). Together, our results highlight that approximating threats by deforestation footprints alone (Ocampo-Peñuela et al. 2016, Symes et al. 2018, Gallego-Zamorano et al. 2020, or by using fixed distances from roads (Allan et al. 2019) may underestimate the footprints of threats substantially.
The footprint of both threats increased since 1985 across all Chaco countries, but at varied rates. Habitat destruction expanded the most in Paraguay, which reflects Paraguay's rampant conversion of forests into pastures (Baumann et al. 2017). Habitat destruction expanded less in Bolivia, partly because two large protected areas cover large forested regions, and because the main deforestation frontiers in Bolivia are in the Chiquitano forest, just north of the Chaco (Hansen et al. 2013). Hunting pressure expanded more in Paraguay and Bolivia, where human population and road construction increased recently, than in Argentina, where   and high hunting pressure (2015), which represent priority areas for threat-specific conservation action. human population density and road density were already high in 1985. In fact, some species, such as Geoffroy's cat, experienced a decreasing hunting pressure in some areas in Argentina. This is likely because subsistence ranchers abandoned some areas as agribusiness expanded (Grau et al. 2008), potentially decreasing hunting pressure but increasing habitat destruction. After 2015, deforestation and forest fires have further advanced in all three countries, most worryingly in some of the last remote areas in northernmost Paraguay (Hansen et al. 2013) and in the northern Bolivian Chaco (Romero-Muñoz et al. 2019a). This highlights the urgency for stronger regulation of deforestation and the expansion of roads across all three countries.
The rapid expansion of threats and the massive declines of core areas, predicted for the first time by our maps, signify the defaunation of the larger mammal community across much of the Chaco. Unfortunately, these trends are widespread in deforestation frontiers (Gibson et al. 2011, Barlow et al. 2016. The declines we detected in most species' core areas often contrast with their generally low-threat global conservation status (Supplementary material Appendix 1 Table A3), highlighting the importance of conducting such assessments at the regional level (de la Torre et al. 2018). Given the varied and key ecological roles of larger mammals, their disappearance can disturb ecosystem functioning, including seed dispersal, carbon storage and nutrient cycling (Dirzo et al. 2014, Periago et al. 2014. This also highlights the importance of the few large remaining core areas for the mammal community as a whole, which are likely to be the last places maintaining the original species assemblage and ecosystem functioning in the Chaco. This reinforces the recognition of the irreplaceable role of 'wilderness' and Indigenous territories in maintaining biodiversity (Ricketts et al. 2010, Gibson et al. 2011). Further, these results underline the importance of halting further agricultural and road expansion into remaining core areas, which could otherwise disappear quickly across the entire Chaco.
A key result of our study was that areas where both threats act together cover increasingly larger portions of the Chaco. This is highly worrying because biodiversity declines even faster where threats synergise (Brook et al. 2008). Such synergistic effects are particularly likely in the Chaco, because its dense and thorny forests make them very hard to access for hunters unless forests are cleared for roads and agriculture. Hunters often kill mammals crossing such clearings; and workers cutting the forest, building fences and producing charcoal actively hunt animals in the remaining forest patches (Altrichter 2006;unpubl.). Accordingly, large mammals tend to disappear from forest strips and smaller forest patches soon after the surrounding areas are deforested (Núñez-Regueiro et al. 2015, Semper-Pascual et al. 2019. Furthermore, in areas already converted to agriculture, ranchers and farmers often persecute carnivores and herbivores thought to cause livestock or crop losses, respectively (Quiroga et al. 2016, Camino et al. 2018. While synergistic effects have been described through non-spatial methods in other deforestation frontiers (Peres 2001), here we provide an approach to map out the individual and combined effect of threats, and thus to track synergistic effects that may be common in deforestation frontiers around the world over time.
Our approach can also provide spatial templates for conservation planning. Our priority areas represent the most important areas for proactive conservation action, such as establishing protected areas. The protected area network currently covers only 9.1% of the Chaco. Extensive priority areas remain unprotected, particularly in northern Paraguay, and northern Argentina, and most are surrounded by threat hotspots. These areas are excellent candidate regions for expanding the existing protected area network and our analyses can serve to update previous prioritization exercises (TNC et al. 2005, Nori et al. 2016). Further, efforts should be directed to ensure Indigenous Peoples' land rights as many of these lands harbour priority areas and are thus important for Chacoan biodiversity.
Our threat hotspots overlapped extensively with previously prioritized areas (TNC et al. 2005, Nori et al. 2016, particularly in the central Chaco. This highlights the need for swift reactive threat management. In hotspots of habitat destruction, potential actions include 1) stopping further agricultural expansion and enforcing existing regulations, 2) securing Indigenous People's rights to land, 3) promoting culturally acceptable livelihoods that encourage sustainable land use and 4) fostering forest recovery. In hotspots of hunting pressure, specific actions include 1) careful planning of new roads and other land changes that foster access for hunters; 2) educational programs and improved management to lower or avoid conflicts with wildlife; 3) enforcing bans on recreational and commercial hunting and 4) ensure the sustainability of Indigenous People subsistence hunting. Several of these recommendations are in agreement with different Indigenous and smallholders visions in the Bolivian and Argentinean Chaco (Noss andCuellar 2001, Camino et al. 2016). Where both threats co-occur, they must be managed simultaneously. Implementing such complementary management actions is more likely to produce conservation gains than addressing single threats alone (Wilson et al. 2007).
Our work represents, to our knowledge, the first spatially explicit and high-resolution mapping of the footprints of multiple threat at the community level. Thereby it advances previous analyses assessing single threats (Ocampo-Peñuela et al. 2016, Benítez-López et al. 2019, threat interactions for individual species (Bleyhl et al. 2015, Romero-Muñoz et al. 2019b, and coarse-grained overlays of multiple threats based on species range maps (Symes et al. 2018, Allan et al. 2019, Gallego-Zamorano et al. 2020). Our study is also the first to reconstruct changes in multiple threats over long time periods, by combining satellite-based land-cover change maps with longitudinal datasets of road networks and over 27 000 subsistence ranches. Still, our work contains some limitations. First, although we gathered the largest occurrence dataset ever collected for the larger mammal community of the Chaco, presence points were scarce for some species in some regions, particularly the northern and southern Chaco for the 1980s. Second, our maps depend on thresholds for classifying threat levels, and we applied common criteria to define them. Still, other thresholds would change our maps. Finally, we used the human population density layers for 2000 also for 1985, because a comparable dataset for 1985 was missing. Although human population has likely not changed markedly in 1985-2000, we may have underestimated changes in hunting pressure for this period. This highlights the importance of long-term human population timeseries to transfer approaches such as ours to other regions (Lloyd et al. 2017).
Mapping the spatial footprints of habitat destruction and overexploitation has been hard, constituting a real barrier towards better understanding their individual versus combined impacts, and for targeting threat-specific conservation planning. Here, we pioneer a new approach to reconstruct the changing footprints of main threats to biodiversity (Fig. 1). Applying this approach to the 1.1 million km 2 Gran Chaco, a global deforestation hotspot, we find that the footprints of habitat destruction, hunting pressure and the areas where they synergize, are rapidly expanding. Such trends are likely common across other deforestation frontiers in Latin America, Africa and south-east Asia and our approach should therefore be broadly applicable to assess the geographies of threat in these regions. Our approach also allows to identify the remaining priority areas for biodiversity and to pinpoint to where threat-specific conservation actions to halt biodiversity declines should be implemented. Overall, our study highlights the importance of understanding and addressing the combined effects of major threats to biodiversity in order to better tackle biodiversity loss.