Environmentally driven sexual segregation in a marine top predator

Sexual segregation in foraging occurs in many animal species, resulting in the partitioning of resources and reduction of competition between males and females, yet the patterns and drivers of such segregation are still poorly understood. We studied the foraging movements (GPS-tracking), habitat use (habitat modelling) and trophic ecology (stable isotope analysis) of female and male Cory’s shearwaters Calonectris borealis during the mid chick-rearing period of six consecutive breeding seasons (2010–2015). We found a clear sexual segregation in foraging in years of greater environmental stochasticity, likely years of lower food availability. When food became scarce, females undertook much longer foraging trips, exploited more homogeneous water masses, had a larger isotopic niche, fed on lower trophic level prey and exhibited a lower body condition, when compared to males. Sexual competition for trophic resources may be stronger when environmental conditions are poor. A greater foraging success of one sex may result in differential body condition of pair mates when enduring parental effort, and ultimately, in an increased probability of breeding failure.

During 2010, 2011 and 2013 there was a significantly smaller overlap between sexes for the 50%, 95% and 99% Kernel UD, compared to 2012, 2014 and 2015, with a mean increase by ~26% in the overlap of the entire UDs (Table 3).
Overall, generalized additive mixed models (GAMMs) showed a good predictable capacity, explaining >22.4% and >22.7% of the deviance in the first-passage time (FPT) duration (proxy of foraging activity) of females and males respectively (according to Table 4). In females, FPT duration increased with increasing bathymetry (BAT) and BAT gradient (BATG) and decreasing SST and distance to colony (DCOL) (Fig. 2). In males, FPT duration increased with increasing BAT, BATG and CHL gradient (CHLG) and decreasing SST and DCOL (Table 4, Fig. 2).

Discussion
Our study documented a clear sexual segregation in foraging by Cory's shearwaters in years of increased environmental stochasticity, i.e. likely years of low food availability. During years of good environmental conditions, both sexes foraged mostly within the colony surroundings, although females tended to exploit in more pelagic waters. In years of greater stochasticity, females undertook much longer foraging trips, similar to those performed during the pre-laying 25 and incubation 30 periods, which decreased their body condition. They also enlarged their isotopic niche, but fed mostly on lower trophic level prey, when compared to males. Conversely, males kept foraging in the colony surroundings, maybe outcompeting females to access comparatively higher trophic level preys, which translated into a higher body condition index and daily mass gain during foraging trips.  28 . Even though at-sea patterns were generally similar between sexes, sexual segregation was observed in years of greater environmental stochasticity, when marine productivity was lowest. Birds travelled greater distances from their colony, and spent more time in foraging areas than they did during years of higher marine productivity. Therefore, previous studies finding no evidence of sexual segregation in foraging strategies of breeding Cory's shearwater 20,31 are partly explained by the fact that they were conducted solely during the incubation period, and presumably during years of relatively good environmental conditions. In fact, Ramos et al. 20 suggested that more subtle sexual differences might exist, for instance in diving behaviour, or be detected in resource partitioning (e.g., sexual differences in the size of consumed prey). While we did not examine such factors in this study, we could show significant sexual differences in Scientific RepoRts | 7: 2590 | DOI:10.1038/s41598-017-02854-2 foraging distribution and oceanographic features within the birds' foraging areas, as well as in several aspects of their trophic ecology. Foraging trip characteristics, features of the foraging areas (e.g. SST) and body condition of females during years of poor environmental conditions were significantly affected, and thus differed from those of males, (and also from females during years of good environmental conditions). Therefore, our study suggests that sexual differences in the foraging ecology of Cory's shearwaters are likely to be perceptible under poor environmental conditions. Our results make sense in the light of the competitive exclusion by the dominant sex, i.e. males over females 3 , which will occur mostly when environmental conditions are poor and resources are scarce. Furthermore, previous studies regarding such effects were conducted mostly during a single breeding season, and thus failed to detect any important influence of environmental conditions on sexual segregation in foraging (e.g. ref. 22).

Trophic ecology
Plasma

Body condition
Adults' body condition index (BCI) Differential effects of foraging choices on the trophic ecology of both sexes. Overall, the carbon isotopic values suggest that males used to feed on prey from coastal environments (i.e. higher δ 13 C value) when compared to females, which exhibited a more pelagic foraging behaviour 32 . While foraging over more pelagic environments, mostly during 2010, 2011 and 2013, females showed lower δ 15 N values, which suggests they were feeding on lower trophic level prey than males. There might be at least two reasons for this; 1) females were more explorative in 2010, 2011 and 2013 (i.e. years of poor environmental conditions, likely with lower food availability), feeding on prey of comparatively lower δ 15 N values, such as small pelagic fish species (e.g. Scomber sp; ref. 23) and enlarging their isotopic niche or, 2) in years of food scarcity birds tend to attend more to fishery discards 33 where males might outcompete females for offal and discards and gain access to higher trophic level prey. This has been reported to occur in several Procellariform species (e.g. northern giant petrels Macronectes halli 34 ), potentially leading to sex-biased bycatch rates (see review by Gianuca et al. 35 ). Even though Cory's shearwaters are considered to interact much less with fisheries and/or to be much less frequently victim of by-catch than their Mediterranean congeneric, they have been observed feeding on fishery discards 36,37 . Although both hypotheses can be true, only the collection of more tracking and blood data along with vessel monitoring system data during subsequent years will allow disentangling the major driver of this sexual isotopic segregation pattern. Nevertheless, both groups were isotopically segregated even in years where both sexes were foraging in similar Final remarks. Our results are compatible with the interpretation that sexual segregation in Cory's shearwaters might be mediated by habitat segregation from neritic to oceanic regimes (by males and females, respectively), with such spatial segregation occurring mostly in years of increased environmental stochasticity.
Multi-year tracking studies are thus crucial (6 breeding seasons in our study) to successfully detect even relatively small sexual differences in the foraging ecology and behaviour of marine predators. Such small differences may be ecologically relevant because sexual competition for trophic resources may be stronger when environmental conditions are poor. A greater foraging success of one sex over the other may result in differential body condition of pair mates to endure on reproductive duties, thus resulting in an increased probability of breeding failure.  Table 1). The devices weighed 17 g, which represented between 2.2% and 2.9% (median = 2.6%) of the birds' mass and was below the 3% threshold advised by ref. 39.
Devices were set to record locations each 5 minutes. GPS loggers were then attached using TESA ® tape to the contour feathers along and in between both scapulas. The whole process took less than 10 minutes, thus minimizing    Table 3. Observed and randomized overlap between female and male Cory's shearwaters Kernel Utilization Distributions (Kernel UD). P represents the proportion of randomized overlaps that were smaller than the observed overlap. Significant differences are shown in bold.
the overall stress to the animal. Upon logger retrieval, a blood sample of about 0.5 ml was collected from the tarsal vein of each individual for stable isotope analysis (SIA). The tracked birds were measured (gonys height and the length of wing, tarsus and culmen) and weighed at both capture and recapture. Body measurements (except body mass) were included on a Principal Component Analysis (PCA) and the PC1 scores were used as a measure of structural body size to compute a body condition index (BCI). This index was obtained from the residuals of the linear regression of body mass on PC 1 scores. We used BCI as an indicator of energetic reserves (i.e. fitness parameter). All animals were handled in strict accordance with good animal practice as defined by the current European legislation. All

Area -Restricted Search (ARS) zones. Fauchald and Tveraa 40 developed a technique, named first pas-
sage time (FPT) to assess the spatial scale that animals use to encounter their prey. FPT is, by definition, the time required for an animal to pass through a circle with a given radius r. By moving this circle along the path of the animal, we will obtain a scale-dependent measure of search effort and therefore the behavioural response of an individual in the environment. Because top marine predators usually forage in a patchy and hierarchical environment 41 , increases in the turning rate and/or decreases in speed of its foraging path should be related to the so-called area-restricted search (ARS) behaviour. ARS will then appear as an individual reaction to changes in resource availability and distribution, by increasing the residence time in the productive patch 40 .
Zones of area-restricted search (ARS) were estimated applying FPT analysis, following 40 and using software R 3.0 (R Development Core Team 2014). Locations were first projected onto a Lambert equal-area projection. Usually, in water positions result in very small-scale ARS zones (<100 m diameter), which considerably increases the variance in FPT and can camouflage larger-scale ARS zone 42 . To address this problem, we removed bouts on the water and interpolated locations to obtain a distance interval of 0.1 km for FPT analysis 43 . We considered positions with speed <3 km as resting or preening behaviours on the water or inland, after inspecting the frequency distribution of speeds. Following the recommendations of ref. 43, FPT analysis was performed in two steps: (1) to detect large-scale ARS we ran the analysis on the whole path, estimating the FPT every 1 km for a radius r from 1 to 50 km; (2) to detect small spatial scale events we run again FPT analysis every 0.1 km for an r varying between 0.1 and 10 km. The plot representing variance in log (FPT) as a function of r allowed us to identify the ARS scales by peaks in the variance. In this calculation, FPT was log transformed to make the variance independent of the Scientific RepoRts | 7: 2590 | DOI:10.1038/s41598-017-02854-2 magnitude of the mean FPT 40 . It is also possible to locate where the bird entered an ARS zone and the time spent in this area by plotting FPT values where a peak of variance occurred as a function of time since departure from the colony. ARS locations were also used to feed the Generalized Additive Mixed Models (GAMMs).
Habitat use. GPS locations of each bird where ARS behaviour was detected (ARS zones) were examined under the adehabitatHR R package 44 generating Kernel Utilization Distribution (Kernel UD) estimates within the R environment 45 . The most appropriate smoothing parameter (h) was chosen via least squares cross-validation for the unsmoothed GPS data, and then applied as standard for the other datasets and grid size was set at 0.05° (to match the grid of environmental predictors). We considered the 50% and 95% kernel UD contours to represent the core foraging areas (FR) and the home range (HR), respectively.
The extent of within-year overlap between male and female home ranges was estimated using using a kernel UD overlap index, which is considered the most appropriate measure of overlapping space use 46 . We used a randomization technique (1000 randomizations of our dataset) to test the null hypothesis that there was no difference in the spatial distribution of males and females in each study year. If the null hypothesis is true, overlap between males and female 50% and 95% kernel UDs should not differ significantly from that calculated if sex were randomly assigned. P-values were determined by the proportion of random overlaps that were smaller than the observed overlap (see refs 18, 47 and 48 for a similar approaches). More details on measures of spatial overlap can be found in ref. 46. Environmental data. The extended winter North Atlantic Oscillation (NAO) index was used as a large-scale environmental predictor for the North Atlantic area, and specifically for the Western Iberia Upwelling Ecosystem (WIUE). The NAO index refers to a north-south alternation in atmospheric mass between the subtropical Atlantic and the Arctic, and thus involves out-of-phase behaviour between the climatological low-pressure centre near Iceland and the high-pressure centre near the Azores. In other words, the index is the result of the difference of normalised sea-level pressures (Pa) between the former two locations (https://climatedataguide. ucar.edu/climate-data/hurrell-north-atlantic-oscillation-nao-index-station-based). In general, a low NAO index value would be depicting an intense upwelling and low SST due to stronger winds (during the previous winter).
Under these conditions, we should expect a lower abundance of plankton 49 and consequently lower abundance and availability of fish prey to top marine predators, such as seabirds. In some upwelling regions, such as the Portuguese coast, an overly intense upwelling may be responsible for a very low recruitment of small pelagic fish (such as sardine Sardina spp.) and low abundance of plankton because, under these conditions, the fish larvae and plankton are driven offshore, which consequently increases larvae mortality and creates a spatial mismatch of plankton and juvenile planktonic fishes a few months later 50 . In addition, the effect of the NAO is regionally dependent, and a reverse effect (i.e. weaker wind fields and higher SST) should occur in the northern and western Atlantic regions 51 .
We used small-scale environmental predictors, such as chlorophyll a concentration (CHL) and sea surface temperature (SST) data, downloaded from http://oceanocolor.gsfc.nasa.gov/, as daily night-time products with a resolution of 0.04° (approx. 4 km) in the SMI-HDF format. Bathymetric data (BAT), taken as water depth, was downloaded from the ETOPO2v2 database at a spatial resolution of 0.03° (approximately 3 km; http://www.ngdc. noaa.gov/mgg/fliers/01mgg04.html). HDF files were converted to raster using the Marine Geospatial Ecology Tools in ArcGIS v10.1 52 , and then to ASCII to create composites. All composites were constructed using the freeware R environment and different functions of the raster package. Spatial gradients of SST, CHL and BAT (SSTG, CHLG and BATG, respectively) were obtained by estimating the proportional change (PC) within a surrounding 3 × 3 cell grid using a moving window as follows: PC = [(maximum value − minimum value) × 100/maximum value] 53 . SSTG and CHLG are believed to be good indicators of oceanic fronts, while the BATG was used as a proxy for slope. Additionally, two static variables were generated. Distance to colony (DCOL) was calculated using the distance tool (spatial analyst toolbox) in ArcGIS v10.1.

Stable isotope analysis (SIA).
Stable-nitrogen isotope ratios ( 15 N: 14 N, expressed as δ 15 N) and stable-carbon isotope ratios ( 13 C: 12 C, expressed as δ 13 C) in the plasma of Cory's shearwater were determined to investigate the trophic choices of each sex during each year. Plasma has a half-life of about 3.5 days 54 (i.e. high turnover rate), therefore it represents prey ingestion and trophic ecology of tracked individuals during the last trips before sampling 55 . The δ 15 N is mainly used to define the trophic position of the consumer 56 , while δ 13 C reflects the foraging habitat of the consumer 57 . There is a gradient of high to low values of δ 13 C from benthic and inshore to pelagic and offshore food webs, because the organic enrichment at the coast is gradually diluted towards the open ocean 58 .
Each of the tracked birds was sampled upon return from a foraging trip. Blood samples (around 0.5 ml) were collected from the tarsal or brachial vein using insulin-syringes with 27 G needles. Blood samples were then separated into plasma and RBC by centrifugation at 12000 rpm for 5 min, within 2-4 hours of sampling and stored frozen at −20 °C until preparation for analysis. Successive rinses with a 2:1 chloroform-methanol solution were performed on the plasma for delipidation 55 . Isotope ratios of carbon and nitrogen of plasma were then determined by continuous-flow isotope ratio mass spectrometry, using an EA-IRMS (Isoprime, Micromass, UK). The analytical precision for the measurement was 0.2‰ for both carbon and nitrogen. All values presented are means ± 1 SD unless otherwise stated.

Statistical analysis. Generalized Linear
Mixed Models (GLMMs) tested the effect of (1) year (2010-2015), (2) sex and, (3) the interaction between year and sex (i.e. independent variables) on the mean values of (1) NAO index (Jun-Sep), (2) CHL, SST and SST anomaly (within 200 km of the breeding colony), (3) trip duration, number of long trips number of short trips −1 (with long trips, ≥5 days of duration and short trips, ≤4 days of duration as defined by ref. 24), max. distance from colony, time spent flying trip −1 day −1 , % of time spent in foraging areas (ARS zones), (4) CHL, SST and SST anomaly within ARS zones, (5) δ 13 C, δ 15 N and SEA B of plasma and (6) adults' body condition index (BCI) and mass gain trip duration −1 (i.e. dependent variables). Trip identity was nested within the individual as a random term to avoid potential pseudo-replication problems, since all individual birds performed multiple trips. Gaussian distribution of error terms and a log-link function were used in the modelling. Post-hoc multiple comparisons with Bonferroni correction were used to identify significant differences between categories of each independent variable. R packages used in the GLMMs were lme4 59 and lmerTest 60 .
When modelling the occurrence of ARS behaviour (First Passage Time -FPT -duration) in male and female Cory's shearwaters we used GAMMs to (1) select the most parsimonious models explaining FPT and (2) estimate smoothers for each of the environmental parameters for the top-ranked models (ΔAICc < 2). GAMMs combine the utilities of linear mixed models 61 and generalized additive models 62 so that random factors, fixed factors and nonlinear predictor variables can all be estimated in the same statistical model. Separate models were developed for male and female birds aiming at simpler interpretations of their outputs (i.e. interactions with the variable sex would be difficult to interpret in complex models). Such models included (1) year and (2) all different environmental predictors of productivity (e.g. SST) as fixed factors, trip identity within bird identity as a random term (to account for pseudoreplication issues).
As part of the GAM functions within the mgcv R package ref. 63 the smoothing parameter is chosen automatically using generalised cross-validation (GCV). In order to model spatial auto-correlation we included an isotropic thin plate spline which is set up as a two-dimensional smoother based on both x and y coordinates (i.e. s(x,y)). Incorporating a spatial smoother is one means of modelling a spatial trend within a model, more details on this approach can be found in ref. 63. Prior to modelling we examined the correlations between all environmental variables in order to ascertain whether collinearity may have occurred. We assumed that a Spearman correlation coefficient higher than 0.5 was problematic, and thus the environmental predictor (from the pair of highly correlated ones) which produced the highest Akaike information criteria (AIC) value on a univariate analysis was excluded. Initially we restricted GAMMs to a maximum of 5 knots to prevent over-fitting, however if GAMMs failed diagnostic checks we increased the number of knots until these checks were satisfactory. For the spatial smoothers in the models we used the default settings in the mgcv package 63 to estimate the number of knots required. When performing GAMMs, minimum adequate models were selected by backwards selection, using K-fold cross-validation, following 18,62 .
To establish the isotopic niche among periods with the stable isotope data we applied the recent metric SIBER (Stable Isotope Bayesian Ellipses in R), which is based on a Bayesian framework that confers a robust comparison to be made among data sets concerning different sample sizes 64 . The area of the standard ellipse (SEAc, an ellipse having a 40% probability of containing a subsequently sampled datum) was adopted to compare female and male isotopic values and their overlap in relation to the total niche width (both groups combined), and a Bayesian estimate of the standard ellipse and its area (SEA B ) was used to test whether females' isotopic niche is narrower than males' isotopic niche (i.e. p, the proportion of ellipses in female birds that were smaller than in male individuals; see ref. 64 for more details). All the metrics were calculated using standard.ellipse and convexhull functions from SIBER implemented in the package SIAR (Stable Isotope Analysis in R 65 ); under R 45 . All data are presented as mean ± SD, unless otherwise stated. Results were considered significant at P ≤ 0.05.