General Approach

Our general approach was to obtain Black Oystercatcher vital rates from published articles and unpublished reports (Appendix) and use them to build a stage-based matrix population model with a model structure informed by the availability of vital rate data (Morris and Doak 2002). Upon constructing a population model, we conducted Monte Carlo simulations, where model inputs were drawn at random from vital rate distributions (Morris and Doak 2002). These simulations produced distributions of model outcomes with probabilistic uncertainties that propagated from uncertainties in the inputs. Our primary outcome of interest from the population model was the finite annual rate of increase (λ), an estimate of potential population growth based on input rates of survival and reproduction. Following model simulations, we regressed λ values against input vital rate values to determine which of the rates were most closely tied to potential population growth (Morris and Doak 2002). We conducted this exercise two times: once using a detailed version of the model with several subcomponents of reproductive output (hereafter, detailed model), and once using a simpler version, with a more integrative measure of reproduction (hereafter, simplified model).

Literature Survey of Vital Rates

We searched both peer-reviewed and gray literature for Black Oystercatcher survival and reproduction data. Information published before 1995 was identified in a comprehensive species account for Black Oystercatcher (Andres and Falxa 1995). With this foundation, we searched the literature to identify newer information published between 1995 and 2017 (Appendix). We accessed original publications and reports whenever possible to verify variable definitions and values. In most cases, vital rates were pulled directly from manuscript tables. In fewer cases, values were pulled from the text of the manuscript or were extracted from figures.

Here, survival estimates were defined as the proportion of individuals in a particular stage that survived a full year. Note that the reported rates are apparent survival rates, which are biased low by an unknown amount due to imperfect detection and emigration (Sandercock 2006; Murphy et al. 2017). The logical consequence of this bias is that population growth rates derived from apparent survival will be biased low, also by an unknown amount. Reproduction data used in this analysis included clutch size, hatching success, fledging success, and pair productivity. Clutch size was defined as the number of eggs laid per pair per year, and generally reflected the size of the first clutch of the season. Hatching success was defined as the proportion of eggs in a clutch that hatched into nestlings and, again, generally reflected the results of the first nesting attempt of the season. Fledging success was defined as the proportion of nestlings (i.e., hatched eggs) that fledged or were assumed to fledge after approximately 35 days. Pair productivity was defined as the number of fledglings, or assumed fledglings, per breeding pair per year.

Stage-based Population Model

We built a stage-based matrix population model for Black Oystercatcher (Fig. 1), where all calculations were restricted to the abundances and vital rates of females (Caswell 2001). We note that vital rates reported in the literature were seldom sex specific; therefore, we aggregated rates reported for either or, usually, both sexes together. Given the granularity of information on survival and reproduction from the literature, we settled on a model with two stages, a subadult stage (S) with individuals ≥ 1 and < 5 years old, and a breeding adult stage (A) with individuals ≥ 5 years old (Andres and Falxa 1995). The abundance (N) per stage (i) at a future time (t + 1) was calculated from current abundance and a transition matrix as

The transition matrix contained parameters reflecting survival within (P_i) and transition between (G_s) stages where

Here, p_i is stage-specific apparent survival rate and d_i is the stage duration in years (Crouse et al. 1987). For this study, we set the subadult duration d_s = 3 and the adult duration d_A = l - 4, where l was the lifespan for the species. In addition to survival information, the transition matrix also contained information on the reproductive output of adults, F_A, where

Here r is the typical sex ratio of chicks (assumed to be 0.5), c is clutch size, e is hatching success, n is fledging success, ƒ is pair productivity, and s is first-year survival (the apparent survival rate of birds after fledging and until they are 1-year old). Note in Eq. 3 that we consider pair productivity as the product of clutch size, hatching success, and fledging success. In the detailed model, we constructed F_A using the more detailed formulation (middle part of Eq. 3), and in the simplified model, we used the simpler formulation (right side of Eq. 3).

Figure 1.

Life stage diagram corresponding with the stage-based matrix population model for Black Oystercatchers used in this analysis. Circles labeled with S and A represent pools of individuals in subadult and adult stages. Vectors labeled with P represent survival within sub adult and adult stages. The vector labeled with G represents survival and transition from sub adult to adult stages. The vector labeled with F represents reproduction by breeding adults and survival of young through their first year.

Simulated Vital Rates and Population Projections

Vital rates collected from the literature (Appendix) were estimates from across the species range, and incorporate variation due to space, time, and imperfect sampling. We used these estimates to produce a probability distribution for each vital rate. Probability distributions were constructed using Bayesian regression where, in all cases, a collection of estimates was modeled with a global intercept (Wade 2002). In many cases, where multiple estimates were available from each of multiple regions, a random intercept for each U.S. state or Canadian province was also included in the model. We used the posterior predictive distributions from these models, which incorporate uncertainty in both the center and scale of vital rate distributions, as probability distributions to be used as inputs for the matrix population model (Wade 2002).

For each Bayesian vital rate model, we specified an appropriate vague prior distribution that depended on the nature of the vital rate. Clutch size at a Black Oystercatcher nest is a discrete variable that ranges from 1 to 5, though it is typically reported as a population average on a continuous scale (Andres and Falxa 1995). Similarly, the number of fledglings produced per pair per year is a discrete variable that ranges from 0 to 5 in Black Oystercatcher, though pair productivity is typically reported as a population average on a continuous scale (Andres and Falxa 1995). Given that clutch size and pair productivity are usually reported as continuous variables with values greater than zero, we modeled them using gamma distributions. Hatching success, fledging success, first-year survival, and adult survival rates are typically reported as continuous proportions (Andres and Falxa 1995). Given theoretical bounds of 0 and 1, we modeled these rates using beta distributions. Bayesian models of vital rates were constructed using the brms package (Bürkner 2017), a library for interfacing the Stan probabilistic programming language (Carpenter et al. 2017) with the statistical program R (R Development Core Team 2016).

Random draws from posterior predictive distributions for vital rates were used as input values over 50,000 simulations of the otherwise deterministic matrix population model. With this many draws from distributions that, in some cases, extend to infinity, an occasional and highly unusual observation is likely to enter the model, resulting in a highly improbable outcome. These improbable outcomes can obscure the qualitative trends that the simulations are meant to uncover (Morris and Doak 2002). We avoided this situation by truncating input distributions to the middle 99% of values, at the 0.5th and 99.5th percentiles of posterior predictive distributions (Morris and Doak 2002). The bounds of the truncated distributions consistently bracketed observed values (Table 1; Appendix).

There were a few parameters explicitly or implicitly included in the population model for which we had little information. First, we assumed that females begin breeding at age 5 and attempt to breed each year thereafter, an assumption supported by observations reported in Hazlitt and Gaston (2002) and Andres and Falxa (1995). Second, there is little information on the maximum lifespan of the species. Anecdotal observations from California, USA, suggest that individuals live at least 16 years (Andres and Falxa 1995). Band recovery data for the American Oystercatcher (Haematopus palliatus), a closely related species, demonstrates that individuals can live until the age of 23 years and 10 months (Lutmerding and Love 2017). Given the lack of information on lifespan, we modeled it as a uniform distribution with a minimum value of 16 and a maximum value of 24. Other model parameters, such as first-year and adult survival, were based on small sample sizes. In these cases, we note that estimates for Black Oystercatchers were well within the ranges reported for other oystercatcher species (Ens and Underhill 2014; Felton et al. 2017; Murphy et al. 2017; Wilke et al. 2017). Also underlying our model were assumptions that vital rates were not strongly correlated with one another and were not correlated with population density. Both of these assumptions are likely incorrect, adding unknown bias and imprecision to our results (Morris and Doak 2002), but data limitations prevented us from evaluating these relationships in a robust way and incorporating relationships into simulations.

Table 1.

Sensitivity and elasticity of Black Oystercatcher potential population growth rate in relation to population model parameters and their ranges used in simulations.

Population projections of the matrix population model were conducted using the popbio package (Stubben and Milligan 2007) written for the statistical program R. For each Monte Carlo iteration, the matrix model was initialized with an initial abundance of 3,000 subadult females and 3,000 breeding adult females, roughly corresponding with estimates of the global population size (Tessler et al. 2014). The model was projected for 30 years, during which time the model quickly converged on a stable stage distribution and λ value. These values, along with total abundance at 30 years and vital rate inputs, were saved for subsequent sensitivity and elasticity analyses.

Sensitivity and Elasticity of Potential Population Growth to Vital Rates

Sensitivity and elasticity of λ to the different vital rates included in the matrix population model were evaluated using two different methods. First, we assessed these properties visually, using bivariate plots of λ outputs vs. vital rate inputs. Second, we used multiple regression, where an output λ was the dependent variable, modeled as an additive linear function of multiple vital rate inputs. In the context of the detailed model, inputs included clutch size, hatching success, fledging success, first-year survival, adult survival, and maximum lifespan. For the simplified model, inputs included pair productivity, first-year survival, adult survival, and maximum lifespan. For sensitivity analyses, variables were centered on the mean before analysis, retaining their original scale. For elasticity, variables were centered and then scaled by the standard deviation before analysis to attain standardized coefficients. In all cases, an additive linear model explained most of the variation in λ (R² > 0.88), indicating that interaction terms were not necessary (Morris and Doak 2002).

Posterior Predictive Distributions of Vital Rates

We were able to find 75 estimates of nine Black Oystercatcher vital rates, from five U.S. States or Canadian Provinces, in the published or gray literature (Appendix). We identified 20 estimates of mean clutch size from Alaska, California, and Washington, USA, and British Columbia, Canada. Percentiles from the posterior predictive distribution for clutch size were 1.84, 2.26, 2.41, 2.56, and 3.09 for the 0th (minimum value), 25th, 50th, 75th, and 100th (maximium value) percentiles, respectively (Fig. 2). Posterior mean and 95% credible intervals for parameter values of the gamma distribution fit to clutch size data were shape, α, = 128.98 (60.45-223.30) and rate, β, = 53.39 (26.07-88.67).

We found 14 estimates for average hatching success, from Alaska, Oregon, and Washington, USA, and British Columbia, Canada. Percentiles (as above) from the posterior predictive distribution for hatching success were 0.05, 0.33, 0.46, 0.61, and 0.93, respectively (Fig. 2). Posterior mean and 95%. credible intervals for parameter values of the beta distribution fit to hatching success data were shape parameter 1, α, = 3.12 (1.09-6.84) and shape parameter 2, β, = 3.53 (1.84-5.18).

There were 12 estimates of average fledging success, from Alaska and Washington, USA, and British Columbia, Canada. Percentiles (as above) from the posterior predictive distribution were 0.09, 0.38, 0.49, 0.60, and 0.90, respectively (Fig. 2). Posterior mean and 95%. credible intervals for parameter values of the beta distribution fit to fledging success data were α = 5.15 (1.65-11.64) and β = 5.36 (2.46-8.48)

Figure 2.

Distributions for Black Oystercatcher clutch size, hatching success, fledging success, first-year survival, annual adult survival, maximum lifespan, and annual pair productivity. These distributions were derived from posterior predictive distributions of Bayesian regression models trained on vital rate data from the literature. Random draws from these distributions served as input values for Monte Carlo simulations of stage-based population models.

We were only able to identify four individual observations of first-year apparent survival rate, from British Columbia, Canada. Percentiles (as above) describing the posterior predictive distribution of first-year survival were 0.09, 0.50, 0.60, 0.70, and 0.98, respectively (Fig. 2). Posterior mean and 95%. credible intervals for parameter values of the beta distribution fit to first-year survival data were α = 10.10 (1.01-35.61) and β = 6.75 (1.27-12.92).

We found four published adult apparent survival rates, from Alaska, USA, and British Columbia, Canada. Percentiles for the posterior predictive distribution for adult survival were 0.47, 0.88, 0.92, 0.95, and 1.00, respectively (Fig. 2). Mean and 95%. credible intervals for parameter values of the beta distribulion fit to adult survival data were α = 33.33 (3.14-102.57) and β = 3.26 (0.72-4.97).

Finally, we identified 15 estimates of average annual pair productivity, from Alaska, Oregon, and Washington, USA, and British Columbia, Canada. Percentiles (as above) for the posterior predictive distribution for pair productivity were 0.11, 0.39, 0.52, 0.68, and 1.43, respectively (Fig. 2). Mean and 95% credible intervals for parameter values of the gamma distribution fit to pair productivity data were α = 6.62 (2.82-12.15) and β = 12.01 (6.29-17.73).

Simulations of Population Growth

Posterior predictive distributions for vital rates were used to generate random inputs to the matrix population model during Monte Carlo simulations. The Monte Carlo distribution for λ resulting from the detailed model, where F_A λ (c)(e)(n), was centered at 1.00. The 2.5th, 25th, 50th, 75th, and 97.5th percentiles for Monte Carlo distribution were 0.85, 0.95, 1.00, 1.04, and 1.13, respectively (Fig. 3). The Monte Carlo distribution for λ resulting from the simplified model, where F_A ∝ƒ was also centered at 1.00, with percentiles of 0.86, 0.96, 1.00, 1.04, and 1.12, respectively (Fig. 3). In both cases, distributions covered values for both positive and negative potential population growth.

Figure 3.

Histograms depicting Monte Carlo output distributions for Black Oystercatcher potential population growth (Lambda) from the detailed and simplified models. Distributions are centered near 1.00, which indicates a stable population. Spread is determined by spatial and temporal variation in vital rates as well as estimate uncertainty.

Sensitivity and Elasticity of Potential Population Growth to Vital Rates

Vital rates appeared to differ in their impact on potential population growth of Black Oystercatcher. Plots of λ vs. vital rates demonstrate that modeled λ was not strongly related to clutch size or maximum age of adults (Fig. 4). In contrast, first-year survival had an intermediate effect and adult survival rate had a strong effect on modeled λ (Fig. 4). In the detailed model, hatching success and fledging success had strong impacts on modeled λ, and in the simplified model pair productivity, a metric that integrates clutch size, hatching success, and fledging success, had a strong effect on modeled λ (Fig. 4).

Sensitivity and elasticity were also assessed via multiple regression with λ regressed against a linear combination of vital rates. The coefficients for sensitivity and elasticity in Table 1 described a similar pattern as the plots in Figure 4, that modeled λ was most strongly related to hatching and fledging success, pair productivity, and adult survival. The effect of adult survival appeared particularly strong when vital rates were regressed on their original scale for the sensitivity analysis (Table 1). However, when rates were standardized for the elasticity analysis, it appeared that hatching and fledging success, pair productivity, and adult survival were of similar importance (Table 1).

Figure 4.

Pairwise relationships between input Black Oystercatcher vital rates and output potential population growth rate (Lambda) from 50,000 Monte Carlo simulations of the stage-based matrix population model. Gray lines indicate a value of 1.00, which indicates a stable population. Values above and below gray lines indicate increasing and decreasing populations, respectively.