Count data suggest that Black Oystercatcher (*Haematopus bachmani*) has locally variable but globally stable populations. A simple, stage-based matrix population model for Black Oystercatcher was built and Monte Carlo simulations of the model were conducted using vital rates from peer-reviewed and gray literature. Simulations yielded a distribution of potential population growth rates that extended from 0.87 to 1.14 and was centered at 1.00, supporting the hypothesis of a globally stable Black Oystercatcher population. Sensitivity and elasticity analyses of the population model showed that potential population growth is particularly sensitive to changes in hatching success (i.e., proportion of eggs hatched), fledging success (i.e., proportion of chicks fledged), and breeding adult annual survival. These rates could be possible targets for population management should it become necessary given future changes in sea temperature, sea level, and coastal development. Pair productivity, which integrates hatching and fledging success, is suggested as a simple and valuable metric for monitoring population growth potential of Black Oystercatcher.

Black Oystercatcher (*Haematopus bachmani)* is a species of management interest because of its limited geographic range, small global population, and reliance on intertidal habitat, which is experiencing ongoing threats related to human development and sea level rise (Tessier *et al.* 2014; U.S. Shorebird Conservation Plan Partnership 2016). The global population of Black Oystercatcher is thought to fall between 12,000 and 17,000 individuals (Tessier *et al.* 2014; Weinstein *et al.* 2014). Ongoing seabird counts (Tessier *et al.* 2014) and Audubon Christmas Bird Counts (Weinstein *et* *al.* 2014) suggest that Black Oystercatcher abundance is increasing in some places and decreasing in others, but is globally stable. Count data, however, can be a misleading indicator of population viability (Van Horne 1983; Pulliam 1988) because counts are affected by movement of individuals, independent of reproductive dynamics. Vital rate data can complement count data when assessing population status (Fryxell *et al.* 2014) and are especially useful when interpreted in the context of population models that allow understanding the net results of integrated rates on potential population growth (Hitchcock and Gratto-Trevor 1997).

The goal of this study was to survey the Black Oystercatcher literature and use the available data on Black Oystercatcher vital rates to build and populate a stage-based matrix population model. The model was intended to be simple and generic, built for heuristic purposes (Hitchcock and Gratto Trevor 1997). Our objectives were to use the model to explore: 1) if vital rate observations reported in the literature are consistent with the hypothesis that Black Oystercatcher populations are globally stable; and 2) how different life stages influence potential for population growth should management actions to increase Black Oystercatcher populations become necessary.

## Methods

### 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*) stages where

_{s}Here, *p _{i}* is stage-specific apparent survival rate and

*d*is the stage duration in years (Crouse

_{i}*et al.*1987). For this study, we set the subadult duration

*d*= 3 and the adult duration

_{s}*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*where

_{A},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).

### 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*^{2} > 0.88), indicating that interaction terms were not necessary (Morris and Doak 2002).

## Results

### 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)

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.

### 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).

## Discussion

We built a stage-based matrix population model for Black Oystercatcher and populated the model with vital rates from the Black Oystercatcher literature to see if observed vital rates yielded population growth rates consistent with the hypothesis that Black Oystercatcher populations are globally stable. We found that, indeed, average empirical rates from across the species range produce potential population growth rates that are not significantly different from λ = 1.00, which would indicate stable populations. Thus, the hypothesis of globally stable Black Oystercatcher populations based on abundance observations is defensible from a vital rates perspective.

Several factors governing vital rates and abundance of Black Oystercatcher will likely change over time. For example, ongoing sea level and sea temperature rise (National Research Council 2012) will likely change the availability of intertidal habitat and prey (Galbraith *et al.* 2002). Coastal development, and subsequent changes in predator communities, will also likely have impacts. If, at some point, it becomes necessary to actively manage Black Oystercatcher populations, it will be useful to know which life stages are likely to have the greatest impact on population growth rates. Sensitivity and elasticity analyses of the population model indicated that activities that increase the likelihood of hatching and fledging success will have a relatively large positive impact on population growth rates. Similar activities have been adopted in the management of other shorebird species (Isaksson *et al.* 2007), such as the Piping Plover (*Charadrius melodus*), where it is common to erect physical barriers to nest predators (Mabee and Estelle 2000; Maslo and Lockwood 2009). The other vital rate with a notable effect on population growth rate was adult survival, a common finding for long-lived vertebrates (Saether and Bakke 2000). Managing adult survival is likely to be more difficult than nesting success given that adults are highly mobile while eggs and chicks are not. Further, survival estimates for Black Oystercatcher are relatively high already, so there is little room improvement (Felton *et al.* 2017).

Shorebird monitoring programs occur across the Americas, and many of these programs yield estimates of hatching success, fledging success, or pair productivity. The results from this analysis provide some insight into how vital rate observations for Black Oystercatchers can be interpreted. For example, Figure 4 suggests that pair productivity observations greater than about 0.65 are likely a good sign for a local population under the assumption that first-year and adult survival are not too far below average for that population. Similarly, pair productivity estimates below 0.35 are likely a sign of caution. These suggestions come from a qualitative reading of Figure 4, which is most appropriate given the nature of our analysis. However, it is interesting to note that these patterns are very similar to those recently described for a closely related species. For instance, Felton *et al.* (2017) suggested that pair productivity of 0.63 indicated stable populations of American Oystercatchers along the Atlantic coast of North America.

We recommend qualitative interpretation of our results given the simplifications made in constructing this model and running the simulations. Given the limited availability of data, we were forced to assume that vital rates were independent of one another. This is likely not true, as conditions that favor hatching success might also favor fledging success (Morris and Doak 2002). It is encouraging that model results were so similar between the detailed model, populated by independent hatching and fledging success inputs, and the simplified model, populated by, in some cases, independent pair productivity estimates that naturally incorporate potential correlations between hatching and fledging success. Model simulations also implicitly assumed that Black Oystercatcher vital rates are not strongly density dependent. Again, this assumption is not likely to be true (Morris and Doak 2002), but estimating density dependence for vital rates requires data on multiple vital rates from single time periods and locations, which were not widely available.

Another modeling decision that impacts conclusions from this analysis is our treatment of uncertainty. There are different schools of thought on how to deal, or not deal, with parametric uncertainty in population simulation (Heard *et al.* 2013). At one end of the spectrum, sampling error is estimated and removed from the total uncertainty associated with simulation results (Morris and Doak 2002). This is sometimes recommended when precise extinction probabilities are a priority outcome of an analysis. At the other end of the spectrum, uncertainty in model inputs is incorporated at every time step in a simulation (in this case, year is the time step), leading to very large uncertainty in results as the simulation proceeds across years (McGowan *et al.* 2011). We took an intermediate approach by incorporating the natural variation and methodological uncertainty in vital rates using Bayesian estimation and posterior predictive distributions for model inputs (Wade 2002), but not incorporating uncertainty at each time step within a single iteration. Our approach leads to uncertainty estimates that fall between those of the other methods. Regardless of which method is chosen, McGowan *et al.* (2011) have shown that each of these three approaches leads to similar mean results.

A relatively simple, stage-based matrix population model for Black Oystercatcher, when populated with empirical vital rates from the literature, yielded results that were consistent with the hypothesis that local Black Oystercatcher populations are increasing in some areas and decreasing in others, and are globally stable. Sensitivity and elasticity analyses suggested that population growth is particularly sensitive to changes in hatching success, fledging success, and their product, pair productivity, and indicate these stages as possible targets for population management should it become necessary. We suggest that pair productivity is a useful measure for monitoring the population growth potential of Black Oystercatcher as it has a straightforward definition and naturally integrates potential correlations between other vital rates. We recommend that modeling activities such as these are updated as more detailed, site-specific information on Black Oystercatcher vital rates becomes available.

## Acknowledgments

We would like to thank the Audubon conservation science team, particularly C. Wilsey, B. Pickens, and C. Burkhalter, and two anonymous reviewers for feedback on early versions of this work. All input data and computing code is available upon request from the corresponding author.

## Literature Cited

*Haematopus bachmani*). No. 155

*in*The Birds of North America ( A. Poole and F. Gill, Eds.). Academy of Natural Sciences, Philadelphia, Pennsylvania; American Ornithologists' Union, Washington, D.C. Google Scholar

*brms:*an R package for Bayesian multilevel models using Stan. Journal of Statistical Software 80: 1. Google Scholar

*Haematopus palliatus*) population growth by targeting nesting season vital rates. Waterbirds 40: 44–54. Google Scholar

*Haematopus palliatus*) populations. Waterbirds 40: 32–43. Google Scholar

*Haematopus bachmani*. M.S. Thesis, University of Victoria, Victoria, British Columbia. Google Scholar

*popbio*package in R. Journal of Statistical Software 22: 11. Google Scholar

*bachmani*) conservation action plan, v. 1.1. Unpublished report, International Black Oystercatcher Working Group, Manomet Center for Conservation Sciences, Manomet, Massachusetts. Google Scholar

*Haematopus*

*bachmani*. International Wader Studies 20: 83–96. Google Scholar

*in*The Ecology, Status, and Conservation of Marine and Shoreline Birds on the West Coast of Vancouver Island ( K. Vermeer, R. W. Butler and K. H. Morgan, Eds.). Occasional Paper Number 75, Canadian Wildlife Service, Sydney, British Columbia. Google Scholar

*in*Population Viability Analysis ( S. R. Beissinger and D. McCullough, Eds.). University of Chicago Press, Chicago, Illinois. Google Scholar

*Haematopus bachmani*in California. Marine Ornithology 42: 49–56. Google Scholar

*Haematopus palliatus*) in Virginia, USA. Waterbirds 40: 55–71. Google Scholar

*(Haematopus bachmani)*Population Dynamics," Waterbirds 41(2), 115-221, (1 June 2018). https://doi.org/10.1675/063.041.0202