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.1⁷¹. 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 Database⁷² (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)⁴⁵. 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 abundance⁶³. 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 relationships⁷³. 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.org⁴, for bird species from https://birdtree.org⁸, and for reptile species from https://vertlife.org⁶. 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 confidence^74,75.

Rarity metrics, IUCN Red List categories and threat types

We defined rarity following a simplified version of the ‘seven forms of rarity’ model⁷⁶, 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 Facility⁷⁷ (https://www.gbif.org) and range data from BirdLife⁷⁸ (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 Conservation⁷⁹ (https://www.iucnredlist.org). The populations in the Living Planet Database⁷² 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 profiles⁷⁹.

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 (σ²) and observation error (τ²) 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 Database^41,81,82. Previous studies have found that not accounting for process noise and measurement error could lead to over-estimation of population declines⁸³, 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 (σ²) and observation or measurement error (τ²) and population trends (μ):

where X_t and X_t−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, σ²). We included measurement error following:

where Yt is the estimate of the true (unobserved) population abundance with measurement error:

We substituted the estimate of population abundance (Y_t) into Eq. 1:

Given

then:

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 km² for all bird species in the Living Planet Database that had records in the BirdLife Data Zone⁷⁸. For mammal species’ geographic range, we used the PanTHERIA database⁸⁴ (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 km² occurrence area based on species occurrence data from GBIF⁷⁷. 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 ranges⁸⁵.

For the geographic ranges of species monitored in the UK, we calculated range extent using a minimal convex hull approach based on GBIF occurrence data⁷⁷. We filtered the GBIF data to remove invalid records and outliers using the CoordinateCleaner package⁸⁶. 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 package⁸⁷.

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 rredlist⁸⁸. 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 profiles⁷⁹.

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 MCMCglmm⁸⁹. 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):

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 (σ²), 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’⁹⁰ 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:

where taxa_i,j,k is the taxa of the ith time series from the jth species; β₀ and β₁ are the global intercept (in categorical models, β₀ = 1) and the slope estimate for the categorical taxa effect (fixed effect), β_0j is the species-level departure from β₀ (species-level random effect); y_i,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 (σ²) 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 MCMCglmm⁸⁹. 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:

where β₀ is the global intercept (β₀ = 1), β_0l is the phylogeny-level departure from β₀ (phylogeny random effect); y_i,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 (σ²) 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 (σ² 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:

where geographic range_i,k,n is the logged geographic range of the kth species in the ith time series; β₀ and β₁ are the global intercept and slope estimate for the geographic range effect (fixed effect), β_0j is the species-level departure from β₀ (species-level random effect); y_i,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 (σ²) 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 Database⁷². 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 species⁶². 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 Database¹⁷. 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 methods^41,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 models³⁶. 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 abundance⁶³. 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 site⁴⁵. 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 datasets⁶⁵. 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 (<1% of the 9286 time series) had very little variance (see Supplementary Table 7 for full references for those studies). The majority of those studies are for bird species and come from the North American breeding bird survey with a measurement unit of an index⁹¹. We have also observed some time series that appear to show logistic relationships with little natural variance (e.g., time series 468, 10193, 17803, see Supplementary Table 8 for full references). Inspecting the raw data showed that some populations have abundances which follow an almost perfect linear or logarithmic increase over time, as could be the case for modelled, versus raw field data. We provide the references for these studies and cannot definitely attribute the low variance to a particular cause across all studies. Some of these studies are reported in units that are an index which may not capture variation in the same way as other raw units of population data. Some of these time series may represent modelled population data based on demographic information rather than only direct observations of populations (e.g., time series 1355⁹²). We chose to not remove studies that may not be raw observation time series based on visual inspection of trends to avoid introducing bias against populations with naturally low variation into our analysis.

Clustering in the values of population trends and fluctuations

We found a clustering of population trend and fluctuation values in some parts of the population change spectrum. For example, we found two peaks—in small increases and in small decreases over time—which were most prevalent in terrestrial bird studies and species, monitored using an index (Fig. 2, Supplementary Fig. 12). Overall 11.4% of time series had trend values between 0.02 and 0.03 and 11.6% of time series had trend values between −0.03 and −0.02. There was also a similar, but smaller, clustering around trends of 0.25 and −0.25. All reported population trends are from models that converged successfully, and visual inspection indicated to us that the μ values are appropriate estimates for the individual time series (Supplementary Fig. 6e). We investigated the population time series where the value of the population trends over time were estimated to be the same value and found that they came from a variety of taxa, locations and survey methods (Supplementary Fig. 6e). We hypothesise that there might be a publication bias against publishing no net change studies, which could explain the trough in μ values of around zero in long-term studies. The clustering of values for some time series may sometimes be associated with the same time series that also have low variance (Supplementary Fig. 6e, see discussion above). With the information available in the Living Planet Database metadata, we cannot fully explain the clustering in population trends. We advocate for more detailed metadata in future versions of the Living Planet Database to allow researchers to filter the database appropriately for individual analyses.

Challenges in estimating geographic range

Estimating geographic range across taxa, and specifically for species that are not birds or mammals, remains challenging owing to data limitations. We used a static measure of geographic range, which does not account for changes in species distributions over time. Furthermore, species could naturally have a small range or the small range size could be due to historic habitat loss⁹³. The UK populations included in the Living Planet Database are predominantly from species with wide geographic ranges (Supplementary Table 5), and our global-scale analysis of the relationship between-population change and geographic range is based on mammal and bird data. As data availability improves, future research will allow us to test the effect of geographic range on the trends of other taxa, such as amphibians and sharks.

Trends relative to null expectation

We tested whether the number of increasing and decreasing populations trends differed from a null expectation using a data randomisation approach (Supplementary Fig. 5b). We used linear models to estimate trends in the data and randomised data with identical structure to the Living Planet Database. We found that there were over 10 times more population declines and increases in the real data relative to the randomised data (2.29% of trends were declining and 2.30% were increasing in the randomised data, versus 28.9% and 32.5% of time series, which had significant negative and positive slopes in the real data, respectively).

Monitoring duration, sampling methods and site-selection bias

To assess the influence of monitoring duration on population trends, we used a Bayesian linear model. We modelled population trend (μ) as a function of monitoring duration (years) for each population, fitted with a zero intercept, as when duration is zero, no population change has occurred. Monitoring duration was weakly positively related to vertebrate population trends, with slightly greater population increases found for longer duration studies (Supplementary Fig. 6, Supplementary Table 2). There was a similar weakly positive effect of number of time points within time series (Supplementary Table 2). In addition, we tested if monitoring duration influenced the relationships between population trends across systems, and population trends across taxa. We found that duration did not influence those relationships, with the exception of reptiles, where declines were more frequent as monitoring duration increased (Supplementary Table 2). Variation in population trends was not explained by sampling method across the five most commonly used abundance metrics (population index, number of individuals, number of pairs, number of nests and population estimate, Supplementary Fig. 12). Following Fournier et al.⁶⁴, we tested the time series that we analysed for site-selection bias. Removing the first five survey points reduces the bias stemming from starting population surveys at points when individual density is high, whereas removing the last five years reduces the bias of starting surveys when species are very rare. The distribution of population trend values across time series was not sensitive to the omission of the first five (left-truncation) or the last five years (right-truncation) of population records (Supplementary Fig. 5a). In addition, we used a data randomisation approach to compare the distribution of trends from the real data to a null distribution and found different patterns (Supplementary Fig. 5b). Overall, our sensitivity analyses suggest that our findings are robust to the potential confounding effects of differences in monitoring duration, sampling method and site-selection.

Comparison of modelling approaches

We conducted the following supplementary analyses: in the second-stage Bayesian models estimating population trends across systems, biomes, taxa and rarity metrics, (1) we weighed μ values by the square of τ², the observation error estimate derived from the state-space models⁴⁰, (2) we used slopes of linear model fits of abundance (scaled at the population level, centred on zero and with a standard deviation of one)⁷³ instead of the μ estimates from state-space models, (3) we modelled the standard error around the slope values of the linear models, the error around μ (half of the 95% confidence interval) and the standard deviation of the raw population data for each time series as additional metrics of population variability. To allow comparison, we scaled the different metrics of population variability to be centred on zero and with a standard deviation of one before they were used as response variables in models. All different analytical approaches yielded very similar results (see main text and Supplementary Figs. 5, 6 and 16, Supplementary Table 2).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com