The Black Oystercatcher (Haematopus bachmani) is a long-lived shorebird that breeds annually from May–August along the west coast of North America, from Baja California, Mexico, to the Aleutian Islands, Alaska, USA. Breeding pairs establish and aggressively defend territories, generally along rocky coastlines, and forage exclusively in the intertidal zone. Black Oystercatchers prey on a wide variety of invertebrate species, although limpets (Lottiidae spp.) and mussels (Mytilus spp.) make up the majority of their diet (Andres and Falxa 1995).
There is concern over the conservation status of the Black Oystercatcher (Tessler et al. 2014) because of relatively restricted habitat requirements (Tessler et al. 2014), small population size (8300–12,500; Andres et al. 2012), vulnerability to a number of coastal environmental threats, e.g., oil spill, introduced predators, human disturbance, and ocean acidification (Sloan and Bartier 2006), and expected contraction of breeding range due to climate change (31% by 2080; Langham et al. 2015). There is also interest in monitoring Black Oystercatcher breeding populations because of their role as indicators of overall intertidal ecosystem health (Sloan and Bartier 2006, Carlson-Bremer et al. 2010) and recovery of the endangered Northern Abalone (Haliotis kamtschatkana; Bergman et al. 2013).
To predict vulnerability to climate change impacts, to identify important areas for conservation activities, i.e., monitoring or protection, and to improve understanding of the conservation status of the Black Oystercatcher, there is a need for improved breeding distribution models (“geospatial habitat models,” “geospatial description of good-to-optimal breeding habitat”) as a basis for prediction (Tessler et al. 2014, Weinstein et al. 2014). These models may be referred to broadly as species distribution models (SDMs). SDMs statistically or physiologically relate species occurrence data to environmental covariates (predictors) to predict species distributions (Guisan and Zimmermann 2000) and are increasingly being used as tools for conservation management in coastal environments (e.g., Lindegarth et al. 2014). To date, no SDM has been developed to predict distribution of Black Oystercatcher breeding pairs at a landscape scale.
Several studies have compared habitat features of Black Oystercatchers measured in situ at known breeding sites vs. random sites (e.g., Vermeer et al. 1992), or have investigated the relationship between reproductive success and habitat quality (e.g., Hazlitt and Butler 2001). However, lack of environmental data covering unsampled shorelines in these studies precludes prediction. McFarland (2010) compared remotely sensed habitat features at known breeding sites vs. random sites in Alaska, but did not generate predictions at unsurveyed shoreline. In California (Weinstein et al. 2014) and Alaska (A. Poe and B. A. Andres, 1999, unpublished manuscript), mean breeding density at different shoreline types was calculated to estimate population size and identify suitable habitat within a GIS. The British Columbia Marine Conservation Analysis program collated and mapped observed breeding pairs from various historical and current surveys at a regional scale, but these maps (available online) do not provide information at unsurveyed shorelines. Finally, probability of observation models developed by the Breeding Bird Atlas of British Columbia predict breeding Black Oystercatcher distribution from a set of distal predictors (latitude, longitude, aspect, elevation, slope) across B.C. at 10 km resolution, a scale generally not useful for landscape-level conservation management (Richmond et al. 2015).
Here, we present a SDM for predicting Black Oystercatcher breeding pair occurrence at a landscape scale (i.e., 10–200 km extent, 100 m resolution; Wisz et al. 2013) in southeastern Haida Gwaii, British Columbia, Canada. British Columbia is estimated to hold 10–20% of the global breeding population (Tessler et al. 2014), approximately one-third of which are thought to inhabit Haida Gwaii (Harfenist et al. 2002). The study area serves as a useful location to test various aspects of model and predictor performance because of availability of a rich set of occurrence and environmental data. Also, because the habitat is relatively undisturbed and in the core of the Black Oystercatcher range, the SDM should provide insights into its breeding habitat requirements throughout much of its range (within the limits of sampled environmental space). Prior to model fitting, we developed a conceptual model for breeding pair occurrence based on a review of relevant literature on breeding biology and habitat; feedback was received from several experts and local practitioners (Dalgarno 2016). This feedback informed our predictor selection, model interpretation, and improved our understanding of the model’s limitations. In general we expect that important exogenous factors influencing breeding pair occurrence are either directly or indirectly related to (1) prey abundance and availability, and (2) reproductive failure and mortality of adults. Endogenous and observer processes influencing observed occurrence are discussed, although not modeled explicitly.
The main objectives of this study are the following: (1) to test the use of several GIS-derived and remotely sensed predictors to predict breeding pair occurrence at a landscape scale; (2) to identify important species-environment associations through interpretation of species-response curves and relative influence of individual predictors on model performance; (3) to evaluate performance (discrimination) and spatial and temporal transferability of the model; and (4) to predict distribution at unsurveyed shoreline in southeastern Haida Gwaii.
The distribution of Black Oystercatcher breeding pairs was predicted on 1423 km of shorelines in southeastern Haida Gwaii, an archipelago located approximately 100 km west of the British Columbia mainland (Fig. 1). Most shoreline lay within Gwaii Haanas National Park Reserve, National Marine Conservation Area Reserve, and Haida Heritage Site (Gwaii Haanas) or K'uuna Gwaay Conservancy and Haida Heritage Site (K'uuna Gwaay; also known as Laskeek Bay). The study area was confined to the Hecate Strait ecosection (Zacharias et al. 1998), as far north as Cumshewa Inlet. These boundaries were chosen to avoid extrapolation into unsampled geographic and environmental space: most notably, higher wave exposures and different oceanographic and climatic conditions west of the study area and increasing human disturbance north of the study area. The Hecate Strait ecosection is characterized by shallow to photic water depth (over 20% 0–20 m), high wave exposure, high subsurface relief, low current speeds, and a mixture of sand and hard bottom substrate (Zacharias et al. 1998). Within the study area, the shoreline is primarily bedrock (38%) and mixed bedrock/sediment (26%; Howes et al. 1994). Of 846 islands in the study area, 80% are less than 1 ha, 32% are less than 0.1 ha, and 28% are forested. Dominant wind directions recorded at Cumshewa Island, the nearest weather station, are southeasterly and, to a lesser extent, northwesterly.
We developed a loosely coupled, vector-based (polyline) modeling approach. Methods used to create model units were informed by the need to integrate ShoreZone attributes and other GIS data sets. ShoreZone used 1:5000 aerial video captured in 1991–1992 to create 1:50,000 scale digital physical shoreline maps delineating discrete, systematically described “units” of homogenous shoreline type (34 coastal classes). Units were also assigned a biological wave exposure class (six ordinal classes), and a designation of “absent,” “patchy,” or “continuous” for “biobands” of highly visible species or species assemblages (for details on how qualitative ShoreZone classes were defined see Sloan and Bartier 2006). To maintain the boundaries of ShoreZone units whilst creating a model unit at higher and more consistent resolution (ShoreZone units range from 10 m to 13,000 m in length), we split ShoreZone units into “segments” with a length as close to 100 m as possible using an automated procedure in a GIS (e.g., a “unit” of 321 m was split into 3 segments of 107 m). The model unit size was chosen to approximate mean Black Oystercatcher breeding territory size (Nysewander 1977, Hazlitt 2001). All occurrence and environmental data, which were of various data types (Table 1), were aggregated to each model unit in a GIS.
Occurrence data were compiled from two available surveys completed 22–29 June 2005 and 6–10 June 2010 in two separate regions, by Canadian Wildlife Service (CWS) and Laskeek Bay Conservation Society (LBCS), respectively (Fig. 1). This timing coincides with a period of relatively minimal permanent breeding pair ingress or egress as territories are usually established mid-May and fledging occurs July–August; despite difference in survey date, the use of both data sets was therefore justified (Vermeer et al. 1992, Hazlitt 1999, Hipfner et al. 2012). Of 1423 km of shoreline within the study area, LBCS and CWS surveys covered approximately 164 km (12%) and 381 km (27%) in the northern (Juan Perez sound and Laskeek Bay) and southern (south of Juan Perez) parts of the study area, respectively. Note that survey routes covered both poor and good quality habitat and were completed prior to commencement of rat eradication efforts in Gwaii Haanas.
In both CWS and LBCS surveys, at least two observers scanned the shoreline by boat for territorial Black Oystercatcher breeding pairs. Boats remained approximately 50 m from land and travelled at a speed of less than 11 km/hr. Surveys were not conducted in rainy (except in light rain) or windy/choppy (> 15 knots) conditions to improve detectability. Birds were deemed territorial if they were observed to exhibit two or more key behavioral criteria: sentry behavior, e.g., one or both birds sitting in a location within the territory where they have a good view of surroundings; alarm behavior, e.g., crouching, scurrying, to draw attention away from nest; group size, e.g., commonly as a pair, particularly when alarmed; location, e.g., clearly associated with potential breeding habitat. Upon finding a territorial pair, shoreline was surveyed by foot to search for evidence of a scrape, chicks, or eggs. Any areas previously occupied by a breeding pair were also searched on foot. In CWS surveys, more poor quality habitat was surveyed, e.g., sand beaches, although these sections were scanned from a farther distance (up to 200 m) and at greater boat speed (Moira Lemon, July 2016, personal communication).
During surveys, GPS coordinates were taken at the location of a territorial pair or nest-site. GPS coordinates that were not available for 11 of 91 sites in CWS surveys were determined post hoc from satellite imagery and field notes. Coordinates were imported into a GIS and assigned to the nearest shoreline segment (see below) using a tool within ArcGIS 10.2. Any shoreline segment with at least one territorial breeding pair was deemed present. Any shoreline segment that was surveyed but did not have a territorial breeding pair present was deemed absent. In 12 cases, breeding pairs were “moved” to an incorrect segment. This occurred when ShoreZone shoreline resolution was too low to accurately represent small islets. The actual nest site location, i.e., small islet, was visible from high-resolution satellite imagery, but was absent in the ShoreZone data. These cases were removed prior to analysis.
Fourteen potential predictors were selected a priori (Table 1), following development of a conceptual model (Dalgarno 2016). Predictors were (1) linked as directly as possible to some known or hypothesized biological mechanism; (2) appropriate to the scale of study; and (3) available for the study area at sufficient quality, coverage, and resolution.
Several predictors were available as attributes within the ShoreZone shoreline mapping system (ShoreZone): biological wave exposure index (BioExp) comprised 5-level ordinal estimates of wave exposure (protected, semiprotected, semiexposed, exposed, very exposed), validated with in situ biological data; shoreline type (ShoreType) comprised 10 classes, which were shoreline types reclassified from 34 coastal classes (Harper et al. 1994; Appendix 1); slope of the intertidal area (Slope) comprised 3-level ordinal estimates (5–20, 20–50, > 50); occurrence of mussel beds (Mytilus californianus) on a 3-level ordinal scale (absent, patchy, continuous), visible from aerial video (Mussel); occurrence of rockweed beds (Fucus spp.) on a 3-level ordinal scale (absent, patchy, continuous), visible from aerial video (Fucus); Area of islands (IslandArea) was calculated from island polygons.
A quantitative wave exposure index (Fetch) was available from the haidawave R package, which was codeveloped by the lead author of this study (Dalgarno and Thorley 2017). This index, which uses algorithms developed by Murtojärvi et al. (2007), is the sum of fetch calculated at 10 radiating lines for each point spaced 10 m along the shoreline. Distance to forest cover (TreeDist) was determined from treeline, i.e., transition from trees to bare rock /shrub/grass, digitized at 1:5000 scale from 2010 cloud-free Digitalglobe satellite imagery freely available in ArcGIS 10.2 at resolution ranging from 0.5 m to 2.5 m. Distance to human disturbance (HumanDist) was calculated by identifying sites of heavy boat traffic and human habitation, i.e., from tourism, seasonal accommodations, and monitoring or research activities. Intertidal area within a circular buffer with a radius of 50 m (IT50), 200 m (IT200), 500 m (IT500), and 1000 m (IT1000) were calculated from intertidal area polygons available from Canadian Hydrographic Service. Different circular buffer sizes were intended to account for various hypotheses regarding Black Oystercatcher foraging distance. Data on occurrence of rats (Rattus rattus or Rattus norvegicus) at islands in 2010 (RatStatus) were obtained from Gwaii Haanas and LBCS. Occurrence was determined through use of camera-trap surveys. Islands with unknown presence/absence were left as NA in the data set, and were handled appropriately by the statistical model (see below).
We hypothesized that IslandArea, TreeDist, HumanDist, and RatStatus are related to nest failure or mortality; wave exposure indices (Fetch and BioExp) are related to both prey abundance/availability and reproductive failure/mortality (via nest washout); Slope, intertidal area (IT50, IT200, IT500, IT1000), Mussel, Fucus, and ShoreType are related to prey abundance/availability and in the case of ShoreType, nesting material.
All statistical analyses were conducted in R version 3.4.1 (R Development Core Team 2017). Results may be reproduced from publicly available code archived within a Github repository (Dalgarno 2017). We fit models using boosted regression trees (BRTs), a machine-learning ensemble method that creates and averages many classification trees (Breiman 2001) using a boosting algorithm. After the initial tree is built, subsequent trees model the residuals, with diminishing overall contribution. The main strengths of BRT are the following: (1) high predictive performance relative to other commonly used SDM methods (Elith et al. 2006, Leathwick et al. 2006, De'ath 2007, Palialexis et al. 2011, Valle et al. 2013); (2) ability to incorporate large number of predictors, including categorical predictors; (3) ability to model nonlinear species response curves and complex interactions; and (4) ability to handle missing data in predictors, i.e., by using surrogate splits, whereby missing values are filled-in with values of a correlated variable. These advantages are balanced by several drawbacks that are generally true for machine-learning methods: (1) tendency toward model overfitting, which can lead to poor model transferability (but see Heikkinen et al. 2012, Bahn and McGill 2013); and (2) limited model interpretability. In general, this trade-off is justified when the goal of the study is prediction, not testing biological hypotheses.
Model fitting and prediction was performed using the gbm package (Ridgeway 2006) and with several functions in the dismo package (Hijmans et al. 2012). Three model parameters are user-defined: learning rate specifies the rate at which the contribution of each subsequent tree to the overall model shrinks; tree complexity specifies the number of nodes per tree, which in effect controls the maximum level of interaction between predictors; and bag fraction determines the proportion of observed data randomly selected for each tree. Altogether, these parameters determine the number of trees used to fit the model and can be adjusted to optimize model performance or, if necessary, limit computation time. Parameters were generally set according to recommendations in Elith et al. (2008). For all models, we used a bag fraction of 0.5 and maximum tree complexity of 3. Because any bag fraction less than one introduces model stochasticity (De'ath 2007), all model summary statistics, e.g., evaluation or interpretation, were calculated from the mean of 10 model iterations. Learning rate was adjusted until the optimal number of trees exceeded 1000, as calculated through a 10-fold cross validation procedure implemented with the gbm.step function in the dismo package.
Prevalence of breeding pairs was low, which led to an unequal balance of presence to absence segments (166:5220; 3%). Low prevalence (especially < 10%) reduces BRT model performance unless corrected (Barbet-Massin et al. 2012). We down-weighted absence data according to prevalence in the training set so that presence and absence data had equal weight in model fitting.
Species-response curves and relative influence scores were used to interpret the association between individual predictors and occurrence. Response curves were determined using partial dependence plots, which plot the effect of different predictor values on the response, while all other predictors are held at their mean (Elith et al. 2008). Relative influence (RI) of each predictor was determined following a method outlined in Friedman (2001). Each time a predictor is selected for splitting, the squared improvement on the model is summed, and averaged across all trees. This value is then normalized so that all predictor RI scores sum to 100, with higher numbers indicating higher influence (Elith et al. 2008). Because model unit size was variable, segment length (SegLength) was included as a model predictor to test its influence.
Several procedures were implemented to reduce risk of model overfitting. First, the aforementioned 10-fold cross validation procedure in the gbm.step function determined the optimal number of trees to use before predictive performance on out-of-bag data, i.e., 10% of data withheld for testing, declined. Second, if any predictors had pairwise Spearman’s correlation coefficient (rho) greater than 0.7, the predictor deemed least ecologically relevant or accurate was removed (Dormann et al. 2013; Table 2). This resulted in the removal of three predictors (BioExp, IT200, IT500); all retained predictors had rho less than 0.6. Finally, noninformative predictors were removed using the gbm.simplify function, within the dismo R package. In this procedure, the model is progressively simplified by dropping the least important predictor and refitting models until predictive performance (average CV error from 10-fold CV) is less than the original model (see Elith et al. 2008 for details). This procedure was run for each of 10 model iterations and the final predictor set excluded any predictors that were removed in five or more iterations.
We assessed the ability of models to discriminate between observed presence and absence using Area Under the receiver operator characteristic Curve (AUC: Fielding and Bell 1997). AUC does not require selection of a threshold to reclassify continuous predicted values into binary presence/absence. In models with good discrimination, higher predicted values will tend to be associated with observed presence and lower predicted values will tend to be associated with observed absence, with minimal overlap (Lawson et al. 2014). AUC value of 1 indicates perfect discrimination and AUC value of 0.5 indicates that the model performed no better than random. AUC of 0.9 indicates excellent, 0.8–0.9 very good, 0.7–0.8 satisfactory, and < 0.7 indicates poor discrimination ability (Hosmer and Lemeshow 2000). All reported AUC values were derived from 10-fold cross-validation.
Evaluation of a model’s ability to predict independent data, i.e., data not used for model training, is generally recommended when sufficient occurrence data are available (Araújo and Guisan 2006). We used three methods to assess predictive ability (indicated by AUC): (1) we trained the model on all data and tested predictive ability using 10-fold cross-validation; (2) we trained the model on 2005 occurrence data (~75% of segments) and tested predictive ability on 2010 occurrence data; and (3) we trained the model on the northernmost 75% of presence and northernmost 75% of absence segments and tested predictive ability on the remaining 25% of segments in each category. Finally, spatial autocorrelation of model residuals was tested using Moran’s I, calculated with the lctools R package (Kalogirou 2017).
A total of 221 territorial breeding pairs were located in the surveys: 118 breeding pairs in 2005 and 103 in 2010. A nest with eggs or chicks was located for 81% of breeding pairs in 2005 (n = 96) and 87% of breeding pairs in 2010 (n = 90). Of the 5386 total segments surveyed, 3% (n = 166) were occupied by at least one territorial breeding pair. Breeding pair density on surveyed shoreline was 0.4 pairs/km.
Approximately 50% of breeding pairs (n = 80) were observed on islands < 10,000 m² (1 hectare), despite these islands composing only 15% of surveyed segments. Only 21 (13%) and 54 (33%) breeding pairs were observed within 20 m and 50 m of a treeline, respectively; 121 breeding pairs (73%) occurred on bedrock (shore types 1 and 2) and 151 (91%) occurred on bedrock or mixed bedrock (shore types 1–6). In general, the range of environmental space within the study area was well covered by surveyed segments, although the 2005 survey sampled a broader range of habitats than the 2010 survey (Table 3 and Table 4).
The final model included eight predictors: TreeDist, IslandArea, Fetch, ShoreType, IT50, IT1000, SegLength, RatStatus. Four predictors (Slope, Fucus, Mussel, HumanDist) were removed from the final predictor set during the model simplification process. The final model had high discrimination ability when evaluated on all data through 10-fold cross-validation (AUC = 0.886). Discrimination ability was reduced when tested on independent data, especially on temporally independent data (2005/2010 temporal partition AUC = 0.829; north/south spatial partition AUC = 0.858; Table 5). There was evidence of moderate positive spatial autocorrelation (Moran’s I = 0.34, p-value = 0).
Distance to treeline (TreeDist) had the greatest influence on the model (RI = 41.5%), followed by Island area (IslandArea; RI = 36.7%), wave exposure (Fetch; RI = 5.1%), shoreline type (ShoreType; RI = 5%), intertidal area within 50 m (IT50; RI = 3.8%), rat occurrence (RatStatus; RI = 3%), and intertidal area within 1000 m (IT1000; RI = 1.6%). Model unit size (SegLength) was retained in the final model, although its influence was relatively low (RI = 3.1%; Table 6).
Partial dependence plots revealed complex nonlinear species response curves (Figs. 2 and 3). Breeding pair likelihood of occurrence increased with distance to the treeline (especially > 15 m), although this effect plateaued beyond ~100 m (Fig. 2, Panel A). Probability of occurrence tended to be higher on smaller islands, although not on islands less than ~2000 m². Islands greater than 10 hectares (100,000 m²) had a negative association with occurrence (Fig. 2, Panel B). Shoreline types that breeding pairs tended to prefer included wide rock, wide rock with sand and gravel, narrow rock, and gravel (Fig. 3, Panel A). Probability of occurrence generally increased with wave exposure, although both very low and very high exposures had a negative influence (Fig. 2, Panel C). Intertidal area within 50 m had a positive effect on breeding pair occurrence up to values of ~6000 m², i.e., > 75% of area covered in intertidal, beyond which the trend was reversed (Fig. 2, Panel D). Effect of intertidal area within 1000 m was generally negative, although this plateaued beyond ~50,000 m² (Fig. 2, Panel F). Breeding pairs tended to occur on islands without rats present (Fig. 3, Panel B). Finally the effect of model unit size (SegLength) on occurrence was generally positive, although the strength of the effect decreased with segment lengths around 100 m (Fig. 2, Panel E).
In Figure 4, we present a map of predicted probability of occurrence, with values generated from the final model (Panel A). As an example, we also show predicted probability of occurrence on Bischof Islands (Fig. 4, Panel B), Murchison/Faraday Island group (Fig. 4, Panel C), and Swan/Bolkus Island group (Fig. 4, Panel D), with known occurrence of breeding pairs displayed as white dots. From these enlarged areas, it is clear that the model correctly predicted low density of breeding pairs on Murchison and Faraday mainland and high density on Bischoff’s Islands. As well, the model generally correctly predicted high occurrence on small, rocky offshore islands with higher wave exposure.
The main output of this study was a map of predicted probability of occurrence of Black Oystercatcher breeding pairs at ~100 m shoreline segments in southeastern Haida Gwaii. Segments approximate mean territory size. Probability of occurrence, which is on a continuous scale (bounded by 0 and 1), can be converted to predicted presence/absence using a threshold, although this may lead to an unnecessary loss of information (Guillera-Arroita et al. 2015). If predictions of presence/absence are required, choice of threshold should depend on the application (Fielding and Bell 1997, Lawson et al. 2014) because different thresholds will maximize different aspects of model performance (Liu et al. 2005).
In this study we aimed to test the ability of several GIS data sets, hypothesized to be proxies for some underlying biological mechanisms, to predict breeding pair occurrence. Although causation cannot be determined from correlative models, species-response curves may inform biological hypotheses that can be tested through experimentation (Dormann et al. 2012). Species-response curves identified here closely matched what was expected based on literature reviewed.
In general, predictors hypothesized to be associated with nest failure/mortality (IslandArea, TreeDist, RatStatus), not prey abundance/availability (Fetch, IT50, IT1000, ShoreType), had higher influence on the models. Tendency for breeding pairs to occur at small islands, far from treeline has been observed previously in Haida Gwaii (Vermeer et al. 1992) and elsewhere (Hazlitt 2001, McFarland 2010). We suggest that these predictors act as proxies for occurrence of a wide variety of predators known to cause nest failure or mortality, including Peregrine Falcon (Falco peregrinus; Hipfner et al. 2012); Gulls (Larus glaucescens; Hartwick 1973, Nysewander 1977), Crows (Corvus caurinus; Hartwick 1973, Nysewander 1977, Hipfner et al. 2012), black bears (Ursus americanus), Common Ravens (Corvus corax), wolverines (Gulo gulo), and Bald Eagles (Haliaeetus leucocephalus; Morse et al. 2006). We expect that breeding pairs avoid nesting in areas that increase their vulnerability to nest failure through the process of natural selection. They are also able to adapt habitat use within several years of environmental change, as revealed by experimental removal of predators and human disturbance from islands (Ainley and Lewis 1974, Nysewander 1977, Byrd et al. 1997, Croll et al. 2015).
Nest failure may also be attributed to washout, i.e., from waves (Vermeer et al. 1992), extreme high tides (Morse et al. 2006), and severe storms (Hartwick 1973, Nysewander 1977). We found some evidence to support the hypothesis that these factors may be driving breeding pair occurrence. First, breeding pairs avoided very small islands (< 1500 m²), likely because these were inundated at high tide, or were more vulnerable to inundation from extreme tides or waves. However, it should be noted that very small islets were not represented by the shoreline data set and any breeding pairs that occurred on these were removed from the occurrence data set. Thus, the apparent avoidance of very small islets should be further examined. Second, breeding pairs tended to occur at moderate wave exposures (Fetch), with a sharp decline in predicted probability of occurrence at extreme exposures. Response to wave exposure appears to reflect a trade-off between prey, which is more abundant (Dalgarno 2016) and available, i.e., from gaping mussels (Falxa 1992) at high exposures, and vulnerability to nest washout or risk of injury at the most extreme exposures. For example, Falxa (1992) observed foraging adults avoiding sites at the highest exposures, which he attributed to avoidance of injury by pounding waves (Falxa 1992). To our knowledge, no previous studies have previously identified these nonlinear relationships, i.e., tendency to avoid very small islands and very high exposures.
Breeding pairs generally preferred rock or gravel shoreline types that were wide, and generally avoided sand, mixed sand/gravel, and narrow shoreline types. A tendency for breeding pairs to occur on rocky or gravel shoreline types has been observed throughout its range (Andres 1998, Poe et al. 2009, Weinstein et al. 2014). However, avoidance of wide rock with gravel and narrow rock with sand and gravel shoreline types observed in this study did not follow this general pattern. Prey abundance is likely to be higher at rocky shore types and at low exposure gravel sites (Dalgarno 2016), which we hypothesize to be the reason for preference of these shoreline types. Wide shoreline types are also likely to provide more habitat for prey and facilitate foraging, because Black Oystercatchers cannot forage on vertical surfaces (Lindberg et al. 1998), and breeding pairs are able to deliver more prey to chicks (Hazlitt et al. 2002).
Similarly, we hypothesize that increasing probability of occurrence with increasing IT50 values (intertidal area within 50 m circular buffer) may be explained by prey abundance and prey availability (Lindberg et al. 1998, Hazlitt et al. 2002). The observed decline in probability of occurrence at extreme values of IT50, i.e., very flat areas with > 75% intertidal area, may be attributed to correlation with shoreline type, i.e., estuary, lagoon, or wide sandy beach. Better performance of the IT50 predictor over IT1000 (intertidal area within 1000 m) provides some evidence that foraging occurs primarily within close proximity to nest-sites, i.e., within a breeding territory, despite observations of foraging adults travelling up to a kilometer from breeding territories (Hartwick 1978).
There is particular interest in the effect of introduced rats on breeding Black Oystercatchers in Haida Gwaii, which has seen eradication of rats from Langara Island (Taylor et al. 2000), and more recently, several islands within Gwaii Haanas (Night Birds Returning Project, unpublished data). Although it is generally assumed that rats depredate Black Oystercatcher chicks and eggs (Kurle et al. 2008, Gruman 2013), the effect of rats on Black Oystercatcher breeding pair occurrence is not well known. Eradication of rats on Hawadax Island, Alaska, led to an increase in breeding population from one pair to six pairs over five years (Croll et al. 2015). Within our study area, Gruman (2013) demonstrated that breeding pair density is lower on islands with rats, although that study did not take into consideration other habitat qualities.
In this study, the influence of rat occurrence (RatStatus) was relatively low. There may be several reasons for the apparent lack of influence of rat occurrence on the model: at the scale of this study, TreeDist may act as a general proxy for occurrence of all predators, including rats; Black Oystercatchers may still occur on rat-infested islands, i.e., out of necessity or ignorance, but because our model does not account for reproductive success, we cannot observe and model the detrimental effect of rats, i.e., mismatch between reproductive success and habitat selection (Arlt and Pärt 2007). It is also possible that the quality of our rat occurrence data set was poor. Further studies may look at change in breeding population following rat eradication in our study area, although these data are not yet published.
Models that explicitly account for the observation process are generally called occupancy models. In this study, we could not assume perfect detection probability (Lyons et al. 2012), nor that detection probability was not correlated with a predictor of occurrence. Thus, model predictions may be biased (MacKenzie et al. 2002, Guillera-Arroita et al. 2015). In repeated boat-based count surveys of Black Oystercatcher breeding pairs in Washington, detection probability was 0.75 (95% BCI: 0.42–0.91; Lyons et al. 2012). Thus, two or more repeated surveys on each shoreline segment would increase detection probability. However, repeating surveys is costly given minimal resources typically available for monitoring activities, and thus the required data for occupancy modeling is rarely available.
In this study, only 20% of shoreline was surveyed twice because of logistical constraints. However, even if all shoreline had been surveyed twice, the length of time between each survey (about a month) may have violated assumptions of occupancy modeling: that repeated surveys occur within a closed population (MacKenzie et al. 2002). Thus, although we acknowledge potential bias in our results given imperfect detection, we also suggest that research be undertaken to weigh the costs and benefits of collecting occurrence data suitable for occupancy modeling of Black Oystercatcher breeding pairs.
Another source of error may have been failure to adequately account for territoriality, which may explain moderate spatial autocorrelation in our model. Breeding pairs aggressively defend territories (Hazlitt 2001) and the study scale was chosen to approximately reflect mean territory size (Nysewander 1977, Hazlitt 2001). However, variability in territory size is large (0.79–1.84 ha; Nysewander 1977; 25–245 m length; Hazlitt 2001). Thus, on shorelines where model unit size is smaller than territory size, segments surrounding an observed breeding pair may have high predicted probability of occurrence, i.e., suitable habitat, but observed absence. An alternative approach to take in future studies may be to model breeding density in relation to features of islands (e.g., Heinänen and Von Numers 2009). As well, a sensitivity analysis could be conducted to determine the effect that changing model unit size has on model performance. However, it should be noted that modeling at the island scale or using larger average model unit size (e.g., ~200m) would lead to reduced accuracy of predictors.
Model unit size (SegLength) was retained in the final model, with probability of occurrence increasing on larger segments. However, we suggest that alternative modeling approaches would lead to reduced predictor accuracy (e.g., highest proportion of shore type rather than exact shore type) and that given its relatively low influence on the model (3%), this trade-off may not be justified.
We developed this model with the intention of providing a modeling framework for predicting Black Oystercatcher breeding habitat in Haida Gwaii, B.C., Canada, and to test a set of predictors that may be used in other regions, or at broader extents. Of the predictors retained in the final model, only ShoreType is currently available with B.C.-wide coverage. The Fetch predictor (and additional fetch-based wave exposure indices) is available for Haida Gwaii, and all remaining predictors, which were generated for use in this study, are only currently available for the study area. However, using methods developed here, and given good quality satellite imagery, shoreline and intertidal area data, predictors could be generated for other regions or broader extents at relatively low cost. We do not recommend collecting data on rat occurrence for the purpose of a Black Oystercatcher model because it is expensive to collect and influence on the model was low.
A regional model will require compilation of occurrence data from various surveys that have been conducted (historically or recently), typically from monitoring activities. There are several challenges associated with this requirement. First, it will be necessary to ensure that occurrence data adequately cover geographic and environmental space of prediction area. Although we show that model predictive ability on independent data was reasonably good (AUC > 0.85), model transferability may be an issue, especially in regions that have different environmental conditions. Moreover, our study demonstrated that predictive ability on temporally split data, i.e., trained on 2005 occurrence data, tested on 2010 occurrence data, was relatively poor. We suspect that this may have been caused by differences in survey methods/detection probability between surveys, or less likely, differences in habitat use between years. Temporal partition of data was also spatially partitioned (2005 occurrence data was collected in southernmost 75% of study area), so it was difficult to determine whether poor transferability could be attributed to different habitat use in the southern part of the study area, which is dominated by the presence of Moresby Island, and is generally less accessible. In general, we demonstrate that transferability may be an issue when occurrence data sets are compiled from multiple years/organizations.
Second, compiling occurrence data from historical monitoring activities may be challenging. In this study, we compiled survey data collected by two organizations: CWS and LBCS. For both data sets, nest locations in the form of GPS coordinates were readily available (as well as some absence data in the form of GPS coordinates from CWS). However, the survey route, which was needed to determine absence, was not immediately available and was determined by browsing field notes, reports, and through direct communication with the organizations. The use of historical survey data sets may be difficult if data or metadata are lost, especially in the event that data are used for purposes that they were not originally intended for (as in this study). Our study generally supports calls to digitize and make freely available historical and current ecological data sets, which are collected carefully and often at great cost (Hampton et al. 2012).
Finally, model predictors will have to be selected based on the study scale or region of interest. For example, water temperature and salinity may have greater influence on intertidal species distributions at broader scales (Raffaelli and Hawkins 1996, Nyström Sandman et al. 2013). As well, tendency of breeding pairs to occur at different shoreline types, e.g., rocky vs. gravel, may vary geographically (Nysewander 1977, Byrd et al. 1997, Andres 1998, Weinstein et al. 2014) and different predators may have more or less influence on occurrence dependent on the region, e.g. foxes (Byrd et al. 1997), rats (Croll et al. 2015), raccoons (Vermeer et al. 1992), etc. Finally, although our study, Morse et al. (2006), and Poe et al. (2009) found no influence of low-level human disturbance on occurrence, human disturbance may be an important predictor of occurrence in regions outside of protected areas or on islands permanently inhabited by humans (Ainley and Lewis 1974, Vermeer et al. 1989, Andres and Christensen 2009).
The conservation status and landscape-level distribution of the Black Oystercatcher is largely unknown because of sparse survey data. The need to develop breeding pair distribution models has been identified as a key research priority by leading experts (Tessler et al. 2014). To our knowledge, this study is the first to develop a SDM for predicting probability of occurrence of breeding pairs at unsurveyed shoreline, at a scale relevant to landscape-level management. This modeling approach may be applied to other regions where adequate environmental and occurrence data exist. Further, our coastal SDM framework, which integrates ShoreZone attributes and novel GIS data sets, may be generally useful for predicting distribution of other coastal species in Haida Gwaii and B.C. more broadly.
This research was funded by Gwaii Haanas National Park Reserve, National Marine Conservation Area Reserve, and Haida Heritage Site, and by Janet Mersey, Department of Geography, University of Guelph. Occurrence data were collected by volunteers and staff at Laskeek Bay Conservation Society (northern area) and by Canadian Wildlife Service (southern area). Funding was provided to Laskeek Bay Conservation Society by Gwaii Haanas. We are grateful to Mark Hipfner for helping to supply occurrence data from Canadian Wildlife Service; Jake Pattison and Vivian Pattison from Laskeek Bay Conservation Society for supplying occurrence data and information on field methods; Patrick Bartier for supplying various GIS data sets; Mika Murtojarvi, Marji Puotinen, and Joe Thorley for help developing the wave exposure model; Laurie Wein and Carita Bergman for supplying rat occurrence data; Thomas Nudds, Anthony Gaston, and Jake Pattison for help with research design; Lorne Bennett for comments on an early draft; Brad Andres, Mark Hipfner, Vivian Pattison, Jake Pattison, and Stephanie Hazlitt for comments on the conceptual model. Finally, we thank the Haida Nation for permitting access to its traditional lands and waters.
Ainley, D. G., and J. T. Lewis. 1974. The history of Farallon Island marine bird populations, 1854-1972. Condor 76:432-446. http://dx.doi.org/10.2307/1365816
Andres, B. A. 1998. Shoreline habitat use of Black Oystercatchers breeding in Prince William Sound, Alaska. Journal of Field Ornithology 69:626-634.
Andres, B. A., and R. E. Christensen. 2009. Dramatic changes in the number of Black Oystercatchers nesting in Sitka Sound, Alaska. Wader Study Group Bulletin 116:181-184.
Andres, B. A., and G. A. Falxa. 1995. Black Oystercatcher (Haematopus bachmani). In P. G. Rodewald, editor. The Birds of North America. Cornell Lab of Ornithology, Ithaca, New York, USA. http://dx.doi.org/10.2173/bna.155
Andres, B. A., P. A. Smith, R. I. G. Morrison, C. L. Gratto-Trevor, S. C. Brown, and C. A. Friis. 2012. Population estimates of North American shorebirds, 2012. Wader Study Group Bulletin 119:178-194.
Araújo, M. B., and A. Guisan. 2006. Five (or so) challenges for species distribution modeling. Journal of Biogeography 33:1677-1688. http://dx.doi.org/10.1111/j.1365-2699.2006.01584.x
Arlt, D., and T. Pärt. 2007. Nonideal breeding habitat selection: a mismatch between preference and fitness. Ecology 88:792-801. http://dx.doi.org/10.1890/06-0574
Bahn, V., and B. J. McGill. 2013. Testing the predictive performance of distribution models. Oikos 122:321-331. http://dx.doi.org/10.1111/j.1600-0706.2012.00299.x
Barbet-Massin, M., F. Jiguet, C. H. Albert, and W. Thuiller. 2012. Selecting pseudo-absences for species distribution models: how, where and how many? Methods in Ecology and Evolution 3:327-338. http://dx.doi.org/10.1111/j.2041-210X.2011.00172.x
Bergman, C. M., J. Pattison, and E. Price. 2013. The Black Oystercatcher as a sentinel species in the recovery of the Northern Abalone. Condor 115:800-807. http://dx.doi.org/10.1525/cond.2013.120182
Breiman, L. 2001. Statistical modeling: the two cultures. Statistical Science 16:199-231. http://dx.doi.org/10.1214/ss/1009213726
Byrd, G. V., E. P. Bailey, and W. Stahl. 1997. Restoration of island populations of Black Oystercatchers and Pigeon Guillemots by removing introduced foxes. Colonial Waterbirds 20:253-260. http://dx.doi.org/10.2307/1521691
Carlson-Bremer, D., T. M. Norton, K. V Gilardi, E. S. Dierenfeld, B. Winn, F. J. Sanders, C. Cray, M. Oliva, T. C. Chen, S. E. Gibbs, M. S. Sepu, and C. K. Johnson. 2010. Health assessment of American Oystercatchers (Haemotopus palliatus) in Georgia and South Carolina. Journal of Wildlife Diseases 46:772-780. http://dx.doi.org/10.7589/0090-3558-46.3.772
Croll, D. A., K. M. Newton, M. McKown, N. Holmes, J. C. Williams, H. S. Young, S. Buckelew, C. A. Wolf, G. Howald, M. F. Bock, J. A. Curl, and B. R. Tershy. 2015. Passive recovery of an island bird community after rodent eradication. Biological Invasions 18:1-13.
Dalgarno, S. 2016. Predictive modelling of Black Oystercatcher (Haemotopus Bachmani) breeding pair occurrence and prey abundance in Haida Gwaii, British Columbia. Thesis. University of Guelph, Guelph, Ontario, Canada.
Dalgarno, S. 2017. Bloy-distribution-model. Github repository. [online] URL: https://doi.org/10.5281/zenodo.997576
Dalgarno, S., and J. Thorley. 2017. haidawave: Wave exposure from fetch and wind data. R package version 0.0.1. R Foundation for Statistical Computing, Vienna, Austria. [online] URL: https://github.com/sebdalgarno/haidawave
De'ath, G. 2007. Boosted trees for ecological modeling and prediction. Ecology 88:243-251. http://dx.doi.org/10.1890/0012-9658(2007)88[243:BTFEMA]2.0.CO;2
Dormann, C. F., J. Elith, S. Bacher, C. Buchmann, G. Carl, G. Carré, J. R. G. Marquéz, B. Gruber, B. Lafourcade, P. J. Leitão, T. Münkemüller, C. McClean, P. E. Osborne, B. Reineking, B. Schröder, A. K. Skidmore, D. Zurell, and S. Lautenbach. 2013. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography 36:27-46. http://dx.doi.org/10.1111/j.1600-0587.2012.07348.x
Dormann, C. F., S. J. Schymanski, J. Cabral, I. Chuine, C. Graham, F. Hartig, M. Kearney, X. Morin, C. Römermann, B. Schröder, and A. Singer. 2012. Correlation and process in species distribution models: bridging a dichotomy. Journal of Biogeography 39:2119-2131. http://dx.doi.org/10.1111/j.1365-2699.2011.02659.x
Elith, J., C. H. Graham, R. P. Anderson, M. Dudík, S. Ferrier, A. Guisan, R. J. Hijmans, F. Huettmann, J. R. Leathwick, A. Lehmann, J. Li, L. G. Lohmann, B. A. Loiselle, G. Manion, C. Moritz, M. Nakamura, Y. Nakazawa, J. M. M. Overton, A. Townsend Peterson, S. J. Phillips, K. Richardson, R. Scachetti-Pereira, R. E. Schapire, J. Soberón, S. Williams, M. S. Wisz, and N. E. Zimmermann. 2006. Novel methods improve prediction of species’ distributions from occurrence data. Ecography 29:129-151. http://dx.doi.org/10.1111/j.2006.0906-7590.04596.x
Elith, J., J. R. Leathwick, and T. Hastie. 2008. A working guide to boosted regression trees. Journal of Animal Ecology 77:802-813. http://dx.doi.org/10.1111/j.1365-2656.2008.01390.x
Falxa, G. A. 1992. Prey choice and habitat use by foraging Black Oystercatchers: interactions between prey quality, habitat availability, and age of bird. Dissertation. University of California, Davis, California, USA.
Fielding, A. H., and J. F. Bell. 1997. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environmental Conservation 24:38-49. http://dx.doi.org/10.1017/S0376892997000088
Friedman, J. H. 2001. Greedy function approximation: a gradient boosting machine. Annals of Statistics 29:1189-1232 http://dx.doi.org/10.1214/aos/1013203451
Gruman, C. A. 2013. Context-dependence of a cross-system trophic cascade in Gwaii Haanas, British Columbia. Thesis. Simon Fraser University, Burnaby, British Columbia, Canada.
Guillera-Arroita, G., J. J. Lahoz-Monfort, J. Elith, A. Gordon, H. Kujala, P. E. Lentini, M. A. McCarthy, R. Tingley, and B. A. Wintle. 2015. Is my species distribution model fit for purpose? Matching data and models to applications. Global Ecology and Biogeography 24:276-292. http://dx.doi.org/10.1111/geb.12268
Guisan, A., and N. E. Zimmermann. 2000. Predictive habitat distribution models in ecology. Ecological Modelling 135:147-186. http://dx.doi.org/10.1016/S0304-3800(00)00354-9
Hampton, S. E., J. J. Tewksbury, and C. A. Strasser. 2012. Ecological data in the Information Age. Frontiers in Ecology and the Environment 10:59. http://dx.doi.org/10.1890/1540-9295-10.2.59
Harfenist, A., N. A. Sloan, and P. M. Bartier. 2002. Living marine legacy of Gwaii Haanas. IV: marine mammal baseline to 2003 and marine mammal-related management issues throughout the Haida Gwaii region. Parks Canada Technical Reports in Ecosystem Science 36:164.
Harper, J. R., W. T. Austin, M. Moms, P. D. Reimer, and R. Reitmeier. 1994. Ecological classification of Gwaii Haanas: a biophysical inventory of the coastal resources. Report prepared for Parks Canada, Calgary, Alberta, Canada by Coastal and Ocean Resources Inc., Sidney, British Columbia, Canada.
Hartwick, E. B. 1973. Foraging strategy of the Black Oystercatcher. Thesis. University of British Columbia, Vancouver, British Columbia, Canada.
Hartwick, E. B. 1978. The use of feeding areas outside of the territory of breeding Black Oystercatchers. Wilson Bulletin 90:650-652.
Hazlitt, S. L. 1999. Territory quality and parental behaviour of the black oystercatcher in the Strait of Georgia, British Columbia. Thesis. Simon Fraser University, Burnaby, British Columbia, Canada.
Hazlitt, S. L. 2001. Territory quality and reproductive success of Black Oystercatchers in British Columbia. Wilson Bulletin 113:404-409. http://dx.doi.org/10.1676/0043-5643(2001)113[0404:TQARSO]2.0.CO;2
Hazlitt, S. L., and R. W. Butler. 2001. Site fidelity and reproductive success of Black Oystercatchers in British Columbia. Waterbirds 24:203-207. http://dx.doi.org/10.2307/1522031
Hazlitt, S. L., R. C. Ydenberg, and D. B. Lank. 2002. Territory structure, parental provisioning, and chick growth in the American Black Oystercatcher (Haematopus bachmani). Ardea 90:219-228.
Heikkinen, R. K., M. Marmion, and M. Luoto. 2012. Does the interpolation accuracy of species distribution models come at the expense of transferability? Ecography 35:276-288. http://dx.doi.org/10.1111/j.1600-0587.2011.06999.x
Heinänen, S., and M. Von Numers. 2009. Modelling species distribution in complex environments: an evaluation of predictive ability and reliability in five shorebird species. Diversity and Distributions 15:266-279. http://dx.doi.org/10.1111/j.1472-4642.2008.00532.x
Hijmans, R. J., S. Phillips, J. R. Leathwick, and J. Elith. 2012. dismo: species distribution modeling. R package "dismo" version 1.1-4. R Foundation for Statistical Computing, Vienna, Austria. [online] URL:http://CRAN.R-project.org/package=dismo
Hipfner, M. J., K. W. Morrison, and A.-L. Kouwenberg. 2012. Biology of Black Oystercatchers breeding on Triangle Island, British Columbia, 2003–2011. Northwestern Naturalist 93:145-153. http://dx.doi.org/10.1898/nwn12-02.1
Hosmer, D. W., and S. Lemeshow. 2000. Applied logistic regression. Second edition. John Wiley and Sons, New York, New York, USA. http://dx.doi.org/10.1002/0471722146
Howes, D., J. Harper, and E. Owens. 1994. Physical shore-zone mapping system for British Columbia. Resources Inventory Committee Publication 8:71.
Kalogirou, S. 2017. lctools: Local correlation, spatial inequalities, geographically weighted regression and other tools. R package version 0.2-6. R Foundation for Statistical Computing, Vienna, Austria. [online] URL:http://CRAN.R-project.org/package=lctools
Kurle, C. M., D. A. Croll, and B. R. Tershy. 2008. Introduced rats indirectly change marine rocky intertidal communities from algae- to invertebrate-dominated. Proceedings of the National Academy of Sciences of the United States of America 105:3800-3804. http://dx.doi.org/10.1073/pnas.0800570105
Langham, G. M., J. G. Schuetz, T. Distler, C. U. Soykan, and C. Wilsey. 2015. Conservation status of North American birds in the face of future climate change. PLoS ONE 10:e0135350. http://dx.doi.org/10.1371/journal.pone.0135350
Lawson, C. R., J. A. Hodgson, R. J. Wilson, and S. A. Richards. 2014. Prevalence, thresholds and the performance of presence-absence models. Methods in Ecology and Evolution 5:54-64. http://dx.doi.org/10.1111/2041-210X.12123
Leathwick, J. R., J. Elith, M. P. Francis, T. Hastie, and P. Taylor. 2006. Variation in demersal fish species richness in the oceans surrounding New Zealand: an analysis using boosted regression trees. Marine Ecology Progress Series 321:267-281. http://dx.doi.org/10.3354/meps321267
Lindberg, D. R., J. A. Estes, and K. I. Warheit. 1998. Human influences on trophic cascades along rocky shores. Ecological Applications 8:880-890. http://dx.doi.org/10.1890/1051-0761(1998)008[0880:HIOTCA]2.0.CO;2
Lindegarth, M., U. Bergström, J. Mattila, S. Olenin, M. Ollikainen, A.-L. Downie, G. Sundblad, M. Bučas, M. Gullström, M. Snickars, M. von Numers, J. R. Svensson, and A. K. Kosenius. 2014. Testing the potential for predictive modeling and mapping and extending its use as a tool for evaluating management scenarios and economic valuation in the Baltic Sea (PREHAB). Ambio 43:82-93. http://dx.doi.org/10.1007/s13280-013-0479-2
Liu, C., P. M. Berry, T. P. Dawson, and R. G. Pearson. 2005. Selecting thresholds of occurrence in the prediction of species distributions. Ecography 28:385-393. http://dx.doi.org/10.1111/j.0906-7590.2005.03957.x
Lyons, J. E., J. A. Royle, S. M. Thomas, E. Elliott-Smith, J. R. Evenson, E. G. Kelly, R. L. Milner, D. R. Nysewander, and B. A. Andres. 2012. Large-scale monitoring of shorebird populations using count data and N-mixture models: Black Oystercatcher (Haematopus bachmani) surveys by land and sea. Auk 129:645-652. http://dx.doi.org/10.1525/auk.2012.11253
MacKenzie, D. L., J. D. Nichols, G. B. Lachman, S. Droege, J. A. Royle, and C. A. Langtimm. 2002. Estimating site occupancy rates when detection probabilities are less than one. Ecology 83:2248-2255. http://dx.doi.org/10.1890/0012-9658(2002)083[2248:ESORWD]2.0.CO;2
McFarland, B. A. 2010. Habitat characteristics of Black Oystercatcher breeding territories. Thesis. University of Alaska, Fairbanks, Alaska, USA.
Morse, J. A., A. N. Powell, and M. D. Tetreau. 2006. Productivity of Black Oystercatchers: effects of recreational disturbance in a national park. Condor 108:623-633. http://dx.doi.org/10.1650/0010-5422(2006)108[623:POBOEO]2.0.CO;2
Murtojärvi, M., T. Suominen, H. Tolvanen, V. Leppänen, and O. S. Nevalainen. 2007. Quantifying distances from points to polygons: applications in determining fetch in coastal environments. Computers & Geosciences 33:843-852. http://dx.doi.org/10.1016/j.cageo.2006.10.006
Nysewander, D. R. 1977. Reproductive success of the Black Oystercatcher in Washington State. Thesis. University of Washington, Seattle, Washington, USA.
Nyström Sandman, A., S. A. Wikström, M. Blomqvist, H. Kautsky, and M. Isaeus. 2013. Scale-dependent influence of environmental variables on species distribution: a case study on five coastal benthic species in the Baltic Sea. Ecography 36:354-363. http://dx.doi.org/10.1111/j.1600-0587.2012.07053.x
Palialexis, A., S. Georgakarakos, I. Karakassis, K. Lika, and V. D. Valavanis. 2011. Prediction of marine species distribution from presence-absence acoustic data: comparing the fitting efficiency and the predictive capacity of conventional and novel distribution models. Hydrobiologia 670:241-266. http://dx.doi.org/10.1007/s10750-011-0673-9
Poe, A. J., M. I. Goldstein, B. A. Brown, and B. A. Andres. 2009. Black Oystercatchers and campsites in Western Prince William Sound, Alaska. Waterbirds 32:423-429. http://dx.doi.org/10.1675/063.032.0307
R Development Core Team. 2017. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. [online] URL: http://www.R-project.org/
Raffaelli, D., and S. Hawkins. 1996. Intertidal ecology. First edition. Chapman & Hall, London, UK. http://dx.doi.org/10.1007/978-94-009-1489-6
Richmond, S., A. R. Couturier, P. D. Taylor, and D. Lepage. 2015. Analyses and mapping. In P. J. A. Davidson, R. J. Cannings, A. R. Couturier, D. Lepage, and C. M. Di Corrado, editors. The atlas of the breeding birds of British Columbia, 2008-2012. Bird Studies Canada, Delta, British Columbia, Canada.
Ridgeway, G. 2006. Generalized boosted regression models. R Package “gbm,” version 2.1.3. R Foundation for Statistical Computing, Vienna, Austria. [online] URL: http://CRAN.R-project.org/package=gbm
Sloan, N. A., and P. M. Bartier. 2006. Living marine legacy of Gwaii Haanas. V: coastal zone values and management around Haida Gwaii. Parks Canada Technical Reports in Ecosystem Science 42:1-413.
Taylor, R. H., G. W. Kaiser, and M. C. Drever. 2000. Eradication of Norway rats for recovery of seabird habitat on Langara Island, British Columbia. Restoration Ecology 8:151-160. http://dx.doi.org/10.1046/j.1526-100x.2000.80022.x
Tessler, D. F., J. A. Johnson, B. A. Andres, S. Thomas, and R. B. Lanctot. 2014. A global assessment of the conservation status of the Black Oystercatcher (Haematopus bachmani). International Wader Study Group 20:83-96.
Valle, M., M. M. van Katwijk, D. J. de Jong, T. J. Bouma, A. M. Schipper, G. Chust, B. M. Benito, J. M. Garmendia, and Á. Borja. 2013. Comparing the performance of species distribution models of Zostera marina: implications for conservation. Journal of Sea Research 83:56-64. http://dx.doi.org/10.1016/j.seares.2013.03.002
Vermeer, K., K. H. Morgan, and G. E. J. Smith. 1989. Population and nesting habitat of American Black Oystercatcher in the Strait of Georgia. Pages 118-122 in K. Vermeer and R. W. Butler, editors. The ecology and status of marine and shoreline birds in the Strait of Georgia, British Columbia, Canada. Canadian Wildlife Service Special Publication, Ottawa, Ontario, Canada.
Vermeer, K., K. H. Morgan, and G. E. J. Smith. 1992. Black Oystercatcher habitat selection, reproductive success, and their relationship with Glaucous-winged Gulls. Colonial Waterbirds 15:14-23. http://dx.doi.org/10.2307/1521350
Weinstein, A., L. Trocki, R. Levalley, R. H. Doster, T. Distler, and K. Krieger. 2014. A first population assessment of Black Oystercatcher (Haematopus bachmani) in California. Marine Ornithology 1972:49-56.
Wisz, M. S., J. Pottier, W. D. Kissling, L. Pellissier, J. Lenoir, C. F. Damgaard, C. F. Dormann, M. C. Forchhammer, J. A. Grytnes, A. Guisan, R. K. Heikkinen, T. T. Høye, I. Kühn, M. Luoto, L. Maiorano, M. C. Nilsson, S. Normand, E. Öckinger, N. M. Schmidt, M. Termansen, A. Timmermann, D. A. Wardle, P. Aastrup, and J. C. Svenning. 2013. The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling. Biological Reviews of the Cambridge Philosophical Society 88:15-30. http://dx.doi.org/10.1111/j.1469-185X.2012.00235.x
Zacharias, M. A., D. E. Howes, J. R. Harper, and P. Wainwright. 1998. The British Columbia marine ecosystem classification: rationale, development, and verification. Coastal Management 26:105-124. http://dx.doi.org/10.1080/08920759809362347