Field data collection.

The data presented in this paper come from camera trap images collected by a global tropical forest monitoring network—TEAM. TEAM is a consortium of three core partners—Conservation International, The Wildlife Conservation Society, and The Smithsonian Institution—and over 80 academic and local partners spanning 15 countries and 17 sites in tropical forests around the world. Data were collected between 2007 and 2014 from sites throughout Africa, Asia, and Latin America. We analyzed data from 15 sites, each of which was surveyed for three or more years (Fig 1 and S1 Table). TEAM monitors mammal and bird communities using camera trap arrays following a standardized protocol [44–46]. At each site, 60–90 camera traps are deployed at a density of 1 camera per 1 or 2 km², depending on the forested extent of the site. Each camera trap point is deployed for at least 30 consecutive days during the dry season, which is defined as months with <100 mm average rainfall, or the drier part of the year at sites with no dry season. More details on the field methods can be found in [47].

TEAM monitoring sites.

Sites were chosen to span moist tropical forest rainfall and seasonality gradients [48], land use change gradients [49], and projected climate change gradients [50]. The sites are distributed in proportion to tropical forest cover on each continent; half of the sites are in South and Central America, roughly one quarter are in Africa, and the rest are in Southeast Asia.

Criteria for species inclusion in TEAM monitoring.

Mammal and bird species captured on camera trap images were only monitored and reported in this study if (1) body mass was greater than or equal to 100 g [51,52], (2) data [53–55] indicated the species spends a large proportion of its time on or near the ground, or (3) data suggested a species was arboreal and there was at least one TEAM site at which the species was detected in five or more photographic events (collection of images separated by at least 1 min) for every year that data were collected. A complete list of all populations and species with relevant metadata is in S2 Table.

Data preparation.

Observational records were condensed into presence–absence matrices, one for each site, species, and year, where the rows corresponded to sampling points and the columns correspond to time periods (days). For each sampling point and day, the cells in these matrices were occupied by a “1” if the species was photographed, a “0” if the species was not photographed, or “NA” if the point was not actively sampled on this day. To reduce model computation time and increase efficiency, we grouped observations into 15 time periods for each site, species, and year so that each time period equaled approximately 7–8 d of sampling. After grouping, observations remained as 0s and 1s. Grouping does not affect model estimates, but only the units of the estimated detection probabilities for each species [15]. All data processing was conducted with scripts programmed in the language R [56].

High detection (Case 1).

Common populations (N = 105 populations). For populations that had an observed average annual detection (number of times species i was detected across all sampling points j and time periods k averaged over t years (∑y_ijkt/t ∑j_tk_t) that was greater than 8%, we used a hierarchical dynamic occupancy modeling approach with covariates [57]. Here, the ecological process that influences occupancy is modeled separately from the observation process (i.e., detection), which is considered imperfect. In the ecological process model, the incidence of a species i in the first year of observation is an unobserved latent variable Z_i1, resulting from a Bernoulli process with expected probability ψ_i1: Z_ij1 ~ Bernoulli(ψ_i1). For subsequent years (year > 1), the site remains occupied by the species with probability ɸ_i, or goes locally extinct from that site with probability (1−ɸ_i), where ɸ_i is the apparent survival of the species from one year to the next and (1−ɸ_i) is the local extinction probability. If the species did not occupy the site in year t (Z_it = 0), it can colonize this site by the following year (t +1) with probability γ_i [57]. The dynamics of a species i can be described recursively for any two contiguous years as: (1)

The observation process of the model assumes the observations for each species i at each sampling point j, year t, and observation period k, y_ijtk, as realizations of a Bernoulli process with mean Z_it p_i, where p_i is the detection probability of species i at year t: y_ijtk ~ Bernoulli(Z_it p_i). This takes into account imperfect detection at the sampling point (i.e., false negatives) and provides an unbiased estimate of occupancy. Under this model, closure is assumed within years but not between years; thus, colonization and extinction only occur on an annual basis. Covariates were selected to test the effects of climate (i.e., rain, temperature) and anthropogenic factors (i.e., people, edges) on occupancy, survival, and colonization. The hypotheses, stated as equations, describe the relationship between state variables and covariates:

For first year occupancy ψ_ij1: (2)

For survival probability ɸ_ijt: (3) and colonization γ_ijt: (4)

We used three types of covariates for each site: sampling point covariates (elev, edge), time-dependent covariates (rain, maxTemp, minTemp), and a covariate that was sampling-point- and time-dependent (peop).

Sampling point covariates were camera trap, elevation above sea level (elev), and camera trap distance to the nearest edge (edge) that could be a road, river, forest, protected area boundary, or settlement. Each camera trap location was determined in the field using a handheld GPS device. Elevation (meters above sea level) of each camera trap was obtained from a 90 m spatial resolution void-filled Shuttle Radar Topography Mission (SRTM) digital elevation model from the Consortium for Spatial Information [58]. Distance to edge was calculated using a geoprocessing tool in ArcGIS that determines the distance from each feature in the “input features” to the nearest feature in the “near features.” The camera trap points were the “input features,” and a polygon or polyline shapefile of the edge type (e.g., roads, rivers) were the “near features.” Local data source shapefiles were used for each edge type. When a local spatial data layer was not available for a particular edge type, we used the Global Roads Open Access Data Set, Version 1 [59], forest classes from Globcover V2.3 [60], and the World Database on Protected Areas protected area boundaries [61].

Time-dependent covariates changed only between sampling years but not between sampling points. These included cumulative rainfall (mm) during the 6 mo prior to the start of the camera trap deployment (rain), average maximum daily temperature (°C) at the site during the 6 months prior to the start of the camera trap deployment (maxTemp), and average minimum daily temperature (°C) at the site during the 6 mo prior to the start of the camera trap deployment (minTemp). Climate data were extracted from local TEAM climate stations installed at each site that follow a standard protocol for data collection [45]. Rainfall data were aggregated from five-minute resolution measurements and summed over six months. Maximum and minimum temperatures were calculated daily and averaged over six months.

The sampling-point- and time-dependent covariate was the presence or absence of people, as detected by the camera trap at that point in a particular year (peop). Presence or absence of people for each site and each year was derived directly from camera trap data at that site; local site managers identified people as another species captured in the images. The people covariate includes all people regardless of their potential benefit or threat to wildlife (e.g., tourists, park visitors, poachers, field personnel, park patrol). All covariate values (except presence or absence of people) were normalized by subtracting the site mean and dividing by the site standard deviation to allow for easier interpretation of effects and to help with model convergence.

We fit the models using Markov chain Monte Carlo (MCMC) methods and the software JAGS [62,63], running through package R2jags in R [64]. Uniformative priors were used for all state variables, latent variables, and covariate parameters (S1 Text). To identify potentially important covariates in the model, we used an indicator variable selection approach where an indicator variable (w_i) multiplies each model term i (w_i*β_i*cov_i) and is given a prior of 0.5 from a Bernoulli distribution (w_i ~ Bernoulli (0.5)) [65,66]. As the MCMC chain runs, the value of w_i reflects the weight of that covariate in the model, or the probability that the variable should be included in the model. If the covariate is very important, E [w_i] is close to 1; if the covariate does not have any influence in model performance, E [w_i] is close to 0. By running long MCMC chains relative to the number of potential model subsets (for this case 2¹¹ = 2,048), it is ensured that most potential models are visited. We used three MCMC chains with 250,000 iterations each, and examined the last 1,000 iterations. At the end of the process, the posterior distribution of the effects (β_i), where w_i > 0.5, were extracted, and 95% highest posterior density (HPD) were calculated. However, the posterior distribution of occupancy (or any of the other state variables) was estimated from the last 1,000 iterations of the chain sampling the entire model ensemble, representing a model average set over all possible models [16]. The approach is similar to stochastic variable selection [63] with the difference being that the indicator variable does not come from a normal mixture process.

Medium detection (Case 2).

Populations with five or more detections per year but less than 8% annual detection (N = 128 populations). For populations that were detected at a site on five or more occasions per year, but that had an annual average observed proportional detection less than 8%, we used the dynamic occupancy modeling approach described above, omitting all covariates. We fitted models using MCMC methods following the approach in [57], using software JAGS [62,63] running through package R2jags in R [64]. JAGS models were based on existing code [57] and were run with three chains of 30,000 iterations each, a thinning rate of three (every third iteration discarded), and a burn-in rate of 29,000 iterations (ignored for calculation of posterior density distributions). Models were checked for convergence by ensuring that the Gelman-Rubin statistic for each parameter was below 1.03 [67]. Model parameters were recovered using the mode of the distribution because many of the posterior distributions were highly skewed. HPD intervals were extracted from the posterior distributions at a 95% level for inference purposes.

Low detection (Case 3).

Populations with fewer than five detections per year (N = 278 populations). For populations i that were detected an average of fewer than five camera trap points j per year t at a given site, we calculated the observed “naïve” occupancy (number of points j at which a population was detected at least once [Z_ijt = 1/number of points sampled]). The error associated with the naïve occupancy was estimated through a Bayesian approach using a beta prior for the occupancy on a binomial distribution [68]. A random sample of 1,000 realizations was drawn from the posterior, and 95% HPD intervals were calculated. Calculations were performed in R using package binom [69].

WPI calculation.

We use the WPI [20] to assess changes in terrestrial vertebrate biodiversity. The WPI is an indicator specifically developed for camera-trap data, designed for use with the TEAM survey, and created to meet the requirements of an indicator for the CBD. Specifically, the WPI tracks changes over time [21] and is sensitive to variation in species richness, species evenness, and changes in species occupancy or abundance. The WPI_s,t is the geometric mean of the occupancies of all the species (i…n) in the community s scaled by the baseline occupancy for each in year t [15]: (5)

This index is unique, because the underlying occupancy models account for variability in detection [72,73] within and among species using unbiased population trends.

The WPI is an ideal metric for evaluating changes in biodiversity because it is sensitive to changes in richness, relative abundance or occupancy, dominance, and other measures of community diversity [21]. We used the approach proposed by Ahumada et al. [74], which allowed us to compute the WPI as a derived quantity from the posterior occupancy distributions (geometric mean of the relative occupancies at each model iteration after model burn-in), resulting in a full posterior distribution for the WPI each year. From these distributions, we extracted the 95% HPD intervals and used the median or mode (whichever was lowest) as the measure of central tendency of the distribution.

The WPI can be aggregated based on spatial extent or species attributes. For example, it can be calculated at the level of a TEAM site (local), of sites on the same continent (regional), and of all sites (global), or for species grouped according to IUCN conservation status or dietary guilds. When aggregating the WPI for sites that started in different years, the WPI distributions are convolved yearly in one distribution, ignoring the baseline of the site that started later in time, because the WPI in the baseline year for a site is always 1 with no error distribution.

WPI Analytics System.

The WPI Analytics System uses the HPE Vertica Analytic Platform. The system integrates existing data capture systems for TEAM camera trap and covariate data. Any new data or changes to the data are automatically migrated from the TEAM Network’s database to an HPE Vertica system using a versioning scheme to maintain unique versions for each migrated data set, thereby maintaining data integrity. After each migration, the WPI system autotriggers a series of computations for new occupancy values of each species for migrated site data. The compute logic is implemented as HPE Vertica User Defined Transform Functions (UDTF), invoked by a set of dedicated compute thread pools in the WPI application. The Compute Engine leverages HPE Vertica-R interoperability and JAGS library for running MCMC simulations [64]. The WPI system provides complete traceability of input data and simulation results for each run. The system provides alerts and feedback on the progress of compute-intensive simulation jobs during compute/recompute and publish/unpublish tasks in the WPI workflow.

The system model implemented on HPE Vertica runs a large number of MCMC simulations on multiple species—more than 33 million iterations for the existing TEAM data from 15 sites. The last 1,000 iterations for each species at a site per year, yielding approximately 2 million data rows, are used for multilevel WPI computation and analysis. The WPI system is capable of processing much larger volumes of data, running hundreds of millions of iterations and enabling TEAM and others to add new sites and new data sources.

Occupancy status categorization.

The dynamic occupancy models were mechanistic and captured year-to-year variation in occupancy due to changes in survival and colonization probabilities. To estimate long-term occupancy trends in these populations, we removed temporal variation in the annual occupancy estimates by fitting logistic regression models. Specifically, we fit logistic regression models to the detection-corrected yearly occupancy estimate distributions (Cases 1 and 2) and the raw occupancy distributions (Case 3) for each population i as a function of time t: (6)

A separate logistic regression through time was fitted for each modeled realization (1,000 realizations for each species population at a site), resulting in a distribution of slopes and intercepts for each population at each site. The median and two-tailed 80th percentile intervals were calculated to determine whether the slopes were significant at the 80th percentile (Fig 2B and S3 Fig). Logistic regressions were fitted in R using the glm function in package stats [56].

Populations were considered to be significantly declining, increasing, or stable based on whether the 80th percentile interval of the slope distribution was entirely below zero, entirely above zero, or neither, respectively (but see below for “unknown” category). In the context of conservation, false positives (i.e., species categorized as declining when actually stable) are preferable to false negatives (i.e., species categorized as stable when actually declining). We therefore used a conservative cutoff (80% credible intervals) for classifying occupancy trends, because it is better to provide an earlier warning signal of occupancy declines that can prompt conservation action [75].

We used log-likelihood goodness of fit tests (G-tests) to assess whether the proportion of populations with increasing, decreasing, or stable occupancy (i.e., demographic status) was independent of case. A significant result indicates significantly different frequency distributions of demographic status among cases [76]. There were no significant differences in the proportion of populations with decreasing, stable, or increasing occupancy between Cases 1 and 2 (i.e., high and medium detection; G-test, G = 0.61, df = 2, p = 0.738), but there were significant differences in occupancy for populations modeled as Case 3 (i.e., low detection) (G-test, G = 35.44, df = 4, p < 0.001) with a higher proportion of populations classified as stable. This is likely due to reduced ability to detect changes because of very few observations (i.e., fewer than five detections annually) rather than to actual occupancy stability. For this reason, and the potentially negative conservation repercussions from failing to detect change, we classified the occupancy status of all Case 3 “stable” populations as “unknown.”

Analysis of occupancy status.

We used log-likelihood goodness of fit tests [76,77]) to examine variation in the proportion of populations with increasing, decreasing, unknown, or stable occupancy based on six categorical variables (i.e., class, IUCN status, body size, guild, landscape type, hunting; see below). We accounted for multiple comparisons and controlled for family-wise error rate by applying a Bonferonni correction (α = 0.05/6 = 0.008).

We created S1 and S4 Figs with relative multiple bar (rmb) plots using the R package extracat [78]. Relative multiple bar plots display the relative frequencies of a response variable (i.e., occupancy status) for each combination of explanatory variables (i.e., class, IUCN status, body mass category, guild, landscape type, hunting status). The total sample is divided over a grid-like graphical display, which illustrates the corresponding sample weights for each variable combination.

We used generalized linear models with a binomial distribution to model the proportion of populations per site as a function of continuous predictors (S1 Table). We examined the proportion of populations with (A) increasing or (B) decreasing occupancy within each site using each of the following continuous variables in a separate model: (1) log of protected area size, (2) human population density, (3) proportion forested, (4) edge density, (5) proportion of hunted populations, (6) proportion of not hunted populations, and (7) years of camera trap data. For example: glm(Site Increasing Proportion ~ Site Human Population Density, family = binomial). Pairwise relationships between these continuous variables are shown in S5 Fig.

Potential explanatory variables.

(1) Body size and guild. We used species-level body mass (g) and guild data compiled from published sources [51–55]. Body mass ranged from 107 g to 39,400 kg (S2 Table), with a median of 3,250 g. We evaluated three body mass categories, defined by orders of magnitude as 100–1,000 g (39 bird and 25 mammal species), 1,001–10,000 g (21 bird and 94 mammal species) and >10,000 g (65 mammal species). Each species was categorized as one of four dietary guilds: carnivore (49 mammal species), herbivore (eight bird and 59 mammal species), insectivore (19 mammal species), or omnivore (52 bird and 57 mammal species) (S2 Table).

(2) Structural landscape connectivity in the zone of interaction (ZOI). Processes outside of the protected area boundaries may affect mammal and bird populations within the camera trap sampling area. We therefore examined the broader landscapes outside of the protected area boundaries of each TEAM site using the zone of interaction (ZOI) to define spatial extent. The ZOI is the area that has the potential to strongly influence biodiversity at the site based on systematic quantification of surrounding watersheds, migration corridors, and human settlements [32].

To map forest cover at each TEAM site, we use the Global Forest Change (GFC) product [79]. After conducting a sensitivity analysis across a range of thresholds, we applied a 75% threshold to the 2000 forest cover layer in the GFC product, taking into account the gain and loss layers in the product to calculate a forest–nonforest map for each site for 2012. We verified the accuracy of these maps through comparison with independently derived land cover classifications. We then filtered the data to set a minimum patch size of 990 m² (S6A Fig). We calculated two measures of landscape connectivity (edge density and the proportion of forested landscape [80]) from the forest cover data for each TEAM ZOI using the “ClassStat” function in the SDMTools library in R [81]. We square root transformed edge density to normalize its distribution and calculated a pairwise distance matrix using a Mahalanobis distance to remove any correlation between the connectivity metrics [82]. We then ran an unweighted pair group method with arithmetic mean (UPGMA) cluster analysis on the distance matrix and identified the landscape types based on the three clusters produced at a height cutoff of 1.5 in the resulting dendrogram [83] (S6B Fig). Based on visual inspection of forest cover (S6A Fig), we classified the three clusters as: (1) intact landscapes in which the protected areas were indistinguishable from surrounding forest (i.e., CAX, COU, NNN, CSN, and YAS), (2) isolated landscapes in which the protected areas were surrounded by non-forested landscape (i.e., BIF, RNF, and UDZ), and (3) patchy landscapes in which the protected area was embedded in a patchwork of forested and nonforested areas (i.e., KRP, PSH, VB, BBS, YAN, BCI, and NAK).

The dichotomy between forest and nonforest cover in our analyses does not differentiate between natural and anthropogenic land uses. For example, some forested areas may be agricultural landscapes with closed canopies (e.g., palm oil surrounding PSH in Malaysia; banana plantations near VB in Costa Rica) that do not provide benefits to wildlife comparable to native forest [84]. Additionally, some nonforested area may represent natural vegetation that arises from breaks in forest cover (e.g., high elevation areas above tree line at YAN in Peru).

(3) Protected area size. To estimate protected area size, we extracted the polygon of the protected area for each TEAM site from the World Database on Protected Areas (WDPA) [85] and verified each polygon with the appropriate local site manager. We calculated the area of each protected area in hectares after reprojecting the polygons to the appropriate local (UTM) coordinate system.

(4) Hunting. Field managers at each TEAM site were surveyed about hunting within the protected area where the core TEAM camera trap sampling occurs. We distributed species lists to managers and asked them to mark whether each species at their site list was hunted (yes), not hunted (no), or whether they did not know (unknown). Our measure of hunting was simply the proportion of populations reported as hunted, which relies on expert opinion and does not take the level of hunting pressure into account. Primary data on hunting warrant further investigation.