Marbled Murrelets prefer stratified waters close to freshwater inputs in Haida Gwaii, British Columbia, Canada

ABSTRACT The Marbled Murrelet (Brachyramphus marmoratus) is a small seabird that is currently listed as threatened in Canada. Understanding this species' marine habitat preferences plays a vital role in our ability to focus conservation planning. We used the longest-running at-sea survey dataset available in British Columbia to examine hotspot persistence and habitat use at Laskeek Bay, Haida Gwaii, BC. The Laskeek Bay Conservation Society has been conducting spring and summer surveys along fixed transect routes in open and shoreline waters from 1997 to 2018. Along with analyzing this long-term dataset, we conducted surveys to measure oceanographic variables (2018–2019) and tested whether Marbled Murrelets in the same area used prey and oceanographic information to select marine habitat in conjunction with physical habitat features. Our hotspot persistence map, defined as areas that repeatedly had counts above a 75% threshold relative to other areas during a given survey, showed that murrelets consistently preferred shoreline transects. Murrelets also preferred shallow marine areas closer to streams, above higher proportions of sandy substrate and closer proximity to abundant nesting habitat. Modeling weather and time variables contributed little additional predictive power. Nonetheless, models that included physical environmental, oceanographic, and prey variables outperformed those with only physical environmental variables. Stratified water was the oceanographic variable most strongly related to higher counts. Our study suggests that stratified waters could work with stream systems to create productive zones for foraging murrelets, and highlights the importance of murrelets having access to marine areas with the preferred physical features. LAY SUMMARY Marine habitat preference studies for the threatened Marbled Murrelet are critical in management and conservation decisions. We used at-sea surveys conducted in Haida Gwaii, BC, by the Laskeek Bay Conservation Society (1997–2018) and performed additional oceanographic variables surveys (2018–2019). Persistence hotspot mapping showed that murrelet distributions were substantially consistent across 22 years of surveys. Murrelet usage was strongly associated with being closer to streams, shallower waters, higher proportions of sandy sediment, proximity to abundant potential nesting habitat, and stratified water conditions.


INTRODUCTION
Many seabird populations around the globe are declining, with marine threats, such as overfishing, bycatch, and warming waters, playing a major role (Birdlife International 2018, Dias et al. 2019). Marine habitat studies that gather baseline information are vital for creating effective management plans. Habitat preference studies aim to describe the behavioral responses that individuals use to select habitat that influence their survival and fitness (Hutto 1985, Block andBrennan 1993). One approach to identifying marine habitat preferences is locating high and low use areas through hotspot mapping of survey counts, which depicts areas that have counts above a threshold relative to other areas within a given study site (Veech 2000, Sussman et al. 2019. Because seabird surveys typically exhibit high variability (Piatt et al. 2007), hotspot maps based on long-term data provide more useful descriptions of probable relative usage at a given location than single surveys, and are robust to temporal variability (Sussman et al. 2019). Once such baseline insight of spatial use is established, understanding the processes behind these patterns can be achieved by quantifying the patterns' relationships with underlying environmental factors (Block and Brennan 1993).
Marbled Murrelets (Brachyramphus marmoratus, hereafter "murrelets") are pursuit diving seabirds that nest widely dispersed in old-growth forests. The harvest of oldgrowth forests has resulted in murrelets having a "threatened" status in Canada (Environment Canada 2014) and in the USA (Lynch et al. 2019). However, changes to marine habitat and prey availability also affect murrelet abundance in a given area, potentially influencing fluctuations in their population (Yen et al. 2004, Bertram et al. 2015. Murrelets' use of marine foraging habitat during the breeding seasons, and their relationships to physical marine characteristics, may vary among geographic regions (Yen et al. 2004, Haynes et al. 2011, Raphael et al. 2015. Marine habitat preferences of this seabird have been described in the southern coastal regions of British Columbia (Yen et al. 2004, Ronconi andBurger 2008); in the waters around Washington, Oregon and Central California , Raphael et al. 2015, Lorenz et al. 2016; and southern regions of Alaska (Kuletz et al. 2008, Haynes et al. 2011, Barbaree et al. 2015. No studies have been conducted in the northern coastal islands of Haida Gwaii, British Columbia, to determine marine spatial patterns and the variables influencing the distribution in these waters. Developing such knowledge will facilitate local and regional conservation planning. Our goal was to identify fine-scale murrelet habitat use in Laskeek Bay, Haida Gwaii, British Columbia. We tested the relationships between usage consistency and coastal, bathymetric, oceanographic, prey, and nesting habitat distance variables. Specifically, we used the longest-running at-sea fixed transect dataset in British Columbia, within Laskeek Bay, Haida Gwaii (1997Gwaii ( -2018, to create a hotspot persistence map identifying locations on a scale of ~0.1 km 2 where birds have been repeatedly seen or are absent throughout the years. We then used this long-term dataset to explore how physical habitat features, weather, and time of day were associated with murrelet distribution. Thereafter, we tested variables collected during surveys in 2018 and 2019 to investigate how murrelet distributions relate to prey and oceanographic features relative to physical features in this bay. We predicted that (1) cooler sea surface temperatures (SST) and more mixed thermal waters would be the oceanographic factors murrelets favored, and (2) that incorporating the number of fish schools available would strengthen the association between murrelets and physical environmental variables.

Study Area
Seabird surveys were conducted in Laskeek Bay, situated on the east side of Louise Island (52.940525°N, 131.663917°W), in the southern portion of Haida Gwaii, British Columbia, Canada ( Figure 1). The study area encompasses a surface area of ~130 km² that includes a mixture of shallow areas and deep zones exceeding 200 m. Twenty-seven kilometers of coastline lies adjacent to the study area, with 10 streams of stream order 2 or higher (Gray 2010) that input freshwater into marine waters. During the breeding season, murrelets often hold prey in their bills for long periods of time until low light hours, and where topography is steep, use streams as flyways to carry food to their offspring on nesting platforms in old-growth trees (Ralph et al. 1995, Haynes et al. 2011.
British Columbia supports breeding murrelet populations estimated most recently in 2002 as 99,100 (72,600-125,600) breeding individuals (Environment Canada 2014, Bertram et al. 2015, comprising ~28% of the world population. Haida Gwaii supports ~16% of the British Columbia total. The species' breeding season in the province extends from late March through early September, but dates vary by region and among individual pairs (Lougheed et al. 2002, Tranquilla et al. 2005. Elsewhere, murrelets preferentially utilize shallow depths and sandy substrates , Yen et al. 2004, Ronconi 2008, which likely contain a higher concentration of forage fish, such as the Pacific sand lance (Ammodytes hexapterus) (Ostrand et al. 2005). In addition to Pacific sand lance, murrelets in Haida Gwaii eat a mix of other fish during the breeding season (Sealy 1975), including northern anchovies (Engraulis mordax), capelin (Mallotus villosus), shiner perch (Cymatogaster aggregata), Pacific sandfish (Trichodon trichodon), and smelt (Osmeridae) (Vermeer and Morgan 1997), as well as Euphausiid crustacea. Oceanographic variables may heavily influence murrelet foraging choices in these waters. Strong tidal currents and rapids occur throughout the waters of Haida Gwaii (Vermeer and Morgan 1997). Upwelling mixes vertical water columns, upturning nutrients, attracting more life, and causing prey that are too weak to swim downward to be more accessible to predators (Hunt et al. 1999). Though not as dominant a force in Haida Gwaii as in more southerly waters, upwelling can still occur (Peterson et al. 2007). Cooler temperature zones are also associated with higher productivity and murrelet presence (Becker and Beissinger 2003, Chavez et al. 2003, Ronconi 2008.

Sea-Survey Data Collection
The Laskeek Bay Conservation Society (LBCS) has been conducting annual seabird surveys during spring and summer since 1997 along fixed offshore linear and shoreline transects (Figure 2A). Biologists completed 90 surveys, mostly in May and June. Each survey consisted of 18 transects, 8 shoreline (~100-300 m offshore) and 10 offshore (~300-9000 m offshore), with a mean length of 3.8 km, ranging from 1.8 to 6.3 km. Offshore transects ran from island to island to create visual points that a boat driver could use to navigate on a straight trajectory. LBCS conducted surveys over a 4-month period (April-July) from 1997 to 2003, and over 3 months (May-July) from 2004 onwards. Surveys were only conducted in fair weather (Beaufort Sea State 3 or less) and included all 18 transects in one day unless the weather turned, in which case a set of surveys might be conducted over 2, usually consecutive, days. A Beaufort Sea State 3 is characterized by small wavelets, crests that do not break, and a light breeze (Canada 2017).
Surveys were conducted by 2 ̶ 4 participants traveling in a small aluminum skiff. Start and end times were recorded for each transect. Using a voice recorder, the primary observer identified all seabirds and dictated the number and time birds were seen on the water, while the secondary observer drove the boat. Any additional surveyors helped with timing, GPS waypoint recordings, and observations. Observations were made out to 50 m on both sides of the boat, producing a summed transect width of 100 m. Because the transect width was narrow, we assumed a 100% detectability of birds. Birds seen on or FIGURE 2. Comparison of survey routes set up by Laskeek Bay Conservation Society (1997-2018) (A) to the oceanographic variables surveys run in 2018 and 2019 (B) to how transects were segmented into ~1-km segments (C) for analysis. For the oceanographic variables surveys (B), dark green represents part one routes and blue represents part two, with both parts run within the same week. Segments used for the analysis (C) are depicted at 5 times the survey width of 100 m.
Ornithological Applications 123:1-17 © 2021 American Ornithological Society just taking off from the water were recorded at the location of their initial sighting. Birds landing on the water while a transect was being conducted were not included as sightings. Between 1997 and 2008, murrelet sighting locations along transects were calculated based on observation time, assuming a constant boat speed, while in the latter years, GPS locations were determined. Further details on historical surveys methods and data digitalization are described in Pastran (2020).
In 2018 and 2019, surveys to measure oceanographic variables along the same transects and murrelet observation protocols were conducted to explore how murrelets are influenced by finer-scale prey and oceanographic features. Because we added stops to take measurements of fish schools and water conditions at ~1.5 km intervals, we conducted these surveys over two days within the same week. Oceanographic variables survey "part one" consisted of 18.7 km length of shoreline and 18.6 km of outer transects, and "part two" consisted of 5.1 km of shoreline and 26.3 km length of outer transects ( Figure 2B). Surveys started at 06:30-07:30 and went until 12:00-13:00. We completed 10 part one and 8 part two surveys between May-July of 2018 and 2019.

Segmenting Data
We binned transects into 100 m × 1.0 km grid rectangles, producing 83 segments ( Figure 2C). This segment length enabled analyses at a fine spatial scale but were long enough to result in measurable aggregations of murrelets. Because shoreline transects were not perfectly linear, most of these transect segments had small deviations from the standard 0.1 km 2 area and rectangular shape ( Figure 2C). We therefore accounted for segment area in our analyses.

Hotspot Persistence Analysis
We examined spatio-temporal variation in murrelet distributions with a hotspot persistence method (Sussman et al. 2019). The method creates a map that defines hotspots for each survey, then calculates the percentage of surveys in which each segment was a hotspot. Since May and June surveys had been run consistently from 1997 to 2018, we only used these months to build the hotspot map. For each survey, 3 steps were taken to classify segments as hotspots. First, we calculated an effort-corrected count to correct for small deviations in segment size resulting from the nonlinearity of shoreline transects, and the exact transect lengths, by dividing each segment's count by the area of the segment. Second, we fit a two-parameter gamma distribution to the effort corrected counts for all segments with each survey day (fitdistrplus package in R; Delignette-Muller and Dutang 2015). Finally, a segment was classified as a hotspot if its effort-corrected value was above the 75th percentile of the gamma distribution for a given survey day.
Using the 75th percentile as the threshold enabled us to illustrate important marine areas without overestimating or underestimating an areas importance. This procedure effectively standardizes surveys for the total number of murrelets present and weighs each survey equally regardless of the total counts. Two surveys in which murrelets were only counted in a single segment violated the assumptions behind this assignment process and were excluded from the analysis. The final number of surveys used was 73. After applying the above steps to each survey event over the 22-year period, we calculated the percent of surveys during which each segment was identified as a hotspot.

Variables Considered
We assembled habitat features we hypothesized would be associated with murrelet use and weather and time were recorded as surveys took place (Table 1; detailed maps in Supplementary Material Figures S1-S4). The physical environmental variables were: distance to shoreline (all shoreline types), distance to sandy shoreline, distance to streams, an index measuring the proximity and abundance of potential nesting habitat, water depth and percent sand bottom substrate. The weather and time variables from LBCS long-term surveys were time of day, percent cloud cover, precipitation and wind speed. The dynamic oceanographic variables surveys in 2018 and 2019 measured SST, thermal mixing, and prey abundance in-situ with at-sea surveys (Table 1).

Physical Environmental Variables
Physical environmental variables were collected in the field or constructed from online sources to evaluate their relationships to murrelet distribution across years (Table  1). We collected depth and seafloor sediment data in the summer of 2019. Sediment collections were used to quantify murrelet associations with percent sandy bottom. We grouped the categories fine and coarse sand together into the more general "sand bottom" term to account for small shifts of grain size that may have occurred over the 22-year study period. Pacific sand lance have predominantly been found in sandy sediment in waters 60 m depth or less (Ostrand et al. 2005), so collections were made down to 60 m as the maximum depth. Points that exceeded 60 m were classified as "Deep" and assigned zero percent sand. For other points, we attached 60 m of crab line to a Petite Ponar grab to obtain sediment. At each collection site, the grab was dropped 3 times. If no sand or only rock was collected after the 3rd drop, we assumed zero percent sand. After collection, we dried the samples on a wood-burning stove (low heat) for 24-48 hr. We then shook the samples through a sieve series (4 mm, 2 mm,1 mm, 0.5 mm, 0.25 mm, 0.125 mm, and 0.063 mm) for ~15 min, then weighed and recorded each layer's mass. We categorized
Cooler SST has been linked to nutrient enhancement and prey aggregations (Chavez et al. 2003).
Thermal mixing (MIX) Mixed or stratified Categorical variable; difference taken between temperature readings at 5m and 10m. sorted into "mixed" and "stratified. " A mixed thermal layer can indicate nutrient mixing, which promotes productivity (Behrenfeld et al. 2006 Higher occurrence of fish schools have correlated to murrelet distribution in reflecting productivity of waters (Haynes et al. 2011).
sediment as sand if it was ≥0.063 mm and ≤2 mm and calculated the percent of the total sample by dry weight, that fell within that range. For water depth, a Lawrence Elite Yi 7 sonar "Fish Finder" was used to record continuously along transects. Potential nesting habitat data adjacent to the transects were taken from a habitat suitability map provided by the B.C. Ministry of Environment (Mather et al. 2010). Nesting habitat was defined as any area where murrelets could potentially nest based on Mather at al. 's (2010) criteria. Following Ronconi (2008, S.A.P. personal communication), we created an index testing the relationship between murrelet counts and potential nesting habitat proximity and abundance, using inverse distance weighting (IDW) in ArcGIS Pro 2.3.0 (details in Pastran (2020)). We screened 3 potential commuting distances to identify the most appropriate spatial scale to calculate potential nesting habitat index values: radii of 5, 10, and 30 km from each segment centroid (Hull et al. 2001, Lorenz et al. 2017. To find which spatial scale was most informative, we plotted the mean relationship of murrelet counts per segment to the nesting index of each given radius. The nesting index using a 5-km maximum distance had the strongest relationship with murrelet counts and therefore was used in the subsequent candidate models. We treated this layer as static because, after inspection from Google Earth Pro images, <4 km 2 of forest had been harvested within the 5-km buffer zone between 1997 and 2003, and no harvesting was detected after 2003. The remaining environmental variables were collected from online sources. Shoreline type was mapped using the physical shore-zone polygon from the GeoBC database (https://catalogue.data.gov.bc.ca/dataset/shore-unitclassifications-line). Details on classification are given in Howes et al. (2005). The distances to shoreline and to shoreline type that contained sandy substrate were calculated from the segments center using the Near tool in ArcGIS Pro 2.3.0. Stream data were taken from the British Columbia Stream Atlas Network (https://catalogue.data. gov.bc.ca/dataset/freshwater-atlas-stream-network). Distance to streams was also calculated with the Near tool as the distance from the center of a given segment to the closest stream head of order 2 or higher (Gray 2010).

Oceanographic and Prey Variables
The oceanographic variables surveys in 2018 and 2019 were taken to measure the influence of oceanographic and prey variables on murrelet distributions (Table 1). We conducted surveys during the morning and early afternoon, with no systematic differences in prey availability expected during this temporal window. A temperature/salinity probe ±0.1 (YSI Pro 30) was used to record temperature at two depths. ArcGIS Pro 2.3.0 was used to interpolate temperature points applying the spline tool with the tension setting to create continuous surface layers for temperature values. The temperature reading at 5-m depth was treated as the SST, as temperature readings closer to the surface represent local heating rather than reflecting vertical mixing conditions (Sakuma et al. 2000). Each transect segment's center point was spatially joined to the corresponding temperature value for a given survey date. To examine the effect of thermal mixing (MIX), the difference between 5-m and 10-m temperature values was calculated, plotted, and interpolated in the same manner as the SST layer (Becker and Beissinger 2003). In correspondence with the temperature probe's accuracy, difference of 0.1 or higher was classified as "stratified" and smaller values as "mixed".
The distributions of potential prey along transects were recorded simultaneously with sea-surveys. We used a Lowrance HST-DFSBL Transom-Mount Skimmer Transducer attached to the skiff 's stern submerged 25 cm below the waterline. This transducer was set to 200 kHz for higher resolution and had a beam angle of 12°. Sonar videos were recorded along each transect as surveys were conducted, with the file stored for later processing. Processing sonar recordings along transect lines was done using Reefmaster 2.0 software, which allows the viewing of sonar videos and the vessel's location at any given time. Prey occurrence was recorded as the number of fish schools observed at a given location down to a depth of 60 m (Haynes et al. 2011), binned into ~1-km segments. We defined a fish school as a free-floating cloud on the screen, or 10 or more individuals counted within 100 m of one another. Because schooling is a visual phenomenon, we set 10 individuals as a threshold value for scoring a school (Gautrais et al. 2008). Two attributes of these tabulations should be kept in mind. First, the transect width of fish recorded underwater by the sonar was smaller than that the 100-m transect width for murrelets, giving the possibility of inflated bird counts relative to fish schools. However, this is a constant difference in all surveys. The second limitation is that the number of occurrences of fish schools does not account for the size of each fish school recorded, thereby does not represent the actual density of prey in the water on a given survey. To assess repeatability of fish school tabulations, two observers analyzed the same 5 transect videos of sonar records. We quantified inter-observer agreement using the intraclass correlation coefficient, calculated with irr in R 4.0.3 (Wolak et al. 2012).

Habitat Preference Analyses
For the murrelet habitat preference analyses, we used the counts per segment as the response variable. We completed sets of candidate generalized linear mixed models (GLMM) predicting murrelet counts per transect segment. The GLMM framework handles non-normal response data and can account for nested, non-independent sampling (Brooks et al. 2017). Because count data were over-dispersed, all models used a negative binomial error distribution (log linked) fit to a "nbinom2" family, which assumes that variance s 2 increases quadratically with the mean m (s 2 = m[1 + m/q] with q >0). Models were fit in R 3.61 using the glmmTMB 1.0.2 function in the TMB package, which uses the Laplace approximation to integrate over random effects. Each candidate model included random effects of year, Julian day, and transect segment nested within transects. The survey length of each segment was included as an offset to adjust for the minor variations in survey area.
Potential environmental predictor variables (Table 1) were checked for collinearity by calculating all pairwise Pearson's correlation coefficients. Significant correlations of r ≥ 0.7 were found between distance to streams and distance to sandy shoreline, as well as distance to shore and depth. Consequently, percent sand bottom (SAND bottom ), depth (DEPTH), distance to streams (STREAM dist ), and the nesting habitat index (NEST index ) were kept for subsequent analyses along with weather and time variables. We standardized and centered predictor variables by subtracting the mean and dividing by the standard deviation to directly compare the magnitude of the effect size of the variables.
To account for spatial autocorrelation, we included the spatial "hierarchical" structure into the GLMM that specified that segments were nested within transects. This method assumes that the dependence of segments within their given transect is constant. We also tested for evidence of spatial autocorrelation after model construction using a correlogram test, which calculates a Moran's I value over increasing spatial lags (Fortin et al. 2002). To compare how well the model's spatial variables accounted for spatial autocorrelation, we first summed all counts across years within their specified segment, ran the correlogram test on the raw murrelet counts, and then ran a second correlogram test on the residuals from the spatial model (Fletcher and Fortin 2018). We compared the results to see what changes in spatial relatedness occurred. This was done separately for both the long-term and oceanographic variables survey datasets.
For the long-term dataset, two sets of candidate models were assembled a priori, consisting of all combinations of physical environmental variables, and all combinations of weather and time variables (Table 2). This was done to compare the relative effects of the two variable types on murrelet counts. Once the top models from both candidate sets were selected from the long-term dataset, a combined model that took variables from both of top scoring models, as well as the random effect coefficients previously listed, was run to test their relative effects on the murrelet counts. The analyses of oceanographic and prey variables based on oceanographic variables surveys (2018-2019) compared a priori candidate models of physical environmental, oceanographic, and physical environmental and oceanographic groupings (Table 2).
Top models within the model set considered were selected using the lowest Akaike Information Criterion corrected for small sample size (AIC c ). Models with ΔAIC c < 2 were considered to have substantial support from the data relative to other candidate models (Anderson et al. 1998, Richards 2005. ΔAIC c refers to the difference in AIC c scores between a given candidate model and the top candidate model (Anderson et al. 1998). We also assessed the statistical significance of independent variables from top models, the incidence rate ratios (IRR), and their 95% confidence intervals (CIs). IRR values are analogous to the odds ratios usually reported to assess results from logistic regressions but applied to negative binomial distribution. The IRR indicates the change in the dependent variable in terms of a percentage increase or decrease of counts with respect to the standardized unit change of the predictor variable (Cummings 2019). Values below 1 indicate a negative relationship and counts above 1 indicate a positive relationship, with 95% CI ranges overlapping 1.0 indicating a null effect.
To evaluate model performance, we calculated the conditional R 2 GLMM , which describes the variance explained by both the fixed and random effects, and the marginal R 2 GLMM , which describes variance explained by fixed effects alone (Nakagawa and Schielzeth 2013). To describe the marginal effects from the given GLMM model, we plotted the predicted values of each response variable with their associated 95% confidence intervals to evaluate the support for each variable at different numeric or categorical values of the given independent variable, when all other independent variables were set to zero (Lüdecke 2018).

Evaluating the Direct Influence of Prey Occurrence
We examined direct differences in the relationships between physical environmental variables and murrelet counts when controlling for fish abundance by fitting predictive counts of murrelet from our model at two different fish school count levels. We used the ggpredict function from the ggeffects 1.0.1 package in R 4.0.3 (Lüdecke 2018) to run predictive models of murrelet counts for each independent variable, showing their conditional relationships with the other variables set to zero, at both the upper and lower quartile of fish school counts. We visually inspected the figures to look for changes in predicted values between fish school levels.

Habitat Preference Persistence Heatmap
Counts within segments ranged from 0 to 92 during the 22 years of May and June surveys used for heatmap construction (Figure 3). A mean of 1.16 murrelets were counted per transect segment. Individual segments were classified as hotspots during 0-63% of surveys over the 22 years. Shoreline segments had a higher percent of hotspots (31%, 14-63%), than offshore segments (13%, 0-43%). The majority of segments classified as hotspots were on the southern part of the shoreline. The furthest northwest transects had two persistent hotspot segments (>50% surveys), which were close to Cumshewa Head, a major landmass farther north. The segments farthest offshore from Louise Island had the lowest percentage hotspots. Murrelet usage of Laskeek Bay was thus clearly strongly biased toward inshore areas, but there was also substantial variation among the inshore segments.

Habitat Preferences from Long-Term Surveys
Correlograms of raw murrelet counts among transect segments provided strong evidence of spatial autocorrelation ( Figure 4A), with positive autocorrelation linearly decreasing until 6,000 m. This pattern disappeared when model residual values were plotted ( Figure 4B), indicating that spatial variables accounted for spatial autocorrelation in the counts. The results of all models are listed in Supplementary Material Table S1. For the long-term analysis, models of physical environmental variables showed strong support for a single model (w i = 0.83; Table 3), which included distance to streams, depth, percent sand bottom and the nesting habitat proximity index. This model TABLE 2. Candidate models used to examine the association between physical environmental and dynamic (weather and time or oceanographic) variables and counts of Marbled murrelets (Brachyramphus marmoratus). In addition to variables listed below, each model also included random effects of year, Julian day, transect segment nested within transects, and an offset of segment length. Combined model from long-term dataset assembled post-hoc therefore not included as candidate model.
Candidate models created from the weather and time variables produced 3 candidate models with similar support ( Table 3). The top-ranked model included rain and time (w i = 0.40), with rain, time, and cloud cover present in the second-ranked model (∆AIC c = 1.71, w i = 0.15) and rain, time and wind speed as the third top model (∆AIC c = 1.93, w i = 0.13). Counts were significantly higher in the morning (IRR = 1.32, 95% CI: 1.05-1.65), and when it rained (IRR = 1.74, 95% CI: 1.05-2.82). No models consisting of weather and time variables alone provided strong predictive power in the absence of random effects, with marginal and conditional R 2 GLMM values accounting for around 1% and 55% of the variance, respectively. We then analyzed combined models to include all the fixed variables from the top physical environmental and weather and time candidate models. The combined top model (Table 4) provides similar results as the separate models. Murrelet counts significantly increased as: water depth decreased (IRR = 0.67, 95% CI: 0.53-0.83; Figure 5A), the distance to streams decreased (IRR = 0.49, 95% CI: 0.38-0.62; Figure 5B), the percent sand along the ocean bottom increased (IRR = 1.16, 95% CI: 1.04-1.30; Figure 5C), locations were more proximal to abundant potential nesting habitat (IRR = 1.23, 95% CI: 1.07-1.42; Figure 5D), counts were made in the morning (IRR = 1.29, 95% CI: 1.03 -1.61; Figure 5E) or when it rained (IRR = 1.72, 95% CI: 1.05-2.82; Figure 5F). For the combined model, the marginal R 2 GLMM explained ~24% of the overall variance, and the conditional R 2 GLMM explained ~52% of the variance, with little difference from the physical variables alone.

Habitat Preferences from Oceanographic Variables Surveys
The oceanographic variables survey data showed that overall number of murrelets in 2018 was almost 4 times lower (n = 246) than in 2019 (n = 926). Similarly, the   -term (1997-2018) and oceanographic variable surveys (2018-2019). The top models of bird counts from at-sea surveys (∆AIC c < 2.0) are reported. Models predict the average number of counts per segment and incorporate the for the long-term survey year (n = 22), Julian day (n = 60), and segment (n = 83) nested within transects (n = 18) as random effects, with oceanographic variable surveys having a similar stucture, only with the number of years (n = 2) and Julian days (n = 18) reduced. K is the number of parameters estimated, AIC c is the Akaike's Information Criterion adjusted to small sample size, ΔAIC c is the differences between the AIC c of each model to the lowest AIC c score, w i is the relative likelihood of each model in relation to all other models in the candidate set, marginal R 2 GLMM is the variation explained by fixed factors, and conditional R 2 GLMM is the variation explained by the fixed and random factors combined. ) for the top model explained ~44% of the variation, and the fixed and random effects (conditional R 2 GLMM ) explain ~64% of the variation. The second top model explained 43% of the fixed effect variation and ~66% of the fixed and random effect variation. The top model includes all oceanographic variables (Table 3) and has 3.9 times more support, explaining ~2% more of the fixed variable variation than the candidate model that includes only the 4 physical environmental variables. From the top model, thermal mixing (IRR = 1.70, 95% CI 0.45-2.60) and fish schools (IRR = 1.20, 95% CI 1.01-1.43) were found to be significant oceanographic variables (Table 4). Contrary to our initial expectations, murrelet counts were significantly higher when water was stratified rather than mixed. We plotted the number of stratified recordings by location for the two field seasons (Supplementary Material Figure  S5). A high number of stratified recordings occurred in the southern bay waters. As expected from previous modeling, there was also a significant relationship showing higher murrelet counts with shallower water depths (IRR = 0.51, 95% CI 0.33-0.80; Table 4) and shorter distances to stream heads (IRR = 0.37, 95% CI 0.24-0.57; Table 4).

Prey Occurrence
There was a high repeatability of fish school the recordings by different observers (ICC score of 0.80 (95% CI: 0.63-0.90). To test whether the number of fish schools influenced the strength of association with physical environmental variables, we simulated counts using only segments with the upper or lower quartiles of fish school counts. We did not find substantial differences in the strength of association when fish schools were high vs. low, with considerable overlap occurring between the 95% confidence intervals (Figure 6).

DISCUSSION
The aim of this study was to better describe the marine distribution of Marbled Murrelets in Haida Gwaii and quantify factors responsible for their habitat preferences. We used a long-term marine survey dataset  and conducted oceanographic variables surveys (2018)(2019) to explore and test relationships between usage consistency and coastal, bathymetric, oceanographic, prey, and nesting habitat features. The hotspot persistence map shows strong consistency in terms of which areas were classified as hotspots across the 22 years of surveys. Higher numbers of murrelets are found adjacent to Louise Island compared with ~0.5 km offshore or farther. Modeling captured this pattern with higher counts being closer to streams, in shallower waters, in marine areas that contain sandy substrate and in transect segments with higher habitat nesting indices. Overall, the physical environmental variables were far better predictors than the concurrent environmental variables, but adding the oceanographic and prey factors increased model performance. Adding in predictor variables to a model can increase R 2 values even if variables are irrelevant and can lead to overfitting Anderson 2002, Hyndman andAthanasopoulos 2018). We therefore cross examined the results of the ∆AIC c , IRR, and R 2 GLMM when considering the impact of variables on murrelet counts. In the oceanographic variables survey results, top AIC c models included the 4 physical environmental variables (distance to streams, depth, percent sandy bottom and the nesting habitat proximity index). However, the IRR scores from the oceanographic variable surveys indicated that only distance to streams and depth significantly affected murrelet counts. This apparent inconsistency is likely due to the lower predictive power we had when modeling the oceanographic variables surveys due to smaller sample sizes, resulting in larger confidence intervals. Because all 4 physical environmental variables were included in the top AIC c candidate model for the long-term dataset and for the top oceanographic models, with IRR scores supporting the significance of these variables in the long-term model, further investigation is warranted into their effects.

Freshwater Runoffs and Stratified Water
A strong correlation between murrelet counts and proximity to streams was repeatedly seen in our results. This relationship has been found in a number of other murrelet studies , Haynes et al. 2011. Two hypotheses may account for the relationship. First, streams are often used as flyways to nesting sites. The reasoning is that murrelets follow streamways when commuting to feed young to avoid unnecessarily expensive climbing flight over watershed boundaries (Barbaree et al. 2015). This hypothesis may be more applicable to sites with more dramatic topography than is present around Laskeek Bay. The second hypothesis is that areas with freshwater and saltwater mixing have higher productivity than areas that do not (Yen et al. 2004), and thus, provide better foraging opportunities. This second hypothesis is a plausible explanation for our results, especially when taking into Stratified water and freshwater runoffs can work together to create productive zones. In Kachemak Bay and Cook Inlet, Alaska, there is a strong association between sheltered stratified waters with an inflow of freshwater from rivers and streams and the abundance of pelagic schooling fish such as Pacific sand lance and juvenile herring (Abookire et al. 2000). The authors believed that areas around river outflows have higher inputs of nutrients, coupled with the fact that stratified waters can create stability and promote primary productivity by keeping nutrients at the surface. Areas that contain both these components are more prone to an abundance of life. A similar phenomenon is likely occurring in the nearshore southern portion of Laskeek Bay.

Correlations with Pacific Sand Lance Habitat Features
Pacific sand lance are an important food source for murrelets and are often found in coarse grain sand in shallow areas (Ostrand et al. 2005). Therefore, it is not surprising that various marine habitat studies have found that sandy shorelines and underwater substrate predict murrelet presence , Yen et al. 2004, Ronconi 2008. As expected, we found a significant positive association between percent sandy sediment and murrelet counts. We also collected a Pacific sand lance from a sediment grab, and saw murrelets holding Pacific sand lance in their bills in Laskeek Bay during the 2018 and 2019 field seasons. This contributes to growing evidence that Pacific sand lance is an important element in murrelet diet.

Connection of Marine Distribution to Potential Nesting Areas
Murrelet densities were higher with greater proximity to and abundance of potential nesting habitat. This type of relationship has been documented a number of times (Yen et al. 2004, Ronconi 2008, S.A.P. personal communication, Raphael et al. 2015, Lorenz et al. 2016, but this FIGURE 6. Conditional relationships between (A) water depth, (B) distance to streams, (C) percent sand bottom, and (D) nesting habitat index to the predicted Marbled Murrelet (Brachyramphus marmoratus) counts within segments using the top-ranked model from the habitat preference analysis. The blue line represents 6 fish school counts held constant, and the red line represents 0 Fish school counts held constant within the model. Bands indicate 95% confidence intervals. study showcases the relationship at a finer geographical scale. Most murrelet nests occur within 30 km of shorelines (Environment Canada 2014, Barbaree et al. 2015), but birds have been documented nesting as far as 145 km inland (Lorenz et al. 2017). Here we show a relationship at a scale of 5 km or less, although we have no information on the specific nesting locations of surveyed birds. If commuting flights expose murrelets to greater predation risk than being on the water, or can substantially increase daily energy expenditure, it remains adaptive to minimize their distances (Hull et al. 2001).

Influence of Fluctuations in Fish Schools
The two years with oceanographic variables surveys had large parallel differences in the number of murrelets and fish schools recorded. These data suggest that ocean productivity in a given year may directly affect murrelet local population abundance. Becker and Beissinger (2003) noticed that murrelets were distributed farther from the two primary breeding area flyways in a year when fewer prey were available at their California site. However, we found no evidence that the strength of association between counts and physical environmental variables was higher when more prey was available. It is possible murrelets forage closer to the stream heads when fish school counts are higher, though this trend was not significant. However, sonar does not detect Pacific Sand Lance schools (Robards et al. 1999) and, this limitation should be taken into account when interpreting model outputs. Despite this, the information on the number of occurrences of fish schools provides a snapshot of how productive the transects and overall waters were at a given time.

Management Implications
The Laskeek Bay at-sea surveys provide the only long-term series available for Haida Gwaii, and this analysis has provided novel information for the area. The hotspot persistence map identifies high and low use areas over 22 years. The consistency in use throughout the years highlights the importance maintaining breeding habitat to support local populations. Inshore waters that are prone to stratification, and are in close proximity to freshwater inputs, have the highest foraging use by Marbled Murrelets and should be considered high-quality marine habitats in management planning. Understanding how murrelets respond to changing marine conditions can help pinpoint explanations for distributional shifts or population declines. This work can aid in the creation of a coastwide marine habitat suitability map for murrelets, facilitating effective policy decisions.

SUPPLEMENTARY MATERIAL
Supplementary material is available at Ornithological Applications online.