Halting predicted vertebrate declines requires tackling multiple drivers of biodiversity loss

TL;DR

Vertebrate populations decline faster when facing multiple threats like disease and climate change, not just habitat loss or exploitation.

Contribution

This study provides a global, population-level analysis of how multiple threats interact to drive biodiversity loss.

Findings

Disease, invasive species, pollution, and climate change are linked to faster vertebrate population declines.

Interactive threats contribute more to population declines than spatial or temporal variation.

Mitigating multiple threats is essential to halt biodiversity loss and achieve nonnegative population trends.

Abstract

Conservation policies aiming to halt biodiversity loss often focus on globally prevalent threats like habitat loss and exploitation, yet direct and interactive effects of multiple threats remain poorly quantified. Here, we go beyond prior meta-analyses or species-level studies by providing a global, population-level empirical analysis of threat interactions by examining 3129 vertebrate population time series worldwide with documented exposure to single and multiple threats. Populations affected solely by habitat loss or exploitation do not exhibit the steepest declines; instead, disease, invasive species, pollution, and climate change are associated with faster declines. Interactive threats contribute more to population declines than temporal or spatial variation. Counterfactual analyses reveal that mitigating multiple threats is essential to achieving nonnegative vertebrate population trends and halting biodiversity loss.

Global vertebrate populations decline faster in the presence of multiple threats compared to single threats.

Figures

Habitat loss and exploitation are the most prevalent threats affecting vertebrate populations.(A) Spatial distribution of the population trends from 3129 population time series from 1281 species contained in the LPD. (B) Proportion of populations affected by the different threats. A total of 20.00% of the populations are not exposed to threats, and 3.04% are exposed to disease, 5.69% to invasive species, 6.07% to climate change, 7.16% to pollution, 27.20% to habitat loss, and 30.60% to exploitation.

Invasive species, disease, and climate change accelerate the decline of vertebrate populations.(A) Effect of threats on vertebrate population trends of 3129 vertebrate population time series from 1281 species contained in the LPD. “None” represents the estimated trend in the absence of threats. Single threats indicate the effect of the threats when acting in isolation. Interactive threats indicate those acting in conjunction with other (one or two) threats. (B) Distribution of the influence of the random effects and the presence of threat interactions on the population trends (see Materials and Methods). (C) Distribution of estimated trends given the realization of the combined model coefficients presented in (A). Trends are calculated as the first derivative of the estimated model slope(s). All density plots are based on 1000 samples from the posterior distribution of the coefficients (A) or trend estimates (B and C). The reported values are the posterior density median values (circles) with 50% (thickest bars), 80%, and 95% (thinnest bars) uncertainty intervals.

Synergies are not the most common interactive effect across threats.Interactive effects of the different threats. Most interactive effects are additive (purple) or antagonistic. The interactive effects were classified by comparing the estimated population trends for interactive threats versus their additive combinations (additions of the population trends exposed to single threats). They were defined as additive when the trends for populations exposed to the combined threats were equal to the summation of the trends for populations exposed to the individual threats. They were defined as synergistic when the trends for populations exposed to the combined threats were larger than the summation of the trends for populations exposed to the individual threats. They were defined as antagonistic when the trends for populations exposed to the combined threats were smaller than the summation of the trends for populations exposed to the individual threats. No interactions were synergistic in the global model, while 84.90% were additive for climate change, 80.60% for disease, 89.20% for exploitation, 85.80% for habitat loss, 78.40% for invasive species, and 85.30% for pollution.

Counterfactuals suggest that lowering exploitation would be the most effective action to slow down the decline of vertebrate populations.The counterfactuals represent the changes in the population trends of the 1740 vertebrate time series affected by threats, had there been no pollution, no disease, no invasive species, no habitat loss, no climate change, no exploitation, or no threats. The continuous, gray line represents when the population trend is 0. The dashed gray line is the median trend of the threatened populations without any intervention. Box plots depict the distribution of population trends, where the thick line in the middle represents the median value; the lower and upper bounds of the box represent the 25th to 75th quartiles, respectively; and the straight lines represent the minimum to the lower quartile (the start of the box) and then from the upper quartile (the end of the box) to the maximum.

Tables

Table 1.. Definition and identification method of singular and interactive threats.

Term	Definition	Method of identification
Singular	The effect of a single stressor	Predicted trend in the presence of a single threat
Interactive	The combined effect of multiple stressors; this classification is made before estimating whether an interaction is additive, antagonistic, or synergistic	Predicted trend in the presence of two or more threats including a statistical interaction term
Additive	The net impact of multiple threats is equal to their single effects	Predicted trend in the presence of two or more threats but the statistical interaction term is “turned off” (see Materials and Methods)
Antagonistic	The net impact of multiple threats is less than their singular effects	The difference between the additively predicted trend and the interactively predicted trend
Synergistic	The net impact of multiple threats is greater than their singular effects	The difference between the additively predicted trend and the interactively predicted trend

Full text

INTRODUCTION

Anthropogenic effects on Earth’s natural ecosystems are now more pervasive than ever (1–3), with threats such as habitat loss, climatic change, invasive species, pollution, and exploitation threatening more than 1 million species with extinction (3–5). Many international conservation objectives, including the Kunming-Montreal Global Biodiversity Framework (6) and the Sustainable Development Goals (7), have been set to halt the loss of biodiversity (8, 9) and protect the services it provides to humanity (10, 11). To achieve these targets, conservation policies often rely on global rankings of threats derived from species assessments [e.g., International Union for Conservation of Nature (IUCN) Red List (1, 2)] or pressure maps (12, 13) that identify habitat loss and exploitation as the leading causes of biodiversity loss (1–3), while climate change is expected to become increasingly important in the coming decades (2, 14). However, such rankings typically reflect the reported prevalence rather than the actual ecological impact of each threat [but see (15)]. Furthermore, these datasets, while invaluable for global synthesis, are not free from biases, such as geographic and taxonomic gaps in data, skewing threat prevalence and impact estimates toward better-studied regions (e.g., Europe and North America) and taxa (e.g., birds and mammals) (16–18). Moreover, global assessments typically evaluate threats at the species level rather than the population scale (19, 20), potentially underestimating the true magnitude of global change impacts (21).

Explicitly accounting for the effects of threats on biodiversity is further complicated by the fact that most species globally [about 80% (2)] are exposed to more than one threat simultaneously (5). Multiple threats can interact in complex and nonlinear ways (22–24), and their cumulative effects can be higher (synergistic effects) or lower (antagonistic effects) than the sum of their individual effects (additive effects) (22, 23). Synergistic interactions are of particular concern for ecology and conservation (23, 25), because they have the potential to further accelerate biodiversity loss by amplifying the impacts of multiple threats (24, 26). In the same vein, antagonistic interactions risk increasing the impact of the other threat(s) if management actions target the mitigation of a single threat involved in the antagonism (27). Identifying the prevalence of synergistic, antagonistic, or additive interactions is therefore critical for prioritizing conservation actions (25), as well as gauging the true impacts of global change on biodiversity. While some studies emphasize the prevalence of synergies as a major cause of biodiversity loss (24), others suggest that additive effects are more common (23, 25). However, the difficulties of accounting for multiple stressors in the analysis of real-world data mean that such studies are largely based on meta-analyses of experimental studies (23, 25), with limited temporal and spatial scales. These constraints might render an incomplete picture of the impacts of multiple threats and preclude identifying which threats are more likely to cause synergistic interactions. Conservation efforts are known to improve biodiversity but are often insufficient to halt biodiversity decline (28), possibly due to many efforts targeting single threats. To date, no study has used empirical data to examine the prevalence of the different types of threat interactions at a global scale and their association with population trends (25).

Here, we show how information on threats to wildlife populations can both provide a more informative estimate of biodiversity trends and guide approaches to reverse the declines. We quantify the trends of populations exposed to single and multiple interacting threats to vertebrate populations at a global scale using a subset (3129) of time series from the Living Planet Database (LPD hereafter; Fig. 1) (29). This dataset contains information about which threats are affecting each population according to the data sources (e.g., scientific manuscripts and technical reports). Threats were classified into six broad categories—exploitation, habitat loss, disease, pollution, climate change, or invasive species—matching those of the IUCN threat categories (for comparability with previous studies) (5), as well as a category of “no threat” exposure. To determine which threats and threat combinations are associated with the most marked declines, we estimated the population trends exposed to single and multiple threats using multilevel Bayesian models, accounting for spatial and temporal autocorrelation (30). To determine the prevalence of synergistic, antagonistic, and additive interactions, we compared the estimated population trends for interactive threats versus their additive combinations (fig. S1 and Table 1). Briefly, we calculated the difference between predicted interactive trends (from our model) and the additive trends (the addition between the two individual trends), with three possible classes: additive (no difference predictions made by additive and interacting threats), antagonistic (interacting threat predictions are less negative than the additive predictions), or synergistic (interacting threat predictions are more negative than the additive predictions). Classes were defined using the 80% credible interval (CI). We consider that a synergy “worsens” trends—i.e., becomes more negative. For a more detailed explanation, please refer to the “Estimation of interactive vs additive effects of threats” section and fig. S1. Last, to assess the potential benefits of management actions, we developed counterfactuals describing how the removal of single or multiple threats would have altered the predicted global trends of vertebrate populations.

Habitat loss and exploitation are the most prevalent threats affecting vertebrate populations.(A) Spatial distribution of the population trends from 3129 population time series from 1281 species contained in the LPD. (B) Proportion of populations affected by the different threats. A total of 20.00% of the populations are not exposed to threats, and 3.04% are exposed to disease, 5.69% to invasive species, 6.07% to climate change, 7.16% to pollution, 27.20% to habitat loss, and 30.60% to exploitation.

RESULTS AND DISCUSSION

Habitat loss and exploitation are not associated with the fastest population declines

Although habitat loss and exploitation are the most prevalent threats (1, 2) and are often the main targets of conservation (7, 31), our results show that they are not associated with the fastest population declines. Our multilevel Bayesian models show that being exposed to any threat accelerates the decline of vertebrate populations worldwide (Fig. 2A and table S1). For instance, vertebrate populations exposed to invasive species decline five times faster (−6.01% year^−1^; 95% CI = [−30.76% year^−1^, 10.10% year^−1^]) than those not exposed to any threat (1.70% year^−1^, CI = [−13.64% year^−1^, 20.00% year^−1^]; Fig. 2A). Similarly, despite being much less prevalent threats (Fig. 1B), pollution (−1.79% year^−1^; CI = [−17.51% year^−1^, 17.05% year^−1^]), climate change (−3.14% year^−1^; CI = [−17.21% year^−1^, 17.05% year^−1^]), and disease (−4.63% year^−1^; CI = [−14.89% year^−1^, 6.49% year^−1^]) also associate with population declines (Fig. 2A). Last, contrary to our initial expectations, exploitation (−0.85% year^−1^; CI = [−15.47% year^−1^, 15.43% year^−1^]) and habitat loss (−0.42% year^−1^; CI = [−17.04% year^−1^, 18.71% year^−1^]) are the threats associated with the slowest average population declines but with large uncertainties (Fig. 2A).

Invasive species, disease, and climate change accelerate the decline of vertebrate populations.(A) Effect of threats on vertebrate population trends of 3129 vertebrate population time series from 1281 species contained in the LPD. “None” represents the estimated trend in the absence of threats. Single threats indicate the effect of the threats when acting in isolation. Interactive threats indicate those acting in conjunction with other (one or two) threats. (B) Distribution of the influence of the random effects and the presence of threat interactions on the population trends (see Materials and Methods). (C) Distribution of estimated trends given the realization of the combined model coefficients presented in (A). Trends are calculated as the first derivative of the estimated model slope(s). All density plots are based on 1000 samples from the posterior distribution of the coefficients (A) or trend estimates (B and C). The reported values are the posterior density median values (circles) with 50% (thickest bars), 80%, and 95% (thinnest bars) uncertainty intervals.

The interactive threats displayed a much wider range of trajectories than single ones, showing the most positive and negative trends (Fig. 2A). Pollution (−2.83% year^−1^; CI = [−44.96% year^−1^, 50.77% year^−1^]) and climate change (−1.42% year^−1^; CI = [−30.06% year^−1^, 33.70% year^−1^]) are the interactive threats displaying the largest variability (Fig. 2A), while the population trends displaying the fastest declines (table S1) are those experiencing the combined influence of pollution, climate change, and disease (−16.74% year^−1^; CI = [−31.55% year^−1^, 1.00% year^−1^]) and the interaction between pollution and invasive species (−10.27% year^−1^; CI = [−67.11% year^−1^, 134.32% year^−1^]). However, the latter is extremely uncertain given its limited representation in our dataset (Fig. 1). In terms of average trends, populations experiencing pollution (−2.83% year^−1^; CI = [−44.96% year^−1^, 50.77% year^−1^]) and disease (−2.97% year^−1^; CI = [−23.61% year^−1^, 18.56% year^−1^]) interactions with any other threats show the fastest average decline (Fig. 2A). Overall, these results underscore the complexity of global conservation, emphasizing the context-dependent nature of these threats and how their interactions add to the uncertainty of biodiversity trends.

Although vertebrate population trends exhibit considerable uncertainty when assessing the effects of interactive threats (particularly between taxa and systems: figs. S2 to S5 and tables S2 to S3), some studies suggest that this uncertainty may arise from correlative factors, such as spatial and temporal autocorrelation (30). Therefore, we compared the impact of spatial and temporal autocorrelation (random effects) with the overall effect of interactive threats on population trends by estimating their differences. Specifically, we contrasted trend predictions without random effects or threat interactions with those including these factors and then subtracted the latter from the former (see Fig. 2B and “Influence of threats on population trends” section for further details). Our findings indicate that, despite the variability in population trends, interactive threats are associated with stronger declines than random effects (Fig. 2B). Therefore, these results provide evidence that interactive threats are contributing to estimated biodiversity loss, more than temporal or spatial sources of variation.

Synergies are the least frequent interaction type

Contrary to the common perception that synergies are a frequent consequence of global change (24), our findings suggest that most threats instead interact additively (61 to 100%; Fig. 3). Here, we tested the prevalence of the additive, synergistic, and antagonistic interactions across the global model, as well as system- and taxon-specific models (see Materials and Methods). Across the different models and threats, most multiple threat effects are classified as additive (Fig. 3, figs. S6 to S8, and tables S4 to S6). For instance, when using the global model, additive threats represent 84.90% of interactions involving climate change, 80.60% for disease, 89.20% for exploitation, 85.83% for habitat loss, 78.47% for invasive species, and 85.29% for pollution (Fig. 3 and table S4).

Synergies are not the most common interactive effect across threats.Interactive effects of the different threats. Most interactive effects are additive (purple) or antagonistic. The interactive effects were classified by comparing the estimated population trends for interactive threats versus their additive combinations (additions of the population trends exposed to single threats). They were defined as additive when the trends for populations exposed to the combined threats were equal to the summation of the trends for populations exposed to the individual threats. They were defined as synergistic when the trends for populations exposed to the combined threats were larger than the summation of the trends for populations exposed to the individual threats. They were defined as antagonistic when the trends for populations exposed to the combined threats were smaller than the summation of the trends for populations exposed to the individual threats. No interactions were synergistic in the global model, while 84.90% were additive for climate change, 80.60% for disease, 89.20% for exploitation, 85.80% for habitat loss, 78.40% for invasive species, and 85.30% for pollution.

All nonadditive effects in the global model are antagonistic, accounting for a proportion varying between 10% (exploitation) and 22% (invasive species) of the interactions (Fig. 3, fig. S6, and table S4). Synergies do not appear in this global model, although the amphibian-specific model does display synergistic threats (fig. S7). These findings therefore suggest that, at the macroecological level, synergies are not as prevalent as perceived in the literature (23, 25, 32).

The predominance of antagonistic interactions mirrors previous findings (33), and their prevalence over synergies could be due to multiple reasons [reviewed in (25)]. For example, a single threat may decrease the sensitivity of a species to a second threat [e.g., a threat decreasing the size of a population might make it easier to avoid predation by invasive species (34)]. However, while a high proportion of antagonistic effects could be perceived as a positive outcome, it demonstrates the difficulties of predicting population responses to multiple threats (22–25). The prevalence of the different types of interactive effects has been shown to be scale-dependent (33). The apparent scarcity of synergies at macroecological scales may thus emerge from nonlinear averaging across heterogeneous local conditions, where antagonisms dominate because of compensatory dynamics among populations and stressors.

The observed predominance of antagonisms contrasts with conclusions drawn from meta-analyses of experimental systems, which generally report synergies as the dominant form of nonadditive interaction (23, 25, 32). Our population-level approach, therefore, extends global threat assessments [e.g., IUCN Red List; (3)] by empirically quantifying how interacting pressures manifest in real-world population trajectories rather than inferred extinction risk. This distinction is critical, as it suggests that large-scale biodiversity responses to concurrent threats may be more buffered, or more context-dependent, than previously assumed on the basis of small-scale experiments.

Still, it is worth highlighting that systems and taxa are idiosyncratic in their threat interactions. For example, with amphibians being the only taxonomic group exhibiting synergies as the main nonadditive interactions (fig. S7), this is in line with their increasingly threatened status (1, 3, 35). Disease, exploitation, and invasive species are the threats associated with the most synergies in amphibians (fig. S7), reflecting their vulnerability to these threats (35). Terrestrial systems are also the only system to exhibit synergies (fig. S8), although this is likely the contribution of terrestrial amphibians. Despite these results, our estimated dominance of nonadditive effects is most likely conservative, given the uncertainty in the population trend estimates and our choice of an 80% CI. Varying this interval threshold alters the classification of threats (see Materials and Methods and figs. S6 and S9), but we have selected 80% as a trade-off between the conservation importance of identifying synergies and our confidence in the effect.

Only the mitigation of multiple threats will reverse predicted vertebrate population declines

Although exploitation and habitat loss are not the threats associated with the fastest population declines, our counterfactuals show that management actions aimed at tackling them may be more effective at slowing the predicted rate of decline in vertebrate populations than targeting other individual threats (Fig. 4 and fig. S10). However, addressing any one threat individually is insufficient to reverse vertebrate population declines, whereas the removal of all threats does, accounting for a 226.44% proportional increase in predicted population trends. This is a key indication that single-target interventions are insufficient for reversing global biodiversity declines.

Counterfactuals suggest that lowering exploitation would be the most effective action to slow down the decline of vertebrate populations.The counterfactuals represent the changes in the population trends of the 1740 vertebrate time series affected by threats, had there been no pollution, no disease, no invasive species, no habitat loss, no climate change, no exploitation, or no threats. The continuous, gray line represents when the population trend is 0. The dashed gray line is the median trend of the threatened populations without any intervention. Box plots depict the distribution of population trends, where the thick line in the middle represents the median value; the lower and upper bounds of the box represent the 25th to 75th quartiles, respectively; and the straight lines represent the minimum to the lower quartile (the start of the box) and then from the upper quartile (the end of the box) to the maximum.

Specifically, these counterfactuals adjust for both the frequency and potential impact of the threats on vertebrate populations by simulating counterfactuals where the effects of threats are removed for each of the populations affected by single or multiple threats (see Materials and Methods). Removing the effects of exploitation results in the largest increase in predicted population trends (81.62% relative increase), although the median trend is still negative, followed by habitat loss (16.87%; Fig. 4). The high prevalence of exploitation and habitat loss makes them more likely to interact with other threats. Consequently, despite their lower impact as single threats, mitigating exploitation and habitat loss should remain a conservation priority at a global scale. In contrast, when the invasive species, disease, and pollution threats are removed, there is a smaller effect on the overall global population trends (14.19, 10.49, and 6.94%, respectively; Fig. 4) because of either their low prevalence or low impact. In addition, we show that managing climate change could substantially improve predicted vertebrate populations, highlighting it as a key component of the current biodiversity crisis (14, 36). These findings contrast with previous studies suggesting that the effects of climate change are still modest (1, 2, 37), highlighting the importance that climate policies could have in reversing biodiversity loss (38).

The discrepancy between our results and previous assessments is due to the different methodologies used to rank global threats and highlights the risk of such approaches to prioritize conservation actions (37). Most global threat analyses are based on information from IUCN Red List assessments [e.g. (1, 2, 16)] and are estimated using the presence/absence of threats to species rather than the magnitude of their impacts (5). Therefore, it is expected that global analyses often find habitat loss and exploitation as the most prevalent threats, given the millennia-long history of human’s exploitation and modification of natural ecosystems (39, 40). On the contrary, here, we incorporate both the potential impacts (population trends) and the frequency (their prevalence in our dataset) of threats on natural systems. The IUCN Red List has recently incorporated threat intensity in their criteria (5, 15), which will be key to characterizing the impacts of threats and guiding conservation priorities.

It is worth noting, however, that population trend estimates show large uncertainties, so our findings should be interpreted with caution. In addition, the counterfactuals presented here assumed that the cause of the population trends are threats. However, we are interpreting the counterfactuals as the “influence of threats on predicted trends” only—i.e., they represent tentative conservation efforts. We do not intend these results to provide a ranking of threats to guide conservation but rather showcase the complexity of setting global threat priorities and the importance of accounting for both the impact and frequency of threats.

Our results show that, to achieve global nonnegative vertebrate population trends, we need to mitigate the effects of multiple threats (Fig. 4 and fig. S10). These findings offer evidence of “death by a thousand cuts” (41), where multiple threats over any one single feature are predictive of biodiversity declines. This is notable, as it implies that, to reverse current vertebrate declines, conservation actions and management must target multiple threats simultaneously. For instance, according to our counterfactual tests, the only actions that could reverse global vertebrate declines would be mitigating habitat loss, climate change, and exploitation or habitat loss, invasive species, and exploitation at the same time (fig. S10). Managing these threats is challenging, given their direct links with the development of human societies (42–44). Therefore, halting vertebrate declines requires ambitious conservation actions tackling multiple drivers of biodiversity loss (8, 45).

From a practical perspective, these results call for integrated threat management approaches that coordinate conservation action across sectors, scales, and governance levels. For instance, land-use planning and habitat restoration can reduce habitat loss while also buffering climate impacts (46), fishery reform and supply-chain traceability can curb exploitation pressures while supporting local livelihoods (47), and biosecurity and invasive species management require international coordination across trade and transport networks (48). At the policy level, embedding biodiversity targets within climate and development frameworks can help integrate responses to multiple stressors simultaneously (8). Collectively, these examples illustrate how multithreat mitigation can be operationalized: through comanagement strategies that link local restoration and sustainable use with global-scale emissions and land-use policy. Our results therefore advance the understanding of not only which threats matter the most but also how their management must be coordinated to achieve measurable biodiversity recovery.

There is much polarizing debate as to whether wildlife populations are truly in decline (29, 49), with recent studies suggesting no net loss in vertebrate populations (19, 50) and that conservation efforts have been positive for biodiversity (28). Here, using population-specific threat data and globally distributed time series, we demonstrate that most vertebrate populations exposed to single and multiple threats are in decline across most systems and taxa. For instance, most vertebrate populations exposed to no threats show positive trends, while those exposed to at least one threat either are declining or show no net change (Fig. 2 and table S1). These results highlight the importance of accounting for the impacts of multiple threats on populations to avoid underestimating the magnitude of biodiversity loss (21).

Our analyses are based on the LPD, which is collected from peer-reviewed literature, gray literature, online databases, and data holders (29), so these data come with inherent biases that need to be considered when interpreting global trends. Of the 25,054 population time series in the LPD, only 7827 have information on threats, and only 3129 were monitored for 10 years or more with at least 5 years of data. Our analyses mitigate these biases in our data by accounting for the temporal and spatial correlative effects, as well as disaggregating models by systems and taxa. It is also important to note that our data do not include the strength and/or recurrency of the threats, which can also have a strong influence on population trends (51) and limit previous work that focuses solely on qualitative changes in IUCN classifications [e.g. (15, 43)].

While we could not account for these different properties of threats, we quantify local population trends and how threats could explain the wide variety of trajectories observed in our study. For example, the differences between models fit to time series monitored for 10 and 20 years of data (fig. S11) may result from the predominantly seabird and fishery representation (fig. S12). These populations often display varying threat strengths through time because of fishing effort changes, which in turn influences time series dynamics (47, 52). Such nonlinear dynamics can drive the overall linear estimate of threats to zero (53). Future research targeting the impacts of the different strengths and recurrences of threats on populations will be crucial to further understand the drivers of population declines (54, 55) and will help to develop more accurate biodiversity predictions (56).

Overall, our results highlight the following: (i) Understanding the impacts of threats on the population dynamics of species provides key knowledge to fully characterize and predict the impacts of global change on biodiversity, and (ii) mitigating multiple global threats should remain a conservation priority to halt biodiversity loss. Because of the limited resources allocated to conservation, rankings of the frequency of global threats are often used to prioritize conservation actions (1–3). However, we show that the main threats affecting vertebrate populations are context-dependent (57, 58) and that accounting for both the impacts and frequency of threats is crucial to developing effective management actions. In particular, multiple threats interacting at the local scale can accelerate biodiversity loss (compared to single threats), and focusing efforts on the mitigation of single threats might not be enough to halt population declines.

MATERIALS AND METHODS

Data

To measure the association between population trends and the different threats, both in isolation and with other stressors, we used one of the largest global population-monitoring databases currently available, the LPD (29). The LPD includes 25,054 population time series of 4392 species, and each time series has repeated monitoring surveys of the population abundance in a given area. These data are collated from a variety of sources, primarily from peer-reviewed literature but also from gray literature, online databases, and data holders (29). This information was digitized by the LPD team, following a standardized protocol (see details below).

Data sources containing population time series for vertebrate species were collated by either scanning articles from journal issues within the conservation biology, wildlife management, and ecology disciplines (this yielded most data sources); conducting keyword searches within academic and generic search engines; or contacting individual data holders to share data directly with the team (this yielded the least data sources). Data sources were selected if they met the criteria for inclusion in the LPD: single vertebrate species abundance estimates from a multiannual survey, monitored using a consistent method and effort (59) in the same location. Supporting texts in the publication, especially sections on limitations, were read to ensure that the abundance estimates were reliable and not an artifact, for example, an increasing trend resulting from increased effort would not be included unless a correction factor was applied. Ancillary information on the survey location, ecology, and human activity (conservation interventions and threats) was collected from the data source alongside the abundance data. To ensure consistency, trained personnel recorded information from the original data sources using a set of guidelines. Once a person had been trained, all data entered were subject to an additional check by an experienced member of the LPD team to ensure quality control and promote consistency in decisions made.

Of the 25,054 population time series making up the LPD (including confidential records), 7827 population time series contain information on whether the populations were exposed to threats. On the basis of information from the data source, for each publication, it was first identified whether the population was threatened or not threatened or whether its threat status was unknown. Threats were identified as direct or indirect human activities or processes that affected the populations for at least 50% of the surveyed years according to the original source of the time series. Information on threat severity was not recorded as it was rarely available, and any mention of threats that was speculative in nature or only identified potential threats was not included. If the population was threatened, the number of threats the population was exposed to was recorded, from one to three. The information within the data sources was sometimes quantitative, e.g., stating the number of individuals hunted annually, but most often it was reported in a qualitative way, e.g., describing a general pattern of hunting that affects the populations. For this reason, and because the impact of the threat was rarely quantified in the data sources, broad categories describing the threat to the population were recorded: climate change, invasive species, habitat loss/degradation, exploitation, pollution, and diseases, following the Red List threat classification (60). As a guide, the subdivisions of the Red List threat classification were used to assign a threat to a broad category. We decided to use these broad categories to be able to compare the results of our study with those used in previous global biodiversity threat assessments (1, 37). Where there was more than one category option for a threat, the threat that better described the impact on the population was selected. For example, for amphibians, we would list the threat of chytridiomycosis as a disease rather than as an invasive species (because of Batrachochytrium dendrobatidis). For species such as Panthera tigris, which are threatened by a loss of prey base from hunting or competition with domestic species, we considered that habitat loss.

To ensure that our data could be associated with the potential impacts of threats on population trends, we only included time series containing information on whether the population was exposed to threats at the time of the study (according to the original source). The number of threats populations were exposed to ranges from zero to three, as populations exposed to four or more threats were rarely found in the literature. To have sufficient information to appropriately capture the directional trends in abundance, we only included populations with at least 10 years of monitoring data and a minimum of five data points (20, 61). We tested the implications of this choice of data quality by running additional analyses with ≥10 and ≥20 years of data. The influence of threats tended toward zero as the time series length increased (fig. S11), with the coverage of threats, systems, and taxa strongly affected (fig. S12). This is likely due to threats and disturbances not acting linearly through time (62), leading to a net zero linear trend (53). The longest time series were also heavily biased toward marine birds and fish, where ~33% of time series are unthreatened in our dataset (fig. S11). We therefore present the results of the analyses run on data with ≥10 years and with a minimum of five data points to minimize system or taxon biases and include time series with persistent threat effects.

Overall, our selection process resulted in 3129 time series of 1281 species, including amphibians (83), birds (1195), fishes and elasmobranchs (947), mammals (698), and reptiles (206). These time series covered all the continents (Fig. 1) and freshwater (616), marine (1271), and terrestrial (1242) systems. The duration of datasets varied between 10 and 65 years, with a median duration of 25 years (Fig. 1), covering a period between 1950 and 2019. Although the original dataset included data up to 2020, none of our selected time series reached that year.

Estimates of population trends

To estimate population trends, we fit a multilevel Bayesian linear model between the natural logarithm of abundance and year and allowed year to interact with each threat/possible combination of threats. To identify the association between abundance and independent/interacting threats, we coded threats and threat combinations as binary factors. This resulted in a model matrix consisting of 36 columns representing the unique combinations of climate change, invasive species, habitat loss/degradation, exploitation, pollution, and disease. For example, if a population experienced exploitation and disease but not pollution, its model matrix would be as follows: exploitation = 1, disease = 1, pollution = 0, exploitation.disease = 1, exploitation.pollution = 0, disease.pollution = 0, exploitation.disease.pollution = 0. We consequently can decompose the additive effects of exploitation and disease (exploitation = 1, disease = 1) from their interaction (exploitation.disease = 1). Defining threats in this way is necessitated by the qualitative threat categorization described above and allows counterfactual analyses.

To account for temporal nonindependence, we modeled the population-level time series with a discrete autoregressive-1 (ar1) temporal process, which assumes that neighboring abundance observations within a time series will be more similar. Furthermore, to capture spatial and phylogenetic correlative nonindependence, we fit covariance structures between species identities and sites (30). Phylogeny was modeled via uncorrelated slopes between species and populations, while shared correlations between sites were fit via a shared slope covariance structure. This site covariance matrix was derived using the Haversine (spherical) distance between each site. To improve run time in the largest datasets, we rounded site coordinates to the nearest integer, i.e., a latitude of 10.65 was set to 11. Rounding in this way has a limited effect on the LPD (fig. S13), given the large geographical distances between sites. We normalized this spatial matrix between 0 and 1, with values close to 1 indicating neighboring sites, while values approaching 0 indicate distant sites.

Year and abundance were centered at zero by subtracting the mean time series year/abundance from each time series’ year/abundance. This parameterization effectively fixes intercepts at 0 for each slope and accounts for random intercept variation without increasing the number of estimated parameters. Consequently, no random intercepts were necessary.

The final model structure was as follows

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Specifically, we assumed that abundance (y) is a Gaussian normally distributed variable with mean μ and variance σ, consisting of n observations. X represents the fixed-effect model matrix containing year and threat variables, with β being the associated fixed-effect coefficients. γ is therefore the random effect coefficients, which are controlled by the independent random slopes for species (r^S^, index i) and population (r^P^, index j) and the correlated site (r^L^, index k) random slopes. Because of challenges fitting the random hyperparameters, the model was reparametrized into a “noncentered” form to prevent correlations between hyperparameters (63). Consequently, rather than modeling slopes varying around the overall year coefficient, we introduce a standard normal parameter for each random effect (z^S^, z^P^, and z^L^), which is then multiplied by each random effect’s respective sigma hyperprior (σ*^S^, σ^P^, and σ^L^*). r^L^ (the correlated site slopes) is also dependent upon the additional variance-covariance Σ, which specifies that covariance is present in neighboring sites. Σ is populated by the Haversine distances between each site. Next, noncentered parameterization allows temporal nonindependence to be modeled by introducing the correlation (ρ) between neighboring residuals (ε) in the linear predictor rather than the observational error term (σ) as is commonly parameterized. Last, we specify that all time series, and in turn species, share the same autocorrelation parameter (ρ). This is necessary given the size of the global dataset, but future work should consider varying ρ by a relevant cluster (such as phylogeny or location). Each realm- and taxon-specific model (see below) does estimate a unique ρ, which partially mitigates such a dataset limitation.

This model structure was then fit across subsets of the LPD to disentangle differential responses of taxa and ecosystems to threats. First, a global model was fit across all appropriate time series to give insight into general trends and relationships between abundance and threats. Second, system-level models were fit upon freshwater, marine, and terrestrial time series alone. Last, taxon-specific models were fit upon time series classified as amphibians, birds, fish, mammals, or reptiles. The category “fishes” incorporated the taxonomic groups Holocephali, Elasmobranchii, Myxini, Cephalaspidomorphi, Actinopterygii, and Sarcopterygii. It was necessary to fit models by system and/or taxon because of the large number of parameters estimated in each model. For example, the inclusion of a taxon variable with five levels (amphibian, bird, fish, mammal, and reptile) would increase the number of fixed-effect coefficients in the global model from 74 (the 36 threat combinations, year, the interaction between year and each threat combination, and global intercept) to 370. The model structure was identical across models, with the spatial variance-covariance matrix Σ recalculated for each data subset.

All models were fit using the brms package version 2.22.0 (64) in R version 4.4.0 (65). Models were run for 5000 iterations over four chains, with a warmup of 2500 iterations. This resulted in a total of 20,000 draws per model. All distributions reported in this study were constructed from 1000 random draws of the posterior distribution.

Multilevel Bayesian model diagnostics

To check the validity of our multilevel Bayesian models, we ran a set of diagnostics. We inspected model convergence by visually examining trace plots and using Rhat values (the ratio of the effective sample size to the overall number of iterations, with values close to one indicating convergence).

Furthermore, we evaluated the model fit by exploring the distribution of the residuals, their variance, their autocorrelation, and the posterior predictive checks (figs. S14 to S17). When the model fits the data, we would expect the residuals to follow a Gaussian distribution and to show constant variance, as shown in our models (figs. S14 to S17). In addition, the posterior predictive checks compare the distribution of the data with the predictions from the model, so if the model is well fitted, the predictions should overlap the data, as shown in our models (fig. S16). Last, our models displayed no evidence of autocorrelation in their residuals (fig. S17).

Influence of threats on population trends

To quantify the trends of populations exposed to single and multiple threats on population trends, we examined coefficients and extracted predictions made by the fitted multilevel Bayesian models. As each threat factor was coded as binary, the year coefficient represents trends in the absence of threats, and each threat/threat interaction coefficient is therefore the moderation of this no-threat trend.

In addition, to estimate the ultimate association of threats with population trends, we estimated the conditional effect for each combination of threats for each year contained within the LPD. This generates a predictive posterior distribution of year-on-year log abundances, which we then converted to a trend by taking the mean derivative (also known as the slope) using the following equation on a draw-by-draw basis

$$[eqn]$$

Henceforth, the term “trend” refers to the derivative of the predicted slope rather than the raw coefficients.

The inclusion of both random and fixed effects in our analysis has the potential to confound predicted trends. We therefore validated our predictions by comparing the variability of trends predicted using random effects alone versus those made with the inclusion of threat fixed effects. To achieve this, we generated null model predictions by estimating the conditional trends (i.e., random effects set to 0) for populations where all threats were “removed” and set to 0. The null trend posterior was then compared to two separate scenarios, estimating the contribution of random versus fixed effects to the observed trends. Random effect predictions were made upon the null dataset but included all random effect identities (species, population, site, etc.). We then averaged these predictions across all time series draws to estimate a marginal trend posterior. Conversely, fixed-effect contributions were estimated via conditional trends with all threats “present” and set to 1. Each of these predictions was then differenced from the null model posterior to create a distribution of trend “contributions.” If the trend difference is zero, then that effect is not contributing to the direction of the trend. We anticipate that random effects will not influence the direction of trends but increase the uncertainty of estimates (30).

Estimation of interactive versus additive effects of threats

We also classified the form of interactive threat effects in a separate analysis using the abovementioned multilevel Bayesian model. By including or excluding threat data, we estimated conditional trends given independent and/or interacting threats and used the difference in trends to classify interaction type. We considered threat combinations at the initial stages as “singular” (they are the only threat acting upon a time series) or “interactive” (they are one of two or more threats acting on a time series simultaneously). We do not know at this stage what form of interaction an interactive threat combination (i.e., a multiple stressor scenario) is. Only once we have estimated our model, which fits a statistical interaction, can we classify what form of interaction those interacting threats represent: a noninteraction (i.e., additive), negative (antagonistic), or positive (synergistic). We define the different types of interactions in Table 1.

In our model, the additive contribution of threats to the trend can be predicted by including the independent threat columns in the threat model matrix solely. For example, if we were interested in the additive effect of exploitation and disease, we would code the model matrix as follows: exploitation = 1, disease = 1, pollution = 0, exploitation.disease = 0, exploitation.pollution = 0, disease.pollution = 0, exploitation.disease.pollution = 0. Exploitation and disease are therefore the only contributing threats to the conditional trend predictions, with their interactions not considered. Consequently, we can alter the model matrix to include exploitation and disease’s interaction by “turning on” the interactive column: exploitation = 1, disease = 1, pollution = 0, exploitation.disease = 1, exploitation.pollution = 0, disease.pollution = 0, exploitation.disease.pollution = 0. We therefore have two posteriors representing the additive and interactive contributions of threats to the population trend.

There is still an ongoing debate on the best approaches to classify interactive effects of multiple threats (25, 66, 67). In this study, we considered that the null model was the simple additive effect (sum of the individual effects), because it was the most conservative approach on the basis of the scale of our study (67). We included organisms from different taxonomic groups and ecosystems, and given that we lacked information about the mechanistic effects of each threat, we decided to maintain a simple additive null effect (67). Consequently, we differenced the interacting threat posterior from the additive null and used the direction of difference to classify threat contributions as additive, antagonistic, or synergistic (fig. S1). If the 80% CI for the differenced posterior transgressed zero, then the threat combination was classified as additive. We assume that threats contribute negatively to trends, and so, if the difference in trend is negative, we classified that threat combination as synergistic—the interaction of threats generates a trend posterior more negative than the additive threat posterior. Conversely, if the difference in trend is positive, we classified that threat combination as antagonistic as the interacting threat posterior is more positive than the additive equivalent. See fig. S1 for a visualization of this process.

It is worth mentioning that our approach to quantify the interactive effects assumes that the population trends measured are predicted by the threats. Nevertheless, in this study, we did not perform any causal analysis, which could establish causal links between the threats and the observed population trends. It is possible that other confounding factors may have played a role in shaping these population trends, but it was out of the scope of this manuscript to establish these causal links. As a result, we advise caution when interpreting the results derived from our analyses of interactive effects, given that we assume causal links that were not formally tested.

Counterfactuals

To explore what would have happened to vertebrate population predictions if threat contributions were removed, we developed six counterfactuals. Each counterfactual represents the predicted population trends where the potential effect of the different single threats (exploitation, climate change, habitat loss, disease, invasive species, and pollution) was removed from the populations exposed to any number of threats. To achieve this, we took the following steps:

1. We first dropped all populations in our dataset that were unthreatened.

2. We then modified the model matrix for each threat in turn by setting all columns where that threat is present to 0. For example, a counterfactual removing exploitation from a population experiencing exploitation and disease would be coded as follows: exploitation = 0, disease = 1, pollution = 0, exploitation.disease = 0, exploitation.pollution = 0, disease.pollution = 0, exploitation.disease.pollution = 0 (exploitation is removed, and disease remains, but the interaction is also removed).

3. We then estimated the trend for each population, given the counterfactual model matrix, to create a distribution of population trends given the individual removal of each threat.