We obtained unbiased estimation of trends in amphibian occurrence probabilities (i.e., changes in the number of populations that occur when local extinctions exceed (re-)colonizations within a metapopulation) from 61 study areas across North America (Fig. 1; Table S1), incorporating data from 389 time-series of 83 species of amphibians (Table S2) using a hierarchical occupancy model30. Study areas were mainly Federal and State protected areas, and we had no a priori expectation that amphibian populations or resource conditions were in decline or deteriorating when surveys were initiated. Data came from field surveys of multiple habitat units (e.g., wetlands, stream segments, forest patches) within each protected area and included habitats where amphibian species were expected to have a chance of occurring. Site selection was unconditional on initial amphibian presence, and sites at each study area were selected under a probabilistic design (e.g., simple random sampling, stratified random sampling); this means that though site selection differed among study areas, statistical inference is appropriate at the scale of each study area. Survey methods differed among habitats and species, but were generally chosen to maximize the probability of detecting the target species and sometimes multiple standard methods for detecting amphibians were used. Multiple species may have been detected during a single survey occasion. Surveys were repeated on multiple occasions within a time period when each species had the possibility of occupying a site (i.e., sites were assumed to be demographically ‘closed’ during the duration of surveys). These detection-nondetection data were analyzed in our occupancy model, and used to estimate detection probabilities, thereby accounting for unequal effort and other sources of observation error so that inference is on the (true) state of a population (i.e., present or absent) and is unconditional on whether or not a species was detected at a site during sampling occasions31.

Statistical model

Using Markov Chain Monte Carlo methods30,32, we estimated trends of amphibian populations for each metapopulation (i.e., each combination of species and study area), and the effects of four of the most widely cited threats [land-use, contaminants, climate change and disease risk5,6,8,13], as there is evidence that these threats have contributed to the extinction of local populations [e.g., (26–29)]. Trends were estimated for each metapopulation, and represent the balance of local extinctions and colonizations over each time-series; a negative trend occurs when extinctions exceed colonizations. We defined 5 regions as geographic clusters of protected areas.

For all studies, multiple observations were collected within a year, which allowed us to estimate unbiased trends across years. We denote year by y and time-series by k, so that ψ yk is the probability a site in the kth time-series in the yth year is occupied by a given species and p yk is the probability of observing a given species during a visit to a site, given that site is occupied by the species. We denote the number of visits to the ith site in a given year for a given time series as N iyk and the number of times the species was observed at the site during those visits as Y iyk . The true, but unobserved, state of a site is denoted by z iyk , where z = 1 if a site is occupied and z = 0 if not.

Observations Y iyk are modeled with a binomial distribution as:

Note the binomial probability is conditional on the true occupancy state of the site, which is a Bernoulli trial with probability of occurrence for the specific combination of year and time-series. To accommodate variation of detection probabilities among different years and time-series, we assumed a random distribution for the detection parameters:

We estimated occupancy trends as a derived parameter within the MCMC algorithm, calculating the average rate of change in number of occupied sites (λ) in each time-series as a function of the estimated trend parameters in the logit-linear model as

where ψ i and ψ f are occupancy at the start and end of each time series.

Our interest was in understanding the factors associated with annual variation in occupancy probabilities. We expected that relationships would not be evident among threats and continental-scale population trends if declines are driven by different factors among populations, but that relationships may be detectable at regional scales if these threats influence populations more consistently at a smaller spatial scale.

We fit a generalized linear mixed model, where the model included two random components, α k , a random-effect for time-series intercept, and δ yk , a random-error component to account for sub-sampling of sites within each combination of year and time-series, where:

Distributions for random-effects were specified as:

A fixed model, βX, was used to specify the effect of year (continuous), threat, and interaction between threat and year. Covariates for threats were scaled and centered to have a mean of 0 and standard deviation of 1 and year was centered to have a mean of 0 at the level of the time-series. To account for model uncertainty in the inclusion of parameters and minimize potential effects of correlation among threat covariates, we fit parameters using a Bayesian LASSO (least absolute shrinkage and selection operator) approach33,34, by specifying the β parameters to come from a double-exponential distribution, where:

and τ, the regularization parameter was estimated from the data used to fit the model (example code to implement a LASSO provided in Supporting Online Material). To allow for a single value of τ to be estimated, we specified year covariates and interactions between year and threats to have a mean of 0 and standard deviation of 1 when fitting the model. We then back-transformed estimates so the estimated effect of year represented the average annual change in occupancy across all time-series and the estimated interaction effects represented the average annual effect of a 1-standard deviation increase in the value of a threat. The Markov Chain Monte Carlo-estimated posterior probabilities of the parameters (i.e., the standardized partial regression coefficients) can be interpreted as the relative contribution of each threat to the observed population trends.

Priors were chosen as follows: logit−1(μ α ) and logit−1(μ p ) were uniform (0,1) on the real scale; σ p , σ α , and σ δ were uniform (0, 3); and τ was uniform (0, 20). Models were fit using JAGS run via R (v.2.15.2; R Core Development Team 2012). We used a burn-in sample of 5000 iterations and then estimated the posterior distribution of parameters based on the next 25000 samples. Three chains were run and model convergence was assessed based on Gelman-Rubin statistics35.

Threats

Disease risk was defined as the suitability for occurrence of the fungal pathogen B. dendrobatidis (Bd)36, as it is classified as one of the “enigmatic” causes of amphibian decline13,37. We used spatially referenced predictions of Bd suitability for the United States, which were generated under the environmental niche model described in36 (provided by D. Olsen, US Forest Service). In a binomial regression model, the detection of Bd was found to be different among major ecosystems of the US, generally increasing with number of species, and increasing with decreasing minimum temperature (details of model found in36). Predicted suitability under this model is qualitatively similar to the global model38 in8. There are no comparable data on other emerging amphibian pathogens such as Ranavirus or Batrachochytium salamandrivorans.

The threat from climate change relates to an increased risk of extinction in response to changes in precipitation, mainly resulting in increased frequency and intensity of droughts39. This threat from climate change was calculated as the mean difference in mean annual precipitation from (2001–2011) and the 30-year normal (1981–2010; available from: http://www.prism.oregonstate.edu/, accessed 13 March 2014), which was a proxy for water availability. To obtain the difference in annual average precipitation during the period of our surveys (2001–11), we subtracted the 10-year average annual precipitation for each sub-basin over the period of record (monthly data obtained from PRISM for 2001–11; http://www.prism.oregonstate.edu/recent/; accessed 13 March 2014) from the 30-year (1981–2010) normal precipitation (http://www.prism.oregonstate.edu/normals/; accessed 13 March 2014).

The threat from land use change was evaluated because it is the single most important threat to biodiversity in general40, and both human influence (Human Influence Index; HII41) and the application of pesticides42,43 are components of this threat. The human influence index (HII), developed by the Wildlife Conservation Society and the Center for International Earth Science Information Network41, is an index of direct human disturbance which incorporates human population density, land transformation, accessibility, and electrical power infrastructure to represent direct human influence on the landscape. We summarized the 1-km2-pixel HII raster map provided by41 as the mean HII within each sub-basin.

The threat from contaminants was described by the estimated annual pesticide use values for compounds known or hypothesized to be important to amphibians (Table in44 and updated by K. Smalling, USGS pers. comm.). We used the ‘EPest-high’ application estimation method from42,43 for each county from 2001–11, which treats non-reporting of application as missing values which are estimated from nearby data. Values were summarized at the sub-basin scale by weighting the average application by the area of the county falling within each sub-basin, and thus reflect the average pesticide application per unit area.

For each threat, we calculated the average, normalized intensity (the value of each threat) within each sub-basin (i.e., United States Geological Survey (USGS) HUC4-scale sub-basins; http://pubs.usgs.gov/wsp/wsp2294/html/pdf.html) for use in our occupancy models. We investigated the correlations among the normalized threats for all 204 sub-basins in the continental United States using the Band Collection Statistics tool in Spatial Analyst in ArcGIS 10.2. Not all sub-basins had amphibian population data; correlations among threats for the subset of sub-basins (n = 55) with amphibian species data were determined independently. All correlations were <0.47, except a positive correlation between Bd and climate for the 55 sub-basins with amphibian data (ρ = 0.648) results from the inclusion of precipitation in the suitability model of36.