The summer distribution, habitat associations and abundance of seabirds in the sub-polar frontal zone of the Northwest Atlantic

Biological production in the oceanic zone (i.e. waters beyond the continental shelves) is typically spatially patchy and strongly seasonal. In response, seabirds have adapted to move rapidly within and between ocean basins, making them important pelagic consumers. Studies in the Pacific, Southern and Indian Oceans have shown that seabirds are relatively abundant in major frontal systems, with species composition varying by water mass. In contrast, surprisingly little was known about seabird distribution in the oceanic North Atlantic until recent tracking showed that relative abundance and diversity peak in the Sub-polar Frontal Zone, west of the MidAtlantic Ridge, now proposed as a Marine Protected Area. However, absolute seabird abundance, distribution, age and species composition, and their potential environmental drivers in the oceanic temperate NW Atlantic remain largely unknown. Consequently, we systematically surveyed seabirds and environmental conditions across this area by ship in June 2017, then modelled the density of common species as functions of environmental covariates, validating model predictions against independent tracking data. Medium-sized petrels (99.8%), especially Great Shearwaters (Ardenna gravis, 63%), accounted for the majority of total avian biomass, which correlated at the macroscale with net primary production and peaked at the sub-polar front. At the mesoscale, the density of each species was associated with sea surface temperature, indicating zonation by water mass. Most species also exhibited scale-dependent associations with eddies and fronts. Approximately 51, 26, 23, 7 and 1 % of the currently estimated Atlantic populations of Cory’s Shearwaters (Calonectris borealis), Great Shearwaters, Sooty Shearwaters (A. grisea), Northern Fulmars (Fulmarus glacialis) and Leach’s Storm-petrels * Corresponding author. E-mail address: Ewan.Wakefield@glasgow.ac.uk (E.D. Wakefield).


Introduction
The oceanic zone (i.e. waters beyond the continental shelves) is the largest habitat on Earth. It remains relatively poorly understood but is undergoing increasingly rapid human exploitation (Crespo et al., 2018;St. John et al., 2016). For example, although seabirds are among its most conspicuous inhabitants, their abundance and relationships with oceanographic processes in the open ocean remain poorly known (Rodríguez et al., 2019). Reducing this uncertainty is important because seabirds are major consumers in pelagic ecosystems (Barrett et al., 2006;Brooke, 2004b;Hunt and McKinnell, 2006), recycle otherwise limiting nutrients (Savoca, 2018;Shatova et al., 2017), and are hyper-mobile (Dias et al., 2012;Edwards et al., 2013;Hedd et al., 2012;Kopp et al., 2011), linking disparate marine and terrestrial ecosystems with important biological and economic consequences (Plazas-Jiménez & Cianciaruso, 2020). Additionally, not only do seabirds consume a similar amount to human fisheries, they are also bycaught in those fisheries, leading to widespread and unsustainable population declines (Croxall et al., 2012;Cury et al., 2011;Dias et al., 2019;Grémillet et al., 2018). They are also likely to be vulnerable to the large-scale effects of climate change Rodríguez et al., 2019). Given the rapid changes the marine environment is undergoing, which are likely to intensify in the next decades, it is urgent to quantify seabird abundance in the open ocean to establish current baseline levels, at least in major hotspots. Moreover, in order to conserve seabirds effectively, it is necessary to understand relationships with environmental drivers and underlying processes driving seabird abundance and distribution (Grémillet & Boulinier, 2009).
In recent years, due to the burgeoning of tracking studies, interest in seabirds over the deep North Atlantic has renewed. Synthesis of tracking data from multiple populations has shown that at least 21 species, originating from breeding locations as far apart as the Arctic (Gilg et al., 2013) and Antarctica (Kopp et al., 2011), aggregate in a relatively small (0.5 million km 2 ) part of this area, west of the mid-Atlantic ridge (MAR), centred at ~37 • W, 50 • N (Davies et al., 2021). This includes breeding adults that routinely commute from northern hemisphere colonies to forage in the area during the boreal summer (Edwards et al., 2013;Paiva et al., 2010a); North Atlantic breeders making stopovers during migrations to and from southern hemisphere non-breeding areas (Egevang et al., 2010;Freeman et al., 2013); and adults from colonies in the northern (Fayet et al., 2017;Fort et al., 2012;Frederiksen et al., 2012) and southern hemisphere (Hedd et al., 2012;Kopp et al., 2011) that spend some or all of their respective non-breeding periods in the area. Seabird abundance in the area is relatively high year-round (Davies et al., 2021). In the boreal winter, the avifauna is dominated by alcids, especially the Little Auk (Alle alle), and Black-legged Kittiwakes (Rissa tridactyla) Fort et al., 2012;Frederiksen et al., 2012). In in the boreal summer, diversity is higher and medium-sized petrels, especially Great Shearwaters (Ardenna gravis), predominate (Davies et al., 2021). It has also been inferred from biologging data that many species undergo moult in the area (Hedd et al., 2012;Kopp et al., 2011), which is a vital maintenance activity for birds (Ellis & Gabrielsen, 2002). It has been recognised that the area of highest seabird concentration coincides with the sub-polar frontal zone (hereafter, SPFZ) (Boertmann, 2011;Skov et al., 1994;Davies et al., 2021), which in the NW Atlantic is an area of particularly dynamic physical oceanography (see Section 2). The SPFZs of the Pacific and Southern Oceans also sustain high seabird diversity and abundance and seabird-habitat associations in these systems have been studied extensively in situ (Ainley & Boekelheide, 1984;Hyrenbach et al., 2007;Pakhomov & McQuaid, 1996;Wahl et al., 1989). In contrast, at-sea studies of seabirds in the NW Atlantic SPFZ have been limited in extent and restricted mainly to descriptions of occurrence and relative density (Bennison & Jessopp, 2015;Boertmann, 2011;Brown, 1986). Only one previous ship-based study described seabird habitat use, and this covered only a small part of the SPFZ (Skov et al., 1994). Furthermore, despite numerous recent tracking studies showing that seabirds use the deep NW Atlantic, tracking has only been used to quantify habitat use by alcids, which predominantly occur in the area in the winter, and Cory's Shearwaters (Calonectris borealis) (Fort et al., 2012;Merkel et al., 2021;Paiva et al., 2010a;Paiva et al., 2010b;Tranquilla et al., 2015). These studies suggest that seabird distribution in NW Atlantic may vary systematically with water mass. However, relationships between most species and fronts, mesoscale turbulence, primary production, etc., especially in the summer, when medium sized petrels predominate, remain essentially unknown, hindering our understanding of the drivers of high seabird abundance in the area.
The core region of high seabird abundance and diversity identified by Davies et al. (2021) is currently being considered by the OSPAR Commission as a potential Marine Protected Area (the North Atlantic Current and Evlanov Seamount Marine Protected Areahereafter NACES pMPA, Fig. 1). However, until recently, it was largely impracticable to track immature seabirds, so the data used to identify this area pertain only to adults. Approximately half of pelagic seabirds are immatures (Brooke, 2004a;Carneiro et al., 2020), which can have markedly different distributions to adults Carneiro et al., 2020;Powers et al., 2020). Similarly, tracking data are lacking for many populations that potentially use the area (e.g. Northern Fulmars Fulmarus glacialis from NW Atlantic colonies), including those of most smaller species, such as the storm-petrels (Hydrobatidae/Oceanitidae). Ship-based surveys, combined with onboard oceanographic sampling, satellite remote-sensing and habitat modelling, can be used to both quantify seabird-habitat associations and estimate their abundance (Buckland et al., 2016;Miller et al., 2013a;Waggitt et al., 2020), validating estimates inferred indirectly from tracking data (Carroll et al., 2019;Louzao et al., 2009). For some species, at sea surveys also allow age and moult status to be observed or inferred, providing information that is difficult to collect remotely (Brown, 1988;Keijl, 2011;Meier et al., 2017).
Here, we report the findings of the first widescale survey of seabirds in the SPFZ of the NW Atlantic, carried out in June 2017. Our primary aims were to (1) quantify the habitat associations of the commonest species in the area (Northern Fulmars, Cory's Shearwaters, Great Shearwaters, Sooty Shearwaters (A. grisea) and Leach's Storm-petrels (Oceanodroma leucorhoa) -hereafter NOFU, COSH, GRSH, SOSH and LHSP, respectively), and (2) estimate their distribution and abundance (individuals and biomass), both in the area as a whole and within the NACES pMPA. In addition, we report the moult status and age of seabirds in the area and the occurrence of less abundant species.

Physical and biogeographic setting
We defined our study area ( Fig. 1; 43 • 18ʹ to 53 • 12ʹN, 29 • 0ʹ to 42 • 6ʹW; extent 1,174,800 km 2 ) to encompass the seabird diversity and abundance hotspot identified by Davies et al. (2021). This area occupies much of the Newfoundland Basin and is bounded to the east by the MAR, to the north by the Charlie-Gibbs Fracture Zone (CGFZ -a large discontinuity in the MAR/Reykjanes Ridge) and to the west by the North American continental rise. Its mean depth is 3814 m (range 907-5031 m) but many shallower seamounts occur, mainly towards the northern, eastern and southwestern margins. The region's physical oceanography was reviewed by Rossby (Rossby, 1996). In brief, surface transport is dominated by the warm North Atlantic Current (NAC), which branches from the Gulf Stream SE off the Grand Banks. From there, the NAC follows the continental rise around the Flemish Cap, before retroflecting sharply north-westwards then eastwards (thus forming the "Northwest Corner"), passing through the CGFZ, and continuing NE beyond the study area. The cold, southward flowing Labrador Current meets the NAC south of the Flemish Cap, returning north-eastwards with the latter. In the upper ocean (0-1000 m), the SPFZ is characterised by a stepped transition from warm/saline subtropical North Atlantic Central Water (NACW; SST > 15 • C, S > 35.5) in the south to colder/fresher Sub-Arctic Intermediate Water (SAIW; SST < 10 • C, S < 35.0) in the north (definitions based on Cook et al. (2013), Søiland et al. (2008) and observations made during our study). West of the MAR, the NAC forms at least three approximately parallel, zonally aligned jets and associated fronts (Belkin & Levitus, 1996;Bower & von Appen, 2008;Miller et al., 2013b;Søiland et al., 2008). The subpolar front has two components -the North Subpolar Front (NSPF) and South Subpolar Front (SSPF). South of this is the less clearly defined Mid-Atlantic Front (MAF) (Belkin & Levitus, 1996;Miller et al., 2013b). The NSPF typically crosses the MAR at the CGFZ but the other fronts are more variably associated with fracture zones in the MAR (Bower & von Appen, 2008). Mesoscale turbulence within the SPFZ is high (Chelton et al., 2011;Miller et al., 2013b;Tilstone et al., 2014) and upper water masses between the MAF and subpolar fronts comprise a matrix of NACW, SAIW, and intermediate water, which together we refer to as Frontal Water (FW) (Cook et al., 2013;Søiland et al., 2008).
Biogeographically, the area is characterised by high summer primary production associated with lateral stirring across the SPF (Longhurst, 1998;Tilstone et al., 2014) and a complex mosaic of pelagic habitats . The SPFZ acts as both a biogeographical barrier separating communities associated with SAIW (relatively low diversity/high abundance) and NACW (relatively high diversity/low abundance), and an ecotone with higher abundance at all trophic levels and its own distinct faunal assemblage (Beaugrand et al., 2002;Vecchione et al., 2010a;Vecchione et al., 2015). Two major ecoregions have been recognised Sutton et al., 2017): To the north, the Sub-Polar Oceanic (Beaugrand et al's terminology), with relatively low zooplankton diversity but high abundance of the copepod Calanus finmarchicus and Hyperiid amphipods. To the south, the Oceanic Warm-Temperate has a more diverse copepod and diatom community and lower seasonal variation in abundance. There are few data on mesotrophic fauna in the area itself but extensive investigations during and surface currents (red/warm, blue/cold). Dashed red lines indicate transient, eastward flowing branches of the North Atlantic Current (NAC) and light blue shading the subpolar frontal zone. Abbreviations are AZ, Azores; CGFZ, Charlie-Gibbs Fracture Zone; GS, Gulf Stream;LC, Labrador Current;MAR, Mid-Atlantic Ridge;NB, Newfoundland Basin;NWC, Northwest Corner;andRR, Reykjanes Ridge. Isobaths indicate 50, 100, 200, 500, 1000, 2000, 3000, 4000 and 5000 m depth. Based on Sy (1988), Rossby (1996), Belkin and Levitus (1996) and Søiland et al. (2008)  the MAR-ECO project (Priede et al., 2013;Vecchione et al., 2010a), conducted over the MAR immediately to the east suggest that Myctophids dominate above 750 m, with shrimps and cephalopods also relatively abundant (Cook et al., 2013;Sutton et al., 2013;Vecchione et al., 2010b).

Methods
To estimate seabird distribution and abundance and investigate habitat associations within the study area, we surveyed birds by ship using standard methods (Camphuysen et al., 2004;Tasker et al., 1984;Webb & Durinck, 1992), while simultaneously measuring sea surface properties in situ and via satellite. We then modelled seabird density as a function of environmental covariates to infer habitat associations and predict seabird distribution and abundance across the study area. To validate these models against independent observations, we also tracked representative samples of four species (NOFU, COSH, GRSH and SOSH) that were abundant during the at-sea survey.

Ship-based seabird survey
We surveyed seabirds from the RRS Discovery (cruise DY080) between the 13th and 29th of June 2017. A sawtooth survey pattern, comprising five transects running perpendicular to the major fronts within the study area, was planned but this was modified during the cruise in order to sample a transient phytoplankton bloom in the NW of the study area ( Fig. 1). We used the Eastern Canada Seabirds At Sea (ECSAS) protocol (Gjerdrum et al., 2012). This follows de facto standard methods for surveying seabirds from ships (Camphuysen et al., 2004;Tasker et al., 1984;Webb & Durinck, 1992), comprising a simultaneous line transect survey for birds on the water, a strip transect survey for birds in flight (transect width of 300 m) and periodic recording of weather conditions that may affect seabird detectability (see Supplementary Methods). Movement of flying birds relative to the ship can bias density estimates (Gaston et al., 1987). To minimise this effect, we used a 'snapshot' method to record birds first detected in flight, flagging records of flying birds if they were within a 300 × 300 m box at the moment of the snapshot (Tasker et al., 1984). We used only records of birds in flight at snapshots (hereafter, in-flight), plus all sightings of birds on the water (hereafter, on-water) to estimate density. Together, we refer to these as 'in-transect' sightings. We noted if birds were obviously following the ship and excluded these records from the density analysis. In order to obtain an approximate estimate of the moult status of birds in the study area, whenever possible we examined birds using binoculars and recorded whether they were clearly in active primary moult (i.e., one or more primaries missing or not fully grown). We recorded the colour phase of NOFU on the four point scale described by van Franeker and Wattel (1982), where birds with the lightest and darkest plumage are termed LL and DD, respectively, and lighter and darker intermediates are termed L and D, respectively. We also recorded the approximate age (calendar year) of gannets, gulls, terns and skuas when this could be discriminated from plumage.

Seabird tracking
In order to obtain an independent estimate of the distribution of the most abundant seabird species in the study area during the ship-based survey, we used light-based geolocation loggers to track these species from their breeding colonies. We then compared the distributions predicted by models fitted to the at-sea survey data to these data (see Section 5.1 for caveats). It was only practicable to track birds from one colony per species. We selected these colonies based on prior understanding of which populations were likely to use the study area during the survey (Edwards et al., 2013;Hedd et al., 2012;Magalhaes et al., 2008;Paiva et al., 2010a,b;Davies et al., 2021) and logistical constraints. In addition, to account for availability bias (see Section 3.7), we deployed Time-Depth Recorders (TDRs) on a subsample of the tracked shearwaters. We tracked NOFU from Eynhallow, Scotland (59 • 8 ′ N, 3 • 8 ′ W); GRSH from Gough, South Atlantic (40 • 19 ′ S, 9 • 56 ′ W); and SOSH from Kidney Island, Falkland Islands (51 • 38 ′ S, 57 • 45 ′ W) for ≥ one year overlapping with the survey period. For logistical reasons, it was not possible to deploy loggers on COSH prior to the survey. Rather, we tracked this species from Corvo, Azores (39 • 41 ′ N, 31 • 7 ′ W) for one year, commencing immediately after the survey. It was not known that LHSP are one of the most abundant species in the study area until the survey was carried out, so we did not track this species. Loggers were deployed and recovered on NOFU as described by Grissot et al. (2020) and on shearwaters as described by Bonnet-Lebrun et al. (2020) (Supplementary Methods for details).

Modelling detection probability
We used model-based distance sampling (Miller et al., 2013a) to estimate the distribution and abundance of the most abundant species (i. e. those with ≥ 60 detections, which is the minimum required to estimate a robust detection function (Buckland et al., 2001)). We assumed that the probability of detecting birds on the water declines with distance from the transect line and with sea and weather conditions. We used the R package Distance (Miller et al., 2019) to model detectability and estimate the probability of detection p w , considering various standard functional forms (Marques et al., 2007), as well as the following covariates -wind speed, wave height, visibility (all continuous) and precipitation (binary). We recorded distances in one of four bands, so we limited the maximum number of parameters in detection functions to three. We therefore considered models with up to two covariates using half-normal and hazard-rate detection functions, plus models with no covariates but including adjustment terms to improve fit, selecting among these models based on their AIC. As distances occurred in four bins, it was not possible to obtain χ 2 goodness-of-fit statistics, so we assessed model fit visually by plotting the observed frequency of detections vs. distance, overlaid with the fitted detection functions. There were insufficient detections of SOSH and COSH to fit robust detection models for these species separately, but SOSH, COSH and GRSH have a similar, predominantly dark, appearance on the water so we pooled detections of these species and fitted a generic detection function for all shearwaters.
We assumed that the probability p f of detecting the larger species (NOFU and the shearwaters) in flight was uniform within the surveyors' search area (Waggitt et al., 2020). LHSP are, in contrast, much smaller and usually fly closer to the water than these species, making them more difficult to detect. Moreover, they forage predominantly by foot pattering or otherwise remaining in partial contact with the sea (Flood & Fisher, 2011), so it was often ambiguous whether detections should be classed as in flight or on the surface. However, this behaviour makes it practicable to estimate the distance of LHSP from the transect line using the Heinemann (1981) technique. We therefore fitted a detection function to pooled in-flight and on-water detections for this species.

Explanatory environmental covariates
We modelled the density of seabirds (see Section 3.5) as functions of three types of candidate explanatory covariates (Table 1): (1) Accessibility from breeding colonies; (2) indices describing environmental phenomena or properties that could affect prey distribution; and (3) spatial location (see Section 3.7). These covariates were either measured shipboard during the survey (hereafter, locally-sensed), and therefore relatively high-resolution but available only along the cruise track, or remotely-sensed at lower resolution across the whole study area.
We considered the following indices: Accessibility (α), which for species breeding at the time of the survey (NOFU, COSH and LHSP), we define as the inverse of the great circle distance to their nearest breeding colony (Wakefield et al., 2011). For non-breeding species (GRSH and SOSH), we assumed that α was approximately uniform across the study area. Sea surface temperature (SST): Different water masses are characterised by different lower trophic level assemblages (Longhurst, 1998) so the distribution of seabirds is hypothesised to vary with water mass and therefore SST (Ainley & Boekelheide, 1984;Ballance et al., 2006;Hunt, 1997;Hyrenbach et al., 2007;Pakhomov & McQuaid, 1996). We obtained gridded Advanced Very High Resolution Radiometer SST data (spatial resolution, 1 km) from the NERC Earth Observation Data Acquisition and Analysis Service (NEODAAS). The study area was frequently obscured by clouds during the cruise, so to avoid missing data we averaged SST over a 28 day window, centred on the study period. We measured the near-surface SST (hereafter, nSST) throughout the cruise (±0.001 • C) using a Sea-Bird SBE 38 Digital Oceanographic Thermometer pump-supplied from an inlet at 5.5 m depth. We also measured salinity (±0.005) using a Sea-Bird SBE45 MicroTSG Thermosalinograph fed by the same supply. However, this was highly correlated with nSST (Spearman's ρ = 0.95, p < 0.001), so we do not consider salinity further in our analysis. Fronts: Seabirds and their prey are often more abundant in the vicinity of thermohaline fronts, either because elevated nutrient supply at fronts enhances primary production or because convergent currents associated with fronts cause biota to aggregate (Bost et al., 2009;Scales et al., 2014). Following Miller et al. (2015), we quantified the presence and intensity of fronts as a continuous metric, front gradient (FG), defined as the mean front gradient magnitude, spatially smoothed to produce a continuous distribution from discrete contours. In brief, we detected fronts in daily-merged microwave and infrared SST maps (9 km spatial resolution (Gentemann et al., 2009)) using a local histogram algorithm (Cayula & Cornillon, 1992) with a minimum SST step of 0.45 • C, and combined the fronts over 3 days . We then smoothed the resulting composite front map using a Gaussian filter (σ = 3 pixels) to convert discrete front contours into a continuous metric. In addition, to identify frontal regions from locally-measured data, we calculated the smoothed gradient of nSST (ΔnSST) by predicting the first derivative of a local regression model (Loader, 1999) of nSST vs. distance along the transect, specifying a bandwidth of 250 km. Eddy kinetic energy (EKE): Ocean currents, like the NAC, give rise to turbulence in the form of mesoscale eddies, jets and waves, which may enhance or suppress primary production through the transport of nutrients into or out of the mixed layer or aggregate biota through convergent advection (Falkowski et al., 1991;Gruber et al., 2011). The intensity of turbulence may therefore affect seabird prey availability (Bertrand et al., 2014;Tew Kai & Marsac, 2010). We quantified mesoscale turbulence using EKE, where and u a and v a are the zonal and meridional geostrophic current anomalies. Sea level anomaly (SLA): Mesoscale eddies/meanders may affect seabird distributions by enhancing or suppressing primary production; containing prey assemblages distinct from surrounding waters; or concentrating planktonic organisms near the surface, especially at eddy margins (McGillicuddy, 2016;Tew Kai & Marsac, 2010). We used daily the SLA to indicate the presence of either cold-(negative MSLA) or warm-core (positive SLA) eddies. We obtained both daily SLA and geostrophic current data from the Global Ocean Gridded L4 Sea Surface Heights and Derived Variables dataset (spatial resolution 0.25 • , ~28 km), downloaded from http://marine.copernicus. eu/ (accessed October 2, 2019). Net primary production (NPP). We assume that seabird prey abundance increases with NPP via trophic cascades (Shaffer et al., 2006;Wakefield et al., 2014). We downloaded 8 day and monthly NPP data (spatial resolution 1/12 • , ~9 km), estimated using a carbon-based production model (Behrenfeld et al., 2005), from http://sites.science.oregonstate.edu/ocean.productivity/index.php (accessed October 7, 2019). To improve their spread, we logtransformed NPP and square root transformed EKE and FG prior to model fitting.

Density modelling
We had two aims in modelling seabird density: (1) to investigate associations between seabirds and environmental conditions and (2) to predict seabird distribution and abundance across the study area. To meet aim 1, we first modelled the density of each species as a function of remotely-sensed environmental indices, plus, for species breeding at the time of the survey, accessibility (hereafter, global models). In addition, because SST and NPP estimated via remote sensing are prone to signal degradation due to clouds (Becker et al., 2010), we fitted a separate density model in which explanatory covariates comprised locally-sensed nSST and ΔnSST, remotely-sensed indices unaffected by clouds (EKE and SLA), and (when relevant) accessibility (hereafter, the local model). Previous habitat modelling studies suggest that seabird distribution is often poorly explained by environmental covariates alone (Block et al., 2011;e.g., Louzao et al., 2009;Wakefield et al., 2017), for example, because environmental indices are poor proxies for prey or because seabirds are imperfectly informed about the distribution of their prey (Grémillet et al., 2008). For the purposes of predicting distribution and abundance (aim 2), we therefore augmented the global model for each species by adding smooth of spatial location (see below), referring to the result as the spatial smooth model.
To model bird density, we divided transects into 12 min segments (equivalent to ~ 3.4 km at a speed of 17 km/h), indexed by j, this length being chosen to reduce serial autocorrelation whilst retaining sufficient spatial resolution to model habitat relationships (Huettmann & Diamond, 2006). Recalling that for NOFU and shearwaters p f ∕ =p w , we specified two simultaneous sets of segments for these taxa, indexed by s, Table 1 Potential environmental explanatory covariates considered during selection of models of seabird density. Covariates available throughout the study area were considered in global models, while covariates unaffected by cloud were considered in local inference models. one containing in-flight detections (s = 0) and the other on-water detections (s = 1) (Miller et al., in review). We then modelled counts per segment n j,s as Generalised Linear Models (GLMs) or Generalised Additive Models (GAMs) with the form where β 0 is the intercept, β state estimates the mean difference between counts of birds on the water and in flight, and f m are either linear or quadratic functions of the environmental covariates x j , plus, in the case of the spatial smooth model, a smooth of segment location (see Section 3.7). In effect, Eq. (1) treats the data as arising from two surveys conducted simultaneously (Miller et al., in review). While the other terms in the model assume that the spatial pattern is the same for birds in flight and on the water, the β state term allows the mean density of birds in these behavioural states to differ. The term p s,j A j , which enters the model as an offset, is the effective area of the segment (sensu Miller et al., 2013a), where A j = wl j , w is the transect width (300 m) and l j the segment length. Based on exploratory analysis, we assumed that n j,s conformed to a negative binomial distribution for all species. We fitted models using the R package dsm (Miller et al., 2021). We took a similar approach for LHSP but specified a single set of segments, containing both birds detected in flight and on the water, and modified equation (1) by removing β state .

Model selection and validation
Selecting species distribution models from a set of candidates that includes very complex or ecologically unrealistic models can result final models that predict poorly in unsampled regions (Bell & Schlaepfer, 2016). We used the following strategies to avoid this. Firstly, we defined candidate models based on prior understanding of how environmental phenomena (measured by available proxies, such as SST), might affect seabird distribution. Secondly, we assumed that, on the scale of the linear predictor, associations between seabird density and environmental covariates could be adequately approximated by linear or second order polynomial terms. Thirdly, we reduced model complexity by backwards selection based on spatial cross-validation, assuming that a good model should predict accurately in areas of space not included in the training data set (Roberts et al., 2017;Wenger & Olden, 2012). Finally, we compared predicted density to patterns of spatial usage observed independently via tracking.
Starting with the terms shown in Table 1, we built the global and local inference models by backwards selection based on cross-validation in a similar manner to Roberts et al. (2017). We divided the study area into nine blocks (Fig. A1), this number being chosen to provide units that both represented the range of environmental conditions across the study area and held similar amounts of survey data. We then fitted the model under consideration to data from eight of the blocks, predicted for the remaining block and calculated the logarithmic score of these predictions. The logarithmic score is the negative log of the probability of obtaining a given count, and has been advocated for the selection of count models due to its simplicity and propriety (Czado et al., 2009). We repeated this process for all blocks and then calculated the mean logarithmic score S, across spatial blocks as: where blocks are indexed by i and there are m i observations within each block, indexed by j. P − i is the probability mass function for the model fitted to all but block i. S decreases as the predicted counts approach their observed counterparts so, during model selection, we accepted a potential model simplification if it resulted in no increase in S. For comparison, we also fitted a model for each species containing only the intercept and the state covariate (hereafter, state-only models). We examined correlations among explanatory covariates using correlation matrices (Dormann et al., 2013) but following Morrissey and Ruxton (2018), we did not automatically exclude candidate explanatory covariates simply because they were correlated with others. Starting with all explanatory covariates structured as quadratic terms (Table 1), we simplified models by sequentially removing second then first order polynomials in order of their significance. We assessed model fit and conformity to assumptions using quantile-quantile and residual plots and residual serial autocorrelation using correlograms (Wood, 2017). To assess residual spatial correlation, we fitted GAMs with a bivariate smooth of location to each model's residuals, reasoning that residuals randomly distributed in space would result in a smooth with zero effective degrees of freedom.

Predicting distribution and abundance
For the purposes of predicting distribution and abundance across the study area, we refitted the best global inference model for each species with an additional bivariate Duchon spline smooth of segment location with a maximum basis dimension of 30. Using these models, we predicted the abundance of each species (birds on the water, plus birds in flight) across a Lambert Azimuthal Equal Area grid (cell size 4 × 4 km) encompassing, and centred on, the study area. COSH, GRSH and SOSH all dive frequently Paiva et al., 2010a;Ronconi et al., 2010b). In order to account for the resulting availability bias (Buckland et al., 2015) for these species, we multiplied their predicted abundance by the inverse of the mean proportion of time they spend at the surface during daylight hours (Winiarski et al., 2014), estimated using the TDR data (Supplementary Methods). Although Leach's petrels and NOFU are capable of diving, they do so infrequently Ortega-Jiménez et al., 2009), so following previous studies (e.g. Waggitt et al., 2020), we assumed that availability bias is negligible for these species. We estimated uncertainty in predicted abundance by posterior simulation, using a technique adapted from Wood (2017; Section 7.2.6). In brief, we randomly drew model parameters from their estimated multivariate normal distribution, propagating uncertainty due to both the detection function and count model following Bravington et al. (2021). We then predicted the abundance in each cell and across the study area, if relevant, applying an availability bias correction randomly drawn from its posterior distribution. We repeated this process 1000 times and then calculated the celllevel and overall means and their corresponding coefficients of variation and 95% confidence intervals. To compare our results to Davies et al. (2021)'s tracking-based estimates of adult abundance, we also used this method to predict abundance in the NACES pMPA, multiplying this by the assumed proportion of adults in the study area. For NOFU and LHSP, the proportion of adults was inferred from our moult observations. For the remaining species, it was assumed, based on published estimates (Brooke, 2004a;Carneiro et al., 2020), to be ~50%. We converted observed and predicted abundances to biomasses using estimates of mean body mass reported by Brooke (2004a).

Geolocator processing and comparison with model predictions
We used the probGLS package to estimate two locations/bird/day from geolocator data, modifying the methods described by Merkel et al. (2016) for trans-equatorial migrants (see Supplementary Methods). Median location errors using this method are ~185 and 145 km during the equinoxes and solstices, respectively . We compared the model-predicted spatial distribution of each species to the distribution of birds in the study area estimated from tracking data following Carroll et al. (2019). In brief, we estimated the utilisation distribution (UD) of each tracked bird using kernel density analysis, implemented in the adehabitatHR package (Calenge, 2006), specifying a fixed smoothing parameter of 75 km and the grid described in Section 3.7. We then averaged UDs across birds. In order to provide sufficient data to resolve distribution patterns, we calculated UDs using bird locations recorded during the seabird survey period ± 20 days. COSH UDs were necessarily estimated using data from the equivalent period in 2018 (see Section 3.2). We quantified the similarity between the modelpredicted distributions and the tracking-based UDs by first cropping the latter to the study area and then normalising each to sum to unity. We then calculated the Bhattacharya affinity, which ranges from 0 (no co-

Fig. 2.
Environmental conditions in the study area during the survey (see Table 1 for definitions). The solid line indicates the cruise track and the dashed lines the 10.5, 13 and 14.75 • C SST isotherms delineating the approximate locations of the North Subpolar Front (NSPF), South Subpolar Front (SSPF) and Mid-Atlantic Front (MAF), respectively. The yellow triangle indicates the centre of the cold-core eddy mentioned in the text. Geostrophic current anomalies are superimposed on SLA. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) occurrence) to 1 (identical distributions), between these two matrices (Fieberg & Kochanny, 2005).

Species associations
To establish which species typically co-occurred with one another at fine scales (~4 km), we calculated the Chao-Jaccard similarity index (Chao et al., 2005) between per-species segment counts and linked similar species using Ward clustering (Kauffman & Rousseeuw, 2005). We restricted this analysis to species detected in at least 20 segments.

Survey effort and environmental conditions
The seabird survey covered 3265 km in 192 h. Beaufort wind force averaged 4.5 ± 1.6 and wave height 1.7 ± 0.8 m during survey bouts. Visibility averaged 9 ± 7 km but fell as low as 300 m at times in the west of the study area due to fog. SST ranged from ~18 • C in the south of the study area to 8 • C in the north (Fig. 2) and the mean latitudes of the MAF, SSPF and NSPF were 46 • 12 ′ , 49 • 39 ′ and 52 • 15 ′ N, respectively (Fig. 3). Multiple mesoscale eddies were crossed, including an intense cold-core ring (diameter ~ 140 km) centred at 46 • 5 ′ N, 40 • 30 ′ W, ~340 km ESE of the Flemish Cap (Fig. 2). NPP was markedly higher north of the SSPF than to the south (means 815 vs. 590 mg C/m 2 /day; Fig. 3) but isolated patches of high NPP also occurred south of this, associated with negative SLAs, including the eddy just mentioned (Fig. 2).

Observed distribution and characteristics of seabirds
We recorded eighteen species of seabird during the survey (Table 2), totalling 7464 individuals. Of these, 4692 were sighted in transect, the vast majority (96.5 %) being medium-sized petrels, all of which were identified to species. Most abundant were GRSH (64.6 % of in-transect sightings), followed by NOFU (22.0 %), COSH (7.1 %), SOSH (2.4 %) and Manx shearwaters (Puffinus puffinus, 0.4 %). LHSP were also detected relatively frequently (2.8 %) but a small proportion of Hydrobatidae/Oceanitidae sightings (3.5 % of 141 birds) were not identifiable to species. Based on size, the majority were suspected to be Wilson's storm-petrels so these records were not included among the data used to model LHSP density.
At the meso-to macroscale (100s-1000s km), mean avian biomass increased with latitude and NPP and decreased with nSST and salinity, peaking between the NSPF and SSPF (Fig. 3). At the mesoscale (10 s-100 s of km), two notable seabird aggregations occurred -one in the NE of the study area, between the SSPF and NSPF (cf. Fig. A2 and Fig. 2) and another smaller one 350 km east of the Flemish Cap, associated with the cold-core eddy described above. Cluster analysis indicated that at coarse scales (1-10 km) and above, the more common species tended to co-occur in two groups (Fig. 4), the first comprising GRSH, NOFU and SOSH, which were most abundant north of the SSPF (Fig. 3), and the second COSH and LHSP, which were mostly confined to the south. Arctic terns (Sterna paradisaea) and Manx shearwaters, occurred throughout the area. Taxa occurring in smaller numbers included all three jaeger spp. (Stercorarius longicaudus, S. parasiticus and S. pomarinus), recorded throughout the study area, and Catharacta skuas, recorded almost exclusively south of the SSPF (Fig. A2). All of the Catharacta spp. (7 of 19 birds) that could be positively identified (following Newell et al., 2013) were South Polar Skuas (C. maccormicki).
Of those birds whose moult could be assessed, 59 % of 334 NOFU, 3 % of 237 COSH, 70 % of 894 GRSH and 50 % of 10 SOSH were in active primary moult. The likelihood of primary moult among NOFU increased significantly with distance from the nearest colony ( Fig. A3; binomial GLM Z 1, 332 = 3.15, p = 0.002). Of the 1201 NOFU for which colour phase was assessed, 94 % were LL, 4 % L, and 2 % D. Of those for which moult was also assessed (n = 307), a significantly higher proportion (25 Fig. 3. Surface hydrography and observed biomass (corrected for imperfect detection) of the five most abundant seabird species averaged by latitude over the survey track. Near-sea surface temperature (nSST) and salinity were measured ship-board at 5 m depth and net primary production (NPP) was estimated from remotelysensed data (Behrenfeld et al., 2005). ρ indicates Pearson's correlation between these indices and bird biomass. nSST and salinity gradients (dotted lines) were calculated prior to averaging by latitude (see Section 3.4). Dashed vertical lines indicate the approximate mean latitudes of the major fronts defined in Fig. 2 and arrows the mean extents of North Atlantic Central Water (NACW), Frontal Water (FW) and Sub-Arctic Intermediate Water (SAIW).

Detection functions
Slightly more NOFU were detected in the 50 -100 m distance band, than in the 0-50 m band. A half-normal model, with scale parameter covariates visibility (continuous, km) and precipitation (true/false), was selected based on AIC (Table A4, Fig. A5). Precipitation (especially fog) may have affected the detectability of NOFU more than other species due to their relatively light colouration. However, compared to the other taxa, the detectability of NOFU (mean probability of detection p w = 0.51) varied relatively little within the range of environmental conditions experienced during the survey (Fig. A5). A half-normal model with scale parameter covariates for visibility and Beaufort wind speed best described the probability of detecting GRSH, SOSH and COSH on the water (p w = 0.48). A half-normal model, with a scale parameter varying with wind speed, best described the probability of detecting LHSP onwater and in-flight (p w = 0.51). Variation in this probability was relatively large over the range of wind speeds experienced during the survey (Fig. A5), presumably because LHSP are harder to detect as sea conditions deteriorate than the other species considered.

Habitat associations
Correlations among candidate explanatory environmental covariates were < 0.7 (Table A5), with the exception that SST was positively correlated with accessibility for COSH (ρ = 0.83, p < 0.001). However, accessibility was rejected during model selection (Fig. 5, Fig. A6). There was no consistent pattern in whether the global or local inference models best described density for each species (Table 3).
For NOFU, both the global and local inference models indicated similar environmental associations (Fig. 5). Density initially rose with SSTs and nSSTs, peaked around 10-12 • C and declined sharply above ~ 13 • C, indicating that this species associates principally with SAIW and FW. In addition, NOFU were positively associated with fronts and negatively associated with EKE. The local model, which had the best predictive performance (Table 3) and lowest residual autocorrelation (Fig. A7), additionally suggested that NOFU density increased with SLA. The global and local inference models for COSH had very similar structures, but the latter predicted more accurately and exhibited lower residual autocorrelation. Both indicated a positive association between COSH density and SST (Fig. 5), the latter indicating highest densities in NACW. In addition, density increased with absolute SLA, with highest COSH densities occurring coinciding with negative SLAs. The local model also indicated a weak negative association with EKE. Both GRSH inference models performed poorly compared either to the GRSH stateonly model or the inference models for the other species (Table 3). Moreover, the global and local GRSH inference models both exhibited  (Bennison and Jessopp, 2015;Boertmann, 2011;Godø, 2004;Priede, 2007;Skov et al., 1994;Davies et al., 2021). relatively high residual autocorrelation (Fig. A7). The former indicated a positive association between GRSH and FG, whereas the latter implied a negative association with nSST (Fig. 5), with density being highest in SAIW/FW. For SOSH, the global inference model performed better than the local model (Table 3, Fig. A7). Both indicated an association with cooler temperatures (SAIW/FW) but the remaining covariates were inconsistent between models: The global model additionally indicated a positive association with FG and a negative association with NPP, while the local model indicated a positive association with SLA. The local and global models for LHSP were very similar to one another but the latter had the better predictive performance (Table 3, Fig. A7). Both indicated a positive association with SST and especially NACW (Fig. 5).

Predicted distribution and abundance
For all species, inclusion of a spatial smooth reduced residual serial and spatial autocorrelation markedly (Fig. A7). The mean proportion of time TDR-equipped shearwaters spent diving during daytime varied seasonally and between species (Fig. A8). During the survey period, COSH, GRSH and SOSH spent on average 0.13 (95% CI 0.08-0.23), 0.24 (0.07-0.81) and 1.50 (1.01-2.23) % of their time diving, respectively, so availability bias had little effect on our abundance estimates (Fig. A9).
GRSH were predicted to be the most abundant species in the study area (~3.9 million individuals, Table 4) and their predicted density highest along its SW/NE axis, peaking in the NE (Fig. 6). Eighteen out of Fig. 5. Partial associations between bird density and environmental indices, predicted by the global and local inference models (see Fig. A6 for associations on the scale of the linear predictor). Global models contain only remotely-sensed environmental covariates, whereas local models contain remotely-sensed covariates unaffected by cloud and locally-sensed environmental covariates (see Table 1 for definitions). Dashed vertical lines show the approximate latitudes of the major fronts defined in Fig. 2. 23 wintering GRSH tracked from Gough Island used the study area during the survey. Their distribution was very similar to that predicted from the at-sea data (BA = 0.92; cf. columns 1 and 3, Fig. 6), their usage being concentrated in the NE and SW. The latter area coincided with the cold core eddy highlighted in Fig. 2. NOFU were predicted to be the second most abundant species (~1.3 million individuals), with higher densities in the northern half of the area (Fig. 6). Only three of 29 breeding NOFU GLS-tracked from northern Scotland during the survey period used the study area, where they were confined to the NE corner (see also Fig. 7). Similarity between the predicted distribution of NOFU and the distribution of tracked birds was therefore relatively low (BA = 0.36). COSH were predicted to be the third most abundant species (~0.5 million individuals), their density increasing towards the south and, to a lesser extent, east of the study area. Predicted density was particularly high in the SW, associated with warm subtropical Gulf Stream waters. Around half (12 of 23) of the breeding COSH tracked from the NW Azores in the May-July 2017 and 2018 used the study area. Their distribution was similar to the model predictions (BA = 0.68) except that their usage was more concentrated in the SE, near to their breeding colony. LHSP (~0.19 million individuals), were predicted to occur predominantly south of the SSPF, with a notable concentration on the southern margin of the cold core eddy east of the Flemish Cap. SOSH were predicted to be the least abundant of the modelled species (~0.15 million individuals), with density highest in the northern half of the study area, peaking in the NE, where both observed and predicted densities of GRSH and NOFU were also high. Although relatively uncertain, the predicted distribution of SOSH was very similar (BA = 0.78) to that of the subset of wintering SOSH tracked from the Falkland Islands that used the study area during early summer, 2017 (7 out of 26 tracked birds). However, the tracked birds predominantly used the area in May, moving onto the Grand Banks at the end of the month (median 30th May; range 4th April -9th June), just prior to the commencement of the survey (Fig. A10).

Discussion
This study is the first detailed ship-based investigation of seabirdhabitat relationships in the SPFZ of the NW Atlantic, an area of highly dynamic physical oceanography. Our results indicate clear speciesspecific habitat associations, including evidence of zonation by water mass analogous to that observed in other major oceanic frontal systems (Ainley & Boekelheide, 1984;Ballance et al., 2006;Hyrenbach et al., 2007;Pakhomov & McQuaid, 1996;Wahl et al., 1989). Moreover, they confirm, as previously inferred from tracking data (Davies et al., 2021), that in the summer, seabirds from both North and South Atlantic breeding populations aggregate in globally important numbers in the SPFZ. Indeed, our results suggest that due to the presence of populations or demographic groups (e.g. immatures) currenty lacking tracking data, the abundance of several species is greater than hitherto supposed. They also confirm that several species moult in the area. Below, we consider potential sources of bias in these results, before discussing seabird/ habitat interactions; age, provenance and moult; and the conservation implications of our findings, in more detail.

Sources of bias and model criticism
During ship-based surveys, attraction of seabirds to the survey vessel can increase their local density, biasing abundance estimates upwards (Hyrenbach, 2001). Although we avoided recording individuals that For each species, the spatial smooth model (used to predict abundance) was identical to the global model, with an additional explanatory covariate comprising a smooth of spatial location. ‡ Mean logarithmic score calculated by 9-fold spatial cross-validation. * Bhattacharya affinity, quantifying the similarity between model-predicted distributions and the utilisation distributions of tracked birds (Leach's Stormpetrels were not tracked during the study). Not calculated for local inference models, because some explanatory covariates retained in these models were measured only along the survey track.  Fig. 6. Estimated distribution of the five most abundant seabird species in the study area. Column one shows the mean density of birds predicted by models fitted to the ship-based survey data, collected in June 2017 (survey track in red) and column two the corresponding coefficients of variation. Grey polygons delineate proposed and existing MPAs (see Fig. 1). For comparison, column three shows the mean Utilisation Distributions (UDs) of birds tracked via light-based geolocation between the dates indicated (Table A1 for details). Leach's Storm-petrels were not tracked during the study. Cumulative percentage UD volume contours are drawn at 50, 75, 90 and 95 % (see Fig. A10 for complete UDs of the tracked birds).
clearly exhibited this behaviour (Gjerdrum et al., 2012), seabirds can react to vessels at distances beyond the perception of ship-based surveyors (Bodey et al., 2014;Collet et al., 2015). Of the species we modelled, NOFU are considered to be the most prone to vessel attraction, and may therefore have been overestimated (Hunt et al., 2000). In contrast, the diving behaviour of some species can cause a downward bias in apparent abundance (Winiarski et al., 2014). Three of the most abundant species that we encountered (COSH, GRSH and SOSH) routinely forage by diving. However, biologging showed that during the daytime, when our survey was conducted, these species spent < 2 % of their time diving, so although we corrected for availability bias, its effect was small. Another potential source of bias in our results is our assumption that flying shearwaters and NOFU were detected perfectly. While frequently made (Tasker et al., 1984;Waggitt et al., 2020), this assumption is difficult to check because of the impracticality of measuring distances to flying birds at sea. If detection probability actually decayed markedly with distance within the 300 × 300 m 'snapshot' count area, then our abundance estimates would be biased downwards (Fifield et al., 2009). In the absence of data to the contrary, we assume that this is not the case for shearwaters and NOFU due to their relative conspicuousness. It was practicable to measure distances to LHSP (see Section 3.3) and detections of this species did decay with distance, presumably due to their much smaller size. We selected inference models by spatial cross-validation, which has the advantage over selection using information criteria, such as AIC, that it explicitly favours models that predict well in unobserved areas of space and therefore tends to reject biologically irrelevant covariates (Bell & Schlaepfer, 2016;Roberts et al., 2017). Nonetheless, a weakness of some of the resulting inference models was that they exhibited marked residual serial or spatial autocorrelation, especially those for GRSH and NOFU (Fig. A7). This is unlikely to have biased our habitat association parameter point estimates, but may have led to underestimates in their uncertainty. Residual autocorrelation often arises due to the omission of one or more important covariates. An obvious candidate in our study is prey availability, which is difficult to measure directly. Instead, we used environmental proxies (SST, NPP, etc.), but mismatches between these and the true distribution of prey are likely due to spatiotemporal lags induced by trophic dynamics, advection, etc. (Grémillet et al., 2008). In addition, there is growing evidence that seabirds use not only current conditions but also memory as foraging cues (Wakefield et al., 2015;Weimerskirch, 2007), so their distributions likely reflect past as well as current conditions. Moreover, aggregation and segregation due to social effects and competition can cause further divergence from an ideal free distribution (Veit & Harrison, 2017;Wakefield et al., 2013). The smooth of spatial location included in our spatial smooth models accounts for much of this variation, albeit not in an interpretable manner. As a result, residual autocorrelation in the latter models was slight (Fig. A7) and their predictions accorded well with independent tracking data ( Table 3), suggesting that our overall abundance predictions are robust.
Given the improvement in fit afforded by the inclusion of a spatial smooth in the prediction models, it is pertinent to explain further why we did not also include this in the inference models. Firstly, this was because many of the environmental covariates that we considered correlate with latitude (e.g. SST, NPP) or longitude (accessibility for NOFU and LHSP) and would therefore be confounded by a spatial smooth (Wood, 2017). Secondly, as explained above, the manner in which memory, social interactions and competition interact to give rise to spatial structure in seabird distributions remains unclear so a model containing a simple smooth of spatial location would lack ecological interpretability (Hodges & Reich, 2010). An alternative to the spatial smooth approach is to model spatial (or spatiotemporal) autocorrelation explicitly as a random effect, thereby partitioning its effects from those of habitat and accessibility. Techniques for achieving this in a distance sampling context are currently being developed (Bachl et al., 2019;Yuan et al., 2017), but at present they remain conceptually challenging and computationally expensive, placing them beyond the scope of the current study.
Although we validated our models against independent tracking data, there are several important caveats to this procedure. Firstly, it implicitly assumes that the mean density of birds at a given location in space (the response in our models) correlates with the mean proportion of time birds spend at that location (i.e. their utilisation probability, which we estimated from the tracking data). While this seems reasonable (Carroll et al., 2019), the relationship between these indices could vary markedly with the movement characteristics of the birds (Whitehead & Jonsen, 2013). Secondly, for logistical reasons, we tracked one species -COSH -a year after the at-sea survey period. Hence, although we tracked these birds at the same time of year, and their distribution (Fig. 6) was similar to that of COSH tracked from the same colony in previous years (Magalhaes et al., 2008;Paiva et al., 2010b), some dissimilarity might inevitably be expected due to interannual variability in the drivers of prey distribution. Finally, tracking data will only provide a fair validation of the model predictions if the tracked birds are representative of the population in the study area. This is more likely to be so for SOSH and GRSH, because in the Atlantic, these species breed predominantly in a few very large colonies (Brooke, 2004a;Clark et al., 2019). The high similarity of the tracked birds' distributions to the model-predicted distributions for these species is therefore reassuring. In contrast, NOFU breed at numerous sites along boreal and Arctic coastlines (Brooke, 2004a). Based on previous studies (Edwards et al., 2013), we tracked NOFU from Northern Scotland, but it has subsequently become clear that NOFU from other regions also use the study area , potentially explaining some of the disparity between the model-and tracking-based distributions (see Section 5.3).

Habitat associations and potential drivers of seabird abundance
For most species, the respective global and local inference models contained different environmental covariates, so it is important to consider the scales of processes that they can resolve (Wakefield et al., 2009). For all covariates, spatial resolution was truncated by segment length to ~ 3.4 km, so we were unable to detect associations with fine scale (<1 km) dynamics that are increasingly thought to contribute to patchiness in oceanic seabird prey (Bertrand et al., 2014;Levy et al., 2012). In addition, to avoid data gaps due to cloud interference, we time-averaged SST and NPP used in the global models (over 28 and 8 days, respectively; Table 1). Due to the relatively rapid surface currents in the study area (~10-20 cm/s (Miller et al., 2013b;Reverdin et al., 2003)), this may have given rise to some blurring of sub-mesoscale eddies, etc. We therefore regard SST in our models as indicative of meso to macroscale (100s− 1000s km) phenomena, such as the distribution of water masses. FG, calculated from blended SST data, was averaged over 3 days and therefore better able to resolve the mesoscale positions of persistent fronts (Miller et al., 2013b). nSST and ΔnSST, used in the local-inference models, were measured in situ, near simultaneously with the seabird data, and could therefore resolve coarse scale (10s km) processes. Although SLA and EKE are unaffected by clouds, their relatively low resolution (~30 km) means that they can only resolve mesoscale phenomena.
At the macroscale (1000s km), the correlation between total seabird biomass and NPP (Fig. 3) presumably reflects the bottom up limitation, ultimately due to the latitudinal gradient in nutrient supply between the subtropical and subpolar gyres (Longhurst, 1998). At finer scales, NPP was rejected during model section for all species except SOSH (weak negative association - Fig. 5), presumably due to the many processes that modulate energy flow to higher trophic levels in pelagic ecosystems. For example, while some mesoscale eddies in the SPFZ remain essentially stationary for six months or more, others translate at 2-3 cm/s in the direction of the NAC (Shoosmith et al., 2005). Assuming that primary production takes weeks to months to reach the trophic levels of seabirds (Lehodey et al., 2010), spatial lags of the order of 10 -100 s of km could occur between the onset of an eddy-induced phytoplankton growth and prey consumption by seabirds.
In almost all models, bird density was most strongly associated with sea surface temperature (Fig. 5), supporting the supposition that at the meso-to macroscale seabirds have a zoned distribution with respect to water mass in the SPFZ of the NW Atlantic (Boertmann, 2011;Skov et al., 1994). Similar zonation occurs in other major oceanic frontal systems (Ainley & Boekelheide, 1984;Ballance et al., 2006;Hyrenbach et al., 2007;Pakhomov & McQuaid, 1996;Wahl et al., 1989) and in particular, the distribution of NOFU, SOSH and LHSP with respect to water mass was comparable to that in the North Pacific (Wahl et al., 1989) (GRSH and COSH are absent from the Pacific). The apparent lack of an association between GRSH and water mass at this scale may be due to our study design. The study area straddles the range edges of other species, providing strong habitat/density contrast, whereas it is almost entirely within the core range of GRSH (Fig. A10). At finer scales however, the local inference model indicated that GRSH tended to occur in patches of cool SAIW/FW. Several potential reasons for seabird zonation present themselves. Firstly, it may be due to the habitat preferences/thermal tolerances of prey (Hunt, 1997;Pakhomov & McQuaid, 1996;Pocklington, 1979;Sydeman et al., 2010). The fact that the meridional temperature and salinity gradients are stepped (i.e., at the NSPF, SSPF and MAF - Fig. 3) rather than continuous in the study area could accentuate this mechanism. Little is known about the diets of seabirds in the area, but studies in adjacent waters suggest that while the four medium-sized petrels predominantly consume fish, squid and swarming crustaceans, NOFU take a higher proportion of zooplankton (e.g., hyperiid amphipods) (Danielsen et al., 2010;Furness & Todd, 1984;Garthe et al., 2004;Ojowski et al., 2001;Phillips et al., 1999) and COSH more fish (Brown et al., 1981;Granadeiro et al., 1998;Neves et al., 2012;Powers et al., 2020;Ronconi et al., 2010a;Xavier et al., 2011). These taxa all exhibit marked zonation by water mass within the study area (Cook et al., 2013;Letessier et al., 2011;Vecchione et al., 2010a;Vecchione et al., 2010b). Secondly, zonation could be mediated by the thermal tolerances of the seabirds themselves (Fort et al., 2012;Pocklington, 1979). Higher latitude species tend to be adapted to cold conditions via higher metabolic rates, sustained by higher food availability, rather than increased insulation (Weathers et al., 2000). COSH and LHSP primarily occur in subtropical-tropical waters Gonzalez-Solis et al., 2007;Halpin et al., 2018;Pollet et al., 2014) and NOFU in polar/subpolar waters  but is unknown if this is directly facilitated via such thermal adaptation. Moreover, GRSH and SOSH routinely range between subtropical and polar waters (Flood & Fisher, 2020), so it is unlikely that these species are thermally limited within the study area. Finally, apparent associations between species and water masses could be due to correlations between SST and other causal phenomena. For example, COSH breed just south of the study area in the Azores, making it difficult to separate the effects of accessibility and SST for this species.
Our models indicated associations between some species and fronts. At the mesoscale, fronts can enhance nutrient supply and therefore primary production and prey availability (Tilstone et al., 2014). At finer scales, they can also aggregate prey due to convergent currents (Belkin et al., 2014;Haney & McGillivary, 1985;Scales et al., 2014). Our global inference models indicated that the former processes may be important for NOFU, SOSH and GRSH, all of which were positively associated with front gradient (Fig. 5). In contrast, the local models indicated that only NOFU were associated with coarse-scale thermal gradients (ΔnSST). This may be due to interspecific differences in foraging behaviour: NOFU and COSH forage predominantly at or < 3 m from the surface Paiva et al., 2010a), conceivably making them more reliant on processes that aggregate prey near the surface (Haney & McGillivary, 1985), while GRSH and SOSH typically forage by diving (to ≥ 15 and 40 m, respectively) Ronconi et al., 2010b;Shaffer et al., 2006). The high proportion of zooplankton in the diets of NOFU may make them particularly reliant on prey aggregation by fronts. Although COSH associate with fronts in the Gulf Stream further south (Haney & McGillivary, 1985), they did not do so in our study area. This may be because they feed primarily on nekton in oceanic waters (Granadeiro et al., 1998;Neves et al., 2012;Xavier et al., 2011). Notably however, 11 % of the COSH in our study were associated with dolphins (principally Delphinus delphis or Stenella coeruleoalba) or tuna (Thunnus spp.), presumably because these taxa drive prey near to the surface (Clua & Grosvalet, 2001;Martin, 1986). In contrast, no NOFU or SOSH and only 2% of GRSH, were associated with subsurface predators.
In common with previous studies (Camphuysen, 2007;Haney, 1986;Tew Kai & Marsac, 2010), the distribution of some species was associated with mesoscale turbulence. Associations between SLA and NOFU (positive), SOSH (positive), COSH (U-shaped but predominantly negative), are consistent with these species preferentially foraging within warm-core (all species) and cold-core (COSH) mesoscale eddies or meanders. COSH also associate with cold-core eddies in the Gulf Stream/continental slope system further south (Haney, 1986), but unlike in our study area, where they occurred most frequently in eddy interiors, they associated with eddy edge fronts there. Haney (1986) also found that GRSH and LHSP associated with the eddy edges (see also Camphuysen, 2007). Although we did not find any widespread associations between these species and SLA or EKE, we observed relatively high densities of both foraging at the margins of the prominent cold-core eddy in the west of the study area noted above (cf. Fig. 2 & Fig. A2). This is one of several standing cyclonic eddies caused by interactions between the NAC and continental slope east of the Grand Banks (Reverdin et al., 2003;Rossby, 1996). Due to their proximity to the shelf, these eddies are hypothesised to be more productive than transient eddies further east (Browning et al., 2021). This combination of high productivity and predictability may make these eddies particularly important for foraging seabirds and other higher predators, a hypothesis that warrants further investigation.
The previous example illustrates one of the limitations of our study. Model selection by spatial cross-validation is conservative in that only covariates that improve predictive performance throughout the study area are retained (Roberts et al., 2017;Wenger & Olden, 2012). As a corollary, and given that our data provide a relatively narrow spatiotemporal snapshot, our models were unlikely to detect less frequent seabird-habitat associations well. Most notably, associations with water masses and fronts only partially explain the peak in seabird density observed in the NE of the study area (~51 • 6 ′ N, 32 • 36 ′ W, Fig. 6). Additional evidence supports the supposition that higher predators concentrated in this area in early summer, 2017: firstly, three of the four shearwater species that we GLS-tracked clearly spent a disproportionate amount of time there (Fig. 6). Secondly, in May, one of three incubationstage NOFU GPS-tracked during a separate study (PMT, unpub. data) made a ~ 2000 km beeline directly there from its colony in Northern Scotland (Fig. 7). Thirdly, counts of baleen whales (Balaenoptera spp. and Megaptera novaeangliae) recorded during our cruise peaked in the area (Wakefield, 2018). This area coincides with both the NSPF and a poorly-defined cold core eddy/meander (Fig. 2), both of which may have enhanced prey availability via the mechanisms discussed above. However, residuals from the global, local and spatial smooth models in the area were relatively high indicating that additional mechanism(s), beyond those described by the environmental covariates in our models, may have enhanced prey availability in the area. Multi-year/multi-species tracking shows that during summer, seabird density peaks on average further SW than this hotspot, at ~ 47 • 42 ′ N, 36 • 0 ′ W, (BirdLife International, 2019), suggesting that seabird distributions can shift meridionally by at least ~ 390 km interannually. Presumably, this reflects variations in physical drivers in the area, which occur at similar spatial scales. For example, the major fronts can migrate meridionally by ~ 300 km and the NAC vary between 2 and 4 major branches at interannual-decadal timescales (Belkin & Levitus, 1996;Bower & von Appen, 2008;Holliday et al., 2020;Miller et al., 2013b). Ultimatley, more data, collected either from ships or via high-resolution seabird tracking, are would be needed to determine the extent and causes of variations in seabird distribution in the SPFZ.
Our results show that the macroscale distribution of seabirds in the NW Atlantic is comparable to that in equivalent systems in other ocean basins, with abundance being much lower in oceanic than adjacent neritic waters, presumably due to differences in primary productivity. E. g., densities of GRSH, SOSH and COSH were approximately an order of magnitude lower than on the Grand Banks (Carvalho & Davoren, 2019;Fifield et al., 2009;Haney & McGillivary, 1985;Powers & Brown, 1988) and those of NOFU were around half (Fifield et al., 2009). Within the ocean zone of the NW Atlantic, seabird abundance is higher west of the MAR than to the east (Bennison and Jessopp, 2015;Davies et al., 2021). Our results indicate a number of potential explanations for this. Firstly, while seabird abundance was not generally correlated with primary production at the mesoscale (Fig. 5), it was at the macroscale (Fig. 3). At this scale, both primary and secondary production are higher to west of the MAR, due, it is thought, to mesoscale turbulence generated by the NAC being more intense in this region and therefore supplying more nutrients from depth and across the SPF (Druon et al., 2019;Letessier et al., 2009;Longhurst, 1998;Tilstone et al., 2014;Vecchione et al., 2010a). Presumably, eddies and fronts aggregate zooplankton more to the west of the MAR for the same reason. In addition, the locations of the major jets/fronts and associated standing eddies are more constrained by bathymetry west of the MAR (Miller et al., 2013b;Rossby, 1996;Shoosmith et al., 2005;Søiland et al., 2008). Given that many seabirds travel 1000s− 10,000s of km to forage (Edwards et al., 2013;Egevang et al., 2010;Hedd et al., 2012;Kopp et al., 2011), the resulting combination of spatial predictability and high primary production may allow the western mid-North Atlantic to sustain higher seabird biomass than waters east of the MAR.

Moult, age and provenance of birds in the study area
Flight feather moult is a vital maintenance activity for seabirds but it is particularly costly because elevated nutrition requirement due to feather synthesis coincides with impaired flight mobility due to feather loss (Ellis & Gabrielsen, 2002). It is not practicable to assess feather growth accurately at sea without catching or closely photographing birds (Keijl, 2011), so our results represent the minimum proportions of birds that would have been in primary moult. Nonetheless, together with earlier studies (Boertmann, 2011;Wynne-Edwards, 1935), our results show that the study area is an important moult site for multiple seabird populations in early summer. For example, contrary to the prevailing view that GRSH moult mainly in inshore productive waters of the NW Atlantic (reviewed by Huettmann & Diamond, 2000), our observations imply that a large proportion of adults from Gough undergo moult partially or fully off-shelf, in the SPFZ. They also support biologging-based inferences that many SOSH and South Polar Skuas also do so (Hedd et al., 2012;Kopp et al., 2011), and indicate for the first time that large numbers of immature or non-breeding adult NOFU moult in the area (see below). Presumably, birds travel very large distances to moult in the SPFZ due to particularly favourable conditions there. These could include not only high food availability but also possibly low competition or predation relative to the neritic zone.
Several inferences about the life history stages and origins of seabirds in the study area can be drawn from our results. Firstly, while tracking has previously shown that adult Arctic Terns and Long-tailed Jaegers use the SPFZ as a migratory staging area in spring/autumn (Egevang et al., 2010;Gilg et al., 2013), our observations show that immatures of both species also use the area in early summer, presumably during their first or second northerly migrations (see also Boertmann, 2011). They also confirm the inference from tracking that immature gannets, most likely from West Atlantic colonies, use the area (Fifield et al., 2014). Secondly, although our tracking data confirm that breeding NOFU from the Eastern Atlantic use the NW Atlantic in June (see also Edwards et al., 2013), the majority of NOFU that we sighted at sea were probably immatures. This is because ≥ 60 % were moulting and whereas immature NOFU begin wing moult in April/May (Brown, 1988), adults do not do so until July-August (Allard et al., 2008;Grissot et al., 2020;Quinn et al., 2016). Thirdly, the presence of dark phase NOFU indicates that birds from Arctic populations use the study area in June, albeit in unknown numbers (dark phase NOFU orginate only from more northerly colonies (van Franeker & Wattel, 1982)). Tracking has recently confirmed that NOFU from Iceland and Jan Mayen use the area during the summer  and post-breeding adults tracked from Devon Island (Nunavut) used the study area extensively . Fourthly, at-sea observations indicated a more north-westerly range limit for COSH than indicated by tracking breeding adults from the Azores (Magalhaes et al., 2008;Paiva et al., 2010a; this study). This could be because immatures or non-breeding adults are not subject to central-place constraint and can therefore exploit foraging areas further from colonies than adults can Zango et al., 2020). Adult and immature COSH usually commence primary moult after our survey period (Flood & Fisher, 2020), so their moult status during our survey is uninformative in this respect. However, tracking shows that COSH from the world's largest COSH colony, Selvagem Grande (30 • 09 ′ N, 15 • 52 ′ W), use the area in their first to fourth summers (Paulo Catry, unpub. data). Tracking immatures from the other breeding locations, especially the Azores, would further resolve this question. Lastly, the relatively high abundance of LHSP in the study area was somewhat unexpected. Although tracking has shown that post-breeding LHSP from Nova Scotian colonies stopover in the study area in autumn during their south-eastward migration , in June incubation-stage adults from Newfoundland colonies remain mainly to its west (Hedd et al., 2018). Notably, 16 out of 17 LHSP caught on the ship during our survey (EDW, unpub. data) were moulting rectrices and/ or primaries. Newfoundland breeders do not commence tail and primary moult until August and October-November, respectively (Hedd & Montevecchi, 2006) so most LHSP in the study area were probably immatures or failed/non-breeders (see also Boertmann, 2011).

Abundance of seabirds in the study area and conservation implications
Our data show that in summer the avifauna of the NW Atlantic SPFZ is overwhelmingly dominated by medium-sized petrels (98% by biomass), and especially GRSH (~63% by biomass; Table 4). Our estimate of the numerical abundance of GRSH in the NACES pMPA is very similar to that made by Davies et al. (2021) from tracking and colony size data (Table 5). For the remaining species, our estimates differ substantially from the latter study but this may be mostly due to differences in how abundance was defined. Firstly, Davies et al. estimated mean abundance during April-June, whereas we do so for June only. Hence, it was not unexpected that we observed a lower abundance of SOSH, because tracking shows that this species uses the area predominantly in April and May, with many birds moving onto the Grand Banks in June (Hedd et al., 2012). Secondly, while Davies et al.'s results are based on tracking between 1999 and 2015 (years covered vary between species), our results are for a single year, so interannual variability could result in differences in observed abundance. Thirdly, Davies et al.'s abundance estimates include only populations with tracking data coverage (they define populations based on the Large Marine Ecosystem or LME from which birds originate). For example, the NOFU tracking data contributing to their study all came from northern Scotland so their abundance estimate in effect pertains only to birds from the North Sea LME.
Notwithstanding these caveats, it is notable that we estimated > 7 times more COSH in the NACES pMPA than Davies et al. Their estimate is based on tracking data from a wide range of colonies, so omission of an important population seems unlikely. Moreover, COSH rarely follow ships (Flood & Fisher, 2020), so this is unlikely to have biased our estimate upwards. Immatures could have been more abundant in the area than we assume (Table 5), but to account for all of the difference, they would have to outnumber adults 6:1. This ratio higher than is typical for petrels (Carneiro et al., 2020), but could perhaps arise if substantial numbers of immatures from other populations (not just the Azores) occur in the area (see Section 5.3). In addition, the range of immature COSH in the NW Atlantic may have shifted northwards recently (Gjerdrum et al., 2018), further increasing their abundance in the study area. A second, not necessarily mutually exclusive possibility, is that the size of the Azores breeding population (upon which Davies et al.'s estimates are largely based) may be substantially larger than hitherto supposed. The current estimate (~188,000 breeding individuals (BirdLife International, 2020)) is based on counts of birds rafting at sea prior to entering colonies, and very limited colony surveys, neither of which are regarded as accurate (Bolton, 2001;Brooke, 2004a;Monteiro et al., 1996;Oppel et al., 2014). The size of the Azores population may therefore warrant reassessment. If the population is larger than currently supposed, it follows that our estimate of the proportion of the Atlantic population occurring in the study area (Table 4) will be biased upwards.
Despite the density of most species being lower than in adjacent neritic waters, our study supports the supposition that globally important numbers of seabirds occur in the SPFZ of the NW Atlantic due to the vast area over which they aggregate (Davies et al., 2021). For example, in June 2017 ~ 26% of the world's GRSH occurred in the study area (Table 4). This includes multiple life history stages and populations and birds engaged in vital maintenance activities (foraging, migratory staging and moult). Moreover, both seabird abundance and diversity are high year-round compared to the rest of the oceanic zone of the North Atlantic, with 29 species having now been recorded. In particular, the area holds high numbers of alcids and kittiwakes in the winter Fayet et al., 2016;Fort et al., 2012;Frederiksen et al., 2012Frederiksen et al., , 2016. Most of these species are undergoing unsustainable population declines . Our results therefore support the view that the area warrants protection (Davies et al., 2021).
During our survey, birds were concentrated further north than the mean centre of aggregation identified by Davies et al., partially overlapping with the recently-established Charlie-Gibbs North and South High Seas MPAs (Figs. 1 & 6). These sites were designed primarily to protect elasmobranchs, demersal fish stocks associated with seamounts, cetaceans, turtles, and the sub-polar front ecosystem (OSPAR Commission, 2010; OSPAR Commission, 2012). Our results underscore that these sites, plus the proposed NACES MPA and existing Mid-Atlantic Ridge North of the Azores MPA to the south (Fig. 1), sustain important aggregations of seabirds. In order to be effective, seabird-specific management measures should therefore be coordinated across this entire MPA network. Although poorly quantified, there is currently thought to be relatively little fishing activity in the area (ICES, 2020) so seabird bycatch is unlikely to be large. However, due to the global trend is for industrial fisheries to intensify, move further offshore and target lower trophic levels (Crespo et al., 2018), this could become a significant cause of seabird mortality in the area in the future. Given that subsurface predators apparently facilitate foraging by Cory's shearwaters (Clua & Grosvalet, 2001;Martin, 1986), declines in heavily exploited species, such as Atlantic bluefin tuna (T. thynnus), could have an indirect negative impact on this species. In addition, the SPFZ ecosystem may be particularly sensitive to climate change (Beaugrand et al., 2015;Pershing & Stamieszkin, 2020). Since the mid-20th century, the Atlantic Meridional Overturning Circulation has slowed, the Gulf Stream has migrated northward and the subpolar gyre freshened, resulting in rapid Table 5 Estimated mean abundance of seabirds (x 1000 individuals) in the proposed North Atlantic Current and Evlanov Seamount (NACES) Marine Protected Area based on at-sea data (this study) and tracking (Davies et al., 2021 changes in surface water temperatures in the NW Atlantic (Caesar et al., 2018). Under current emission scenarios, these changes are forecast to continue or accelerate (Sgubin et al., 2017). In response, regional NPP is declining (Bopp et al., 2013;Osman et al., 2019;Saba et al., 2016) and higher trophic level community structure altering (Beaugrand et al., 2010;Beaugrand et al., 2015;Pershing & Stamieszkin, 2020;Reygondeau & Beaugrand, 2011;Villarino et al., 2015). Resulting impacts on seabirds and other higher predators are difficult to anticipate but for example, based on the energetic niche analysis, it has been predicted that wintering seabird distributions in the NW Atlantic will shift north and west (Clairbaux et al., 2021). Monitoring the behaviour, demography and distribution of seabird populations that use the SPFZ could allow early detection of these effects and wider ecosystem changes (Brisson-Curadeau et al., 2017), guiding dynamic management and if necessary, realignment of the regional MPA network. Our results provide a baseline against which changes in seabird abundance in the oceanic NW Atlantic can be assessed and indicate which potential drivers of seabird distribution in the area could be most usefully be investigated further.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.