Habitat destruction and overexploitation drive widespread declines in all facets of mammalian diversity in the Gran Chaco

Global biodiversity is under high and rising anthropogenic pressure. Yet, how the taxonomic, phylogenetic, and functional facets of biodiversity are affected by different threats over time is unclear. This is particularly true for the two main drivers of the current biodiversity crisis: habitat destruction and overexploitation. We provide the first long‐term assessment of multifaceted biodiversity changes caused by these threats for any tropical region. Focussing on larger mammals in South America's 1.1 million km2 Gran Chaco region, we assessed changes in multiple biodiversity facets between 1985 and 2015, determined which threats drive those changes, and identified remaining key areas for all biodiversity facets. Using habitat and threat maps, we found, first, that between 1985 and 2015 taxonomic (TD), phylogenetic (PD) and functional (FD) diversity all declined drastically across over half of the area assessed. FD declined about 50% faster than TD and PD, and these declines were mainly driven by species loss, rather than species turnover. Second, habitat destruction, hunting, and both threats together contributed ~57%, ~37%, and ~6% to overall facet declines, respectively. However, hunting pressure increased where TD and PD declined most strongly, whereas habitat destruction disproportionally contributed to FD declines. Third, just 23% of the Chaco would have to be protected to safeguard the top 17% of all three facets. Our findings uncover a widespread impoverishment of mammal species richness, evolutionary history, and ecological functions across broad areas of the Chaco due to increasing habitat destruction and hunting. Moreover, our results pinpoint key areas that should be preserved and managed to maintain all facets of mammalian diversity across the Chaco. More generally, our work highlights how long‐term changes in biodiversity facets can be assessed and attributed to specific threats, to better understand human impacts on biodiversity and to guide conservation planning to mitigate them.


| INTRODUC TI ON
Human activities are driving the global biodiversity crisis, and the two biggest threats are habitat destruction and overexploitation (Díaz et al., 2019;Hooper et al., 2012;IPBES, 2019). Assessing patterns of biodiversity change due to these threats is therefore crucial for conserving biodiversity and achieving sustainability goals Díaz et al., 2019;IPBES, 2019). Most efforts assessing biodiversity change across broad spatial scales have focused on taxonomic diversity and have not connected those changes to multiple threats (e.g. Dornelas et al., 2014;Kerbiriou et al., 2009;Tingley & Beissinger, 2013). Yet, a focus on taxonomic diversity neglects evolutionary history and long-term evolutionary potential (Winter et al., 2013). Likewise, taxonomic diversity overlooks the diverse ecological functions of species in ecosystems (Cadotte et al., 2011;Winter et al., 2013), which maintain ecosystem integrity and functioning and ultimately provide nature's contributions to people (Cadotte et al., 2011;Cardinale et al., 2012;Díaz et al., 2019).
Therefore, assessing how different threats contribute to long-term changes in all three biodiversity facets, taxonomic (TD), phylogenetic (PD), and functional (FD) diversity, is crucial for a more comprehensive understanding of human impacts on nature. Likewise, this understanding is also vital to develop conservation strategies that better account for all biodiversity facets.
The varied and widespread change across biodiversity components is better understood through temporal assessments-rather than space-for-time substitutions-particularly across rapidly changing regions (Damgaard, 2019). Long-term studies have reported decreasing (Tingley & Beissinger, 2013), increasing (Kerbiriou et al., 2009), or no net change of TD (Dornelas et al., 2014). Likewise, the few studies focusing on more than one biodiversity facet have reported both similar (Jarzyna & Jetz, 2017) and different (Monnet et al., 2014;Villéger et al., 2010) temporal trends among facets. Only two studies, focusing on birds across France (Monnet et al., 2014) and the USA (Jarzyna & Jetz, 2017), have simultaneously assessed the long-term changes of all three biodiversity facets at broad scales.
All facets increased in both studies, except for FD in France, which remained stable.
In contrast, the long-term multifaceted changes of biodiversity in the tropics-where the overwhelming majority of Earth's biodiversity resides, and where their main threats are expanding the fastest (Barlow et al., 2018;Bradshaw et al., 2009)-remain largely unexplored. The downward trends of TD and of natural habitats reported in the tropics (Barlow et al., 2016;Hansen et al., 2013;Romero-Muñoz et al., 2020) suggest that biodiversity trends likely differ substantially from those reported from temperate regions.
Recent advances in remote sensing and ecological modelling allow us to reconstruct detailed land-use change histories, as well as the distributions of species and the spatial footprints of threats for multiple species across several decades and large regions (Baumann et al., 2017;Benítez-López et al., 2019;Romero-Muñoz et al., 2020).
Together with increasingly available trait and phylogenetic information, these developments open the opportunity to assess long-term changes in multiple biodiversity facets across rapidly changing regions.
Habitat destruction and overexploitation are the leading drivers of global biodiversity decline (IPBES, 2019;Maxwell et al., 2016).
Both threats are rapidly expanding into previously natural areas across the tropics due to the increasing human demand for agricultural products, such as beef, soy (predominantly used as livestock feed), and palm oil (Kehoe et al., 2017;Laurance et al., 2009). Yet, these threats affect species differently (Ripple et al., 2017;Romero-Muñoz et al., 2020), and therefore may affect biodiversity facets differently. For instance, habitat destruction and degradation through land-use change may disproportionately affect species within specific phylogenetic lineages (Frishkoff et al., 2014;Nowakowski et al., 2018) or with specific traits (Newbold et al., 2020;Wordley et al., 2017). Likewise, species from some lineages (D'agata et al., 2014;Davis et al., 2018), or with certain traits, like large body size (Benítez-López et al., 2019;Ripple et al., 2017), are more vulnerable to overexploitation.
Habitat destruction and hunting pressure are widespread in tropical regions (Gallego-Zamorano et al., 2020;Romero-Muñoz et al., 2020) and, where these threats co-occur, they exacerbate biodiversity loss even more than either threat alone (Brook et al., 2008;Mouillot et al., 2013;Peres, 2001). Despite the importance of assessing and mapping the combined impact of these major threats, previous studies have focused either on individual threats, usually only habitat modification, when assessing changes in several biodiversity facets (e.g. Chapman et al., 2018;Frishkoff et al., 2014;Wordley et al., 2017), or on multiple threats for single species (Romero-Muñoz et al., 2019) or only on TD (Romero-Muñoz et al., 2020). The contribution of the individual versus combined effects of threats to changes in multiple biodiversity facets remains so far unexplored. This is unfortunate, as learning how the three biodiversity facets are impacted by threats in space and time would enable more effective conservation planning, through threat-specific targeting of conservation actions (Devictor et al., 2010;Pollock et al., 2017).
Although conservation planning often assumes that one facet also represents others, recent studies found considerable spatial mismatches among facets (Devictor et al., 2010;Mazel et al., 2018;Safi et al., 2011). Therefore, conservation planning could benefit from identifying the most important areas for each biodiversity facet, as well as where those areas overlap. However, current methods to map facets are not ideal, because most rely on expert-based species range maps, which represent single snapshots in time, contain errors, and are built at a coarse scale (Ficetola et al., 2014). As threats are expanding and intensifying (Allan et al., 2019;Benítez-López et al., 2017;Kehoe et al., 2017), there is an urgent need to map the spatial congruence of all three biodiversity facets at resolutions fine enough to inform conservation planning. This is particularly urgent in tropical deforestation frontiers, which are global hotspots of biodiversity loss (Barlow et al., 2018;Bradshaw et al., 2009;Hoekstra et al., 2005). Many such frontiers, particularly those in tropical dry forests, are weakly protected (Hoekstra et al., 2005;Kuemmerle et al., 2017). The Gran Chaco (hereafter 'Chaco') in South America is one of the most at at-risk regions globally, due to rapid expansion of cattle ranching and soy cultivation WWF, 2015). The region is a global hotspot of habitat conversion and defaunation (Baumann et al., 2017;Romero-Muñoz et al., 2020), yet despite calls for assessing the facets of biodiversity in this region (Periago et al., 2014), no such assessment exists.
Here, we provide the first assessment of TD, FD, and PD for the large-and medium-sized mammals (>~1 kg, hereafter 'larger mammals') in the Chaco, a global deforestation hotspot. These species represent a phylogenetically diverse group, including some lineages endemic to the Chaco, such as those represented by the Chacoan peccary (Catagonus wagneri), the Chacoan mara (Dolichotis salinicola), and several armadillo species (Nori et al., 2016). The enormous variation in size and morphology among larger mammals in the Chaco also translates into a high diversity of ecological roles and resource uses.
Several ecosystem functions are unique to larger mammals, such as dispersing the seeds of the largest trees or regulating the populations of other large animals (Lacher et al., 2019). Such roles have significant effects on nutrient cycling and energy flows, and in structuring ecological communities and thus promoting high biodiversity and ecosystem stability (Lacher et al., 2019;Terborgh, 2015). In turn, larger mammals provide various contributions to people directly, such as by being sources of protein for local communities, and indirectly, such as by enhancing forest regeneration and carbon storage capacity (Bello et al., 2015;Noss et al., 2004). Many larger mammals in the Chaco are highly vulnerable to habitat destruction and hunting (Romero-Muñoz et al., 2020;Semper-Pascual et al., 2018) and their declines threaten to erase unique evolutionary histories, affect ecosystem integrity, and negatively impact nature's contributions to people.
Here we aim to assess 30 years of change in the three facets of mammalian diversity in the Chaco, and explore how habitat destruction and overexploitation contributed to these changes. To our knowledge, this represents the first assessment of this kind for (a) mammals, (b) the tropics, and (c) in relation to multiple, interacting threats. Specifically, we aim to answer:

| Study region
The Gran Chaco region extends across 1.1 million km 2 in Argentina, Paraguay, and Bolivia, and is the largest tropical and subtropical dry forest in the world. Xeric forests are the dominant vegetation formation, with interspersed mosaics of natural savannas and gallery forests. The climate ranges from tropical in the north to subtropical in the centre and south. Precipitation ranges from 1,400 mm in the east to 400 mm in the west and south. The Chaco is rich in biodiversity, with over 150 mammal species, 500 birds, and over 3,000 plant species (TNC, FVS, FDSC, & WCS, 2005). Over the last decades, the Chaco has become a global deforestation hotspot, losing 20% of its forests since 1985 due to the expanding croplands, mainly in Argentina, and livestock ranching, mainly in Paraguay and Bolivia (Baumann et al., 2017). These pressures are likely impacting ecosystem functioning over large scales (Periago et al., 2014), although this has not been yet quantified. Despite these pressures, only about 9% of the Chaco is currently protected (Nori et al., 2016).

| Datasets used
We produced maps of the habitats and the footprints of habitat destruction and hunting pressure separately for 48 larger mammals between 1985 and 2015 at a 1 km 2 resolution in an earlier study (Romero-Muñoz et al., 2020). This was done using habitat suitability models to track habitat suitability across space and time, and hunting pressure models to do the same for the hunting risk. We performed separate multi-temporal habitat suitability models and hunting pressure models for each species.
To assess habitat suitability, we used the largest database of presence records of larger mammals ever collected for the Chaco, containing occurrences from 1985 to 2015, which we analysed using maximum entropy modelling (Maxent v3.4.1; Phillips et al., 2017).
Maxent predicts a species' occurrence' by comparing the locations of recorded presences to the overall distribution of environmental predictors for a study region, which are sampled through background points (Phillips et al., 2017). We generated seven predictors related to habitat suitability for mammals for 1985 and 2015, at a 1 km 2 resolution: % forest, % pastures, % cropland, % forest edge, distance to water, mean annual temperature, and mean annual precipitation. All Maxent models were parameterized using only hinge features to avoid overfitting, a regularization multiplier of 1, and a prevalence value of 0.5, and we controlled for sampling bias in our occurrence and background datasets. We cross-validated all models using aver- We refitted the original global model to Neotropical mammals only (n = 1,974 abundance changes). We included the distance to hunter's access points and human population density as spatial predictors of hunting risk, and species' body mass as predictor of species-specific vulnerability to population decline due to hunting (Benítez-López et al., 2019). The result of our hunting pressure model is a hunting pressure index ranging from 0 (no decline in abundance) to 1 (local extirpation; see Romero-Muñoz et al., 2020).
We projected the habitat suitability and hunting pressure models to 1985(Romero-Muñoz et al., 2020. For each year, we classified the maps resulting from the habitat suitability model into good and poor habitat suitability using the maximum sensitivity plus specificity threshold (Liu et al., 2013). Similarly, we classified the hunting pressure maps for each species into low and high hunting pressure. We used a 30% abundance decline as threshold to classify a species as threatened (here representing high hunting pressure), following the IUCN Red List criteria (IUCN, 2012). Overlapping these classified, binary habitat suitability and hunting pressure maps per species highlighted four areas differently affected by threats at each time step: (a) areas with good habitat suitability and low hunting pressure (hereafter: core areas), (b) areas with good habitat suitability but high hunting pressure (hunting pressure), (c) poor habitat suitability but low hunting pressure (poor habitat), and (d) poor habitat and high hunting pressure (co-occurring threats). Furthermore, we considered the change from core area in 1985 to 'poor habitat' in 2015 as 'habitat destruction' and from core to 'hunting pressure' as increasing hunting pressure (hereafter simply 'hunting pressure'), because such increasing threats can be solely attributed to anthropogenic impacts during this period.
We then calculated the three biodiversity facets at a 5 × 5 km 2 resolution, where each gridcell represents a community of mammals.
The 25 km 2 gridcell size is meaningful in our case as it allows for integrating across species with a wide range of home range sizes (mean = 9.7 km 2 , SD = ±22.2; gathered from Jones et al., 2009), while being fine enough for regional conservation planning. To aggregate our species level, 1 km 2 resolution maps to the 5 km grid, we assigned the most frequent habitat or threat category at the 1 km grids. We considered a species present in a grid cell if it had core area in it.

| Depicting biodiversity facets
We calculated metrics for each of the three facets for the years 1985 and 2015 individually (Figure 1). We derived TD as species richness (i.e. number of species per gridcell). In addition, we assessed the change in community composition over time, (i.e. temporal community dissimilarity) and the contribution of its components (Baselga, 2010). Changes in the species composition of a community result from changes in species richness and the replacement of some species by others (i.e. turnover), or a combination of both (Baselga, 2010). Measures of total dissimilarity, such as the Sørensen dissimilarity index, can thus be decomposed into its turnover and species richness change components (Baselga, 2010). While the Sørensen index measures total dissimilarity, the Simpson dissimilarity index accounts only for the turnover component. Thus, the difference between both indices accounts for the species richness change component of dissimilarity (Baselga, 2010). We assessed the temporal community dissimilarity only for TD, as equivalent methods have not yet been developed for the other biodiversity facets (Baselga & F I G U R E 1 Framework to quantify and map changes in the three facets of mammalian diversity across our study region. FD, functional diversity; PD, phylogenetic diversity; TD, taxonomic diversity Orme, 2012). We determined the contribution of species richness change and turnover to total dissimilarity change between 1985 and 2015 for each community using the beta.temp function in the R package betapart v1.5.1 (Baselga & Orme, 2012).
For measuring phylogenetic diversity, we used Faith's PD index, which represents the minimum total length of the phylogenetic tree's branches of the species within each community (Faith, 1992).
To account for phylogenetic uncertainty, we extracted an average tree for the entire set of species from a set of 1,000 trees available in the PHYLACINE database  using the function averageTree in the phytools package (Revell, 2012). Based on this phylogenetic tree, we calculated the PD index using the Picante package in R (Kembel et al., 2010).
As our measure of FD, we calculated functional richness (also known as 'FRic'; Villéger et al., 2008). This index represents the amount of multidimensional functional space occupied by all the species in the community (Villéger et al., 2008). We chose this index instead of the commonly used dendrogram-based index (Petchey & Gaston, 2002), because the latter has been shown to produce biased estimates of FD, leading to inaccurate biogeographical patterns (Maire et al., 2015). We first gathered a database of traits related to resource use (Table S1), assessed the collinearity among traits through a Pearson's correlation test, to ensure that traits had r < .5 to and therefore non-redundant (Villéger et al., 2008; Figure S1). Our final list included seven traits: diet, use of forest strata, use of day/night, home range size, body mass, generation length, and number of offspring per year, gathered from several sources (Table S1; Jones et al., 2009;Myers et al., 2019;Tacutu et al., 2012;Wilman et al., 2014). To calculate FD, we used a distance-based framework, based on the Gower distance among traits (Laliberté & Legendre, 2010). We weighted traits equally and applied a principal coordinates analysis (PCoA) ordination based on the distance matrix to build a multidimensional functional space. We calculated FD using the FD package v1.0-12 in R (Laliberté & Legendre, 2010).
The quality of the functional space used to calculate FD (i.e. how well it represents the initial trait values) depends on the number of dimensions of the multidimensional functional space spanned by the PCoA axis (Maire et al., 2015). We compared the quality of functional spaces produced by two to seven dimensions by calculating the mean squared deviation (mSD) metric (Maire et al., 2015). This metric assesses the degree of consistency between the initial and final functional distances (the closer mSD is to 0, the higher the quality of the functional space). As can be expected (Maire et al., 2015), functional space quality increased with the number of dimensions ( Figure S2a). However, there is trade-off between functional space quality and the spatial comprehensiveness of our analyses, as FD can only be calculated for communities with more species than the number of dimensions (Villéger et al., 2008). Therefore, the more dimensions are used to calculate FD, the higher the number of gridcells that will be dropped from the analyses ( Figure S2b). Aiming to produce a high-quality functional space that still allows us to estimate FD across a large portion of our regions, we opted for a four-dimensional functional space in our case (mSD = 0.023). About

| Threat effects on biodiversity facets
A challenge in assessing the impacts of multiple threats on multiple biodiversity facets is the attribution of community-level declines in facets to particular threats, because different threats can affect species differently in different areas. We addressed this challenge using newly developed measures for assessing species' functional and phylogenetic distinctiveness within communities Violle et al., 2017). Specifically, we assessed the relative importance of habitat destruction and hunting pressure on the decline of each facet by first assessing which species lost core areas in a given gridcell (i.e. one or both threats became prevalent for that species in that gridcell), and second summing up threats responsible for that loss (habitat destruction, hunting, or co-occurring threats) while weighting species according to their distinctiveness.
Distinctiveness is defined as the mean distance in the diversity measure used to the N other species within each community . For TD, these distances between species are always 1. For FD and PD, we calculated distinctiveness per species for the community in 1985-here considered as the baseline year-using the 'distinctiveness' function in the funrar package v1.4.1 in R (Grenié et al., 2017). A species' taxonomic, phylogenetic, or functional distinctiveness is calculated according to: where N is the number of species within the community, and d ij is the taxonomic, functional, or phylogenetic distance between species i and j. We scaled d ij between 0 and 1 using a minmax transformation .
We used this approach to calculate the relative importance of habitat destruction, hunting pressure, or co-occurring threats for biodiversity loss, separately for each of our three facets. To assess whether the relative importance of threats changed at higher values of facet loss, we repeated this procedure at different thresholds of loss per facet (e.g. for the top 75%, 50%, and 25% of the gridcells with the highest losses per facet). We also assessed the spatial congruence among the three biodiversity facets for 2015 by measuring the pairwise correlation of facet values across all gridcells. We then mapped the gridcells with the top 5%, 10%, 17%, and 25% values per facet (Brum et al., 2017). We assessed the overlap of the top 17% gridcells among facets-the minimum surface recommended for protection by the Convention on Biological Diversity's Aichi target 11 (Tittensor et al., 2014).

| RE SULTS
We Regarding the overall community composition change over time, the contribution of species richness change was larger (mean = 0.11 ± 0.14 SD; median = 0.06) than that of turnover (0.09 ± 0.18 SD; median = 0.00) to total Sørensen temporal dissimilarity (mean = 0.20 ± 0.21 SD; median = 0.16; Figure S4). These differences were highly significant (Wilcoxon signed ranks test, V = 217,830,000, p < .001). Among communities that changed, species richness change (and here specifically species loss) was a larger contributor to dissimilarity than turnover (61% vs. 39% of communities).
The expansion of habitat destruction and hunting pressure be- and co-occurring threats (6%; Figure 4a). Focussing on those areas experiencing highest loss in a facet revealed some interesting differences compared to the region-wide results. When considering the 25% of gridcells with highest loss per facet, the relative importance of hunting increased (to 38% and 39% for TD and PD, respectively), whereas the relative importance of habitat destruction increased to 64% for FD. In addition, the relative importance of co-occurring threats increased slightly for all facets when focusing on those areas experiencing highest loss in a facet (Figure 4b). The contributions of threats to facet declines were practically identical for the entire Chaco and the areas of the Chaco with communities with five or more species ( Figure S7).

F I G U R E 3
The

| D ISCUSS I ON
To better understand the impacts of people on nature, we need to learn how different facets of biodiversity change in response to anthropogenic threats. Here, we provide the first multi-decadal, broadscale assessment of changes in all three biodiversity facets caused by specific threats. Furthermore, this is to our knowledge the first assessment of changes in multiple biodiversity facets for mammals and in the tropics. Using habitat maps for 48 larger mammals and the spatial footprints of habitat destruction and hunting, we assessed how these threats, individually and jointly, drive changes in mammalian diversity over 30 years across the 1.1 million km 2 Chaco. Our analyses reveal a general biotic impoverishment. This is illustrated by the rapid and widespread changes in mammalian communities, resulting in declines across all biodiversity facets. These changes were mainly driven by defaunation rather than by species turnover.
Habitat destruction was the main threat responsible for declines such as forest regeneration and the regulation of herbivore and mesopredator abundances (Bello et al., 2015;Terborgh, 2015).
Temporal changes in species composition were mainly driven by changes in species richness, whereas species replacement played a smaller role. These findings strongly suggest that the overall dominance of changes in total richness, specifically species loss, also explains the downward trends we found for phylogenetic and functional diversity. Nevertheless, although turnover was less important, the replacement of some species with distinctive functions (e.g. the Azaras's capuchin monkey Sapajus cay) or phylogeny (e.g. the Chacoan naked-tailed armadillo) by less distinctive species may drive PD and FD declines in some communities. Our finding of the dominance of species loss contrasts with reported global trends, where species richness has been on average relatively stable over time across local studies, and turnover has been the main driver of community changes (Blowes et al., 2019;Dornelas et al., 2014). However, it is unclear what proportion of the communities assessed at global scale are under similarly high pressure from habitat destruction and hunting as communities in the Chaco.
The varying rates and spatial patterns of decline in biodiversity facets that we found suggest that trends observed for a single facet may conceal how other facets change. For instance, the higher similarity in the spatial patterns of decline of TD with PD than with FD suggests that TD is an imperfect surrogate for changes in other facets-although this is often assumed (Dornelas et al., 2014;McGill et al., 2015). Despite varying geographical patterns of decline, on average all facets declined. This adds further evidence that change in biodiversity facets is context-specific, as exemplified for birds, where all facets changed in parallel in the United States (Jarzyna & Jetz, 2017), but not in France (Monnet et al., 2014). This also could imply that recently reported observations of long-term stability of local species richness, but high species turnover across the world (Blowes et al., 2019;Dornelas et al., 2014) may conceal important trends in other facets. Our work, the first from the tropical and subtropical biomes, thus reinforces calls based on studies from temperate regions (Jarzyna & Jetz, 2017;Monnet et al., 2014) to assess long-term change of all facets of biodiversity.
Our results advance the previously limited understanding of the downward trends of biodiversity facets in tropical deforestation frontiers. Although the calculation of FD and the comparison of temporal changes among facets were limited to communities with five or more species (see Section 2), these communities, mainly located in the northern half of the Chaco, faced most of the changes in land use and threat levels over the last three decades (Baumann et al., 2017;Romero-Muñoz et al., 2020). In contrast, areas in the southern Chaco for which FD could not be calculated remained largely stable since 1985 (Baumann et al., 2017;Romero-Muñoz et al., 2020). Overall, the changes we report represent considerable declines in one or more facets of mammalian diversity for tens of thousands of mammal communities. These troubling widespread declines across biodiversity facets may be common in deforestation frontiers, highlighting the urgency of assessing them in such regions. Furthermore, our work provides spatially explicit, fine-scale trends of change of biodiversity facets for individual communities in the Chaco.
A major advance of our approach is the ability to attribute declines in biodiversity facets to specific threats. The greater contribution of habitat destruction to declines across all facets is partly due to habitat destruction expanding over a ~41% larger area than hunting pressure. However, in the 25% of gridcells with the highest declines per facet, hunting pressure contributed more to TD and PD declines. This 25% of gridcells includes many remote and highly diverse areas. In total, over half the Chaco is currently impacted by one or the other threat, with co-occurring threats contributing much less to declines in all facets. This is partly because we only considered changes from core areas to areas under threat over time; when Romero-Muñoz et al. (2020) considered changes from one threat to multiple threats for TD, the importance of synergistic threats was substantial. Furthermore, it is important to note that we focus on core areas, whereas species could remain outside them. Yet, outside core areas species are affected by one or more threats, and may in many cases be locally functionally extinct, or committed to local extinction in the near future (Semper-Pascual et al., 2018). Overall, these results uncover the substantial and mutually amplifying importance of habitat destruction and overexploitation in deteriorating all facets of mammalian diversity.
Such large and widespread impacts of habitat destruction and hunting pressure on biodiversity facets may be common across tropical and subtropical deforestation frontiers. This is corroborated by various studies reporting negative effects of habitat modification on biodiversity facets (e.g. Chapman et al., 2018;Frishkoff et al., 2014;Wordley et al., 2017). However, only studies focusing on TD have assessed the relative impact of multiple threats. In such studies, habitat destruction contributed more than overexploitation to TD declines in tropical mammals and birds, but with wide variation among species (Gallego-Zamorano et al., 2020;Romero-Muñoz et al., 2020;Symes et al., 2018), implying that responses into threats may differ between TD and other facets.
Indeed, while habitat destruction was the main contributor to decline across facets in our case as well, it was particularly important for declines in FD. In turn, hunting pressure was particularly prevalent where TD and PD declined the most. Such differences are likely due to different susceptibility of species to different threats as well as the different spatial footprints of threats. For instance, habitat destruction affected many species across much of the Chaco, whereas hunting pressure typically affected larger species more strongly and increased mainly in more remote areas ( Figure S5). Our work therefore highlights the importance of simultaneously assessing the impact of multiple threats on multiple biodiversity facets and provides a novel framework to do so.
The three biodiversity facets we mapped showed moderate spatial congruence across the Chaco in 2015, with FD being more distinct than TD and PD (which were more congruent). This pattern appears to be mainly driven by FD being high in the 'Bañados del Ibera' area in easternmost Chaco, likely because a few species with unique trait combinations, such as the marsh deer Blastoceros dichotomus or river otter Londra longicaudis, concentrate in this savannah wetland (Romero-Muñoz et al., 2020). Our finding that FD had a more distinct spatial distribution than TD and PD at the regional scale corroborates patterns found for terrestrial birds and mammals (Pollock et al., 2017), and for marine mammals (Albouy et al., 2017) at the global scale. This underlines the importance of considering functional diversity in biodiversity assessments.
Our work contributes to conservation planning by detecting and mapping priorities for all facets of biodiversity at fine spatial resolutions. Identifying the most valuable areas per facet can help us to identify protection gaps and to prioritize areas for closing these gaps, which is important given that most conservation planning exercises have focused only on TD, but largely ignored PD and FD (Pressey et al., 2007). Furthermore, our approach allows us to map which threats affect different biodiversity facets and where. This can identify the best locations to implement threat-specific management actions, such as promoting forest recovery, fostering sustainable hunting (including the preferential targeting of hunting-resilient species such as rodents), through culturally appropriate education or sensitization programmes, or ensuring the land rights of Indigenous Peoples (Ripple et al., 2016). Such threat-specific actions are an increasingly important focus of conservation planning, yet have so far not focussed on multiple facets of biodiversity (Pressey et al., 2007;Tulloch et al., 2015). Our work illustrates how this can be achieved using the Chaco, a severely under-protected region in need of conservation planning, as a demonstration case (Periago et al., 2014).
In conclusion, our study reveals a widespread impoverishment in mammalian diversity, including in overall richness, evolutionary history, and ecological functions, across large areas in the Chaco since 1985. Our approach linking changes in biodiversity facets with specific threats allowed us to uncover how habitat destruction and hunting pressure individually and jointly drive declines across all biodiversity facets in this global deforestation hotspot.
Our approach, along with the resulting indicators and maps, can inform conservation planning by governments and conservation organizations by identifying priority areas for facet protection and for threat-specific management interventions. Overall, our study advances the understanding of where and how multiple biodiversity facets change over time in response to different extinction drivers.

ACK N OWLED G EM ENTS
We thank many researchers for sharing species and spatial datasets

DATA AVA I L A B I L I T Y S TAT E M E N T
Data available on request from the authors.