Study area and study species
Our study was carried out at Raso Islet (16°36’ N, 24°35’ W), Cabo Verde, a flat and inhabit islet integrated in the Integral Natural Reserve of Santa Luzia (Vasconcelos et al. 2015). Boyd’s shearwater is a subspecies of little shearwater (c. 5,000 pairs), belonging to the lherminieri complex, breeding in the archipelago of Cabo Verde (BirdLife International 2020; Semedo et al. 2020). It is currently classified as “Least Concern” in the IUCN Red List; however there are some signs of decline owing to the impacts caused by introduced species (BirdLife International 2020). This small-sized pelagic seabird (~160 g) is an endemic subspecies of Cabo Verde, and it is the nearest counterpart of the Macaronesian shearwater, which breeds in Azores, Madeira, Selvagens and Canary Islands (BirdLife International 2020). It is a winter breeder, and like other Procellariiformes, lays a single egg each breeding season. Briefly, adults arrive at the colony in August-September to prospect and defend their breeding borrow, females lay the egg in January-February, which hatches about 50 days later (mid-March) and the chick is fed approximately for 60 days, leaving the nest between the last half of May and the first half of June (Zajková et al. 2017).
GPS deployment and sample collection
From March to April 2018 and 2019, mini-GPS loggers (nanoFixTMGeo & Geo+, PathTrack Ltd., UK) were attached to the four central tail feathers of breeding adults, using TESA© tape (Wilson et al. 1997). Each logger together with the tape did not exceed 4 g weight, representing less than 3% of adults’ body mass (~160 g), the standard and recommended threshold to not compromise individuals’ foraging abilities (Phillips et al. 2003). The body mass of tracked adults did not differ prior (162.0 ± 18.3 g) and after (170.0 ± 25.8 g) carrying the device (paired t-test, t9=1.591, p=0.146). GPS deployment did not last more than 5 minutes, and adults were returned to the respective nests. Each logger was programmed to record each geographical position every 10 minutes (~140 locations per day). During deployment sessions, some breast feathers were collected for molecular sexing (see Table S1 for more details), while during logger retrieval, a blood sample (~0.8 ml) was collected from the brachial vein, centrifuged to separate plasma from red blood cells (RBC), and both blood partitions were kept in ethanol (70%) until preparation for stable isotope analysis (SIA). Four tags were successfully retrieved in 2018 (3 males and 1 female) and 24 devices (12 males and 12 females) in 2019. Prey samples were collected along the breeding seasons of 2018 and 2019 for subsequent SIA. Main prey groups were created according to prey type (squid or fish), life-stage (larval or adult), and distribution in the water column (epipelagic or mesopelagic).
Sample preparation and stable isotope analysis
Plasma was selected for carbon and nitrogen isotopic analysis, because its turnover rate corresponds approximately to the tracking period duration, i.e., around 5-7 days (Inger and Bearhop 2008), while RBC would reflect a larger timeframe of about 3-4 weeks (Bearhop et al. 2002; Cherel et al. 2005b). Nitrogen (δ15N; 15N/14N) isotopic values are commonly used as a proxy of predator’s trophic level, increasing about 2-5 ‰ at each trophic level (Minagawa and Wada 1986), while carbon (δ13C; 13C/12C) values are used as a habitat indicator, because it only suffers a slight increase (ca. 0-1 ‰) at each trophic level (Kelly 2000). Seabird plasma and prey muscle were dried during 24 and 48 hours, respectively, at 60 °C. Next, the samples were rinsed with a 2:1 chloroform: methanol solution to remove the overload of lipids that can deplete 13C values (Cherel et al. 2005c; Post et al. 2007). All samples were ground to a powder, weighted (~0.35 mg) in tin capsules, and analysed through an elemental analyser/isotope ratio mass spectrometry (EA/IRMS). The results were expressed using the standard δ notation, following the equation: δX = [(Rsample/Rstandard) –1] × 1,000, where X is 13C or 15N, and R is the ratio 13C:12C or 15N:14N, respectively. Rstandard values correspond to the Vienna PeeDee Belemnite (V-PDB) and atmospheric N2, for 13C and 15N respectively (Bond and Jones 2009). Replicate measurements of internal laboratory standards (acetanilide) indicate a precision of ± 0.2 ‰ for both carbon and nitrogen isotopic ratios. The C/N ratio was examined to verify if lipid removal was effective in all plasma and muscle samples.
Prey assemblage
Four main prey groups were created: epipelagic fish, mesopelagic fish, squid, and fish larvae. Epipelagic fish was comprised by three fish species which inhabit the epipelagic layer of the ocean (i.e., the upper 200 m of the water column): Tylosurus acus, Sardinella maderensis, Selar crumenophthalmus (mean ± SD: δ13C = –16.95 ± 0.48 ‰, δ15N = 9.72 ± 0.70 ‰, N = 16; C/N = 3.59 ± 0.06); mesopelagic fish was comprised by Myctophum affine and Hygophum sp. (mean ± SD: δ13C = –18.66 ± 0.51 ‰, δ15N = 10.23 ± 0.43 ‰, N = 9; C/N = 3.19 ± 0.06); squid included individuals of two different species: Hyaloteuthis pelagica and Callimachus rancureli (mean ± SD: δ13C = –17.02 ± 1.56 ‰, δ15N = 11.66 ± 2.18 ‰, N = 6; C/N = 2.89 ± 0.11); finally, fish larvae included fingerlings captured near surface during pelagic tours, identified as Ophioblennius sp. and Synodus saurus (mean ± SD: δ13C = –18.47 ± 0.24 ‰, δ15N = 8.38 ± 0.43 ‰, N = 10; C/N = 3.04 ± 0.04). All potential prey were identified to the lowest possible taxonomic level, weighted, and measured the body-length (for fish) or mantle-length (for squid). Prey were initially identified using local guides or catalogues and, specifically, squid were identified using their beaks (Xavier and Cherel 2009). Also, a small piece of muscle tissue of each species was collected to create a DNA reference collection (see Table S1), either to confirm the previous identification or to achieve a lower taxonomic level. These species are among the most abundant of each group within Cabo Verde archipelago and are generally ingested by local breeding seabirds (Carreiro personal communication).
GPS data analysis: behavioural classification and kernel estimation
To estimate missing locations and standardize sampling effort, GPS tracks were resampled by linear interpolation to exactly 15 minutes interval. Individual foraging trips were divided in short (< 1 day) and long (≥ 1 day), after inspecting trip duration frequency using an histogram (Fig. S1). To avoid potential disturbance caused by social interaction and flying movements during landing at the colony, a distance to colony filter of 1 km radii was applied, to discard those locations. Maximum distance to colony, latitude and longitude at the distal point of each foray were computed using several functions within trip R package (Sumner et al. 2020). The classification of at-sea behaviours was carried out using a combination of instantaneous flight speed and path sinuosity (calculated as the ratio of instantaneous flight speed given the speed between every third positions). Histograms of the frequency of these two variables showed adults were drifting on the water, i.e., resting, when flying speed was below 2 km h-1; intensive search, i.e., foraging behaviour, was assigned when the flight speed was between 2-10 km h-1 and path sinuosity was above 7; extensive search, i.e., relocating behaviour, was assigned when the flight speed was equal or above 10 km h-1 and path sinuosity was above 7; traveling behaviour was assigned when the flight speed was simultaneously above 2 km h-1 and below or equal to 7 (Fig. S2). After behaviour classification, the proportion of time spent on each behaviour was calculated for each foray within each individual.
We computed the 50% kernel UD contours to represent adult foraging areas (FA), and the 95% UD contours to describe adult home ranges (HR). Following the methods and R scripts described by (Lascelles et al. 2016) we calculated mean ARS zones radii for short and long foraging trips, and used those values as smoothing parameters (h) in the computation of Kernel UDs. A smoothing parameter of 4 km was used for short foraging trips, and a smoothing parameter of 8 km for long foraging trips. Kernel UD contours (95% and 50%), and respective areas, were calculated using the ‘kernelUD’ and ‘kernel.area’ functions within the adehabitatHR R package (Calenge 2006). The overlap of UD contours was calculated between sexes within short and long foraging trips through the ‘kerneloverlap’ function, using the Bhattacharyya's affinity (BA), under the adehabitatHR R package (Calenge 2006).
Environmental predictors and habitat suitability models
Monthly values of ocean (1) bathymetry (BAT, blended ETOPO1 product, 0.01° spatial resolution, m), (2) chlorophyll a concentration (CHL, 0.04° spatial resolution, mg m–3), (3) ocean mixed layer thickness (OMLT, 0.08° spatial resolution, m), (4) sea surface height above geoid (SSH, 0.08° spatial resolution, cm), and (5) sea surface temperature (SST, 0.08° spatial resolution, °C) were extracted within the foraging range of Boyd’s shearwater for March 2018-2019 and April 2019, and the mean raster was calculated for each environmental predictor. Variable 1 was downloaded from https://www.ngdc.noaa while variables 2-5 were downloaded from http://marine.copernicus.eu. Spatial gradients of all environmental predictors were calculated using an estimating proportional change within a 3 x 3 cell grid following Louzao et al. (2009). Bathymetry gradient (BATG) identifies the presence of oceanic topographic features, such as seamounts or shelf-breaks (i.e., slope areas); the gradient of CHL (CHLG) and SST (SSTG) can be used as a proxy of oceanic fronts, while OMLT gradient (OMLTG) indicate the change level on the mixed layer thickness, and consequently, the depth of the thermocline which drives the abundance and distribution of marine prey; the gradient of SSH (SSHG) could help identify the occurrence of mesoscale eddies. All environmental predictors were rescaled to the coarsest spatial resolution (i.e., 0.08°) and extracted for each GPS location, before running habitat suitability models. All computations were conducted under several functions within raster R package (Hijmans et al. 2020).
Habitat suitability models were computed separately for males and females for long and short foraging trips, i.e., four modelling exercises. Tracking data from 2018 and 2019 was jointly analysed given the low sample size of tags retrieved in 2018 and similar foraging range and distribution between years (Fig. S3). Prior to habitat modelling, all environmental predictors, and respective gradients, were inspected for multicollinearity. Multicollinearity was tested using the variation inflation factor (VIF) and Pearson correlation coefficients (r > 0.6) (Table S2), under the usdm R package (Naimi 2017). Testing for collinearity issues enables to account only with non-redundant variables, avoiding model overfitting and inflated errors (Zurell et al. 2020). Ensemble Species Distribution Models (ESDM; Marmion et al. 2009) were conducted using all GPS locations (presence data; Fig. S4) through the ‘ensemble_modelling’ function from the SSDM R package (Schmitt et al. 2017). We tested 7 modelling techniques: Artificial Neural Network (ANN), Classification Tree Analysis (CTA), Generalized Additive Models (GAM), Generalized Linear Models (GLM), Multiple Adaptive Regression Splines (MARS), Random Forest (RF), and Support Vector Machine (SVM). Each model algorithm was computed ten times using a 10-fold cross-validation procedure, using 70% of all data set for model calibration and the remaining 30% grid squares as random test for model validation (Araujo et al. 2005; Marmion et al. 2009; Zurell et al. 2020). Model goodness of fit was examined using the area under the receiver operating characteristic (ROC) curve (AUC). Models were classified excellent when AUC > 0.90, good when 0.80 < AUC < 0.90, reasonable when 0.70 > AUC < 0.80, and not acceptable when AUC < 0.70 (Araujo et al. 2005). The relative importance of environmental predictors to the probability of occurrence of each sex within short or long foraging trips, was given by the average contribution calculated from all models.
Isotopic niche and mixing models
Stable Isotope Bayesian Ellipses in R (SIBER) package was used to obtain the isotopic niches of females and males. The standard ellipse area corrected for small sample size (SEAC) and considering 40% of all observations, was used to calculate the isotopic niche width of each sex, and to calculate niche overlap between sexes, using the ‘maxLikOverlap’ within the SIBER package (Jackson et al. 2011). Bayesian standard ellipse areas (SEAB) were calculated using 10,000 iterations of Markov-chain Monte Carlo (MCMC) simulation to test for the probability of group 1 (e.g., females) being smaller than that of group 2 (e.g., males), using the rjags R package (Plummer et al. 2019). The contribution of each prey group for Boyd’s shearwater diet was estimated using Bayesian mixing models within the simmr package (Parnell and Inger 2016). This package provides a wide range of new functions for a more accurate and realistic diet estimation, including the incorporation of prior diet information to the model (Parnell and Inger 2016); however, to the best of our knowledge there is no data about Boyd’s shearwater diet, so no prior information was added to the model. Trophic Discrimination Factors (TDFs) are needed to accurately run the isotopic mixing model and these are often tissue-specific, species-specific and diet-specific, meaning that they may vary according to consumer’s species and its diet, and the tissue analysed (Phillips et al. 2014; Jenkins et al. 2020). To our best knowledge there are no TDF available for Boyd’s shearwaters; hence, we opted to use the average values of fractionation between prey and plasma of 3 seabird species, from captive experiments (Barquete et al. 2013; Jenkins et al. 2020). So, we used a TDF of − 0.44 ‰ and + 2.10 ‰ enrichment for carbon and nitrogen, respectively. Although these species have different feeding ecologies than Boyd’s shearwaters, we believe that these average TDFs are the most adequate for our mixing model exercises, considering other values available in the literature for seabirds (Hobson and Clark 1992; Bearhop et al. 2002; Cherel et al. 2005a; Sears et al. 2009; Chiaradia et al. 2014; Ciancio et al. 2016). A standard deviation of ± 1.0 ‰ was used to account for possible differences on enrichment factors between species. Before running the models a simulation method proposed by Smith et al. (2013) was used to inspect the feasibility of the isotopic mixing polygons. The sensitivity analysis (using 1500 iterations) applied to mixing polygons indicated that adult isotopic signatures were within 95% of the simulated mixing regions (probability ranges: 0.21 to 0.61 for males, 0.43 to 0.63 for females), validating our models (Fig. S5). We ran a mixing model for each sex computed using the function ‘simmr_mcmc’ from the simmr R package (Parnell and Inger 2016).
Statistical analysis
Generalized linear mixed models (GLMMs) using the appropriate family or linear mixed models (LMMs) following a normal distribution, were used to test the effect of sex on adult trip parameters and at-sea foraging behaviour, separately for short and long foraging trips: (1) trip duration, (2) maximum distance to colony, (3) total distance travelled, (4) percentage of time spent foraging, (5) relocating, (6) resting, and (7) travelling, (8) latitude and (9) longitude coordinates at the maximum distance to colony. All models were run using sex as a fixed factor, while the bird identity (i.e., individual) was included as a random factor to avoid pseudo-replications. Years were pooled together due to the lower sample size recorded in 2018 (3 males and 1 female). Mixed models were conducted with the lme4 (Bates et al. 2015) and lmerTest (Kuznetsova et al. 2017) R packages. Sexual differences on the size of FA and HR within short and long foraging trips were tested using t-tests for independent samples, when data followed a normal distribution, or Mann-Whitney tests, when data did not follow a normal distribution. All response variables were tested for normality, homoscedasticity, and log (total distance travelled), square root (maximum distance to colony and trip duration), or arcsine (time spent foraging, resting, and relocating) transformed whenever necessary. Throughout the results values are expressed as mean ± SD. All analyses and modelling computations were performed using the R software ver. 4.0.0 (R Core Team 2020) and the significance level was set at p ≤ 0.05.