Rare and common vertebrates span a wide spectrum of population trends
Workflow
All data syntheses, visualisation and statistical analyses were conducted using R version 3.5.171. For conceptual diagrams of our workflow, see Supplementary Figs. 1 and 2. Effect sizes plotted on graphs were standardised by dividing the effect size by the standard deviation of the corresponding input data.
Population data
To quantify vertebrate population change (trends and fluctuations), we extracted the abundance data for 9286 population time series from 2084 species from the publicly available Living Planet Database72 (http://www.livingplanetindex.org/data_portal) that covered the period between 1970 and 2014 (Supplementary Table 1). These time series represent repeated monitoring surveys of the number of individuals in a given area, hereafter, called ‘populations’. Monitoring duration differed among populations, with a mean duration of 23.9 years and a mean sampling frequency of 23.3 time points (Supplementary Fig. 3, see Supplementary Figs. 6 and 7 for effects of monitoring duration on detected trends). In the Living Planet database, 17.9% of populations were sampled annually or in rare cases multiple times per year. The time series we analysed include vertebrate species that span a large variation in age, generation times and other demographic-rate processes. For example, from other work that we have conducted, we have found that when generation time data were available (~50.0% or 484 out of 968 bird species, and 15.6% or 48 out of 306 mammal species), the mean bird generation time is 5.0 years (min = 3.4 years, max = 14.3 years) and the mean mammal generation time is 8.3 years (min = 0.3 years, max = 25 years)45. Thus, we believe that most vertebrate time series within the LPD capture multiple generations.
In our analysis, we omitted populations which had less than five time points of monitoring data, as previous studies of similar population time series to the ones we analysed have found that shorter time series might not capture directional trends in abundance63. Populations were monitored using different metrics of abundance (e.g., population indices vs number of individuals). Before analysis, we scaled the abundance of each population to a common magnitude between zero and one to analyse within-population relationships to prevent conflating within-population relationships and between-population relationships73. Scaling the abundance data allowed us to explore trends among populations relative to the variation experienced across each time series.
Phylogenetic data
We obtained phylogenies for amphibian species from https://vertlife.org4, for bird species from https://birdtree.org8, and for reptile species from https://vertlife.org6. For each of the three classes (Amphibia, Aves and Reptilia), we downloaded 100 trees and randomly chose 10 for analysis (30 trees in total). Species-level phylogenies for the classes Actinopterygii and Mammalia have not yet been resolved with high confidence74,75.
Rarity metrics, IUCN Red List categories and threat types
We defined rarity following a simplified version of the ‘seven forms of rarity’ model76, and thus consider rarity to be the state in which species exist when they have a small geographic range, low population size, or narrow habitat specificity. We combined publicly available data from three sources: (1) population records for vertebrate species from the Living Planet Database to calculate mean population size, (2) occurrence data from the Global Biodiversity Information Facility77 (https://www.gbif.org) and range data from BirdLife78 (http://datazone.birdlife.org) to estimate geographic range size and (3) habitat specificity and Red List Category data for each species from the International Union for Conservation79 (https://www.iucnredlist.org). The populations in the Living Planet Database72 do not include species that have gone extinct on a global scale. We extracted the number and types of threats that each species is exposed to globally from their respective species’ IUCN Red List profiles79.
Quantifying population trends and fluctuations over time
In the first stage of our analysis, we used state-space models that model abundance (scaled to a common magnitude between zero and one) over time to calculate the amount of overall abundance change experienced by each population (μ,40,80). State-space models account for process noise (σ2) and observation error (τ2) and thus deliver robust estimates of population change when working with large data sets where records were collected using different approaches, such as the Living Planet Database41,81,82. Previous studies have found that not accounting for process noise and measurement error could lead to over-estimation of population declines83, but in our analyses, we found that population trends derived from state-space models were similar to those derived from linear models. Positive μ values indicate population increase and negative μ values indicate population decline. State-space models partition the variance in abundance estimates into estimated process noise (σ2) and observation or measurement error (τ2) and population trends (μ):
$$X_t = X_{t-1} + mu + varepsilon _t,$$
(1)
where Xt and Xt−1 are the scaled (observed) abundance estimates (between 0 and 1) in the present and past year, with process noise represented by εt~ gaussian(0, σ2). We included measurement error following:
$$Y_t = X_t + F_t,$$
(2)
where Yt is the estimate of the true (unobserved) population abundance with measurement error:
$$F_tsim gaussianleft( {0,,{it{T}}^2} right)$$
(3)
We substituted the estimate of population abundance (Yt) into Eq. 1:
$$Y_{t} = {it{X}}_{{it{t}} – 1} + mu + varepsilon _{it{t}} + {it{F}}_{it{t}}.$$
(4)
Given
$${it{X}}_{{it{t}} – 1} = {it{Y}}_{{it{t}} – 1} – {it{F}}_{{it{t}} – 1}$$
(5)
then:
$${it{Y}}_{it{t}} = {it{Y}}_{t – 1} + mu + varepsilon _t + F_t – F_{t – 1}$$
(6)
For comparisons of different approaches to modelling population change, see ‘Comparison of modelling approaches section’.
Quantifying rarity metrics
We tested how population change varied across rarity metrics—geographic range, mean population size and habitat specificity – on two different but complementary scales. In the main text, we presented the results of our global-scale analyses, whereas in the SI, we included the results when using only populations from the UK—a country with high monitoring intensity, Thus, we quantified rarity metrics for species monitoring globally and in the UK. The three rarity metrics used in this study were weakly correlated at both UK and global scales (Supplementary Fig. 11).
Geographic range
To estimate geographic range for bird species monitored globally, we extracted the area of occurrence in km2 for all bird species in the Living Planet Database that had records in the BirdLife Data Zone78. For mammal species’ geographic range, we used the PanTHERIA database84 (http://esapubs.org/archive/ecol/E090/184/). To estimate geographic range for bony fish, birds, amphibians, mammals and reptiles monitored in the UK (see Supplementary Table 5 for species list), we calculated a km2 occurrence area based on species occurrence data from GBIF77. Extracting and filtering GBIF data and calculating range was computationally intensive and occurrence data availability was lower for certain species. Thus, we did not estimate geographic range from GBIF data for all species part of the Living Planet Database. Instead, we focused on analysing range effects for birds and mammals globally, as they are a very well-studied taxon and for species monitored in the UK, a country with intensive and detailed biodiversity monitoring of vertebrate species. We did not use IUCN range maps, as they were not available for all of our study species, and previous studies using GBIF occurrences to estimate range have found a positive correlation between GBIF-derived and IUCN-derived geographic ranges85.
For the geographic ranges of species monitored in the UK, we calculated range extent using a minimal convex hull approach based on GBIF occurrence data77. We filtered the GBIF data to remove invalid records and outliers using the CoordinateCleaner package86. We excluded records with no decimal places in the decimal latitude or longitude values, with equal latitude or longitude, within a one-degree radius of the GBIF headquarters in Copenhagen, within 0.0001 degrees of various biodiversity institutions and within 0.1 degrees of capital cities. For each species, we excluded the lower 0.02 and upper 0.98 quantile intervals of the latitude and longitude records to account for outlier points that are records from zoos or other non-wild populations. We drew a convex hull to most parsimoniously encompass all remaining occurrence records using the chull function, and we calculated the area of the resulting polygon using areaPolygon from the geosphere package87.
Mean size of monitored populations
We calculated mean size of the monitored population, referred to as population size, across the monitoring duration using the raw abundance data, and we excluded populations, which were not monitored using population counts (i.e., we excluded indexes).
Habitat specificity
To create an index of habitat specificity, we extracted the number of distinct habitats a species occupies based on the IUCN habitat category for each species’ profile, accessed through the package rredlist88. We also quantified habitat specificity by surveying the number of breeding and non-breeding habitats for each species from its online IUCN species profile (the ‘habitat and ecology’ section). The two approaches yielded similar results (Supplementary Fig. 10, Supplementary Table 2, key for the profiling method is presented in Supplementary Table 6). We obtained global IUCN Red List Categories and threat types for all study species through their IUCN Red List profiles79.
Testing the sources of variation in population trends and fluctuations
In the second stage of our analyses, we modelled the trend and fluctuation estimates from the first stage analyses across latitude, realm, biome, taxa, rarity metrics, phylogenetic relatedness, species’ IUCN Red List Categories and threat type using a Bayesian modelling framework through the package MCMCglmm89. We included a species random intercept effect in the Bayesian models to account for the possible correlation between the trends of populations from the same species (see Supplementary Table 1 for sample sizes). The models ran for 120,000 iterations with a thinning factor of ten, a burn-in period of 20,000 iterations and the default one chain. We assessed model convergence by visually examining trace plots. We used weakly informative priors for all coefficients (an inverse Wishart prior for the variances and a normal prior for the fixed effects):
$$Prleft( mu right) sim Nleft( {0,,10^8} right)$$
(7)
$$Pr(sigma ^2) sim Inverse,Wishart,left( {V = 0,,nu = 0} right)$$
(8)
Population trends and fluctuations across latitude, biomes, realms and taxa
To investigate the geographic and taxonomic patterns of population trends and fluctuations, we modelled population trends (μ) and population fluctuations (σ2), derived from the first stage of our analyses (state-space models), as a function of (1) latitude, (2) realm (freshwater, marine, terrestrial), (3) biome (as defined by the ‘biome’ category in the Living Planet Database, e.g., ‘temperate broadleaf forest’90 and (4) taxa (Actinopterygii, bony fish; Elasmobranchii, sharks and rays; Amphibia, amphibians; Aves, birds; Mammalia, mammals; Reptilia, reptiles). We used separate models for each variable, resulting in four models testing the sources of variation in trends and four additional models focusing on fluctuations. Each categorical model from this second stage of our analyses was fitted with a zero intercept to allow us to determine whether net population trends differed from zero for each of the categories under investigation. The model structures for all models with a categorical fixed effect were identical with the exception of the identity of the fixed effect, and below we describe the taxa model:
$$mu _{i,j,k} = beta _0 + beta _{0,j} + beta _1 ast taxa_{i,j,k},$$
(9)
$$y_{i,j,k}sim gaussianleft( {mu _{i,j,k},sigma ^2} right),$$
(10)
where taxai,j,k is the taxa of the ith time series from the jth species; β0 and β1 are the global intercept (in categorical models, β0 = 1) and the slope estimate for the categorical taxa effect (fixed effect), β0j is the species-level departure from β0 (species-level random effect); yi,j,k is the estimate for change in population abundance for the ith population time series from the jth species from the kth taxa.
Population trends and fluctuations across amphibian, bird and reptile phylogenies
To determine whether there is a phylogenetic signal in the patterning of population change within amphibian, bird and reptile taxa, we modelled population trends (μ) and fluctuations (σ2) across phylogenetic and species-level taxonomic relatedness. We conducted one model per taxa per population change variable—trends or fluctuations using Bayesian linear mixed effects models using the package MCMCglmm89. We included phylogeny and taxa as random effects. The models did not include fixed effects. We assessed the magnitude of the random effects (phylogeny and species) by inspecting their posterior distributions, with a distribution pushed up against zero indicating lack of effect, as these distributions are always bounded by zero and have only positive values. We used parameter-expanded priors, with a variance-covariance structure that allows the slopes of population trend (the μ values from the first stage analysis using state-space models) to covary for each random effect. The prior and model structure were as follows:
$$Prleft( mu right) sim Nleft( {0,,10^8} right),$$
(11)
$$Prleft( {sigma ^2} right) sim Inverse,Wishart,left( {V = 1,,nu = 1} right),$$
(12)
$$mu _{i,k,m} = beta _0 + beta _{0,k} + beta _{0,m},$$
(13)
$$y_{i,k,m} sim gaussianleft( {mu _{i,k,m},,sigma ^2} right),$$
(14)
where β0 is the global intercept (β0 = 1), β0l is the phylogeny-level departure from β0 (phylogeny random effect); yi,k,m is the estimate for change in population abundance for the ith population time series for the kth species with the mth phylogenetic distance.
To account for phylogenetic uncertainty for each class, we ran ten models with identical structures, but based on different randomly selected phylogenetic trees. We report the mean estimate and its range for each class.
Population trends and fluctuations across rarity metrics
To test the influence of rarity metrics (geographic range, mean population size and habitat specificity) on variation in population trends and fluctuations, we modelled population trends (μ) and fluctuations (σ2) across all rarity metrics. We conducted one Bayesian linear model per rarity metric per scale (for both global and UK analyses) per population change variable—trends or fluctuations. The response variable was population trend (μ values from state-space models) or population fluctuation (σ2 values from state-space models), and the fixed effects were geographic range (log transformed), mean population size (log transformed) and habitat specificity (number of distinct habitats occupied). The model structures were identical across the different rarity metrics and below we outline the equations for population trends and geographic range:
$$mu _{i,k,n} = beta _0 + beta _{0,k} + beta _1 ast geographic,range_{i,k,n},$$
(15)
$$y_{i,k,n} sim gaussianleft( {mu _{i,k,n},,sigma ^2} right),$$
(16)
where geographic rangei,k,n is the logged geographic range of the kth species in the ith time series; β0 and β1 are the global intercept and slope estimate for the geographic range effect (fixed effect), β0j is the species-level departure from β0 (species-level random effect); yi,k,n is the estimate for change in population abundance for the ith population time series from the jth species with the nth geographic range.
Population trends across species’ IUCN Red List Categories
To investigate the relationship between-population change and species’ Red List Categories, we modelled population trends (μ) and fluctuations (σ2) as a function of IUCN Red List Categories (categorical variable). We conducted one Bayesian linear model per population change metric per scale (for both global and UK analyses). To test variation in population trends and fluctuations across the types and number of threats to which species are exposed, we conducted a post hoc analysis of trends and fluctuations across threat type (categorical effect) and number of threats that each species is exposed to across its range (in separate models). The model structures were identical to those presented above, except for the fixed effect which was a categorical IUCN Red List Category variable.
The analytical workflow of our analyses is summarised in conceptual diagrams (Supplementary Figs. 1 and 2) and all code is available on GitHub (https://github.com/gndaskalova/PopChangeRarity, DOI 10.5281/zenodo.3817207).
Data limitations: taxonomic and geographic gaps
Our analysis is based on 9286 monitored populations from 2084 species from the largest currently available public database of population time series, the Living Planet Database72. Nevertheless, the data are characterised by both taxonomic and geographic gaps that can influence our findings. For example, there are very few population records from the Amazon and Siberia (Fig. 1b)—two regions currently undergoing rapid environmental changes owing to land-use change and climate change, respectively. In addition, birds represent 63% of all population time series in the Living Planet Database, whilst taxa such as amphibians and sharks where we find declines are included with fewer records (Fig. 2 and Supplementary Fig. 4). On a larger scale, the Living Planet Database under-represents populations outside of Europe and North America and over-represents common and well-studied species62. We found that for the populations and species represented by current monitoring, rarity does not explain variation in population trends, but we note that the relationship between population change and rarity metrics could differ for highly endemic specialist species or species different to the ones included in the Living Planet Database17. As ongoing and future monitoring begins to fill in the taxonomic and geographic gaps in existing datasets, we will be able to reassess and test the generality of the patterns of population change across biomes, taxa, phylogenies, species traits and threats.
Data limitations: monitoring extent and survey techniques
The Living Planet Database combines population time series where survey methods were consistent within time series but varied among time series. Thus, among populations, abundance was measured using different units and over varying spatial extents. There are no estimates of error around the raw population abundance values available and detection probability likely varies among species. Thus, it is challenging to make informed decisions about baseline uncertainty in abundance estimates without prior information. We used state-space models to estimate trends and fluctuations to account for these limitations as this modelling framework is particularly appropriate for analyses of data collected using disparate methods41,81,82. Another approach to partially account for observer error that has been applied to the analysis of population trends is the use of occupancy models36. Because the precise coordinates of the polygons where the individual populations were monitored are not available, we were not able to test for the potential confounding effect of monitoring extent, but our sensitivity analysis indicated that survey units do not explain variation in the detected trends (Supplementary Fig. 12).
Data limitations: temporal gaps
The population time series we studied cover the period between 1970 and 2014, with both duration of monitoring and the frequency of surveys varying across time series. We omitted populations that had less than five time points of monitoring data, as previous studies of similar population time series data have found that shorter time series are less likely to capture directional trends in abundance63. In a separate analysis, we found significant lags in population change following disturbances (forest loss) and that population monitoring often begins decades to centuries after peak forest loss has occurred at a given site45. The findings of this related study suggest that the temporal span of the population monitoring does not always capture the period of intense environmental change and lags suggest that there might be abundance changes that have not yet manifested themselves. Thus, the detected trends and the baseline across which trends are compared might be influenced by when monitoring takes place and at what temporal frequency. Challenges of analysing time series data are present across not just the Living Planet Database that we analysed, but more broadly across population data in general, including invertebrate datasets65. Nevertheless, the Living Planet Database represents the most comprehensive compilation of vertebrate temporal population records to date, allowing for analyses of the patterns of vertebrate trends and fluctuations around the world.
Data limitations: time series with low variation
Eighty populations ( More
