Occupancy trends of overwintering coastal waterbird communities reveal guild‐specific patterns of redistribution and shifting reliance on existing protected areas

Climate change and anthropogenic stressors are redistributing species and altering community composition globally. Protected areas (PAs) may not sufficiently protect populations of species undergoing distributional shifts, necessitating that we evaluate existing PAs and identify areas for future protection to conserve biodiversity across regional and temporal scales. Coastal waterbirds are important indicators of marine ecosystem health, representing mobile, long‐lived, higher trophic‐level consumers. Using a 20‐year citizen science dataset (1999–2019) with a before‐after control‐intervention sampling framework for habitat protection, we applied dynamic occupancy models to assess winter occupancy trends along the Pacific coast of Canada. Specifically, we sought to understand potential drivers of regional declines, spatial commonalities among guilds, and changes in habitat use before and after PA designation, as well as between PAs and non‐PAs. Occupancy trends varied regionally, with greater declines in the south compared to the north. Regional differences underlined potential range shifts, particularly for species with traits linked to temperature tolerance, movement, and high productivity foraging, as cold‐tolerant, migratory benthivores and piscivores wintered farther north relative to 20 years ago or retreated to cold‐water fjords. While 21 of 57 (36.8%) species responded positively to PA designation (before‐after), greater occupancy declines tended to occur in PAs established pre‐1999 relative to non‐PAs (control‐intervention). Since PAs are currently concentrated in the south, negative associations were most apparent for species retreating northward, but existing PAs may have a stabilizing or transitory effect on southern wintering species shifting into the region from farther south. We emphasize that conservation strategies must balance persistence of current communities with preserving the climate‐adapted biodiversity of tomorrow by accounting for community‐level effects of species moving into and out of existing PAs. Incorporating range shifts into PA planning by predicting distributional changes will allow conservation practitioners to identify priority habitats, such as cold‐water refugia, for persistent wildlife communities.

Specifically, we sought to understand potential drivers of regional declines, spatial commonalities among guilds, and changes in habitat use before and after PA designation, as well as between PAs and non-PAs.Occupancy trends varied regionally, with greater declines in the south compared to the north.Regional differences underlined potential range shifts, particularly for species with traits linked to temperature tolerance, movement, and high productivity foraging, as cold-tolerant, migratory benthivores and piscivores wintered farther north relative to 20 years ago or retreated to cold-water fjords.While 21 of 57 (36.8%) species responded positively to PA designation (before-after), greater occupancy declines tended to occur in PAs established pre-1999 relative to non-PAs (control-intervention).Since PAs are currently concentrated in the south, negative associations were most apparent for species retreating northward, but existing PAs may have a stabilizing or transitory effect on southern wintering species shifting into the region from farther south.We emphasize that conservation strategies must balance persistence of current communities with preserving the climate-adapted biodiversity of tomorrow by accounting for communitylevel effects of species moving into and out of existing PAs.Incorporating range shifts into PA planning by predicting distributional changes will allow conservation

| INTRODUC TI ON
Protected areas (PAs) remain a central component of area-based conservation, aimed at preserving biodiversity, ecological processes, and biological resources critical to wildlife populations (Dudley, 2008;Gray et al., 2016).Due to climate change and other anthropogenic stressors, nearly all ecosystems will shift to an altered climatic state, driving the redistribution of biodiversity and individual species (Chen et al., 2011;Elsen et al., 2020;Lenoir & Svenning, 2015;Mason et al., 2015), particularly marine species (Lenoir et al., 2020).As most PAs have fixed boundaries, their presence alone may not sufficiently contribute to the persistence of species undergoing distributional shifts (Araújo et al., 2011;Thomas & Gillingham, 2015), resulting in pronounced changes to community structure and ecosystem functioning (Antão et al., 2022;Blowes et al., 2019;Pecl et al., 2017).
Conservation efforts must therefore combine protecting habitats that promote resistance to climate change effects (e.g., refugia) with identifying and preserving future biodiversity hotspots, thereby protecting species retreating from and expanding into current PAs (Hylander et al., 2022;Lawler et al., 2020;Lehikoinen et al., 2019).
Effectively linking conservation initiatives across regional and temporal scales requires evaluating reliance on current PAs and identifying habitats for future protection within the context of shifting species ranges or community composition (D'Aloia et al., 2019;Scridel et al., 2021).Marine birds, waterfowl, and waders wintering along coastal habitats (hereafter 'coastal waterbirds') can serve as important case studies for area-based conservation under climate change.Specifically, coastal waterbirds are considered indicators of marine ecosystem health, including changes in water mass distribution, ocean productivity, and prey resources, as they are higher trophic-level consumers, relatively long-lived, and highly mobile (Sydeman et al., 2001;Veit et al., 1996).During winter in the northern hemisphere, coastal habitats and estuaries are ecologically important, providing the necessary resources to withstand winter conditions and prepare for the energetically costly activities of migration and reproduction (Pelletier et al., 2020).These coastal habitats are also at risk from pervasive anthropogenic stressors, including but not limited to over-fishing, pollution, hunting, habitat degradation, and recreational disturbances (Harley et al., 2006;Islam & Tanaka, 2004;Madsen, 1995;Mahoney & Bishop, 2017).Coastal waterbird conservation therefore requires measures that address the interacting stressors of climate and habitat use change.
Globally, existing PAs have the potential to conserve coastal waterbird communities, particularly with proper management interventions (Wauchope et al., 2022), yet are also at risk of diminishing effectiveness as species shift poleward (Gaget et al., 2021;Johnston et al., 2013;Pavón-Jordán et al., 2020).In these studies, species-specific abundance or species turnover provided insight into how communities change over time within and beyond existing PAs.
However, to properly evaluate the influence of habitat protection, it is necessary, but relatively uncommon, to assess changes in wildlife communities before and after PA implementation, as well as compare among spatially replicated control (unprotected) and intervention (protected) sites (Green, 1979;Stewart-Oaten & Bence, 2001).
An integrated before-after control-intervention (BACI) design is a powerful approach to separate the effects of PAs on wildlife populations from unrelated spatial and temporal changes, including those driven by range shifts (Christie et al., 2021;Wauchope et al., 2022).
Additionally, explicitly modeling detection probability allows for more accurate estimates of changes in abundance or occupancy, particularly when the probability of detection varies among species or over time (Guélat & Kéry, 2018;Kéry & Schmidt, 2008).It is also necessary to fully partition the parameters underlying metapopulation dynamics, such as colonization and extinction (MacKenzie et al., 2003).A framework that combines a BACI approach with explicitly modeling detection probability is therefore paramount to separating regional declines from range shifts in the context of areabased conservation effects.
In recent decades, the abundance of multiple overwintering coastal waterbird species has declined in the Salish Sea, a marginal sea located along the southern extent of Canada's Pacific coast (Bower, 2009), and an area of global importance for millions of waterbirds migrating along the Pacific Flyway or wintering in the region (Butler et al., 2021).However, spatial patterns in declines vary regionally and are not consistent among groups with similar ecological traits (Ethier et al., 2020), suggesting differential mechanisms are at play.In isolation, local changes in abundance are insufficient to distinguish range-wide species declines from shifting winter ranges as individuals track annual changes in climatic conditions or food resources across larger areas (Hyrenbach & Veit, 2003).Instead, regional patterns in systems where resources are spatially fragmented are best understood in a metapopulation framework that links the local processes that influence occupancy (i.e., colonization and extinction) across larger spatial scales (Hanski, 1999).It is therefore necessary to assess changes in colonization-extinction dynamics and their contribution to range-wide, spatially explicit occupancy trends to understand species reliance on existing PAs and drivers of regional declines.
practitioners to identify priority habitats, such as cold-water refugia, for persistent wildlife communities.

K E Y W O R D S
area-based conservation, before-after control-intervention, citizen science, climate change, cold-water refugia, community ecology, dynamic occupancy model, range shift Standardized, long-term monitoring programs spanning large geographic areas present unique opportunities to monitor occupancy trends and potential range shifts, including within existing PAs.Using a 20-year dataset (1999-2019) collected from the British Columbia Coastal Waterbird Survey (BCCWS), a citizen science program managed by Birds Canada, we sought to (1) resolve regional trends in coastal habitat use for the winter waterbird community; (2) assess general occupancy patterns among families and ecological traits (i.e., dietary niche, migration distance, and temperature tolerance) to better understand spatial commonalities among species groups; (3) examine the extent of community-level range shifts; and (4) understand temporal changes in overwintering coastal waterbird occupancy in protected versus unprotected habitats, as well as before and after PA designation.Bird counts included survey effort and sampled habitats within PAs designated prior to 1999, before and after PA designation during the study period, and unprotected control habitats, providing the rare opportunity to combine a BACI design with species-specific detection probability estimates.In addition, comparing patterns among groups of species with shared ecological traits provides a strong framework in which to understand current community composition and predict climate-driven spatial and temporal changes, as trait-based approaches characterize niche space and thus directly link species to the environment (Pigot et al., 2020;Schleuning et al., 2023).
We hypothesized occupancy trends would be influenced by guild-level characteristics, and that these would vary spatially (Ethier et al., 2020).Specifically, we predicted that occupancy declines would be more prevalent in higher trophic-level consumers due to greater sensitivity to changing or declining prey availability and quality (Hyrenbach & Veit, 2003;Vilchis et al., 2015).We also predicted that occupancy declines would be greater in cold-tolerant species, which may track their climate envelope and shift beyond the study region to the north or to isolated, cold-water refugia, compared to warm-tolerant species increasing throughout the study region as populations shift northward from previous core wintering habitats farther south (Hyrenbach & Veit, 2003).Additionally, we expected occupancy trends of year-round residents to remain relatively stable compared to increases for migrant species, reflecting a tendency to remain farther north.Finally, since existing PAs are largely concentrated in the Salish Sea, we expected that species exhibiting retreating winter range shifts out of the Salish Sea would experience reduced protection over time as species ranges move beyond existing PA boundaries, compared to species expanding into the Salish Sea which may benefit from increased protection coverage.

| Study area
The Pacific coast of Canada represents over 25,000 km of shoreline along the province of British Columbia (BC), consisting of a wide range of habitats including estuaries, rocky intertidal zones, mudflats, and sand beaches (Sebert & Munro, 1972).Mild winter temperatures and the location within the Pacific Flyway make these habitats globally valuable overwintering sites for a diverse and abundant bird community (Butler et al., 2021).
Our study area was bounded by the BCCWS (see Section 2.2), from 48.29° to 54.33° N latitude and −132.10° to −122.79°W longitude (Figure 1a).We categorized the study area into three general regions representing different environments and human impacts: (1) inner coast (hereafter "Salish Sea"); (2) outer coast (hereafter "Pacific coast"), and; (3) the North-Central coast (Figure 1a).The Canadian portion of the Salish Sea occurs between Vancouver Island and the mainland at the southern extent of BC, including the Strait of Georgia and the Strait of Juan de Fuca.This region is characterized by high population density and human impacts stemming from industrial, urban, and rural activities, co-occurring with the greatest diversity of overwintering coastal waterbirds among the three regions (Gaydos & Pearson, 2011;Stocks & Vandeborne, 2017) The BC coast contains 437 PAs with an average size of 266 ± 3999 km 2 (mean ± SD; range = 0.0004 to 82,515 km 2 ; Government of Canada, 2022).The BCCWS covers 104 of these PAs, 31 of which were established during the study period  to preserve critical habitat and stabilize coastal marine ecosystems.Including those established pre-1999, sampled PAs varied in their protection length from 4 to 96 years (Figure 1b).There are no consistent differences between PAs established pre-1999 and during the study period in terms of habitat type or management activity.However, while some of the oldest PAs are migratory bird sanctuaries, many tended to be designated for recreational activities and did not directly prioritize waterbirds, compared to

| Field surveys
The BCCWS is a volunteer-based initiative started in 1999 to monitor population trends of coastal species (Birds Canada, 2021a).
The survey consists of 368 routes, the majority of which are within the Salish Sea, while 15% occur along the Pacific and North-Central coast.Each route covers a similar area, extending from shoreline to 1 km offshore, and surveys are conducted on the second weekend of each month following a standardized protocol (Birds Canada, 2021a).Counts occur within 2 h of high tide, during fair weather days with wind speeds between zero to three on the Beaufort scale, and last between 30 min and 5 h (average = 2 h) depending on the observer and site.All birds are counted, excluding non-native species, vagrants, and hybrids.
Flyovers are not included as they are not considered to be using the surveyed habitat.

| Data processing
For our analysis, we used surveys conducted during peak winter (December-February) from 1999 to 2019.Only routes that were visited during all 3 months for at least 3 years within the 20-year study period were retained, resulting in a final sample of 310 routes (Figure 1a).Survey routes were visited for an average of 13.1 ± 6.3 years (mean ± SD).With respect to protection status, sites protected pre-1999 (n = 73), sites protected during the study period (n = 31), and unprotected sites (n = 206) were visited for an F I G U R E 1 Distribution of coastal waterbird survey routes in British Columbia, Canada, divided into three general regions: Salish Sea (circle), Pacific coast (square), and North-Central coast (diamond), along with their associated habitat protection status (a)."Prior" indicates routes within PAs established pre-1999, "protected" are routes within PAs established during the study period (1999-2019), and "unprotected" are routes that have never been within PA boundaries.Row (b) shows the number of years since the habitat was protected.Both include insets of the Salish Sea region, with the highest density of survey routes.Map lines delineate study areas and do not necessarily depict accepted national boundaries.average of 13.4 ± 6.3, 11.4 ± 6.3, and 13.2 ± 6.4 years, respectively.Survey data for each visit were converted from counts to occurrence data by zero-filling.All species recorded on at least 1% of visits were retained for the analysis resulting in a final list of 57 species (Appendix S1; Table S1).
Each species was also categorized into ecologically similar groups based on phylogenetic relatedness (family), dietary niche, migration distance, and temperature tolerance.Species were classified as benthivores, herbivores, omnivores, and piscivores to reflect their primary winter diet following Bower (2009) and Ethier et al. (2020).Migration distance included resident, short-distance, and long-distance migrants based on the distance between the centroids of the breeding and wintering range of each species (Davidson et al., 2015;Ethier et al., 2020).Temperature tolerance was derived from winter species temperature indices (STI), or the long-term average temperatures experienced by a species over its wintering range (December-February), representing a seasonal climate envelope (Devictor et al., 2008).Temperature tolerance differs from migration distance in that it reflects global range size as well as an affinity for higher latitudes.STI values were accessed from Lehikoinen et al. (2021), where temperatures were averaged from 1950 to 2000 using 30-arc resolution WorldClim data (Hijmans et al., 2005) and BirdLife International range estimates (Santangeli & Lehikoinen, 2017).To evaluate the average response of species with similar winter climate envelopes, we categorized species as cold (≤0°C; n = 15), cool (>0°C and ≤5°C; n = 17), and warm-tolerant (>5°C; n = 25).These categories were chosen because they divided species along an approximately evenly balanced temperature gradient and 5°C was the median winter water surface temperature for the study region (Appendix S1; Figure S1).See Appendix S1; Table S1 for ecological traits of each species.

| S TATIS TIC AL ANALYS IS
We used a dynamic occupancy approach to address changes in the probability of occupancy among years in response to habitat protection (MacKenzie et al., 2003).Dynamic occupancy models relax the assumption of a closed population and allow metapopulation dynamics such as colonization and extinction to be estimated among years (MacKenzie et al., 2009;Weir et al., 2009).The detection of species during a survey is imperfect and can be influenced by multiple factors such as weather and habitat structure (Guélat & Kéry, 2018;Moilanen, 2002).Dynamic occupancy models explicitly estimate differences in detection probability using variation in observations among visits within the same season when populations are assumed to be more stable (i.e., limited partial emigration; MacKenzie et al., 2003).By accounting for imperfect detection, the latent, unobserved variable "true occupancy" can be estimated for each site and year, allowing one to estimate trends over time, as well as metapopulation characteristics such as the proportion of occupied sites even if each site is not visited annually (Weir et al., 2009).
We acknowledge that individual coastal waterbirds are likely somewhat nomadic during the winter, challenging the assumption of a closed population.However, occupancy is less sensitive to fluctuations in abundance driven by individual movements, and we assume that species use the same general area throughout the survey period to access spatially aggregated resources.Therefore, we believe our survey design sufficiently approximates a closed population for the average coastal waterbird species.

| Model covariates
Dynamic occupancy models include sub-models to estimate detection probability (ρ), initial occupancy (ϕ), colonization (γ), and extinction rate (ε).For detection probability, we tested the influence of ordinal date, time of day, survey duration (minutes), and visibility.Time of day was the survey start time rounded to the nearest half hour and visibility was collapsed into two categories: (1) <1 km and (2) ≥1 km.Observer ID was not included, given that with over 400 different observers across the 20-year survey, a model with observer ID would not converge.Instead, we tested a polynomial spatial autocorrelation term to control for unmeasured observer effects shared among adjacent sites that may influence detection probability and defined as latitude + longitude + latitude2 + longitude2 + latitude × longitude (Schuster & Arcese, 2013;Tozer, 2016).
The probability of initial occupancy is influenced by route-level covariates or habitat variables that are consistent across years.
We considered minimum sea surface temperature, average salinity, net primary productivity (NPP), bathymetric slope, wave exposure, coastal habitat type, and the initial protection status in 1999 as route-level covariates.Temperature, salinity, and NPP were extracted from the Bio-ORACLE database, representing long-term averages between 2000 and 2014 of water surface measurements interpolated at a 9.2 km resolution (Assis et al., 2018;Tyberghein et al., 2012).Slope was derived from the MARSPEC geophysical database at a 1-km resolution (Sbrocco & Barber, 2013).Sea surface temperature is the minimum estimated within the 3-month winter period at each location.Average salinity is the mean across the 3 months in practical salinity units (PSU; g/L), and NPP (g/m 3 /day) is the annual minimum, occurring during winter.Slope is the angle of the ocean floor from the shoreline, influencing both near-shore water depth and tide duration.Exposure and coastal habitat type were extracted from the BC physical shore-zone mapping database (Howes et al., 1994).Exposure is derived from topographical measurements categorizing coastal wave exposure, which we collapsed into three representative levels: protected, semi-protected, and exposed.We quantified habitat type into six categories approximating the physical structure and substrate of the shoreline: estuary/marsh, sand and mudflat, sand and gravel beach, gravel beach, rock cliff, and man-made structure.Finally, initial protection status was derived from the Canadian Protected and Conserved Areas Database Canada, 2022), indicating whether a survey route overlapped with a PA in 1999 (1) or not (0).We only considered a route to be within a PA if there was at least 5% overlap.
For colonization and extinction rates, we included a habitat protection covariate, which indicates whether a survey route falls within a PA for a given year, weighted by the proportion of overlap between the survey and PA (0%-100%) to avoid treating survey routes on the periphery of PAs as equal to one that is fully protected.In this instance, habitat protection is a year-route variable such that there is a before and after protection status for each route depending on what year the PA was established (Figure 1b), allowing the resulting change in colonization and extinction rates to be estimated.This differs from the habitat protection status for initial occupancy, which is a static variable that differs between control and intervention sites.For all parameters, the same spatial autocorrelation term described for detection probability was tested.All covariates were standardized prior to model evaluation.

| Model fit and selection
We performed a goodness-of-fit test on the global, fully parameterized model using the MacKenzie-Bailey chi-squared test (MacKenzie & Bailey, 2004), which involved generating χ 2 estimates from 1000 bootstrap iterations and comparing it to the original χ 2 statistic.
A non-significant result (p > .05)indicates a "good fit" or that the predictors adequately explain variation in each of the estimated parameters.From this, we also derived a c-hat statistic, where values less than or greater than 1 indicate under-and over-dispersion, respectively.C-hat values were used to correct model comparisons of all nested sub-models using Quasi Akaike Information Criterion (QAICc) for small sample sizes, which introduces an additional penalty for over-or under-dispersion (Lebreton et al., 1992;MacKenzie & Bailey, 2004).
We used a stepwise model selection approach where we subsequently estimated the top model for detection probability and initial occupancy, followed by colonization and extinction rate simultaneously (Tozer, 2016).All nested sub-models for each of the parameters were fit and compared using QAICc.The top model from the previous parameter was used in model selection for the subsequent parameter.A top model for each parameter and each species was determined by considering all models within 2 QAICc of the top model (Burnham & Anderson, 2002).Any models with uninformative, nested predictors were removed and the model with the highest model likelihood was selected for further model evaluation steps.If multiple unnested models had likelihoods >0.80, then the most parsimonious model was selected.This approach reduces overfitting while also selecting models that describe the most variance.See Table 1 for the global model structure.

| Model inference and post-hoc analysis
Once the top model was selected for each species, we employed Markov chain Monte Carlo (MCMC) sampling in the STAN environment to produce posterior distributions for each parameter.
We fit the top model for each species using four chains of 3000 iterations, each with a 1000 iteration burn-in and thinning rate of 1.The mean and 95% credible interval (CrI) were used to evaluate effect size support.From the posterior distributions of the MCMC model output, we estimated latent occupancy for each route-year combination and then calculated 20-year occupancy trends using a geometric mean approach (Smith et al., 2014(Smith et al., , 2020)).We averaged estimated occupancy within the first and last 5 years and then the difference was converted to an annual proportional change over a 20-year time period using the following equation: In doing so, we smoothed over any large year-to-year fluctuations that may overestimate trends.The annual proportional change in occupancy for each route was then averaged across routes to generate a regional occupancy trend for each species, all species combined (community-wide), and for each guild.We also categorized routes by protection status to test for spatial differences in trends and median occupancy over the 20-year study between habitats within PAs prior to 1999, PAs established during the study period, and unprotected habitats.All median occupancy and trend comparisons are presented with a 95% CrI derived from the posterior distributions.
Finally, we used the top model to predict annual occupancy of each species across the entire study region, and then overlayed the proportional occupancy change between 1999 and 2019 to identify hotspots of increasing or declining occupancy.We estimated core occupancy for each species per year, defined as the range encompassing map cells with greater than 50% of occupancy probabilities (Rushing et al., 2019(Rushing et al., , 2020)).We then calculated the weighted centroid for 1999 and 2019 for each species and estimated the Haversine distance and bearing between the TA B L E 1 Global model structure prior to model selection, which was conducted sequentially.All possible nested models of detection probability (ρ) were fit, holding the other parameters constant as intercept-only models, and then compared using QAICc to determine the top detection probability model for subsequent analyses.The top detection model was then held constant for subsequent model selection steps until a full top model was determined.The quadratic terms for date and time were included to test for nonlinear differences over the winter season or day, respectively.two points to estimate the extent and direction of shifts in core occupancy.As such, we used the relative highest likelihood winter habitats to be comparable among species and avoid out-weighted influence of ephemeral sites.

Parameter
Statistical analyses were conducted in R version 4.2.1 (R Core Team, 2022).Environmental covariates were collated and extracted using the package "sdmpredictors" (Bosch, 2020).

| Model fit
Goodness-of-fit tests were non-significant for all 57 species with an average p-value of .30(range = 0.06 to 0.59), indicating that the tested predictors adequately explained model variance.As well, chat values were close to 1 (median = 1.32, range = 0.77-2.79)demonstrating minimal issues with under-or overdispersed residuals.
We therefore applied the stepwise model selection process to all 57 species.

| Detection probability and initial occupancy
For the average species, detection probability was influenced by four predictors (range = 1-6).Survey duration and visibility were included in most top models, where longer surveys and greater visibility were associated with higher detection probabilities (Table 2).Date and time were also important for over half of all species, but their effects were variable with no clear pattern (Table 2).Support for the quadratic terms of date and time suggests detection cannot be fully accounted for using linear terms in this system.Finally, spatial autocorrelation was identified for all but two species, suggesting systematic spatial biases in detection probability that may reflect shared observer effects among adjacent routes (Table 2).
The influence of habitat and environmental predictors on initial occupancy varied widely, with no single variable included in the top model for more than 24 species (42.1%;Table 2).On average, three predictors (range = 0-7) occurred in the top models.Minimum sea surface temperature, NPP, salinity, bathymetric slope, and exposure were all selected relatively frequently (Table 2).Shore type was infrequently selected but there was a demonstrable preference for mud flats or estuaries.Temperature, NPP, salinity, and exposure generally had positive influences, while slope had a negative effect (Table 2).Routes within PAs established prior to the study period (1999) also tended to have higher initial occupancy (Table 2).

| Species-specific occupancy trends
The median proportion of occupied routes over the 20-year study period ranged from a low of 0.03 for snow goose (Anser caerulescens) to a high of 0.90 for bufflehead (Bucephala albeola) and glaucouswinged gull (Larus glaucescens; Appendix S1; Table S2).Average occupancy across species was 0.37 ± 0.03 (mean ± SE) but tended to be slightly higher along the outer Pacific and North-Central coast (0.46; Appendix S1; Table S2).Uncertainty around predicted occupancy was relatively low and consistent across all coastal regions, suggesting limited influence of geographic bias in survey sites (Appendix S1; Figure S2).
On average, 21 species exhibited significant declines in annual route occupancy, including species like eared grebe (Podiceps nigricollis) which decreased by up to 2.81% per year, while seven species experienced annual increases of up to 2.33% per year (e.g., brant Branta bernicla; Appendix S1; Table S3).These values were largely TA B L E 2 Influential predictors of detection probability and initial occupancy for coastal waterbirds in British Columbia, Canada.The total number refers to the number of species, out of 57, where the specific variable occurred within the top model, regardless of its effect size (i.e., model selection).The (+) and (−) rows refer to the direction of influence and the values are the number of species where the 95% credible interval of the variable estimate did not include zero (i.e., strong effect).For initial occupancy, a positive influence of exposure means greater wave exposure (i.e., semi-protected or exposed), while a positive influence of shore type represents a preference for sandy estuaries or mudflats versus rocky or urbanized shores.driven by trends in the Salish Sea, with 22 species experiencing significant declines compared to seven species increasing.Trends along the North-Central and Pacific coast were more positive, with five and nine species increasing compared to two and six species declining, respectively (Appendix S1; Table S3).Differences among regions could also be seen in the average route occupancy trends across species, with −0.52 ± 0.16 in the Salish Sea compared to 0.23 ± 0.32 and −0.09 ± 0.26 in the North-Central and Pacific coast regions, respectively.

| Guild-specific occupancy trends
The greatest community-wide declines in occupancy occurred in the Salish Sea and along the exposed, outer coastal waters of the Pacific coast, while positive occupancy rates concentrated around the outer coast of Haida Gwaii and fjords of the North-Central coast (Figure 2).Separating trends into dietary guilds, only herbivores were stable or increasing in the Salish Sea compared to omnivores, which were stable with spatial clusters of decline, and both benthivores and piscivores which decreased steeply.Benthivores and piscivores declined precipitously along the entire coast, except for increases in the sheltered fjords of the North-Central coast.Omnivores and herbivores were largely responsible for the increasing occupancy along Haida Gwaii, and in general, both increased or remained stable along much of the North-Central coast (Figure 2).For migration distance, resident species were relatively stable throughout the Salish Sea, with increased occupancy along the central coast but decreasing occupancy in the north.Short-distance migrants declined along most of the coast, except for small pockets of stability in the North-Central fjords, while long-distance migrants exhibited substantial gains in the north and decreased in the south (Figure 2).Finally, coldtolerant species (i.e., ≤0°C) declined in both the north and the south, with slight gains in the central coast and northern fjords.In contrast, cool species (0-5°C) increased widely along the North-Central coast and warm-tolerant species (>5°C) were relatively stable throughout the study region with patches of increased occupancy in the north (Figure 2).

| Winter range shifts
Core winter occupancy across species shifted poleward by an average of 114.1 ± 110.9 km across the 20-year study period or ~5.7 km/ year.Five of ten families (50%) shifted north or northwest by ~60 to 280 km, with the remainder shifting southwest to the outer Pacific Coast or southeast to the Salish Sea (Figure 3).In terms of dietary niches, benthivores were the only species consistently shifting north, with herbivorous species split between northern and southern shifts, and core occupancy moving south for both piscivores and omnivores (Figure 3).Long-distance migrants also shifted to the northwest, compared to pronounced south-eastern shifts by short-distance migrants and residents.Finally, cool species tended to move northwest, while cold-and warm-tolerant species experienced apparent shifts to the south (Figure 3).

| Reliance on protected areas
After controlling for the proportion of overlap between protected and surveyed areas, before-after comparisons suggested 21 of 57 species (36.8%) responded positively to protection established between 1999 and 2019, with seven species exhibiting significantly greater colonization rates and 14 species with significantly lower extinction rates (Figure 4a).However, six species appeared to have a negative response to protection (i.e., increasing extinction probability and/or decreasing colonization probability; Figure 4a).Across ecological traits, protection status had minimal effects.On average, herbivores experienced greater colonization and lower extinction rates after PAs were established (Figure 4b).No consistent association was apparent across migratory guilds (Figure 4c), while warm-tolerant species had higher colonization rates in response to protection (Figure 4d).
Control-intervention comparisons of occupancy trends from 1999 to 2019 indicated that the greatest annual declines tended to occur within PAs established prior to 1999, rather than PAs established during the study period or unprotected areas (Figure 5a).

Across dietary guilds, benthivores and piscivores declined within
PAs, while omnivores decreased regardless of protection status, and herbivores increased in unprotected areas (Figure 5b).Both shortand long-distance migrants declined within but not beyond PAs on average (Figure 5c).Finally, cold-tolerant and cool species declined in established PAs only, compared to warm-tolerant species which had no consistent trends among protection categories (Figure 5d).In fact, at the species level, only three species experienced increased occupancy within PAs established pre-1999 over the 20-year period compared to 27 significant declines and an overall annual rate of −0.70 ± 0.17.In contrast, PAs established during the study period and unprotected habitats included 5 and 11 species with increasing occupancy relative to 8 and 11 declining species, respectively (Appendix S1; Table S4).Similar patterns can be seen in family-wide occupancy trends (Appendix S1; Figure S3).

| DISCUSS ION
Winter coastal waterbird occupancy declined overall since 1999, with pronounced regional variation.More species declined in the Salish Sea compared to the Pacific coast, while only the North-Central coast exhibited an average increasing trend across species.
These patterns generally support previous work comparing waterbird abundance in the Salish Sea and Pacific coast regions (Ethier et al., 2020).While trend direction aligned with Ethier et al. (2020), we often found support for more pronounced changes, both positive and negative.This difference could be because we assessed routes from a broader geographic area, accounted for detection probability, or informed nearby trends using spatial clustering.Occupancy and abundance can also be influenced by environmental characteristics operating at different temporal or spatial scales (Orrock et al., 2000).
Thus, while we would not expect perfect concordance between approaches, the same general pattern in results provides confidence in the respective findings.

F I G U R E 2
Predicted change in the probability of occupancy of coastal waterbirds from 1999 to 2019 in British Columbia, averaged across all species (top left), as well as for each of the ecological traits: dietary niche, migration distance, and temperature tolerance.Temperature tolerance is based on the species temperature index derived from the average temperature within the winter range of each species.
Cold-tolerant species are those associated with temperatures below 0°C, "cool" from 0 to 5°C, and "warm" above 5°C.See Figure S1 for underlying environmental variable distributions.Map lines delineate study areas and do not necessarily depict accepted national boundaries.
Occupancy trends supported our hypothesis that declines would be more prevalent in higher trophic-level marine consum- Figure S1).They may therefore harbor preferred habitat and foraging conditions for cold-tolerant, higher trophic consumers during the winter months, particularly as admixing brackish conditions can promote prey availability (Arimitsu et al., 2016).These results support the potential for cold-water refugia maintained by complex topography along the BC coast, which may contribute to the persistence of certain species in the region (Bonett-Lebrun et al., 2022;Morelli et al., 2020).
Occupancy declines across dietary guilds in the Salish Sea, except for herbivores and some omnivores, suggest that avifauna loss in this region may be indicative of marked changes in ecosystem health within a highly industrialized region subjected to intense anthropogenic pressure.Shifts from "high-productivity" avifaunal assemblages, characterized by diving seabirds inhabiting regions with dense prey aggregations, to "low-productivity" assemblages, including surface feeding and pursuit taxa that subsist on more diverse food resources, are noted in other systems experiencing changes in ocean productivity (e.g., California coastal systems; Hyrenbach & Veit, 2003).Nutrient pollution from rivers, nonpoint source runoff, and wastewater sources are also recognized as primary threats to biodiversity in the Salish Sea (Khangaonkar et al., 2019).Given that the Salish Sea is a center of human population density, these pressures are much greater relative to the Pacific and North-Central coast.
Regional occupancy trends also underline the potential for winter range shifts.While changes in core occupancy suggest the average species winters farther north relative to 20 years ago, there was considerable variation among ecological traits.
Apparent northern and southern core occupancy shifts for benthivores and piscivores, respectively, highlight a retreat of higher trophic-level consumers to fjord habitats along the North-Central coast.In contrast, a mixed, but predominantly northward shift for herbivores indicates winter conditions may be improving for herbivorous species along the entire coast, driving expansion northward as well as into the Salish Sea from farther south.Omnivores appeared most variable, exhibiting the weakest shift in core occupancy and spatially clustered regions of stability along the coast, suggesting diet breadth may promote persistence under changing conditions (Buckley & Kingsolver, 2012).
Cold-tolerant and migrant species had higher rates of occupancy

| Reliance on protected areas
Evaluating the relative importance of colonization-extinction processes within a BACI framework is necessary to properly assess reliance on PAs (Wauchope et al., 2022).Our before-after results indicate that PAs designated in the past 20 years are associated with greater colonization and reduced extinction rates for several species, suggesting their collective value in conserving coastal bird biodiversity.These processes indicate that the winter bird community is experiencing less structural turnover after habitat protection relative to before, both attracting individuals (colonization) and contributing to persistence (reduced extinction).At the species level, reduced extinction was the most common response to habitat protection, highlighting that PA designation criteria in the past 20 years are identifying habitats with greater refugial capacity to promote resistance to climate change and anthropogenic effects (Rojas et al., 2022).While variation in response to F I G U R E 4 Before-after effect sizes for the influence of habitat protection on colonization (yellow) and extinction rates (blue) along the coast of British Columbia between 1999 and 2019.Species are associated with a colonization and extinction effect size if habitat protection was selected in the top model for those parameters.All species-specific rates significantly different from zero (dashed orange line) are denoted with an asterisk (a), similar to averaged rates across ecological traits (b-d).Cold-tolerant species are those associated with temperatures below 0°C, "cool" from 0 to 5°C, and "warm" above 5°C.Points represent standardized effect sizes, while error bars are the 95% CrI.Gray bands are included to aid visualization.
protection was highly variable within ecological traits, positive effects on colonization were most apparent for guilds that are increasing or stable along the south coast, such as herbivores and warm-tolerant species, where PA coverage is the greatest.This indicates the potential for PAs as "stepping-stone" sites, providing suitable winter habitat for species shifting into the Salish Sea from the south (Morelli et al., 2020).
Control-intervention comparisons revealed that occupancy trends in PAs designated prior to 1999 were more negative in comparison with those protected recently (1999-2019) and unprotected sites.This was particularly true for guilds with the greatest range shifts, such as migrants, cool species, and higher trophic-level consumers.A contributing factor is likely the spatial and temporal structure of habitat protection in the region.Most areas protected pre-1999 occur within the Salish Sea, a region experiencing overwhelmingly negative changes in occupancy.Given this spatial bias, we would expect species moving into the Salish Sea from farther south to drive positive trends regardless of when PAs were established.On average, herbivores, year-round resident, and warm-tolerant species, which have experienced a net southern core range shift, did not exhibit decreasing occupancies in PAs established pre-1999, lending some support to this notion.These results should be interpreted with caution, as conservation value cannot be evaluated solely on winter distributions.For example, we did not address breeding or migratory stopover use, which may reveal different occupancy patterns.However, regionspecific trends and apparent range shifts can pose a problem for existing PA networks that were designed for a previous ecosystem state (Elsen et al., 2020).Initial occupancy was positively associated with PAs established pre-1999 for 15.8% (n = 9) of species despite subsequent declines, further supporting a change in occupancy patterns and winter community structure over the 20-year study period.Therefore, the distribution of historical PAs along BC's southern coast may no longer fully capture current or future winter waterbird distributions.

| CON CLUS ION AND CONS ERVATION IMPLIC ATIONS
Using citizen-science data, we document widespread declines in winter waterbird occupancy along the Pacific coast of Canada and fundamental changes in community structure.Declines were most pronounced in the Salish Sea, an area of global importance to waterbirds subjected to often intense anthropogenic stressors (Ban & Alder, 2008;Robb, 2014).We demonstrated that recently established PAs (1999-2019) contribute to species occupancy and F I G U R E 5 Control-intervention differences in occupancy trends between PAs established prior to 1999 (purple), habitats protected during the study period (blue), and unprotected habitats (red) across all species (a), dietary guilds (b), migration distance (c), and temperature tolerance (d).Cold-tolerant species are those associated with temperatures below 0°C, "cool" from 0 to 5°C, and "warm" above 5°C.Combinations where the mean is significantly different from zero (dashed gray line) are marked with an asterisk (α = .05).
community stability, despite potential range shifts and the gradual decline in occupancy observed in PAs established pre-1999.We also identified largely unprotected regions, like fjords, which may act as refugia to aid the persistence of waterbird communities undergoing the greatest change, as well as the importance of PAs in the Salish Sea as potential stepping-stone refugia for species moving into the region from farther south.
. The Pacific region extends along the western coast of Vancouver Island, as well as Johnstone Strait and Queen Charlotte Strait at the northeast of the island.Finally, the North-Central region extends from the northern tip of Vancouver Island to Haida Gwaii inclusive.The Pacific and North-Central coast have substantially lower human populations and consequently are considerably less impacted by anthropogenic activities; however, the west coast of Vancouver Island is a popular tourist destination, and the North-Central coast contains Canada's third largest port in Prince Rupert.
more modern PAs that tend to represent land trusts or Indigenousgovernment collaborations centered around wildlife conservation and cultural preservation.Additionally, older PAs occur in the south, with all sampled PAs established >50 years ago occurring in the Salish Sea, compared to the Pacific and North-Central coast where PAs have been established primarily in the past 15-25 years (Figure 1b).
ers and species restricted to colder climate envelopes.Benthivores and piscivores demonstrated substantial declines across the study area, except for North-Central coastal fjords.Fjord habitats were also the only sites where cold-tolerant and cool species exhibited stable or increasing occupancy.Fjords represent colder, less saline water relative to the exposed coast(Wang et al., 2001; Appendix S1; However, only three species have increased within older PAs compared to 27 species decreasing, suggesting that range shifts and a spatial bias are insufficient to fully explain the declines.Changes in designation criteria from a focus on recreational activities to wildlife-based decisions and collaborations with local land trusts or Indigenous peoples for ecosystem-level solutions may have contributed to the difference between older and newer PAs (Schuster et al., 2019).Future research should investigate how designation criteria, habitat type, or management activities influence colonization and extinction rates among PAs (Le Saout et al., 2013; Wauchope et al., 2022).
Canada's coastal areas currently protected, coastal PAs are integral to reaching this goal for effective biodiversity conservation.The Northern Shelf Bioregion Network Action Plan is a government and Indigenous-led strategic framework for establishing a network of PAs along the North-Central coast of BC (MPA Network BC Northern Shelf Initiative, 2023).Given our results, we recommend incorporating the needs of species retreating poleward and those expanding into the region to facilitate the persistence of existing wintering populations and the transition to more climateadapted communities, respectively (D'Aloia et al., 2019; Hylander et al., 2022).Specifically, we recommend establishing PAs, including Indigenous Protected Areas (IPAs), along BC's North-Central coast, targeting species undergoing distributional shifts and including the protection of cold-water refugia.Additionally, for the Salish Sea region, we recommend reinforcing ongoing conservation efforts that preserve high-quality coastal habitat, particularly focusing on management practices within older PAs that target rapidly declining species and expanding newer PAs into sites with stable occupancy patterns.Coordinating a network of PAs prioritizing community stability and connectivity for shifting distributions along the Pacific coast, especially in partnership with Indigenous peoples and local communities, is crucial for sustaining biodiversity and ecosystem health in an era of global change (