The study was conducted between July and August in 2013, 2018 and 2019, in the Parco Nazionale dell’Arcipelago di La Maddalena (Italy) located between Sardinia and Corsica, and hosting 200-500 breeding pairs of Scopoli’s shearwaters (from our surveys) distributed in about 12 islands. This archipelago is at minimum 7 km far from Lavezzi islands (France), hosting 350-400 pairs (Péron et al., 2018). We sampled individuals breeding on two small islets of this archipelago and located only 7 km apart, a very close distance for such a wide-ranging species (Figure 1). These two sub-colonies, Barrettini (41°17’3.59’’N, 9°24’5.96’’E) and Spargiotto (41°14’59.93’’N, 9°19’25.01’’E), host the higher number of pairs of Scopoli’s shearwaters per island in the archipelago, namely 40-50 and 60-80 pairs respectively. Barrettini and Spargiotto are besides the most eastern and western sub-colonies in this archipelago, respectively (among the ones with more than 10 breeding pairs). There, shearwaters breed in relatively easily accessible burrows, enabling the capture and identification of many breeding individuals (between 24 and 37 each year on each islet), as well as the recapture of individuals equipped with GPS devices. Scopoli’s shearwaters typically mate in late April, lay eggs in late May, and incubate them for about 55 days. They feed mostly on plankton and small pelagic fish (Grémillet et al., 2014) but they can also exploit discards from fishing vessels (Cecere et al., 2015; Cianchetti-Benedetti et al., 2018). During incubation, both partners undertake long lasting foraging trips (1-18 days, 7 days on average; Cecere et al., 2013), while upon chick-rearing (in July – October), they alternate few long trips lasting 6-20 days with shorter trips (1-4 days; Cecere et al., 2014).

During early chick-rearing period, the parent attending the nest was captured during night-time directly in its nest, ringed with a metal ring or identified if already ringed, and equipped with a GPS logger attached on the back feathers (see De Pascalis et al., 2020 for more details). The GPS loggers were set to record a position every 10 min. Similarly to De Pascalis et al. (2020), we removed positions within 5 km from the sub-colony. We relied on manual foraging trip identification for 2013 (from Cecere et al., 2014), considering a foraging trip as an at-sea excursion lasting more than five hours, starting and ending at the colony. For 2018 and 2019, we followed the same method as in De Pascalis et al. (2020), and automatically identified foraging trips following Lascelles et al. (2016) and its tripSplit function, considering only trips lasting at least 5 hours, with inner-and outer buffers of 5 and 20 km respectively. We ensured the correct discrimination by visually inspecting each foraging trip in QGIS (v.3.16.10, QGIS.org, 2022). Overall, we relied on data from 32 individuals tagged in 2013, 44 in 2018, and 15 in 2019; 16 individuals were tracked in two years and one individual in three years. Sex-ratio (proportion of females) was well balanced in all years and sub-colonies (0.50 on average, range: 0.41-0.63). See Table S1 for a detailed description of sample size per year, sex, and sub-colony.

The entire statistical procedure is described in detail in Figures S1, S2. The first step consisted in deriving Kernel Density Estimates (KDEs), as this allowed then (1) to test whether individuals are site-faithful across trips and years and thus for us to retain only independent trips, (2) to test whether there is spatial and temporal heterogeneity in the spatial data, (3) to derive ICAs, i.e., delineate areas within the KDEs that most of the tracked population uses and so is deemed important for the whole population. All analyses were performed in R v4.1.0 (R Development Core Team, 2021), with functions from the R package track2KBA (Beal et al., 2021b).

For all subsequent analyses, we considered both 50% and 75% Kernel Density Estimates (KDEs). The choice of 50% KDEs (i.e., the smallest area on which the probability of finding an individual is 50%) was driven by the fact that this level best represents the core areas of activities for our population of Scopoli’s shearwaters (see Figure S3, following BirdLife International, 2009) and is the usual level used when proposing conservation areas (Beal et al., 2021b). However, we also used 75% KDE as this may lead to more unified (i.e., less divided) ICAs than 50% KDEs (see Discussion), and it might be advisable to focus not only on the core areas of activities (typically related to foraging areas) but also areas used for travelling back-and-forth from the colony and those used for rafting. All KDEs were derived at the trip level, with the smoothing factor href (determined using the tripSummary and findScale functions; see Table S2 for the values href found for each sub-sample), and a resolution of 0.75 km. As already discussed in (Beal et al. 2021b), the choice of smoothing factor and resolution is crucial as these will affect the size and the precision of the KDEs and ICAs. It is thus important to select them based on the biology of the species and the aim of the conservation actions. Here, we chose href (instead of the mag value provided by the findScale function for instance) because our aim was to derive ICAs encompassing not solely foraging sites, but any important area used by shearwaters for foraging, resting, and travelling. This is also why we retained all sampled locations (and not only the ones corresponding to foraging behaviour; except locations within 5 km of the colony, to release some of the emphasis on the colony for this central place forager). We chose a resolution of 0.75 km based on preliminary analyses and graphical assessment (see Figure S4): this resolution precisely enough fit shearwaters’ movements here, while limiting computational time.

We initially relied on 397 foraging trips performed by 73 individuals. We first tested if, for a given year and sub-colony, individuals were faithful to foraging sites as this would lead to autocorrelation among consecutive foraging trips of a same individual (function indEffectTest). In all but two combinations of years and sub-colonies, individuals showed site fidelity across foraging trips (p > 0.1 for Barrettini in 2013 and 2019, with both 50% and 75% KDEs; p < 0.003 for all other cases). We thus down-sampled trips until site fidelity was no longer significant: we retained either all trips (Barrettini in 2013 and 2019), every second trip (Spargiotto in 2019), every third trip (Barrettini in 2018), every fourth trip (Spargiotto in 2018), or only one trip per individual (Spargiotto in 2013). Since some individuals were tracked in more than one breeding season, we ensured that individuals were not site faithful across years. To do so, we followed the same principle as the one proposed in the indEffectsTest function. We estimated, separately for each individual sampled several years (N=17), the overlap between its KDEs derived at the trip level and differentiated between and within year overlaps. We then combined the overlap values of all 17 individuals and compared between and within year overlaps with a bootstrap Kolmogorov-Smirnov test (1000 iterations, p > 0.07). This allowed us to consider the foraging trips performed in different years as independent. In the end, we analysed 199 foraging trips, performed by 91 unique individual-year entities.

To detect spatial heterogeneity, we tested whether individuals breeding in the two sub-colonies spatially segregated while foraging at sea, to assess the need of implementing the rest of the analyses at the sub-colony level. We used the function indEffectTest to test whether within- and between-sub-colony overlaps of individual KDEs were different. We expected here the within-sub-colony overlaps to be significantly greater than between-sub-colony ones, which would indicate that the area used by individuals (i.e., KDEs) from the same sub-colony were more similar than areas used by individuals from different sub-colonies. We performed this test for each year separately. To detect temporal heterogeneity, we used the same method and tested whether within- and between-year overlaps of individual KDEs were different. We performed this test for each sub-colony separately.

We identified ICAs following the methodology initially developed by Lascelles et al. (2016), and recently improved by Beal et al. (2021b). This approach implied three steps (see Figures S1, S2 for a graphical summary). First, based on this restricted dataset (i.e., with down-sampled independent trips) we aimed to identify KDEs at the foraging trip level (function estSpaceUse). Second, we ensured that sampled individuals were representative enough of their sub-colony by assessing the representativeness of the samples from each year and each sub-colony separately using the function repAssess and the KDEs (with 600 iterations, because the stability of estimates was reached from this number onwards for the smallest sample size, see Figure S5). Third, we identified ICAs with the findSite function, which delineates areas that are used by a certain percentage (at least 10%) of the tracked individuals given the location of the KDEs, and so are deemed important for the whole population. This percentage is given by the argument ‘thresh’ in the findSite function. Here, we retained the default value for ‘thresh’, which adapts the percentage threshold depending on the representativeness of the sample (20% threshold for 70-80% representativeness, 12.5% threshold for 80-90% representativeness, and 10% for representativeness above 90%; see Appendix S1 in Lascelles et al., 2016; Beal et al., 2021b). We chose not to provide a total population size for the function findSite, as this parameter is only used for visualizing the number of individuals in the entire population predicted to use the identified ICA (instead of the default percentage of the tracked population). In other terms, the choice of population size does not affect in any way the contour of the ICA itself, which was the object of focus here. We could then extract the contour of the ICAs (Figure S2). This 3-step procedure – deriving KDEs, assessing representativeness, and deriving ICAs – was applied for each year and each sub-colony separately given the presence of spatial and temporal heterogeneity in foraging areas (see previous paragraph and results).

The function indEffectTest applied to individual KDEs (paragraph 2.5) allowed to test for the presence of spatial and temporal heterogeneity in foraging sites. Yet, whether these differences in overlaps of KDE within and between sub-colonies or years (see results) cascade to actual low overlap of ICAs (i.e., the parts of the KDEs supposedly used by the majority of the population), between sub-colonies or years, remained to be explored. Our first aim was then to quantify potential differences in ICAs resulting from accounting for a single sub-colony (spatial heterogeneity) or a single year (temporal heterogeneity). To quantify differences due to spatial heterogeneity, we estimated, for a given year, the overlaps between ICAs of both sub-colonies (e.g., Barrettini vs. Spargiotto in 2013). To quantify differences due to temporal heterogeneity, we estimated, for a given sub-colony, the overlap between ICAs obtained for different years (e.g. 2013 vs. 2018 for Barrettini). The proportion of overlap Overlap_i,j between pairs of ICAs (i, j) was estimated as:

A referring to the areas of the ICAs of i or j, their intersection (i∩j), and their union (i∪j).

Our second aim was to propose a method to derive a single and reliable ICA while accounting for the spatial and temporal heterogeneities in the data (method A) and show why this is better than ignoring these heterogeneities (method B). Once one has sampled several years and/or multiple sub-colonies, one can identify a global ICA by either (A) merging the areas found for each year and sub-colony separately, or (B) directly identifying an ICA from the larger dataset comprising all years and sub-colonies (i.e., not accounting for possible temporal or spatial heterogeneity). We compared then the differences in total areas and the overlap between the ICAs obtained with these two (A, B) alternative approaches, 1) for each year, pooling both sub-colonies, 2) for each sub-colony, pooling the three years, and 3) for the three years and both sub-colonies pooled.

Considering the tracks of each foraging trip (Figures 1 and S6), there was a clear segregation of used areas between individuals breeding in the two sub-colonies: individuals from Spargiotto moved mostly south-west of the colony, while individuals from Barrettini moved towards north-east. Accordingly, overlap of KDEs was greater within than among sub-colonies (p < 0.001 for each year and both for 50% and 75% KDEs). In other words, individuals breeding in the same sub-colony used more similar areas during foraging trips than individuals breeding in different sub-colonies. Overlap of KDEs was also greater within than among years for individuals breeding in Spargiotto (at least for core areas: p < 0.001 for 50% KDE, and marginally for larger areas: p = 0.07 for 75% KDE). Yet, this was not the case for individuals breeding in Barrettini (p > 0.5 for both 50% and 75% KDEs). Given this spatial and (partially) temporal heterogeneity in the spatial data, we implemented the rest of the analyses for each sub-colony and year separately.

Representativeness was generally high (Table S3), and above 72% for all the year- and sub-colony specific sub-samples, except Spargiotto in 2019 when considering 50% KDEs. This latter sample, which was the smallest of all years and sub-colony combinations (Table S1), reached only 45.9% of representativeness and was thus discarded for the subsequent corresponding ICA estimations. Foraging trips from Spargiotto in 2019 were still considered for ICAs estimation based on 75% KDE (with a representativeness of 72.8%).

In accordance with the significant spatial and temporal heterogeneity found in section 3.1, the overlap between ICAs from different sub-colonies or years was rather low (below 0.5; Figure 3). In particular, the overlap between ICAs of the two sub-colonies was low each of the three years, whether considering ICAs based on 75% KDE (Figures 2, 3) or ICAs based on 50% KDE (Figures S7 and 3). ICAs’ shape and location were rather different for each of the three years of the study, and this was particularly striking for 2019 compared to the other years (Figures 2 and S7). Accordingly, the overlap between ICAs from different years were relatively low (Figure 3): from 0.16 to 0.49 for ICAs based on 75% KDE, and from 0.01 to 0.40 for ICAs based on 50% KDE.

Combining data from different sub-colonies and/or years before deriving ICAs (method B) usually led to smaller ICAs than when merging the ICAs obtained for these years and/or sub-colonies (method A; Table 1; Figures 2, 4, S7 and S8). These differences were particularly striking when ignoring both the spatial and temporal heterogeneity: pooling data from both sub-colonies and all years led to ICAs that were more than 70% smaller than when merging the ICAs from each of the sub-colonies each year (Figures 4C, S8C and Table 1). One exception was the comparison of ICAs based on 50% KDEs for both sub-colonies in 2013: the area obtained by merging the ICAs of Barrettini and Spargiotto was 16% smaller than the ICA obtained by pooling initial data (Table 1). Yet, this peculiarity did not hold for ICAs based on 75% KDEs (Table 1). In 2019, the number of foraging trips from individuals of Spargiotto was so low compared to Barrettini (10 versus 64 trips) that the potential ICA derived from data of Spargiotto only was almost not represented at all when combining both sub-colonies (Figure 4C).

We showed that sampling only one year, only one sub-colony, or simply not accounting for spatial and temporal heterogeneity in large multi-scale datasets led to a very partial representation of ICAs, especially so when sample sizes greatly differed across groups. For example, in 2019, we obtained only 10 foraging trips for Spargiotto and 64 for Barrettini after accounting for foraging site fidelity and removing correlated trips. These 10 and 64 foraging trips were representative enough of their sub-colony to derive ICAs (72.8% and 94.3% representativeness respectively when considering 75% KDE). When ignoring the spatial heterogeneity and deriving an ICA on the 2019 dataset (method B), we obtained an ICA (based on 75% KDE) encompassing nearly the entire ICA found for Barrettini alone, and almost none of the ICA found for Spargiotto alone (Figure 2C). In 2013, when sample sizes were more balanced across sub-colonies, we did not find this effect (Figure 2A). When sample sizes differ across groups (sub-colony and year here, but it could be sex for example), we advise to merge ICAs derived for the different groups, as it allows for usually more extended and representative ICAs. Indeed, up to 75% of the ICAs area was lost when not accounting for these sources of heterogeneity, biasing the location of ICAs toward those found in years or sub-colonies with the most foraging trips. In only one occasion, the ICA derived by pooling initial data was larger than the ICA obtained by merging group-specific ICAs (Table 1, Figure S7A). Surprisingly, this larger ICA also included an area just north of Barrettini, which was not included in any of the sub-colony or year-specific ICAs (based on 50% KDE), and in that respect may not be that relevant conservation-wise if only a minority of the tracked breeders used it each year (less than 10%, the chosen minimum threshold for ICA definition here). The methodology we propose, by allowing the derivation of reliable ICAs from unbalanced sample sizes, should be of broad use for most conservation actors, which often do not have access to spatially and temporally balanced samples. Another possibility to account for spatial heterogeneity in the data is to use a predictive modelling approach, when individuals are known to select similar physical habitat features when foraging (e.g., similar bathymetry, Ceia et al., 2015), an hypothesis that would remain to be tested across different sub-colonies.

Testing for within-colony spatial heterogeneity is a fundamental step before planning any conservation action, but this may be far from easy. The main obstacle for such investigation is practical: it can be highly challenging to capture and tag wild animals at their natural breeding sites. For example, not all breeding sites, beaches, or cliffs might be accessible, and/or one might focus on individuals located at the edge of a colony to minimize disturbance to the rest of the group. In that respect, the Scopoli’s shearwater colony studied here is ideal: pairs breed in well separated burrows between rocks on the ground, distant from each other by few meters, in easily accessible islands. As in previous studies (Bogdanova et al., 2014), we advocate whenever possible to collect spatial data from various locations of the colony, to be able to generalize the detected exploited areas. Yet, if collecting spatial data from several parts of a colony is not feasible, we propose several other options. For example, if individuals tend to follow the route taken by conspecifics when departing from the colony or closely located rafts (Weimerskirch et al., 2010; Thiebault et al., 2014), observing the departure and returning angles of individuals from different sub-colonies might inform on the possibility that members of different sub-colonies segregate at sea. If one can only sample one part of a colony but knows its extent, it is possible to use a modelling approach to test for spatial segregation. Indeed, one can test whether the observed distribution corresponds to that of simulated tracks ignoring between-sub-colony competition for example (Ito et al., 2020), or, on the contrary, whether the observed distribution is biased towards one area, away from the rest of the colony (similar to colony-level simulations in Aarts et al., 2021). If spatial segregation of foraging grounds is then suspected, care should be taken when deriving ICAs. The combination of accurate resource availability maps and modelling approaches accounting for competition between- but also within-colonies (e.g. Wakefield et al., 2011; Aarts et al., 2021), with a back-up verification on existing spatial data, may open exciting perspectives for accurately deriving ICAs from partial datasets.

Our results highlight the importance of dynamic placement of ICAs given that environmental changes affect the distribution of the sentinel species and their prey (Perry et al., 2005; Poloczanska et al., 2016; and see Cashion et al., 2020 for a torough discussion on the importance of dymamic ICAs in marine ecosystems). Indeed, we observed among-year differences in the areas exploited by Scopoli’s shearwaters from La Maddalena Archipelago, which is not surprising given their foraging habits. Scopoli’s shearwater feed mostly on plankton and small pelagic fishes (Grémillet et al., 2014), two resources whose distribution can vary between years (Marty et al., 2002; Saraux et al., 2014). Shearwaters also exploit discards from fishing vessels (Cecere et al., 2015; Cianchetti-Benedetti et al., 2018), whose locations are naturally dynamic and rather stochastic (changing within few minutes or hours) compared to the temporal scale of shearwaters foraging trips (1-20 days, Cecere et al., 2014). Yet, Scopoli’s shearwater breeding in the eastern sub-colony showed no significant difference in used areas (KDE) across the three studied years, unlike individuals breeding in the western sub-colony. One hypothesis could be that resources in the northern and eastern part of the colony were more stable across years than those located west, yet we unfortunately do not have any empirical evidence to formally test that. Other populations of the same species also showed both consistency (Grémillet et al., 2014) and differences (Courbin et al., 2018) in the areas exploited across years, which may suggest either a stochastic effect of the sampled year, or geographical differences in resources availability. Such discrepancies among populations also highlight the importance of testing for temporal and spatial heterogeneities when multi-year or sub-colony data are available, and not solely rely on findings from other populations or years. It is to be noted that we exclusively used data from the chick-rearing stage, and not from the incubation stage, because individuals are more spatially constrained when they must combine self-feeding with regular trip to the colony to provision offspring.

Among-year differences in exploited areas may arise from environmental stochasticity and longer-term environmental changes. Long-term investigation of variation in foraging behaviour in seabirds are still scarce, but one striking example, which also investigates sub-colony segregation, is the 15-year study of European shags in Scotland (Bogdanova et al., 2014). This study also provides a methodological framework enabling to test how many years are necessary to accurately represent the population-wide used areas. Our approach here that complements that previous study, shows how to best account for temporal heterogeneity when one does not have access to enough data to fully and certainly represent the population at-sea distribution (yet with representative enough within-year sub-sample). With three years of data and slightly contrasting sample sizes (Table S1), we cannot exclude that the observed temporal heterogeneity is the results of stochasticity. Yet, longer term datasets will certainly become more common in the coming years and will allow for a thorough assessment of longer-term temporal heterogeneity and account for that when deriving ICAs.

Our study focused on the effect of sub-colony and temporal heterogeneity in the implementation of ICAs. Yet, spatial segregation can also occur between sexes (Catry et al., 2006), age classes (Pelletier et al., 2014; Haug et al., 2015; Pettex et al., 2019), or personality types for example (Krüger et al., 2019). Here, the sex-ratio in our sample was well balanced (0.5 on average), and sampling was performed randomly regarding individuals’ age and personality (random selection of individuals among the captured ones). Depending on the population or species, any source of heterogeneity can be accounted for with the same procedure as the one proposed here (see Figure S1 for a summary).

One aspect that should also be considered is the focus of the conservation actions: either on the foraging areas exclusively, or on any exploited areas, even the ones used for travelling between patches or rafting. Here, we chose to rely on any used areas (i.e., all GPS positions obtained, whether they corresponding to foraging, travelling, or resting events), which led to a specific emphasis really close to the colony (Figures 2, 4, emphasis reduced by the withdrawal of locations within 5 km of the colony, see methods). This was motivated by the fact that shearwaters may suffer from disturbance both while foraging and travelling (e.g., Rodríguez et al., 2015; Dierschke et al., 2016; Genovart et al., 2018). Yet, depending on the species and the conservation context, it is possible to select only specific locations (e.g., foraging vs. non-foraging, core-foraging areas vs. home-range). We therefore focused on ICAs based on 75% KDE (larger scale used areas), but nonetheless compared them to ICAs based on 50% KDE (considered as core used areas) which are more classically used. The ICAs based on core used areas (50% KDE) were naturally smaller than ICAs based on larger used areas (75% KDE) but they were also interestingly more divided (see Figure S9 for a graphical representation of these ICAs). ICAs based on larger used areas enabled to integrate corridors of movement between core used areas (as typically seen in Spargiotto in 2018 or Barrettini in 2019, Figure S9). Therefore, for conservation actions aiming at preserving any used sites, 75% KDE might lead to more ecologically relevant ICAs.

Our study system also exemplifies the importance of deriving reliable ICAs in international contexts (here Italy and France, see also for instance Beal et al., 2021a; Davies et al., 2021), which may imply international agreements, and, here, would require a careful assessment of foraging distribution of individuals breeding in the nearby Lavezzi archipelago. Furthermore, we focused here on one species. Yet, using data from other species sharing the same foraging habitats (for instance here Yelkouan shearwaters Puffinus yelkouan) would also constitute a critical step to apply conservation actions benefiting multiple species.