All is not decline across global vertebrate populations

The populations of Earth’s species are changing over time in complex ways, creating a mixture of winners and losers in a time of accelerating global change. A critical research challenge is to test if there are specific biomes, taxa or types of species that are experiencing rapid alterations in abundance over time. We conducted an analysis of the Living Planet Database using a state-space modelling approach including nearly 10,000 abundance time-series from over 2,000 vertebrate species. We integrated the population abundance data with information on geographic range, habitat preference, taxonomic and phylogenetic relationships, conservation status and threats from occurrence, phylogenetic and conservation assessment data. We found that 15% of populations were declining, 18% were increasing, and 67% showed no net changes in abundance over time. Against a backdrop of no biogeographic and phylogenetic patterning in population change, we uncovered a distinct taxonomic signal. Amphibians were the only taxa that experienced net declines in the analyzed data, while birds, mammals and reptiles on average became more abundant over time. The continuum of abundance changes over time was poorly captured by species’ rarity and global-scale threats. Capturing the full spectrum of abundance change of species around the planet will inform conservation efforts and improve projections of biodiversity and ecosystem function change under the accelerating influence of the Anthropocene.

Page 6 of 44 globally is explained by differences in species' rarity and IUCN conservation status. We measured 116 rarity using three separate metrics -geographic range, mean population size (number of animals 117 that were recorded by monitoring) and habitat specificity. In a post-hoc analysis, we categorized   Figure S1). In the second stage, we modelled the population trend and fluctuation estimates from 132 the first stage across latitude, realm, biome, taxa, rarity metrics, phylogenetic relatedness, 133 species' conservation status and threat type using a Bayesian modelling framework ( Figure S2).

134
We included a species random intercept effect to account for the possible correlation between 135 the trends of populations from the same species (see table Table S1 for sample sizes). As     195 dip test, D = 0.04, p < 0.01, Figure 1a). Bony fish population trends were centered on zero (slope 196 = -0.001, CI = -0.004 to 0.002, Figure 1a-b) and sharks and rays showed net declines, but the 197 credible intervals overlapped zero (slope = -0.01, CI = -0.02 to 0.01). Fluctuations were most 198 common for amphibian populations (slope = 0.04, CI = 0.036 to 0.049, Figure 2d), which were 199 monitored for the shortest time period on average (11 years, Figure S1,

248
Critically endangered: slope = 0.035, CI = 0.028 to 0.041), but not more population declines, with 249 increasing conservation status ( Figure 4, Table S2). Populations from species that were exposed    Table S2 for model outputs.

335
Contrary to our initial predictions, we did not find a distinct geographic patterning of population 336 change around the world, nor a consistent trend of increasing declines in population abundance with increasing threat level (Figures 1 and 4). Similar lack of biogeographic signal has been 338 documented in regional studies of population change from the Netherlands 33 and in temperate  we use species' Red List conservation status as a proxy for extinction risk 75 , an approach used 487 recently in other studies 76 . We extracted the number and types of threats that each species is 488 exposed to globally from their respective IUCN profiles 75 .

490
Quantifying population trends and fluctuations over time

491
In the first stage of our analysis, we used state-space models that model abundance (scaled to a 492 common magnitude between zero and one) over time to calculate the amount of overall  Table S5 for species list), we 530 calculated a km 2 occurrence area based on species occurrence data from GBIF 73 . Extracting and 531 filtering GBIF data and calculating range was computationally-intensive and occurrence data 532 availability was lower for certain species; thus, we did not estimate geographic range from GBIF 533 data for all species part of the Living Planet Database. Instead, we focused on analyzing range 534 effects for birds and mammals globally, as they are a very well-studied taxon and for species

535
Page 26 of 44 monitored in the UK, a country with intensive and detailed biodiversity monitoring of vertebrate 536 species. We did not use IUCN range maps, as they were not available for all of our study species, 537 and previous studies using GBIF occurrences to estimate range have found a positive correlation 538 between GBIF-derived and IUCN-derived geographic ranges 82 .

540
For the geographic ranges of species monitored in the UK, we calculated range extent using a 541 minimal convex hull approach based on GBIF 73 occurrence data. We filtered the GBIF data to We calculated mean size of the monitored population, referred to as population size, across the 554 monitoring duration using the raw abundance data, and we excluded populations which were not 555 monitored using population counts (i.e., we excluded indexes).

558
To create an index of habitat specificity, we extracted the number of distinct habitats a species 559 occupies based on the IUCN habitat category for each species, accessed through the package 560 Page 27 of 44 rredlist 85 . We also quantified habitat specificity by surveying the number of breeding and non-561 breeding habitats for each species from its online IUCN profile (the 'habitat and ecology' section).

562
The two approaches yielded similar results ( Figure S3, Table S2, key for the profiling method is 563 presented in Table S4). We obtained global conservation status and threat data for all study 564 species through the IUCN Red List classification 75 .

566
Testing the sources of variation in population trends and fluctuations

567
In the second stage of our analyses, we modelled the trend and fluctuation estimates from the 568 first stage across latitude, realm, biome, taxa, rarity metrics, phylogenetic relatedness, species' 569 conservation status and threat type using a Bayesian modelling framework through the package 570 MCMCglmm 86 . We included a species random intercept effect in the Bayesian models to account 571 for the possible correlation between the trends of populations from the same species (see Table   572 S1 for sample sizes). The models ran for 120 000 iterations with a thinning factor of ten and a  where β0 is the global intercept (β0 = 1), β0l is the phylogeny-level departure from β0 (phylogeny 618 random effect); yi,k,m is the estimate for change in population abundance for the ith population 619 time-series for the kth species with the mth phylogenetic distance. 620

621
To account for phylogenetic uncertainty, for each class, we ran ten models with identical 622 structures but based on different randomly selected phylogenetic trees. We report the mean   types and number of threats to which species are exposed, we conducted a post-hoc analysis of 650 trends and fluctuations across threat type (categorical effect, determined by species IUCN threat 651 categorization) and number of threats that each species is exposed to across its range (in 652 separate models). The model structures were identical to those presented above, except for the fixed effect which was a categorical IUCN conservation status variable. The analytical workflow 654 of our analyses is summarized in conceptual diagrams (Figures S1 and S2) and all code is 655 available on GitHub (link to public repository to be added at time of publication).   (Figures 2 and S5). On a larger scale, the Living Planet Database 667 under-represents populations outside of Europe and North America and over-represents common 668 and well-studied species 59 . We found that for the populations and species represented by current 669 monitoring, rarity does not explain variation in population trends, but we note that the relationship 670 between population change and rarity metrics could differ for highly endemic specialist species 671 that are not included in the Living Planet Database. As ongoing and future monitoring begins to 672 fill in the taxonomic and geographic gaps in existing datasets, we will be able to re-assess and 673 test the generality of the patterns of population change across biomes, taxa, phylogenies, species  has been found to be enough to detect directional trends in 80-90% of cases 68 . In a separate 693 analysis, we found significant lags in population change following disturbances (forest loss) and 694 that population monitoring often begins decades to centuries after peak forest loss has occurred 695 at a given site 41