Assessing year-round habitat use by migratory sea ducks in a multi-species context reveals seasonal variation in habitat selection and partitioning

1 –––––––––––––––––––––––––––––––––––––––– © 2020 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Jean-Michel Gaillard Editor-in-Chief: Miguel Araújo Accepted 24 July 2020 43: 1–18, 2020 doi: 10.1111/ecog.05003 43 1–18


Introduction
Long-distance migration of wildlife involves a delicate balance between the rewards of occupying favorable habitats year-round and the challenges of the extreme journeys required to access those habitats (Alerstam et al. 2003). High fidelity to specific breeding and wintering sites and migration routes, as well as reliance on evolutionary and hereditary cues, make long-distance migrants highly susceptible to mismatches between timing of migratory movements and availability of key resources (Robinson et al. 2009, Liedvogel andLundberg 2014). In addition, many species of migratory wildlife concentrate in favorable habitat sites during migration, making them vulnerable to localized perturbations (Berger 2004, Buehler andPiersma 2008). Migratory birds, which represent about 19% of all avian species, undertake some of the most complex and challenging migrations of any taxa (Alerstam et al. 2003) and are currently experiencing widespread population declines (Kirby et al. 2008). Despite these negative trends, however, habitat conservation for migratory birds often falls short of levels achieved for nonmigratory species (Runge et al. 2015).
Long-distance migrants present a challenge to traditional, site-based models of habitat conservation developed for non-migratory species (Martin et al. 2007, Singh andMilner-Gulland 2011). Some migrants rely on discrete habitat sites with varying characteristics, which often span multiple nations with varying conservation priorities (Berger 2004, Shillinger et al. 2008). In addition, because populations of migratory species can be affected by factors occurring throughout the annual cycle, improvements to habitat used during one period or season may not be sufficient to meet conservation targets if habitat during another part of the annual cycle is limiting (Marra et al. 2015). As a result of this complexity, conservation strategies for migratory species requires a detailed understanding of habitat use throughout the annual cycle in order to develop coordinated management across multiple sites and seasons (Wilcove and Wikelski 2008). Understanding the habitat requirements of migratory birds throughout the annual cycle can help prioritize sites and resources for conservation, as well as identify potential periods of vulnerability (Kirby et al. 2008, Allen andSingh 2016).
As part of improving habitat conservation for migratory species, conservationists may need to account for species co-dependency on the same habitat, which requires sound multi-species approaches to evaluating habitat requirements (Witting and Loeschcke 1995, Mittermeier et al. 1998, Brooks et al. 2006. Evaluating habitat needs across guilds of similar species allows managers to set conservation priorities at the landscape scale, focusing on areas of particular ecological significance to maximize the impact of conservation action (Moilanen et al. 2005, Block et al. 2011, Hindell et al. 2011, Raymond et al. 2015. Multi-species assessments can be useful for evaluating the effects of future habitat loss or environmental change on migratory pathways (Martin et al. 2007), and identifying representative species that can be targeted for monitoring (Lambeck 1997, Bonn andSchröder 2001). While migration strategies vary among birds, the physiological demands imposed by long-distance migration often lead to concentrations of multiple avian migratory species at locations with highly predictable resource availability and favorable geographic features (Roberts et al. 2001, Mehlman et al. 2005, and the resulting interspecific competition for shared resources may have population-level effects (Péron and Koons 2012). Thus, identifying sites or habitat requirements shared by multiple species of long-distance avian migrants at specific periods in the annual cycle can provide important insights into the biogeographic and conservation needs of migratory species that transcend species-specific habitat models (Donovan et al. 2002).
Sea ducks (tribe: Mergini) are a suite of migratory bird species that breed in boreal and arctic habitats and winter along the coasts of large water bodies. Although their nonbreeding ranges broadly overlap and they often congregate in mixed-species flocks during migration and in the winter, different species of sea ducks have evolved to exploit subtly different habitats and resources within shared habitat areas, and these differences may result in interspecific variation in key habitat areas and movement networks (Lamb et al. 2019). Sea ducks also present unique challenges for habitat conservation and management, as their habitat needs vary widely throughout the annual cycle from coastal aquatic non-breeding areas to terrestrial breeding sites far inland (Johnsgard 2010). The long-distance migrations required to travel between those sites must be delicately adjusted to the availability of resources in diverse habitats separated by thousands of kilometers, and often involve migratory bottlenecks in which large portions of the population occupy in shared staging areas (Oppel et al. 2009, Boere andPiersma 2012). Across North America, sea ducks are experiencing widespread population declines, leading to increased interest in identifying key habitats and conservation strategies (Bowman et al. 2015). However, since sea ducks spend much of the year in remote, uninhabited breeding and staging areas or offshore wintering sites, direct observations of their abundance and distribution are expensive, labor-intensive, and may be limited in space and time.
In eastern North America, a continental-scale, partnership-based project designed to track numerous individuals from five high-priority sea duck species (common eider Somateria mollissima, black scoter Melanitta americana, surf scoter Melanitta perspicillata, white-winged scoter Melanitta deglandi and long-tailed duck Clangula hyemalis) using satellite telemetry provides a solid foundation for multi-species assessments of annual-cycle habitat selection (Sea Duck Joint Venture Management Board 2014, Lamb et al. 2019). Data from this project has filled large gaps in understanding of sea duck breeding habitat in North America, revealing previously unknown distributional patterns of species (particularly scoters) at boreal and arctic breeding sites and informing the development of surveys in areas difficult to access using traditional methods (Reed et al. 2017). Contributed data also has shown how sea ducks use coastal wintering habitats, and potential spatial overlap with proposed wind energy installations (Loring et al. 2014, Beuth et al. 2017, Meattey et al. 2019. However, each of these assessments has focused on species-specific, single-season habitat selection. There remains a need to address sea duck habitat requirements throughout the annual cycle in a more comprehensive and integrated manner, both to understand how ecological risks are distributed and to optimize use of limited funding for management and protection (Sea Duck Joint Venture Management Board 2014).
We present a multi-species assessment of year-round habitat use patterns of sea ducks in eastern North America derived from year-round location data on individual sea ducks. Our main goals for the study were 1) to describe single-and multi-species patterns of habitat suitability throughout the annual cycle; 2) to assess between-species habitat partitioning and selection and 3) to identify key biophysical features selected by sea ducks during different periods of the annual cycle. Based on observed ranges of the focal species , we expected to observe a high degree of habitat partitioning during the breeding season, an intermediate level of partitioning during migration, and a relative lack of partitioning during winter with extensive use of similar habitats. Further, based on previous studies, we expected to see static physical habitat features (distance to shorelines, elevation and slope) as the strongest predictors of sea duck occurrence and species counts in both breeding and nonbreeding habitats, with greater reliance on preferred habitat features during migratory staging to facilitate rapid acquisition of energy. Our analyses provide a framework for identifying specific locations, time periods and habitat features for which targeted conservation can provide the greatest shared benefit to species with similar habitat requirements, and for developing coordinated range-wide, multi-species conservation strategies.

Material and methods
We used satellite telemetry to track sea ducks in eastern North America throughout the annual cycle (Lamb et al. 2019). After filtering location data to include only sedentary points, we measured occurrence by season and species across the study area. To relate sea duck occurrence to environmental features, we used a multivariate ordination of available sites based on their ecological characteristics to measure single-and multi-species habitat suitability and evaluate between-species habitat partitioning, and used linear models to examine the relationships of individual ecological covariates to multi-species occurrence.
Study area -we captured sea ducks at multiple areas along the Atlantic coast and Great Lakes of eastern North America during molting, staging and winter (August-March) between 2002 and 2017 (full details of capture locations and species ranges summarized in Supplementary material Appendix 1 Fig. A4). Since sampling the entire range of each species was not possible, we selected sampling locations to represent locations and time periods with particularly dense species concentrations, which occur primarily during nonbreeding. Capture efforts for long-tailed ducks, and surf scoters focused primarily on wintering sites, and transmitters were distributed proportionally based on concentrations of ducks observed during the Atlantic Winter Sea Duck Survey. Capture efforts for white-winged scoters occurred mainly on molting sites and the ones for black scoter focused on pre-breeding migration sites in the St Lawrence River, with additional captures of the three scoter species at other annual stages (Sea Duck Joint Venture Management Board 2014). Additional sampling of long-tailed ducks at Lake Michigan was added late in the project to account for gaps in observed data (Fara 2018). Sampling of common eider was limited to one of three eastern subspecies (S. m. dresseri) during breeding and wintering periods.
Transmitter deployment -we captured both sexes of subadult and adult ducks on water using a combination of decoys and netting. We captured the majority of sea ducks using mist nets (1.3 × 18 m 2 , 127-mm mesh), which we positioned either above the water on floating poles to catch ducks in flight (Brodeur et al. 2008), submerged as gillnets to catch birds during dives (Breault and Cheng 1990), or in the Great Lakes, suspended nets placed horizontally underwater that were raised to a vertical position ducks just before birds flew over the site to capture them in flight (Ware et al. 2013). Additional capture techniques included night-lighting and dip-netting for wintering ducks roosting on the water, netgunning and fishing nets (Lamb et al. 2019). Veterinarians experienced in avian surgery implanted 26-50 g coelomicimplant Platform Transmitter Terminals (PTT) (Microwave Telemetry, Columbia, MD, USA; Telonics, Mesa, AZ, USA; Geotrak, Apex, NC, USA; mention of commercial products does not imply US or Canadian Government endorsement) into the abdominal cavity following implantation techniques described by Korschgen et al. (1996). Individuals were selected for transmitter attachment based on body mass, such that transmitter mass represented less than 5% of overall body mass (Phillips et al. 2003). Transmitters followed varying duty cycles consisting of 2−4 h 'on' periods followed by 10−120 h 'off' periods, resulting in one location every 0.5-5 d (for specific duty cycles by deployment event, see Lamb et al. 2019, Supplementary material Appendix 1). We excluded data collected during the first 14 d following surgeries to minimize potential biases in assessments of habitat use patterns and movement dynamics due to surgery (Esler et al. 2000).
Spatial data processing -unprocessed satellite telemetry data vary in quality of location estimates based on the configuration and number of satellites used to obtain each location, and time intervals between location estimates also may vary. Data were initially processed using a hybrid Douglas-Argos filter to remove redundant and erroneous locations (Sea Duck Joint Venture 2015). We then used a switching state-space model (Jonsen et al. 2005) to simultaneously determine the most probable track for each individual given the observation error associated with each location (i.e. error correction) and produce a regular track from irregular data with varying uncertainty (i.e. interpolation). This modeling approach also allowed us to classify successive locations based on patterns in inter-location distances and turning angles. Although behavior varies at finer temporal scales, our ability to detect state changes was limited by the spatiotemporal scale of data collection, which would not have identified fine-scale changes such as foraging movements. We therefore restricted behavioral classification to two day-to-day states: sedentary, in which movements between successive locations were characterized by short distances and frequent directional change, and transient, in which locations were widely spaced and directional change infrequent, corresponding to multiday directed movements such as migration and dispersal.
To allow the model sufficient information to interpolate individual tracks, we removed all individuals that had fewer than 50 good-quality locations (Argos Location Classes 1-3, or < 1500 m error radius; typically, one month or less of data) prior to analysis, leaving 476 individuals (Supplementary material Appendix 1 Table A1). We did not interpolate over time periods of > 7 d between successive locations, because longer temporal gaps produce unrealistic movement trajectories (Jonsen et al. 2005). Based on the duty cycles of transmitters, the maximum programmed period between locations for a correctly-functioning unit was 120 h, with most units sampling more frequently; thus, 90% of locations were separated by ≤ 4 d, and 78% of locations were separated by ≤ 3 d. Average sampling intervals varied among species from 2.3 d (surf scoters) to 3.6 d (white-winged scoters).
We ran all models in the 'bsam' package (Jonsen et al. 2005, Jonsen 2016) in R (R Core Team) using a switching first difference correlated random walk model with a one-day timestep, 5000 burn-in samples for model training, and 5000 posterior samples for analysis. To reduce autocorrelation, we retained every 5th posterior sample and used a 0.1 smoothing parameter. Model outputs included probable daily locations, as well as a score from 1 to 2 (hereafter, b) indicating the average assignment of the location to either a transient (1) or sedentary (2) behavioral state across all retained posterior samples. We checked model fit and convergence using the 'diag_ssm' function in 'bsam', which includes trace plots, density plots, autocorrelation plots and shrink factor plots, and visually assessed the resulting tracks to ensure that state assignments corresponded to periods of migration and residency. Behavioral state assignments were proportionally similar across Argos location classes (Supplementary material Appendix 1 Fig. A6).
To assess habitat use, we removed all transient locations (b < 1.5) and assigned the remaining sedentary locations to one of four seasons -wintering, pre-breeding staging, breeding and post-breeding staging/molt -based on its timing within the annual cycle and relative to long-distance displacements, thus limiting our analysis to terrestrial and marine habitats. Although aerial habitats used for migrating between breeding, staging and wintering sites are undoubtedly important, measuring flight altitudes and atmospheric covariates was beyond the scope of the present study. Median start dates were 3 November for wintering, 22 April for pre-breeding migration, 1 June for breeding and 25 July for post-breeding migration and molt. Although most flight feather molt takes place following the start of post-breeding migration, breeding females may molt their flight feathers before departing the nest site; thus, although most molt locations are included in post-breeding migration, some may also be included in the breeding season. We then calculated 95% kernel density estimates for all species within each season using the ks package (Duong 2007) in R.
To examine variation in habitat use across the study area, we created a 100 km 2 hexagonal grid within each season's 95% kernel area to match the resolution of our habitat variables (≤ 0.1 degree, or approximately 10 km). Although suitable habitat likely exists outside the 95% kernel area of occurrence points, we chose not to predict habitat suitability beyond observed use areas in order to ensure that background habitat values accurately represented available habitat characteristics (Franklin 2010). For each season, we calculated habitat suitability for each species and overlaid maps to obtain an estimate of multi-species suitability (see 'Habitat suitability' section). We excluded dresseri common eider from our analysis of breeding season habitat selection, because their breeding habitats occur primarily on offshore islands (Goudie et al. 2000). Breeding eiders forage and raise chicks in the marine environment, and therefore rely on different prey resources than breeding scoters and long-tailed ducks occupying freshwater systems. Moreover, habitat used by breeding eiders remains open year-round, whereas access to inland arctic habitats is limited by freeze/thaw dynamics. Thus, the habitat characteristics driving eider occurrence during breeding are likely very different from those governing occurrence of the other four species (Johnsgard 2010).
Environmental covariates -we chose separate sets of predictor variables to predict occurrence in inland habitats (primarily used during breeding) and nearshore aquatic habitats (primarily used during migration and wintering time periods) (full details of environmental datasets are given in Supplementary material Appendix 1 Table A2). The variables we used were not strongly collinear (R 2 < 0.6 for all pairs). In all cases, we used long-term averages to represent dynamic covariates (Supplementary material Appendix 1  Table A2). Since sea ducks have strong interannual fidelity to both breeding (Phillips andPowell 2006, Takekawa et al. 2011) and non-breeding sites (Robertson and Cooke 1999, Lepage et al. 2020, we expect that habitat selection is driven largely by individual life histories and previous years' experience rather than year-to-year conditions. To represent potential predictors of occurrence in inland terrestrial breeding habitats, we used nine variables, including three dynamic climate variables (total annual precipitation, annual maximum temperature and annual minimum temperature [Worldclim2; Fick and Hijmans 2017]) and six static biophysical covariates (distance to marine coastline [World Vector Shoreline; Wessel and Smith 1996]; distance to lake > 50 km 2 and distance to large lake > 300 km 2 [North American Rivers and Lakes; U.S. Geological Survey]; elevation [ETOPO2; National Geophysical Data Center 2006]; slope [ETOPO2]; and landcover type [Landsat 2010; U.S. Geological Survey]). We collapsed eleven landcover types into five general categories: barren, polar (sub-polar/polar grassland-lichen-moss, sub-polar/polar shrubland-lichenmoss and sub-polar/polar barren-lichen-moss), taiga (subpolar taiga needleleaf forest), temperate/sub-polar (mixed forest, temperate/sub-polar broadleaf forest, temperate/subpolar needleleaf forest, temperate/sub-polar grassland and temperate/sub-polar shrubland), and wetland. We selected biophysical habitat variables to represent the nesting habitats of sea ducks, which generally occur in low-lying habitats with extensive wetland cover (Reed et al. 2017), while climate covariates represented factors potentially affecting vegetation structure, nest site availability and foraging conditions at the regional scale. We assessed dynamic climate variables using annual averages, rather than only during the period of occurrence, since precipitation and temperature prior to the breeding season can also affect vegetation structure and nest site conditions (Fu et al. 2014).
We measured environmental characteristics of nearshore aquatic non-breeding habitat using nine variables. Six of these variables were used to assess habitat across all locations, while three were available only for a subset of locations. ). We chose these variables to represent a suite of likely drivers of nearshore habitat variation, particularly the distribution of mollusks and prey populations. Because limited data are available on fine-scale variation in features such as currents and eddies, which have a high degree of shortterm variability in coastal areas (Kaltenberg et al. 2010), we used distance to coastline, bottom slope and depth as proxies for these processes. Net primary production, which integrates chlorophyll concentrations over a range of depths (Behrenfeld and Falkowski 1997), provides an index of aquatic productivity that influences the distribution of consumers at higher trophic levels. Salinity and temperature also influence the distribution of aquatic prey species depending on their osmotic and thermal tolerances. Positive values of depth indicate shallower areas.
Data on the remaining three static variables -bottom substrate, tidal current velocity and aquatic vegetation presence -were available only for some sections of our study area, including the Gulf of St Lawrence (bottom substrate only), Great Lakes (bottom substrate and submerged aquatic vegetation) and U.S. Atlantic coast (bottom substrate, seagrass presence and tidal current velocity). Despite the limited spatial extent of these datasets, we expected that these features would likely play a role in structuring sea duck prey distributions and/or foraging habitat suitability. We therefore excluded these covariates from overall habitat selection and suitability analyses, but separately modeled their relationships to sea duck occurrence and species counts over the measured subset of locations (see 'Relationships of sea duck occurrence and species counts to individual environmental covariates').
We standardized all variables using the seasonal 100 km 2 hexagonal grids used to calculate occurrence. We calculated distance values as the distance from the hexagon centroid to the feature of interest. For all other variables, we resampled the data using the mean value for each hexagon.
Habitat suitability -to map habitat suitability across the study area for individual species and the multi-species assemblage in ecological space, we conducted a multivariate ordination of all habitat variables using a Hill-Smith principal components analysis (PCA) (Hill and Smith 1976), which allows the inclusion of both categorical and continuous variables. For each grid cell, we evaluated species-specific habitat suitability using Ecological Niche Factor Analysis (ENFA: Hirzel et al. 2001), which measures the distance of the cell from the center of the species' distribution (from presenceonly data) in multivariate space. We then identified the top 10% and 25% of cells for each species based on suitability values, and overlaid these sets of cells to map multi-species suitability (i.e. the number of species for which each cell fell within either the top 10% or top 25% of suitability values across all available cells for that season).
Habitat selection and partitioning -to characterize species-level habitat partitioning across the suite of habitat variables available for the full study area, we compared niche position and breadth among species using an Outlying Mean Index (OMI) analysis (Dolédec et al. 2000) using the 'adehabitatHS' package (Calenge 2006). Briefly, OMI calculates the direction of maximum marginality of occupied sites (from presence-only data) relative to available habitat in ordination space, allowing comparisons of niche position and breadth among co-occurring species. The position of each habitat characteristic on the first axis of the OMI represents the marginality of occupied sites on that variable, with positive scores indicating higher-than-average values and negative scores indicating lower values. In ecological terms, greater positive or negative OMI values indicate more specialized niches, while values close to zero indicate more generalist use of available habitats. OMI does not assume specific resource selection functions, and allows differences in individual niche selection to be taken into account when describing the distribution of a group of animals. Using the same PCA ordination we used to map habitat suitability, we conducted OMIs for each season on all individuals, then we averaged the scores of individuals of each species on the first OMI axis to calculate niche position for that species, and calculated the 95% confidence interval of the mean as a measure of niche breadth. To assess habitat partitioning among species, we examined the overlap of individual species' niches, with less overlap indicating greater partitioning on that variable.
Importance of environmental covariates -since our ordination analyses were primarily descriptive, we also quantitatively assessed relationships of species counts to specific environmental covariates using generalized additive models. We modeled the number of species occurring in a given cell (response) as a function of various subsets of smoothed environmental covariates measured across the full study area along with a smoothed latitude-longitude interaction to account for spatial structure in the data (predictors), and compared models in an information theoretic framework. We ran zero-inflated generalized additive models in the 'mgcv' package (Wood 2011, Wood et al. 2016. We selected smoothing parameters using residual maximum likelihood, and used a binomial distribution with a logit link to model zero-inflation and a truncated poisson distribution with log link to model count data. We began with the global model and dropped covariates using backward stepwise selection until dropping additional covariates no longer reduced the Akaike's information criterion (AIC) value of the resulting model. We visually verified the fit of the final models using quantile-quantile plots. Full model selection results are given in Supplementary material Appendix 1 Table A3.
For covariates that were not available across the entire study area, we used forward stepwise selection to test whether these additional data improved the explanatory power of our models. For each additional variable, we ran the top model from the full dataset with and without the covariate on a reduced dataset including only grid cells for which the additional covariates were available. We considered additional variables to be useful if their inclusion improved the AIC value of the resulting model by ≥ 2 points (Burnham and Anderson 2004).
Habitat suitability: breeding -during the breeding season, suitable habitat areas for different species were relatively distinct in space and particularly latitude (Fig. 1a-d). Surf scoters habitat was primarily located to the east of Hudson Bay, white-winged scoter habitat to the south and west, black scoter habitat at higher latitudes on both sides of 7 the bay, and long-tailed duck habitat primarily to the north and northwest of the bay. Areas of high multi-species breeding suitability occurred primarily in the Barrenlands directly to the west of Hudson Bay, as well as in small areas of northern Quebec (Fig. 1e-f ). No grid cells were within the top 10% of highly suitable habitats for all four species (Fig. 1e).
Habitat selection and partitioning: breeding -during breeding, species partitioned habitat primarily according to precipitation, minimum temperature, proximity to marine coastlines, elevation and land cover type (Fig. 2). Longtailed ducks generally preferred areas nearer to coastlines, with lower temperatures, precipitation and elevation, than did scoters (Fig. 2). Among scoter species, surf scoters nested at greater elevations, precipitation levels and temperatures. White-winged scoters nested comparatively further from marine coasts in areas with greater minimum temperatures, while black scoters nested closer to lakes (Fig. 2). In general, most species selected for areas with lesser-than average precipitation, greater minimum and lesser maximum temperatures (i.e. narrower temperature ranges), and lower slopes (Fig. 2). Landcover classes were useful for distinguishing habitat niches: long-tailed ducks nested exclusively in areas with polar landcover types, black scoters occupied areas of polar and taiga cover, white-winged scoters used both taiga and temperate/subpolar cover, and surf scoters used exclusively temperate/subpolar landcover categories.
Habitat suitability: non-breeding -during winter, suitable habitats were relatively similar across all species (Fig. 3a-e). Habitats suitable for black and surf scoters extended over a greater latitudinal range (Fig. 3a-b), while white-winged scoters and long-tailed ducks had more suitable habitat areas in the Great Lakes (Fig. 3c-d) and common eiders were largely 8 limited to the northeastern portion of the multi-species range (Fig. 3e). Areas of multi-species overlap were concentrated along the mid-Atlantic coast of the United States ( Fig. 3f-g). The waters of southern New England, including Cape Cod Bay and Nantucket Shoals, occurred in the top 10% of available habitats for all five species (Fig. 3f ).
During migration, habitat suitability varied by season. Suitable pre-breeding staging habitats were concentrated in the St Lawrence Estuary, Gulf of St Lawrence, Great Lakes and James Bay. Suitable habitats for white-winged scoters and long-tailed ducks were located predominantly in the western part of the study area (Great Lakes and James Bay; Fig. 4c-d), while black and surf scoters and common eiders utilized predominantly eastern habitat areas (St Lawrence Estuary, Gulf of St Lawrence, Nova Scotia and northern New England; Fig. 4a, b, e). Highly suitable habitat for all  species co-occurred in the southern portions of the Gulf of St Lawrence and along the eastern coasts of Nova Scotia and New Brunswick (Fig. 4f-g). In contrast, during post-breeding migration and molt, only common eiders extensively used the same habitat areas as during pre-breeding (Fig. 5e). The other four species (black, surf and white-winged scoters and long-tailed ducks) co-occurred along the shorelines of southern Hudson Bay and northern James Bay, particularly around Sanikiluak, as well as along the Labrador coast and southern Ungava Bay (Fig. 5a-d).
Habitat selection and partitioning: non-breedingthroughout the non-breeding period, sea ducks selected nearshore, shallow-water habitat areas (Fig. 6). Within the occupied habitat area, species showed the greatest amount of habitat partitioning on dynamic habitat variables (net primary production, sea surface temperature and salinity) and bottom slope (Fig. 6).
During winter, species partitioned habitat based on productivity, temperature and bottom slope (Fig. 6a). Black and surf scoters showed strong positive selection for net primary production, weak positive selection for sea surface temperatures, and moderate negative selection for bottom slope. In contrast, white-winged scoter, long-tailed duck and common eider showed weak negative selection for net primary production, strong negative selection for sea surface temperature, and weak positive selection for bottom slope. All species except long-tailed duck selected positively for salinity.
Migratory habitat selection varied between seasons and species. During pre-breeding migration, selection for habitat characteristics was generally weak (i.e. close to zero), with the exception of net primary production, which was strongly positive for all species. Habitat partitioning was limited, with considerable between-species overlap and within-species variation (Fig. 6b). All species showed weak negative selection for sea surface temperature, weak positive selection for salinity (with the exception of long-tailed ducks), and weak or varying positive selection for bottom slope (with the exception of long-tailed ducks). During post-breeding migration, sea ducks displayed relatively strong selection on habitat covariates, with a high degree of partitioning and limited overlap among species (Fig. 6c). All species except long-tailed duck selected strongly for net primary production. For sea surface temperature, salinity and bottom slope, selection varied among species from weak to strong and positive to negative.
Importance of environmental covariates -the best model for species counts during breeding included precipitation, minimum and maximum temperatures, distance to lake, distance to large lake, distance to coast, elevation, landcover and latitude-longitude interaction and explained 90.4% of total deviance (Table 1a). The best model for winter and prebreeding staging species counts across all sites included all covariates (distance to coast, sea surface temperature, salinity, net primary production, depth, slope and latitude-longitude interaction) and explained 78.5 and 76.2% of total deviance respectively (Table 1b-c). The best model for species counts during post-breeding migration included all covariates except for slope and explained 79.7% of total deviance (Table 1d).
For the subset of the study area where tidal current velocity, aquatic vegetation presence and bottom substrate were measured, models of sea duck species counts were improved during all seasons by the inclusion of the additional covariates (Table 1b-d). Tidal current velocity appeared in all highlysupported models for winter and pre-breeding and two of three highly-supported models for post-breeding, with a positive relationship to species counts (Table 2). Aquatic vegetation presence was a positive predictor in the top models for winter and pre-breeding and in one of three highly-supported models for post-breeding (Table 2). Substrate appeared in one of three highly-supported models for winter species counts, with sand, sediment, silt and hard bottom habitats favored over clay and gravel (Table 2).

Discussion
Examining habitat use in a seasonal, multi-species framework allowed us to identify covariates that affected occurrence across a suite of similar species, as well as assess differences in habitat associations among five species of sea ducks. We found that both habitat selection and interspecific habitat partitioning varied by season, with weak selection and strong partitioning during the breeding season, strong selection and partitioning during post-breeding migration and molt, moderate selection and partitioning during winter, and weak selection and partitioning during pre-breeding migration.
Habitat suitability -during the breeding season, most sea ducks in our study occupied nesting habitats in areas of subpolar vegetation near or above the northern limits of boreal forests. Individual species differed substantially in their use of breeding habitats. Some areas, including large portions of the Barrenlands west of Hudson Bay, were among the top 25% of suitable habitats for all four species. Our results expand on previous analysis of black and surf scoters by Reed et al. (2017), which also identified the Barrenlands as an area of high multi-species importance. The areas identified by our models as being of highest multi-species importance, including the Barrenlands and parts of northern Quebec, fall outside the boundaries of most North American breeding bird surveys including the Waterfowl Breeding and Population Habitat Survey (Roy et al. 2019); however, the Barrenlands are currently being targeted for additional pilot surveys to quantify the extent of use by nesting sea ducks (Reed et al. 2017). In our analysis, no sites were in the top 10% of suitable breeding habitats for all species. This suggests that optimally conserving key breeding habitats would require targeting habitats with high single-species suitability as well as multi-species hotspots. Hotspots of multi-species non-breeding habitat suitability were located in southern New England south to the Chesapeake Bay (winter); St Lawrence Estuary, the Gulf of St Lawrence, and southern James Bay (pre-breeding); and Hudson Bay and northern James Bay (post-breeding). These areas correspond to key locations identified in other studies (Silverman et al. 2013, Lamb et al. 2019). Multi-species overlap was high during winter, but decreased during migratory periods as migration routes diverged. Long-tailed ducks and white-winged scoters tended to use more inland habitats during pre-breeding migration. During post-breeding migration and molt, long-tailed ducks remained farther north than other species, and dresseri common eiders used unique coastal staging areas. Thus, while multi-species hotspots are useful for identifying key migratory habitats, single-species models may be necessary to fill gaps in important habitats.
It is important to note that we assessed habitat suitability only within areas used by tracked individuals, and not across the full range of each species. Although non-breeding areas used by tracked birds generally corresponded to known eastern North American non-breeding ranges of their respective species, breeding areas occupied by tracked birds differed in some cases from published ranges (Sea Duck Joint Venture 2015). Since many sea duck species have distinct subspecies or subpopulations that winter in either eastern or western North America, further study would be useful to clarify whether these wintering subpopulations also occupy distinct segments of the breeding range, which might help to explain why some known breeding areas were not used by individuals in our study.
Habitat selection and partitioning -overall, the four species included in our analysis of breeding habitat preferred to nest in relatively flat areas near large lakes. This supports a previous predictive model of scoter breeding habitat (Reed et al. 2017), in which lake proximity and area were key predictors of scoter breeding locations. Our results also suggest that sea ducks select for areas of lesser annual precipitation and smaller ranges of annual temperatures. These mild climate conditions could improve both nest cover and foraging opportunities on invertebrate prey during the breeding season, and may also reduce energy costs of incubation and danger of nest and duckling loss from exposure to cold and wet conditions (Mallory 2015).
During non-breeding, geophysical aquatic habitat features that were consistently important across species and seasons included shallow water depths and nearshore locations. These results are consistent with previous studies based on aerial survey data suggesting that non-breeding sea ducks are closely confined to nearshore waters (Silverman et al. 2013, Winiarski et al. 2014, Smith et al. 2019. The sea ducks we investigated generally selected for waters with relatively high productivity across all seasons. Most species also preferred lower temperatures, higher salinity and greater bottom slopes, although selection varied between seasons and species. Previous studies of survey data from the same five species during winter by Zipkin et al. (2010) and Silverman et al. (2013) showed similar species-specific relationships with sea surface temperature and bottom slope to the patterns we observed; however, their models did not include either salinity or net primary production.
Niche partitioning among species was strong during breeding. Of the four species included in our breeding habitat analysis, long-tailed ducks and surf scoters were the most dissimilar, with black and white-winged scoters occupying intermediate habitats. For several habitat variables, including temperature, precipitation and elevation, long-tailed ducks and surf scoters occupied the maximum and minimum values among the four species included in our analysis. A previous assessment of sea duck breeding habitat suitability by Reed et al. (2017) grouped all three scoter species and did not include long-tailed ducks; our work builds on this analysis by suggesting additional differentiation between scoter species on several variables. We also observed between-species habitat partitioning during winter, with black and surf scoters preferring higher productivity, higher temperatures and flatter bottom slopes than the other three species.
Our analysis suggests a contrast in habitat selection and partitioning between pre-and post-breeding migrations. 13 Niche specialization and habitat partitioning were greater during post-breeding migration and molt than during prebreeding migration. The strong and consistent selection for preferred habitats we observed during post-breeding migration and molt reinforces previous work suggesting the importance of this period of the sea duck annual cycle (Lamb et al. 2019), during which sea ducks molt their flight feathers and are temporarily flightless (Salmonsen 1968). To complete the molt, birds must elevate their rates of energy consumption despite their limited mobility (Scott et al. 1994). For female sea ducks that have successfully raised young, post-breeding migration follows the energetically-demanding nesting period (Mallory 2015); therefore, they are also using staging and molt sites to replenish depleted resources to complete long-distance migration to wintering areas. Thus, access to habitats with abundant resources is particularly crucial during post-breeding migration and molt, restricting the range of suitable habitats available. In contrast, pre-breeding staging areas are used for shorter periods of time and do not include a flightless period (Johnsgard 2010).
Habitat use during migration is also subject to environmental and phenological constraints. During pre-breeding migration, sea ducks time their northward movement to the spring thaw in order to arrive at breeding sites as soon as possible after snowmelt (Takekawa et al. 2011); thus, the timing and spatial extent of sea ice breakup may limit access to Table 1. Model selection values for sea duck occurrence and species counts during a) breeding, b) wintering, c) pre-breeding staging and d) post-breeding staging and molt, 2002-2017. The best model for the full dataset was chosen using backward stepwise selection, and additional non-breeding habitat covariates (tidal current velocity, aquatic vegetation presence and substrate) were added to the top model for the subsets of locations for which they were available. Bold text denotes final models that were highly supported. aquatic habitats. Indeed, factors contributing to the likelihood of ice cover (sea surface temperature, salinity and tidal currents) were more important predictors of habitat selection during pre-breeding migration than during post-breeding and molt. Previous work on several species of eiders has demonstrated that sea ice plays a key role in structuring pre-breeding migration patterns (Mosbech et al. 2006, Petersen 2009, habitat use at stopovers (Solovieva 1999, Oppel et al. 2009) and prey composition and availability along migratory routes (Lovvorn et al. 2015). In addition, while pre-breeding staging sites are typically located along flyways between breeding and wintering sites, some segments of sea duck populations undertake additional migratory movements in order to molt at preferred sites (Salmonsen 1968, Lepage et al. 2020) and show high site fidelity to molt locations (Phillips and Powell 2006, Lepage et al. 2020). The relatively stronger habitat preferences and selection exhibited by all species during the post-breeding migration in our study is consistent with patterns of avian migration in multiple taxa, in which post-breeding migrants prioritize maximizing energy intake , while pre-breeding migrants prioritize speed and timing to ensure early arrival at breeding sites (Morris et al. 1994, La Sorte et al. 2016. While we were able to characterize some important aspects of habitat partitioning in our study, additional differences in habitat use may have occurred below the spatial (100 km 2 ) and temporal (seasonal averages) scales of our analysis (Holm and Burger 2002, Mohd-Azlan et al. 2014, Péron et al. 2019. Individual species may respond differently to changes in density-dependent competition, resource availability and environmental conditions at fine spatiotemporal scales, which would further allow them to partition resources in shared habitat areas. Additional study of the behavioral dynamics of mixed flocks of sea ducks could provide interesting insights into whether and how resource partitioning may occur even in areas where there is no apparent spatial partitioning of habitat. Importance of environmental covariates -our analysis included several variables that were available for only some portions of the study area: tidal current velocity, aquatic vegetation presence and bottom substrate. Tidal current velocity and aquatic vegetation presence were positive predictors of species counts during all seasons. High tidal velocities may be associated with high primary productivity, and some sea duck species have been observed to forage in or near tidal currents (Holm and Burger 2002). Previous studies have shown strong associations between waterfowl and seagrass, which is either frequent or incidental in the diet of several sea duck species and also serves as habitat for other aquatic prey (Kollars et al. 2017). Our results further suggest that aquatic vegetation may be a particularly critical food source during pre-breeding migratory staging. Bottom substrate was also a positive predictor of winter species counts, supporting previous studies suggesting that bottom sediments have an important influence on prey availability and foraging habitat selection for wintering sea ducks (Loring et al. 2013, Žydelis andRichman 2015). Developing these three data layers at a continental scale could help to improve modeling and prediction of sea duck occurrence.
Conservation recommendations -analyzing migration patterns in a multi-species framework can be used to identify key areas for conservation and restoration (Lamb et al. 2019), however, choosing sites within these areas and implementing habitat management activities requires a detailed understanding of habitat use and preferred habitat features. Areas of high multi-species suitability identified in our study could be useful targets for conservation and monitoring. Notably, our analysis of breeding habitat shows multi-species hotspots outside the range of traditional breeding surveys, suggesting that accurate assessment of population and habitat change may require expanding survey areas. The importance of climate characteristics to breeding habitat selection further suggests that these areas may be vulnerable to future changes in climate conditions, which additionally emphasizes the urgency of monitoring in these habitats. In non-breeding habitat areas, our results suggest that management of aquatic habitats for sea ducks likely depends on the specific timing of species occurrence. Given the strong selection we observed during post-breeding migration, habitats preferred by multiple species during this period -shallow, nearshore environments with high productivity -may be high-impact targets for conservation or restoration. Conversely, weaker selection at prebreeding migratory staging sites suggests that factors other than the measured covariates, such as varying sea ice cover or ephemeral resources, may be driving patterns of habitat use. Given that occurrence and intensity of use do not always reflect habitat quality (Winker et al. 1995), further studies of individuals occupying suitable and less-suitable sites could be used to establish whether conserving preferred habitats provides fitness benefits to sea ducks, and during which seasons habitat conservation would be maximally beneficial. Finally, while large-scale telemetry studies such as ours can provide important insights into annual-cycle conservation strategies, they are often too cost-and labor-intensive to be practical for monitoring long-term impacts of future environmental and habitat change or success of conservation measures. Thus, understanding how telemetry data can contribute to population models that incorporate other more readily-available datasets (Oppel andPowell 2010, Arnold et al. 2018) is crucial to developing cost-effective and feasible annual-cycle monitoring strategies for waterfowl.