Research Paper

INTRODUCTION

Despite their presence in huge numbers on coasts and islands across the world, seabird populations are challenging to monitor. Seabird colonies are often large, diffuse, geographically isolated, and in some cases include many non-breeding individuals that confound estimates of breeding totals (Mercker et al. 2021). Many seabird species also nest in burrows under variable substrate and across complex topography, which can require survey methods to vary among and even within colonies (e.g., Buxton et al. 2016, Arneill et al. 2019, Lavers et al. 2019). Because of the sporadic and inconsistent nature of these surveys, it has therefore been difficult to combine infrequent and imprecise estimates of colony size into regional models of population change for status and trend assessment.

Adding urgency to these estimation difficulties, numerous colonial seabird populations are known or suspected to be in decline (Paleczny et al. 2015) from an array of threats including fisheries bycatch, offshore energy production, light attraction, pollution, invasive species, and climate change (e.g., Anderson et al. 2011, Sandvik et al. 2012, Regular et al. 2013, Ronconi et al. 2015, Dias et al. 2019, Gilmour et al. 2023, Phillips et al. 2023). Yet, population dynamics among colonies may be synchronized or idiosyncratic depending on the spatiotemporal patterns of these threats, with subsequent effects on the variability of regional population dynamics (Che-Castaldo et al. 2017). For example, large-scale environmental processes can synchronize population dynamics within a region, potentially reducing regional population viability (Palmqvist and Lundberg 1998) but simultaneously improving the statistical precision of population estimates (e.g., through hierarchical modeling that acknowledges the shared population trajectories among colonies). Monitoring population status therefore requires statistical tools that accommodate partially shared trajectories among colonies, while accounting for the sparseness and uncertainty inherent to most seabird monitoring data.

Bayesian hierarchical generalized additive mixed modeling (GAMM) offers a promising framework for estimating regional trajectories and trends of bird populations (Smith and Edwards 2021). This approach can handle the challenges of infrequent and imprecise colony surveys, incorporates flexible population trajectories without assuming fixed patterns of change over time, and allows uncertainty to be fully propagated to estimates of regional-level abundance. These models also share information among monitored colonies (if supported by the data) to improve statistical precision and can provide insights into the spatio-temporal scales at which environmental drivers influence population dynamics. After fitting the model, interval-specific rates of mean annual population change (i.e., trends) can then be derived from any stretch of the regional trajectory time series, with the length dictated by the required application. For example, three-generation trend estimates are required for national population status assessments by the Committee on the Status of Endangered Wildlife in Canada (COSEWIC) and the International Union on the Conservation of Nature (IUCN), while shorter-term (e.g., 10-year) trend estimates can be more useful for assessing species recovery or extinction risk.

To assess the value of this approach for colonial seabird monitoring, we apply Bayesian GAMMs to colony census data from two priority burrow-nesting species in eastern Canada: Leach’s Storm-Petrel (Hydrobates leucorhous, hereafter “storm-petrel”) and Atlantic Puffin (Fratercula arctica, hereafter “puffin”). Storm-petrels are small-bodied pelagic seabirds that nest in small and often concealed burrows. They have a broad distribution across the Northern hemisphere, including 100+ breeding colonies across eastern Canada ranging in size from a few individuals to several million (COSEWIC 2020; Fig. 1). In Canada, a wide variety of threats (Lieske et al. 2020, Pollet et al. 2023) have led to a population decline of over 50% since the 1970s (COSEWIC 2020, d’Entremont et al. 2020, Wilhelm et al. 2015, 2020, Bond et al. 2023) that led to their designation in 2020 as Threatened by COSEWIC. Similar trends have been observed across the Atlantic basin and the species is listed as Vulnerable by IUCN (BirdLife International 2018a). Puffins also nest in burrows, have a similarly diffuse breeding distribution of 100+ colonies across the coasts of eastern Canada, also ranging in size from a few to hundreds of thousands of individuals (Fig. 1), and can often be found on the same islands as those occupied by storm-petrels (Wilhelm et al. 2015). The IUCN designated the Atlantic Puffin as Vulnerable based on ongoing and projected population declines in Europe (BirdLife International 2018b), though the Canadian puffin population has seemingly grown in recent decades (Lowther et al. 2020), possibly linked to the closure of gill-net fisheries in the region (Regular et al. 2013). It follows that these two species provide ideal case studies to demonstrate how a GAMM-based approach can be informed by time series of colony census data of varying precision and frequency across space and time, and can generate robust and transparent regional trajectories and trends with readily interpretable estimates of uncertainty.

We assess the utility of a Bayesian hierarchical approach for estimating seabird population trajectories and trends by (1) validating the approach through simulations of realistic colony monitoring data for sporadically surveyed burrow-nesting species, and (2) applying the method to estimate colony-level and regional population trends based on empirical monitoring data for storm-petrel and puffin in eastern Canada. In the context of ongoing environmental change and increasing pressures on seabirds in the North Atlantic and oceanwide, this is a timely exercise for these case study species and for demonstrating the power of Bayesian methods to optimize hard-won but limited data to meet pressing seabird conservation needs.

METHODS

Case study species and census methods

Storm-petrels and puffins are burrow nesters that have widespread and overlapping breeding distributions in eastern Canada (Fig. 1) but have different life histories and behaviors that could lead to divergent population trends. At the colony, storm-petrels are nocturnally active, while adult puffins are strictly diurnal. Storm-petrels travel for several days and over many hundreds of kilometers to feed on planktonic crustaceans and deep water mesopelagic lanternfish (Myctophidae; Hedd et al. 2009) and raise their chicks for ~65 days (Hedd et al. 2018, d’Entremont et al. 2020), whereas puffins forage locally to feed their chicks small forage fish several times a day over ~40 days (Pratte et al. 2017, Delord et al. 2020). Storm-petrels show high natal dispersal and can recruit in any colony across the North Atlantic (Bicknell et al. 2012), whereas puffins exhibit higher natal philopatry (reviewed in Kersten et al. 2021).

We focused our analysis on data from colonies in eastern Canada that were surveyed at least twice between 1965 and 2024 and where at least 10 burrows were counted on average, resulting in the omission of one small colony. Data were primarily collected as part of Environment and Climate Change Canada’s Seabird Colony Monitoring Program, supplemented with additional survey data from published literature and technical reports where possible. A total of 13 storm-petrel and 28 puffin colonies were included in our analysis (Fig. 1). In Newfoundland, storm-petrels and puffins can nest on the same island but rarely in the same habitat, showing distinct nesting habitat preferences (e.g., storm-petrels prefer fern or forested habitat while puffins prefer grassy slopes; Wilhelm et al. 2015, 2020, Bond et al. 2023). Elsewhere in eastern Canada, storm-petrel and puffin colonies do not overlap, except for Machias Seal Island where they share similar nesting habitat. Thus, the general survey approach for both species is the same, but depending on the size of the colony, methods can vary as follows.

Complete hole count

For sites where the entire island can be searched, surveyors attempt to count all potential burrow entrances (or holes, defined as a hole in the ground that appears to a standing observer to lead to a burrow). This approach was used for smaller colonies in Newfoundland and Labrador (usually < 3000 pairs) with methods previously published in Robertson and Elliot (2002a,b) and Robertson et al. (2002a). Briefly, an island-wide systematic hole count is done using a transect approach and having multiple observers walk next to each other 1–2 meters apart and count all holes encountered, or in the case of Machias Seal Island, counting all holes along transects within quadrats (Durham et al. 2024).

For smaller islands with fewer than 200 holes (e.g., Bacalhao and Tinker Islands) all holes were assessed by inserting an arm (and using a wooden spoon to extend the reach if necessary) and assigning the hole as either being (1) an extra entrance to a burrow, (2) too short to be a burrow (< 30 cm), (3) an empty burrow, (4) an occupied burrow with either an adult and/or egg and/or chick, or (5) unknown (i.e., the observer could not reach the end of the burrow to confirm it being occupied or not). If an island in Newfoundland and Labrador had more than 200 holes, hole assessments were done using a plot approach by laying a rope grid (either 3 X 3 m or 5 X 5 m) in areas occupied by the species (puffin or storm-petrel) and representative of the island. Depending on the island’s size, a range of 5–20 plots were randomly placed in different areas of the island at both the periphery and center of the colony. Hole occupancy rates were calculated as the number of occupied burrows divided by the total number of holes counted with known occupancy status. Each plot was weighted by the number of holes assessed to calculate island-wide occupancy rates. Standard errors were based on the number of plots assessed on each island.

All monitored puffin colonies in Québec are found within the North Shore Migratory Bird Sanctuaries, seven of which were consistently monitored (Fig. 1B). For these colonies, complete hole counts have been largely phased out to minimize disturbance and replaced by counts of adults attending the colony on land and water. Burrow counts continue at the Betchouane Migratory Bird Sanctuary, while a system of transects and quadrats is used to estimate area of occupancy and mean occupied burrow density at the Baie de Brador Migratory Bird Sanctuary (Rail and Chapdelaine 2002, Rail 2021). Storm-petrel colonies in Québec generally consist of a few scattered pairs, with burrows often hidden in dense vegetation where counting holes is inefficient. Storm-petrel surveys at these sites have used playback at night and automated recording units to confirm the presence of the species, but do not provide reliable estimates of colony size; thus, no Québec storm-petrel colonies were used in trend analyses.

Grid or transect approach

For larger colonies where complete hole counts are not feasible, island-wide grids or transect lines were set 25–100 m apart to determine the area occupied by puffins or storm-petrels, and 20–80 plots were assessed to calculate occupied burrow density as described above (see also Robertson and Elliot 2002a, Robertson et al. 2002b, Pollet and Shutler 2018). This approach is suitable for islands where puffins and storm-petrels are nesting in relatively flat areas and thus where correcting for slope is not a concern.

Habitat delineation

An alternative for the grid or transect approach to estimate occupied area, the use of a Geographic Information System (GIS) approach has proven to increase the efficiency of surveys and provides more accurate estimated sloped areas occupied by puffins and storm-petrels that can then be applied to the estimated occupied burrow density assessed through plots as described in the other two approaches (e.g., Wilhelm et al. 2015, 2020, d’Entremont et al. 2020, Bond et al. 2023). Briefly, estimating the area of habitat occupied by puffins or storm-petrels was done using a hand-held GPS and walking around the boundary of the nesting area and/or using high resolution imagery to delineate the various habitat types utilized by each species. On large and convoluted colonies where puffins and storm-petrels nest on slopes, maps of contour lines were also incorporated into the GIS approach to provide more accurate estimated areas on slopes (e.g., Wilhelm et al. 2015, 2020, Bond et al. 2023).

Description of statistical model

We analyzed survey count data by adapting Bayesian hierarchical state-space GAMMs described in Smith and Edwards (2021). These models decompose spatio-temporal variation in population counts into contributions from biological processes of interest (i.e., temporal changes in colony-level population indices) and observation processes (i.e., variation due to imprecision in survey counts). The full model is described in Figure 2.

In our application to seabird colony counts, the observation process model assumes survey counts (Y_i,t; the estimated number of occupied burrows in a colony i in year t) arise from an over-dispersed Poisson process with mean λ_i,t (Equation 1). The model assumes λ_i,t is a normally distributed random variable with mean equal to an expected count (E_i,t) and a variance term (ε²_i,t) that characterizes the magnitude of sampling error during a survey. This explicitly recognizes that there are inevitable discrepancies between the numbers of birds estimated from surveys and the total number of birds that would be counted with a true colony census (i.e., ε_i,t is the standard error associated with a survey).

For plot-based surveys of burrow-nesting seabirds, ε²_i,t can be calculated via design-based estimators of spatial population totals (Horvitz and Thompson 1952). However, because many surveys lack empirical estimates of sampling error, we included an additional sub-model to estimate the magnitude of sampling variation for cases where ε²_i,t was missing (i.e., where raw plot-level data to calculate this quantity were unavailable). In this dataset, surveys that reported both quantities showed a strong positive relationship between log⁡(ε_i,t) and log(Y_i,t). To leverage this empirical pattern, we embedded a log-linear regression within the Bayesian model to impute missing values of ε_i,t (Equation 3). Reported standard errors were treated as known data, while only surveys lacking error estimates were assigned latent ε_i,t values drawn from the regression model.

The population process model decomposes expected counts E_i,t into contributions from two terms (Equation 4) through a logarithmic link function where μ_i,t represents the general shape of each colony’s trajectory across years, and δ_i,t represents annual departures of population indices from the colony-level smooths (sometimes called “process variance”). The primary goal of our analysis is to describe temporal patterns in μ_i,t and its variation among colonies. We accomplished this using hierarchical GAMMs that fit smoothed temporal trajectories to each colony (Equation 5) using an approach described in Smith and Edwards (2021). The term μ_i,t therefore includes effects of colony-level intercepts (β_0,i) and a semi-parametric “smoothed” temporal process defined by a generalized additive function involving the product of a series of colony-level smoothing coefficients (β_i,k) and a design matrix (X_t,k) that is constructed from a series of K basis functions multiplied by a smoothness penalty. In this study, we used thin plate splines with K=12 basis functions to construct X_t,k, providing flexibility to model a wide range of complex trajectories (Appendix 1).

Following Smith and Edwards (2021), we estimated the smoothing parameters β_i,k as random effects arising from a shared distribution across colonies (Equation 5). This allows the shape of colony-level trajectories to be partially conserved among colonies, if supported by the data, which can improve model predictions for colonies with sparse data. This is also a valuable feature if large-scale environmental processes (e.g., climate oscillations) affect an entire regional population in a similar way, leading to trajectories with similar shapes.

We modeled annual deviations δ_i,t from the temporal smooth as normally distributed (Equation 4). These deviations can be driven by a wide variety of population processes, including random annual variation in breeding propensity or environmental effects affecting recruitment. While δ_i,t describes genuine year-to-year fluctuations in colony abundance around its temporal smooth, we focus inference for status and trend assessment on changes in μ_i,t, which represents longer-term and more persistent changes in the expected annual abundance at colonies.

Estimation of population indices, trajectories, and trends

For each colony, we calculated estimates of relative annual abundance, or population index, using Equation 6. This definition of population index removes the effects of observation error and random annual process variation, yielding an index represented by each colony’s long-term temporal smooth. The addition of the log-normal variance term in Equation 6 re-scales the indices such that they are centered on the mean of the observed counts (see Sauer and Link 2011). We calculated the regional population total as the sum of annual population indices across all the colonies within each year of study, using Equation 7. This constructs a regional population trajectory that weights colonies according to their relative abundance; larger colonies have a stronger impact on the estimated regional population dynamics than smaller colonies.

We defined the population trend from year t_start to t_end as the geometric mean annual rate of population change over that time interval (Sauer and Link 2011), which can be expressed as an annual percent rate of change using Equation 8. We reported trends in two ways. First, to describe how trends have changed across the study, we calculated sequential “rolling” 5-year trends (centered on each year). Second, we reported trend over the last three generations for each species (44 years for storm-petrels, 33 years for puffins; Bird et al. 2020) to illustrate how this approach can inform status assessments (e.g., for COSEWIC).

Bayesian model specification

We fit statistical models in a Bayesian framework using JAGS version 4.3.0 (Plummer 2003) through the jagsUI library version 1.5.2 (Kellner et al. 2019) within the R programming language version 4.0.2 (R Core Team 2024). We used the “jagam” function from the mgcv package (Wood 2017) to prepare a template for constructing hierarchical GAMMs within the JAGS language, which we manually modified to incorporate colony-level random effects and observation error.

We specified vague priors on all model parameters. All variance parameters were priors of Uniform(0,25), except for process variance which we assigned a prior of Lognormal(–1.6, 624); this prior remained highly vague but improved model convergence. After a burn-in of 1,000,000 iterations, we stored every 50,000th sample until we accumulated 1000 posterior samples from each of three MCMC chains. We assessed chain convergence by visual examination of MCMC traceplots and by evaluating that the Gelman-Rubin convergence statistic was close to 1 for all model parameters. We conducted a posterior predictive check on each species’ model by comparing discrepancy measures (root mean square error or RMSE) from simulated datasets based on the posterior distribution of the model parameters to discrepancy measures from the actual empirical data, whereby a well-fitting model will produce simulated data that closely resembles the actual data. We calculated the “Bayesian p-value” as the proportion of RMSEs from the actual data that were greater than RMSEs from the simulated data. Bayesian p-values between 0.3 and 0.7 indicate that the observed data appear to be consistent with “perfect datasets” that were simulated directly from the fitted model, indicating good model fit. Code to repeat these analyses are available in Appendix 4.

Simulations

We conducted a series of simulations to confirm that the statistical model could correctly estimate regional population trends with minimal bias and appropriate credible interval coverage across a range of hypothetical scenarios of population change, using only limited and imprecise survey data. Each simulation assumed there were I = 15 colonies that experienced stochastic population dynamics across a T = 50 year period. Population dynamics at each colony were affected by two autocorrelated environmental covariates: (1) an “unshared” covariate that affected each colony independently, and (2) a “shared” covariate that affected all colonies simultaneously, thereby imposing correlations among all 15 trajectories.

We simulated each environmental driver as a first-order Markov process (i.e., a random walk), starting with a value of 0, resulting in a temporal trajectory for each environmental driver that influenced the temporal dynamics at each colony. In all simulations, ENV_UNSHARED,i,t was modeled as Normal(ENV_{UNSHARED,i,t–1}, 0.004). We considered two scenarios for ENV_SHARED,t. The first scenario omitted ENV_SHARED,t altogether, allowing each colony to have fully independent trajectories described entirely by ENV_UNSHARED,i,t. The second scenario modeled the shared environmental driver as Normal(ENV_SHARED,t–1, 0.01), thereby imposing a moderate correlation among colony-level trajectories because ENV_SHARED,t varied more through time than ENV_UNSHARED,i,t. Annual log-scale expected counts log(E_i,t) were modeled as the sum of log(E_i,t–1), both environmental covariates, and additional non-Markovian annual noise drawn from Normal(0,0.01).

Initial population indices for each colony were drawn from a lognormal distribution according to E_i,1~Lognormal(log(10000), 0.25). We assumed each colony was surveyed in 2–6 randomly selected years across the 50-year period, with at least one survey occurring in the first and last 10 years of the simulation. Observation error was simulated according to Equation 2, with parameters chosen from the posterior mean from the empirical analysis for storm-petrels (Appendix 1). Counts were simulated as Poisson random variables using Equation 1. The magnitude of observation error was then retained for 65% of surveys, to replicate the empirical data where estimates of error were not available for all surveys.

In each simulation, we calculated annual regional abundance using Equation 7 (Fig. 2) and estimated the “true” smoothed regional trajectory by fitting a GAMM to the full simulated series of R_t using the mgcv package in R. The estimated regional long-term trend was calculated across all years using Equation 7. We then fit the Bayesian statistical model to the simulated observed data (i.e., only 2–6 imprecise surveys at each colony), and compared the estimated trend to the true trend. An example of simulated population trajectories for 15 colonies are illustrated in Appendix Figure S1.1, model-estimated trajectories are shown in Figure S1.2, and both the true and simulated regional trajectory are illustrated in Figure S1.3. We repeated this exercise 1000 times for each unique simulation scenario. Finally, we calculated the mean bias in trend estimates as the difference between estimated and true trend, and credible interval coverage as the proportion of simulations where the 95% credible intervals on the regional trend estimates contained the true simulated regional trend.

RESULTS

Simulations

Simulations confirmed that the model could recover regional trends with negligible bias (less than 1%), regardless of whether colonies had correlated or fully independent temporal trajectories (Appendix 1; Fig. S1.4). Credible interval coverage was slightly lower than nominal; 95% credible intervals contained the true regional trend in approximately 91% and 80% of simulations for 50-year and 10-year regional trend estimates, respectively (Fig. S1.4). Estimates were more likely to be negatively biased when true trends were highly positive and vice versa because of shrinkage imposed by random effect structures in the models.

Storm-petrels

Survey counts in at least two years were available from 13 storm-petrel colonies in Newfoundland and Labrador, Nova Scotia, and New Brunswick over the 57-year period from 1966 to 2024 (Fig. 1A). One colony had 15 surveys within this period, two colonies had 6–7 surveys, five had 5 surveys, four had 3 surveys, and one had 2 (Fig. 3), for 67 surveys in total (9% of all possible surveys years at 13 colonies). Of the 67 storm-petrel surveys, 36 (54%) had an estimate of the error of the count (Fig. 3).

Posterior predictive checks yielded a Bayesian p-value of 0.73 and predicted counts were highly correlated with observed counts, indicating adequate model fit (Appendix 2; Fig. S2.1). Population trajectories were highly consistent among storm-petrel colonies (Fig. 3), reflected by a low estimate of (mean = 0.07; 95% CRI = 0.003 to 0.19), thereby allowing the model to share trajectory information among colonies and improve estimates in data-sparse colonies. In general, storm-petrel colony abundances increased until the late 1980s, followed by declines until the mid to late 2000s, with apparent population stability in the most recent decade. The regional storm-petrel population trajectory experienced a similar temporal pattern (Fig. 4A & C), resulting in a trend of -1.03% per year (95% CRI = -2.13 to -0.05) over the most recent three generations, with a 98% chance the population declined (Fig. 4E).

Puffins

For puffins, survey counts in at least two years were available from 29 colonies in Québec, Newfoundland, and Labrador, Nova Scotia, and New Brunswick over the 56-year period from 1965 to 2024 (Fig. 1B). One colony had 20 surveys within this period, four had 11 surveys, three had 6–10 surveys, 14 had 3–5 surveys, and seven had 2 (Fig. 5), for 158 surveys in total. Surveys were conducted for at least one colony in 37 years across the 60-year study period. Of the 158 puffin surveys, 63 (40%) had an estimate of the standard error of the count (Fig. 5).

Posterior predictive and plots of observed versus predicted values indicated adequate model fit (Bayesian p-value = 0.43; Fig. S2.2 in Appendix 2). In contrast to storm-petrels, puffin trajectories varied substantially among colonies (Fig. 5). Several colonies experienced large population increases (e.g., Gull Island, Machias Seal Island) while others remained relatively stable (e.g., Great Island, North Shore Migratory Bird Sanctuaries) or declined (e.g., Gannet Clusters 3 and 4). Accordingly, the random effect structure of our model shared little information on trajectory shape among puffin colonies (mean estimate of = 0.43; 95% CRI = 0.31 to 0.58), leading to large uncertainty in some colony estimates (e.g., Baccalieu Island) and reducing the precision of the regional trajectory estimate in some periods of the time series. The regional puffin trajectory was highly uncertain until the late 1970s, after which there was strong evidence of population growth in the late 1980s, and ongoing growth or stability thereafter (Fig. 4B & D). Over the most recent three generations, the regional puffin trend estimate was +1.66% per year (95% CRI = +0.21 to +4.16), with a 99% chance the population trend was positive (Fig. 4F).

DISCUSSION

Seabird colony monitoring data are hard-earned and critical for many conservation and management purposes, including status and trend assessments, impact assessments, and as indicators of marine ecosystem health (Bird et al. 2021). However, sound inference requires accounting for complex features of their population biology (i.e., non-linear population trajectories, highly variable counts among years, large differences in relative abundance among colonies; Paleczny et al. 2015) as well as observation processes (i.e., substantial observation error, highly intermittent and imbalanced survey data; Bird et al. 2021). Our model accommodates these complexities and fully propagates these multiple sources of error and uncertainty to estimates of regional population size and trend via Bayesian methods.

For colonial seabirds in general, regional population trajectories tend to depend heavily on the dynamics of a few extremely large colonies where routine monitoring is an obvious priority for regional status assessments. For storm-petrels in eastern Canada, more than 70% of the population currently breeds on Baccalieu Island, while more than 70% of the regional puffin population breeds on three islands (Baccalieu, Great, and Gull Islands). However, monitoring smaller colonies remains crucial for a comprehensive understanding of spatio-temporal population dynamics. In our storm-petrel case study, data from smaller colonies helped to constrain estimates of the population trajectory at the large Baccalieu colony in the earlier years of the time series. Our dataset included 21 surveys at 11 storm-petrel colonies prior to 1984 (i.e., prior to the first reliable survey available at the large Baccalieu colony). Our model was able to provide an informed “hindcast” of the trajectory at Baccalieu and therefore the regional trajectory prior to 1984 based on observed population increases at both the larger and smaller colonies during this period, and evidence throughout the remainder of the time series that colonies within the region share similar trajectories. The estimated trajectory at Baccalieu prior to 1984 is also consistent with indirect paleoenvironmental evidence that the Baccalieu colony grew substantially during these earlier years (Duda et al. 2020).

Trajectory estimates from our model are less certain in periods of the time series when fewer surveys were conducted, especially at the beginning and end of each colony’s time series. In addition, although the GAMM structure in our model is effective at describing complex population trajectories across years in which surveys were conducted (i.e., for interpolation), GAMMs can produce unreasonable predictions outside the range of data used for model training (i.e., for extrapolation and forecasting; e.g., Elith et al. 2010). Care should therefore be taken to identify data-rich segments of estimated population trajectories for conservation-based uses. Our model was developed to inform the State of Canada’s Birds Report (Birds Canada and Environment and Climate Change Canada 2024), which summarizes the best available information on national and regional avian population trajectories since 1970. However, we caution that estimates prior to 1980 are based on extremely sparse data and fitted population trajectories are largely extrapolations for many colonies during that period. Similarly, although 10-year trends are commonly reported by the Canadian and U.S. Federal governments, those shorter-term estimates would have much lower reliability if the largest colonies had not been recently surveyed. In this study, we reported trends over the most recent 3-generations (1978–2023 for storm-petrels; 1991–2024 for puffins) to illustrate how model results would be interpreted relative to a conservation objective such as the IUCN or COSEWIC’s species status assessment criteria. Despite higher uncertainty in the most recent years of the estimated trajectories, our 3-generation trends provided high confidence that the regional puffin population in eastern Canada experienced growth (99% probability) while the storm-petrel population declined (97% probability).

Although a strength of this approach is that data from colonies of varying sizes, census frequency, and count accuracy can contribute to estimates of regional trajectories and trends, discretion should still be used if including data from colonies that are very small. Small populations can be subject to high demographic stochasticity and local scale extinction and colonization processes that differ fundamentally from those at larger populations (Lande 1993). Inclusion of these extremely small colonies in our hierarchical model could affect estimates of spatio-temporal random effects (i.e., estimates of temporal process variance and/or shared trajectory patterns), with consequent effects on estimates of larger-scale regional population dynamics. This phenomenon occurred in our dataset, where the storm-petrel colony on Machias Seal Island is over 30 times smaller than the next smallest colony, and appeared to experience fundamentally different temporal dynamics than the 12 other colonies (Fig. 3). Additionally, that colony had a more complete time series than any other colony (particularly from 1965 to 1980) but had no estimates of survey precision to accompany the counts. Omitting that colony from the analysis resulted in stronger sharing of trajectory information among the larger colonies in the early years of the study, resulting in a stronger signal of regional population growth from 1965 to 1980 (Appendix 3, Figs. S1.1, S1.2). Nevertheless, the estimate of 3-generation trend from 1978 to 2023 was nearly identical whether that colony was included or not (Fig. S1.3), owing to the more complete survey counts and highly consistent trajectories from other colonies in more recent years.

By sharing trajectory information among colonies, our analysis provides indirect evidence of the scales at which shared environmental drivers and dispersal among colonies affect our case study populations. Although our analysis did not evaluate synchrony among year-to-year fluctuations (owing to extremely limited time series at individual colonies), 12 of 13 Canadian storm-petrel colonies in our study experienced similar population trajectories over a 58-year period. This is significant because synchrony in population dynamics is a strong predictor of metapopulation extinction risk (Heino et al. 1997). Leach’s storm-petrel populations throughout the North Atlantic Ocean are presumed to act as a single metapopulation linked by high dispersal of immature birds recruiting into non-natal colonies across the breeding range; this theory is supported by a homogeneous genetic structure across colonies (Bicknell et al. 2012). During the breeding season, individuals also forage across vast oceanic areas (Hedd et al. 2018), potentially exposing individuals from all colonies to shared environmental conditions. The apparently coordinated population changes we observed are thus consistent with the influence of large-scale environmental forces (e.g., shifts in ocean productivity, prey availability, or climate patterns) on the regional population.

Simultaneously, demographic studies suggest that adult survival in these colonies is currently between 0.82 and 0.88 (Calvert et al. 2024), which is substantially lower than on Canada’s Pacific coast (Rennie et al. 2020). In contrast, hatching (68–91%) and fledging success (90–91%) at several major colonies remain high (Pollet et al. 2021), suggesting that conditions at breeding sites are not currently limiting. Declines in northeast Atlantic populations (COSEWIC 2020, Deakin et al. 2021, Pollet et al. 2021, Burnell et al. 2023) further support the hypothesis that population trajectories are shaped by broad-scale pressures. Although our correlative, non-mechanistic model cannot isolate the drivers of change, the combination of similar colony trajectories across the region, widespread survival declines, and evidence of shared trans-Atlantic population trends points to the importance of international collaboration in addressing large-scale threats to this species.

Local-scale and spatially variable processes likely play a larger role in determining colony-level dynamics for puffins in eastern Canada, where trajectories varied considerably among colonies (Fig. 5). At the three largest colonies in eastern Newfoundland (Gull, Great, and Baccalieu Islands), population increases are likely due to population recovery after reductions in incidental mortality of breeding adults following a moratorium on groundfish gillnet fishing in the early 1990s (Regular et al. 2013). Interestingly, these increases were predicted by Calvert and Robertson (2002) on the basis of demographic monitoring and population projection analysis. The smaller colonies on the northeast coast of Newfoundland would also have benefited from the fishing moratorium, as this area sustained an important cod fishery for centuries. However, colony-specific trends are more variable in this area, possibly linked to avian and mammalian predation pressures that can have significant impacts on smaller populations (e.g., Russell and Montevecchi 1996, Zabala Belenguer et al. 2024). In Newfoundland and Labrador, Arctic foxes (Alopex lagopus) and polar bears (Ursus maritimus) are known to access seabird islands during the breeding season (Birkhead and Nettleship 1995, Robertson and Elliot 2002a, Burke et al. 2011), which likely has a strong effect on productivity in some years. Similarly, in Québec North Shore Migratory Bird Sanctuaries, seabird declines on some islands have been associated with the presence of mammalian predators (Rail 2021), including Red (Vulpes vulpes) and Arctic foxes, and black bears (Ursus americanus). Puffin reproductive success is also tied to the availability of capelin (Mallotus villosus) and sand lance (Ammodytes spp.). In eastern Newfoundland, capelin biomass has declined and become less predictable (Buren et al. 2014, Pedersen et al. 2017), resulting in variable annual chick growth and/or mass breeding failures (Wilhelm et al. 2013, Fitzsimmons et al. 2017, Wilhelm et al. 2021). However, similar reductions in reproductive success have been reported in the Gulf of Maine near the Machias Seal Island colony, where rapidly warming waters have reduced the quality of prey fed to chicks (Scopel et al. 2019, Diamond et al. 2021, Major et al. 2021) but colony abundance has apparently increased (Major et al. 2024, this study; Fig. 5).

Monitoring wildlife populations comes with many challenges, and survey methods need to be adapted to the biology of the species and logistical realities of data collection. Colonial seabirds are no exception, as their highly clumped distributions in remote locations and dangerous habitats pose challenges for conducting balanced, design-based surveys among colonies. Furthermore, frequentist-based analyses are not amenable to analyzing sparse, inconsistent, and highly skewed survey data, and require significant simplifying assumptions to generate trend estimates. The approach used here appears to be highly robust to the vagaries of seabird monitoring data (Appendix 1) and is capable of synthesizing data with variable precision owing to changes in survey methodologies over time, differences in data sparsity among colonies, huge differences in colony sizes, and potentially shared non-linear trajectories among multiple colonies. Instead of calculating a specific trend over time, our Bayesian GAMM approach models all the data in a holistic way, allowing for non-linear dynamics, appropriate weighting of colony size, and sharing of trajectories across colonies if the data support that sharing. Trend estimates can be extracted at any interval required, and longer-term non-linear dynamics can be examined. Non-linear trends will become increasingly important as time series mature, and especially for seabirds that are influenced by natural oceanographic cycles.

Although our hierarchical Bayesian analysis provides robust estimates of population trends, we do not prescribe specific recommendations for optimal survey intervals at individual colonies. Survey frequency decisions require balancing multiple, often competing considerations: the precision required to inform management decisions and the urgency of detecting ecologically meaningful population changes, the demographic significance and connectivity of each colony within the broader metapopulation, and the logistical and financial constraints inherent to long-term monitoring programs (Nichols and Williams 2006). Given that most seabird colony surveys in Canada are conducted by government agencies operating under resource constraints, survey frequency must be strategically prioritized among competing monitoring objectives. We recommend that survey design decisions be guided by clearly articulated conservation objectives and, where feasible, informed by formal monitoring design frameworks that explicitly evaluate the value of additional information relative to survey costs (e.g., adaptive monitoring designs or value-of-information analyses; McDonald-Madden et al. 2010, Canessa et al. 2015).

Given the demonstrated suitability of this method to seabird populations, a future priority is to assess this modeling approach at larger spatial scales and with a broader suite of species that have different nesting habits and survey methods. That these methods have proven useful to model other bird populations with very different types of data input (Birds Canada and Environment and Climate Change Canada 2024) suggests that this approach is likely to be applicable for many species. Simulation testing with different data types and different underlying colony-level dynamics will also be needed before expanding the method more broadly. Of course, any statistical inference is only as good as the data input, so future applications will depend on continued investment in quality monitoring data. Finally, modifications of this modeling framework are likely to be fruitful. In particular, local extinction and colonization processes could be incorporated into the model to accommodate situations in which colonies are temporarily abandoned, and spatio-temporal covariates could be incorporated to simultaneously improve the precision of trend estimates and test hypotheses about environmental drivers. Taken together, accurate assessments of trends and robust identification of drivers of seabird populations will focus conservation efforts on populations of highest concern and on mitigating the key impacts faced by these birds.

LITERATURE CITED

Anderson, O. R., C. J. Small, J. P. Croxall, E. K. Dunn, B. J. Sullivan, O. Yates, and A. Black. 2011. Global seabird bycatch in longline fisheries. Endangered Species Research 14(2):91-106. https://doi.org/10.3354/esr00347

Arneill, G. E., C. M. Perrins, M. J. Wood, D. Murphy, L. Pisani, M. J. Jessopp, and J. L. Quinn. 2019. Sampling strategies for species with high breeding-site fidelity: a case study in burrow-nesting seabirds. PLoS ONE 14(8):e0221625. https://doi.org/10.1371/journal.pone.0221625

Bicknell, A. W. J., M. E. Knight, D. Bilton, J. B. Reid, T. Burke, and S. C. Votier. 2012. Population genetic structure and long-distance dispersal among seabird populations: implications for colony persistence. Molecular Ecology 21(12):2863-2876. https://doi.org/10.1111/j.1365-294X.2012.05558.x

Bird, J. P., R. Martin, H. R. Akçakaya, J. Gilroy, I. J. Burfield, S. T. Garnett, A, Symes, J. Taylor, Ç. H. Şekercioğlu, and S. H. Butchart. 2020. Generation lengths of the world’s birds and their implications for extinction risk. Conservation Biology 34(5):1252-1261. https://doi.org/10.1111/cobi.13486

Bird, J. P., B. K. Woodworth, R. A. Fuller, and J. D. Shaw. 2021. Uncertainty in population estimates: a meta-analysis for petrels. Ecological Solutions and Evidence 2(3):e12077. https://doi.org/10.1002/2688-8319.12077

BirdLife International. 2018a. Leach’s Storm-Petrel (Hydrobates leucorhous). The IUCN Red List of Threatened Species 2018:e.T132438298A132438484. https://dx.doi.org/10.2305/IUCN.UK.2018-2.RLTS.T132438298A132438484.en

BirdLife International. 2018b. Atlantic Puffin (Fratercula arctica). The IUCN Red List of Threatened Species 2018:e.T22694927A132581443. https://dx.doi.org/10.2305/IUCN.UK.2018-2.RLTS.T22694927A132581443.en

Birds Canada and Environment and Climate Change Canada. 2024. The state of Canada’s birds report. Birds Canada, Port Rowan, Ontario, Canada, Environment and Climate Change Canada, Gatineau, Québec, Canada. https://doi.org/10.71842/8bab-ks08

Birkhead, T. R., and D. N. Nettleship. 1995. Arctic fox influence on a seabird community in Labrador: a natural experiment. Wilson Journal of Ornithology 107: 397-412.

Bond, A. L., S. I. Wilhelm, D. W. Pirie-Hay, G. J. Robertson, I. L. Pollet, and J. Arany. 2023. Quantifying gull predation in a declining Leach’s Storm-Petrel (Hydrobates leucorhous) colony. Avian Conservation and Ecology 18(1):5. https://doi.org/10.5751/ACE-02388-180105

Burke, C. M., A. Hedd, W. A. Montevecchi, and P. Regular. 2011. Effects of an Arctic fox visit to a low Arctic seabird colony. Arctic 64(3):302-306. https://doi.org/10.14430/arctic4120

Burnell, D., A. J. Perkins, S. F. Newton, M. Bolton, D. T. Tierney, and T. E. Dunn. 2023. Seabirds count. A census of breeding seabirds in Britain and Ireland (2015-2021). Lynx Nature Books, Barcelona, Spain.

Buxton, R. T., A. M. Gormley, C. J. Jones, and P. O. B. Lyver. 2016. Monitoring burrowing petrel populations: a sampling scheme for the management of an island keystone species. Journal of Wildlife Management 80(1):149-161. https://doi.org/10.1002/jwmg.994

Calvert, A. M., and G. J. Robertson. 2002. Using multiple abundance estimators to infer population trends in Atlantic Puffins. Canadian Journal of Zoology 80:1014-1021. https://doi.org/10.1139/z02-081

Canessa, S., G. Guillera-Arroita, J. J. Lahoz-Monfort, D. M. Southwell, D. P. Armstrong, I. Chadès, R. C. Lacy, and S. J. Converse. 2015. When do we need more data? A primer on calculating the value of information for applied ecologists. Methods in Ecology and Evolution 6(10):1219-1228. https://doi.org/10.1111/2041-210X.12423

Che-Castaldo, C., S. Jenouvrier, C. Youngflesh, K. T. Shoemaker, G. Humphries, P. McDowall, L. Landrum, M. M. Holland, Y. Li, R. Ji, and H. J. Lynch. 2017. Pan-Antarctic analysis aggregating spatial estimates of Adélie penguin abundance reveals robust dynamics despite stochastic noise. Nature Communications 8(1):832. https://doi.org/10.1038/s41467-017-00890-0

Committee on the Status of Endangered Wildlife in Canada (COSEWIC). 2020. COSEWIC assessment and status report on the Leach’s Storm-Petrel (Oceanodroma leucorhoa) Atlantic population in Canada 2020. COSEWIC, Ottawa, Ontario, Canada. https://www.canada.ca/en/environment-climate-change/services/species-risk-public-registry/cosewic-assessments-status-reports/leachs-storm-petrel-2020.html

Deakin, Z., E. S. Hansen, R. Luxmoore, R. J. Thomas, M. J. Wood, O. Padget, R. Medeiros, R. Aitchison, M. Ausden, R. Barnard, et al. 2021. Decline of Leach’s Storm Petrels Hydrobates leucorhous at the largest colonies in the northeast Atlantic. Seabird 33:74-106. https://doi.org/10.61350/sbj.33.74

Delord, K., C. Barbraud, D. Pinaud, B. Letournel, B. Jaugeon, H. Goraguer, P. Lazure, and H. Lormée. 2020. Movements of three alcid species breeding sympatrically in Saint Pierre and Miquelon, northwestern Atlantic Ocean. Journal of Ornithology 161:359-371. https://doi.org/10.1007/s10336-019-01725-z

d'Entremont, K. J. N., L. Minich Zitske, A. J. Gladwell, N. K. Elliott, R. A. Mauck, and R. A. Ronconi. 2020. Breeding population decline and associations with nest site use of Leach’s Storm-Petrels on Kent Island, New Brunswick from 2001 to 2018. Avian Conservation and Ecology 15(1):11. https://doi.org/10.5751/ACE-01526-150111

Dias, M. P., R. Martin, E. J. Pearmain, I. J. Burfield, C. Small, R. A. Phillips, O. Yates, B. Lascelles, P. Garcia Borboroglu, and J. P. Croxall. 2019. Threats to seabirds: a global assessment. Biological Conservation 237:525-537. https://doi.org/10.1016/j.biocon.2019.06.033

Duda, M. P., G. J. Robertson, J. E. Lim, J. A. Kissinger, D. C. Eickmeyer, C. Grooms, L. E. Kimpe, W. A. Montevecchi, N. Michelutti, J. M. Blais, and J. P. Smol. 2020. Striking centennial-scale changes in the population size of a threatened seabird. Proceedings of the Royal Society B 287(1919):20192234. https://doi.org/10.1098/rspb.2019.2234

Durham, S. E., S. P. Saunders, A. W. Diamond, T. V. Riecke, and H. L. Major. 2024. Divergent population trends of two sympatric auk species in the rapidly warming Gulf of Maine. Ecology and Evolution 14(11):e70495. https://doi.org/10.1002/ece3.70495

Elith, J., M. Kearney, and S. Phillips. 2010. The art of modelling range-shifting species. Methods in Ecology and Evolution 1(4):330-342. https://doi.org/10.1111/j.2041-210X.2010.00036.x

Gilmour, M., S. Borrelle, L. Elliott, R. Okawa, and A. Rodríguez. 2023. Pollution—lights, plastics, oil, and contaminants. Pages 177-216 in L. Young and E. VanderWerf, editors. Conservation of marine birds. Academic, London, UK. https://doi.org/10.1016/B978-0-323-88539-3.00012-1

Hedd, A., W. A. Montevecchi, G. K. Davoren, and D. A. Fifield. 2009. Diets and distributions of Leach’s Storm-Petrel Oceanodroma leucorhoa before and after an ecosystem shift in the Northwest Atlantic. Canadian Journal of Zoology 87:787-801. https://doi.org/10.1139/Z09-060

Hedd, A., I. L. Pollet, R. A. Mauck, C. M. Burke, M. L. Mallory, L. A. McFarlane Tranquilla, W. A. Montevecchi, G. J. Robertson, R. A. Ronconi, D. Shutler, S. I. Wilhelm, and N. M. Burgess. 2018. Foraging areas, offshore habitat use, and colony overlap by incubating Leach’s Storm-Petrels Oceanodroma leucorhoa in the Northwest Atlantic. PLoS ONE 13:e0194389. https://doi.org/10.1371/journal.pone.0194389

Heino, M., V. Kaitala, E. Ranta, and J. Lindström. 1997. Synchronous dynamics and rates of extinction in spatially structured populations. Proceedings of the Royal Society of London Series B: Biological Sciences 264(1381):481-486. https://doi.org/10.1098/rspb.1997.0069

Horvitz, D. G., and D. J. Thompson. 1952. A generalization of sampling without replacement from a finite universe. Journal of the American Statistical Association 47(260):663-685. https://doi.org/10.1080/01621459.1952.10483446

Kellner, K., M. Meredith, and M. K. Kellner. 2019. Package ‘jagsUI’. A wrapper around ‘rjags’ to streamline ‘JAGS’ analyses. R Package Version 1(1).

Kersten, O., B. Star, D. M. Leigh, T. Anker-Nilssen, H. Strøm, J. Danielsen, S. Descamps, K. E. Erikstad, M. G. Fitzsimmons, J. Fort, et al. 2021. Complex population structure of the Atlantic Puffin revealed by whole genome analyses. Communications Biology 4(1):922. https://doi.org/10.1038/s42003-021-02415-4

Lande, R. 1993. Risks of population extinction from demographic and environmental stochasticity and random catastrophes. American Naturalist 142(6):911-927. https://doi.org/10.1086/285580

Lavers, J. L., I. Hutton, and A. L. Bond. 2019. Changes in technology and imperfect detection of nest contents impedes reliable estimates of population trends in burrowing seabirds. Global Ecology and Conservation 17:e00579. https://doi.org/10.1016/j.gecco.2019.e00579

Lieske, D. J., L. M. Tranquilla, R. A. Ronconi, and S. Abbott. 2020. “Seas of risk”: assessing the threats to colonial-nesting seabirds in Eastern Canada. Marine Policy 115:103863. https://doi.org/10.1016/j.marpol.2020.103863

Lowther, P. E., A. W. Diamond, S. W. Kress, G. J. Robertson, K. Russell, D. N. Nettleship, G. M. Kirwan, D. Christie, C. J. Sharpe, E. Garcia, and P. F. D. Boesman. 2020. Atlantic Puffin (Fratercula arctica), version 1.0. in S. M. Billerman, editor. Birds of the world. Cornell Lab of Ornithology, Ithaca, New York, USA. https://doi.org/10.2173/bow.atlpuf.01

McDonald-Madden, E., P. W. Baxter, R. A. Fuller, T. G. Martin, E. T. Game, J. Montambault, and H. P. Possingham. 2010. Monitoring does not always count. Trends in Ecology & Evolution 25(10):547-550. https://doi.org/10.1016/j.tree.2010.07.002

Mercker, M., N. Markones, K. Borkenhagen, H. Schwemmer, J. Wahl, and S. Garthe. 2021. An integrated framework to estimate seabird population numbers and trends. Journal of Wildlife Management 85(4):751-771. https://doi.org/10.1002/jwmg.22026

Nichols, J. D., and B. K. Williams. 2006. Monitoring for conservation. Trends in Ecology & Evolution 21(12):668-673. https://doi.org/10.1016/j.tree.2006.08.007

Paleczny, M., E. Hammill, V. Karpouzi, and D. Pauly. 2015. Population trend of the world’s monitored seabirds, 1950-2010. PLoS ONE 10(6):e0129342. https://doi.org/10.1371/journal.pone.0129342

Palmqvist, E., and P. Lundberg. 1998. Population extinctions in correlated environments. Oikos 83:359-367. https://doi.org/10.2307/3546850

Phillips, R. A., J. Fort, and M. P. Dias. 2023. Conservation status and overview of threats to seabirds. Pages 33-56 in L. Young and E. VanderWerf, editors. Conservation of marine birds. Academic, London, UK. https://doi.org/10.1016/B978-0-323-88539-3.00015-7

Plummer, M. 2003. JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling. In Proceedings of the 3rd International Workshop on Distributed Statistical Computing 124(125.10):1-10.

Pollet, I. L., A. L. Bond, A. Hedd, C. E. Huntington, R. G. Butler, and R. Mauck. 2021. Leach’s Storm-Petrel (Hydrobates leucorhous), version 1.1. in Birds of the world. Cornell Lab of Ornithology, Ithaca, New York, USA. https://doi.org/10.2173/bow.lcspet.01.1

Pollet, I. L., A. K. Lenske, A. N. M. A Ausems, C. Barbraud, Y. Bedolla-Guzmán, A. W. J. Bicknell, M. Bolton, A. L. Bond, K. Delord, A. W. Diamond, D. A. Fifield, C. Gjerdrum, L. R. Halpin, E. S. Hansen, A. Hedd, R. Hoeg, H. L. Major, R. A. Mauck, G. McClelland, L. McFarlane Tranquilla, W. A. Montevecchi, M. Parker, I. Pratte, J.-F. Rail, G. J. Robertson, J. C. Rock, R. A. Ronconi, D. Shutler, I. J. Stenhouse, A. Takahashi, Y. Watanuki, L. J. Welch, S. I. Wilhelm, S. N. P. Wong and M. L. Mallory. 2023. Experts’ opinions on threats to Leach’s Storm-Petrels (Hydrobates leucorhous) across their global range. Avian Conservation and Ecology 18(1):11. https://doi.org/10.5751/ACE-02370-180111

Pollet, I. L., and D. Shutler. 2018. Leach’s Storm Petrel Oceanodroma leucorhoa population trends on Bon Portage Island, Canada. Seabird 31:75-83. https://doi.org/10.61350/sbj.31.75

Pratte, I., G. J. Robertson, and M. L. Mallory. 2017. Four sympatrically nesting auks show clear resource segregation in their foraging environment. Marine Ecology Progress Series 572:243-254. https://doi.org/10.3354/meps12144

R Core Team. 2024. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://doi.org/10.32614/R.manuals

Rail, J.-F. 2021. Eighteenth census of seabirds breeding in the sanctuaries of the North Shore of the Gulf of St. Lawrence, 2015. Canadian Field-Naturalist 135(3):221-233. https://doi.org/10.22621/cfn.v135i3.2675

Rail, J.-F., and G. Chapdelaine. 2002. Quinzième inventaire des oiseaux marins dans les refuges de la Côte-Nord: techniques et résultats détaillés. Technical Report Series No. 392. Environment Canada, Sainte-Foy, Québec, Canada.

Rennie, I. R., D. J. Green, E. A. Krebs, and A. Harfenist. 2020. High apparent survival of Leach’s Storm Petrels Oceanodroma leucorhoa in British Columbia. Marine Ornithology 48:133-140. https://doi.org/10.5038/2074-1235.48.1.1357

Regular, P., W. Montevecchi, A. Hedd, G. Robertson, and S. Wilhelm. 2013. Canadian fishery closures provide a large-scale test of the impact of gillnet bycatch on seabird populations. Biology Letters 9(4):20130088. https://doi.org/10.1098/rsbl.2013.0088

Robertson, G. J., and R. D. Elliot. 2002a. Changes in seabird populations breeding on Small Island, Wadham Islands, Newfoundland. Canadian Wildlife Service Technical Report Series 381. Canadian Wildlife Service, Atlantic Region, Mount Pearl, Newfoundland and Labrador, Canada.

Robertson, G. J., and R. D. Elliot. 2002b. Population size and trends of seabirds breeding in the Gannet Islands, Labrador. Canadian Wildlife Service Technical Report Series 393. Canadian Wildlife Service, Atlantic Region, Mount Pearl, Newfoundland and Labrador, Canada.

Robertson, G. J., R. D. Elliot, and K. G. Chaulk. 2002b. Breeding seabird populations in Groswater Bay, Labrador, 1978 and 2002. Canadian Wildlife Service Technical Report Series 394. Canadian Wildlife Service, Atlantic Region, Mount Pearl, Newfoundland and Labrador, Canada.

Robertson, G. J., J. Russell, R. Bryant, D. A. Fifield, and I. J. Stenhouse. 2006. Size and trends of Leach’s Storm-Petrel Oceanodroma leucorhoa breeding populations in Newfoundland. Atlantic Seabirds 8:41-50.

Robertson, G. J., J. Russell, and D. A. Fifield. 2002a. Breeding population estimates for three Leach’s Storm-Petrel colonies in southeastern Newfoundland, 2001. Canadian Wildlife Service Technical Report Series 380. Canadian Wildlife Service, Atlantic Region, Mount Pearl, Newfoundland and Labrador, Canada.

Ronconi, R. A., K. A. Allard, and P. D. Taylor. 2015. Bird interactions with offshore oil and gas platforms: review of impacts and monitoring techniques. Journal of Environmental Management 147:34-45. https://doi.org/10.1016/j.jenvman.2014.07.031

Sandvik, H., K. E. Erikstad, and B. E. Sæther. 2012. Climate affects seabird population dynamics both via reproduction and adult survival. Marine Ecology Progress Series 454:273-284. https://doi.org/10.3354/meps09558

Sauer, J. R., and W. A. Link. 2011. Analysis of the North American breeding bird survey using hierarchical models. Auk 128(1):87-98. https://doi.org/10.1525/auk.2010.09220

Smith, A. C., and B. P. Edwards. 2021. North American Breeding Bird Survey status and trend estimates to inform a wide range of conservation needs, using a flexible Bayesian hierarchical generalized additive model. Condor 123(1):duaa065. https://doi.org/10.1093/ornithapp/duaa065

Wilhelm, S., A. Hedd, G. J. Robertson, J. Mailhiot, P. M. Regular, P. C. Ryan, and R. D. Elliot. 2020. The world’s largest breeding colony of Leach’s Storm-Petrel Hydrobates leucorhous has declined. Bird Conservation International 30:40-57. https://doi.org/10.1017/S0959270919000248

Wilhelm, S. I., J. Mailhiot, J. Arany, J. W. Chardine, G. J. Robertson, and P. C. Ryan. 2015. Update and trends of three important seabird populations in the western North Atlantic using a geographic information system approach. Marine Ornithology 43:211-222. https://doi.org/10.5038/2074-1235.43.2.1133

Wood, S. N. 2017. Generalized additive models: an introduction with R. Second edition. Chapman and Hall/CRC, Boca Raton, Florida, USA.