Niche segregation mechanisms in marine apex predators inhabiting dynamic environments

Understanding the mechanisms that allow the coexistence of species is key to preserve full ecosystem functioning. In dynamic environments, the study of ecological niches faces the complexity associated to the three dimensionality of the habitat and requires information that reflects such heterogeneity. Within this context, this study intends to identify the segregation mechanisms behind the co‐occurrence of five phylogenetically related pelagic birds by applying a functional perspective based on seabirds' vertical ranges and prey availability features such as depth and body size.


| INTRODUC TI ON
The niche concept has been a major theme in ecology, mostly influenced by the definition stated by Hutchinson (1957), who described it as a "n-dimensional hypervolume of environmental states within which a species is able to survive." However, even the classical definition may create a dichotomy that affects the way in which the entire concept is approached (Chase & Leibold, 2003).
In fact, much of the confusion surrounding the term results because no distinction is made between the responses of organisms to their environment and the effect of organisms on their environment (Chase & Leibold, 2003;Peterson et al., 2011). Under such circumstances, two different niches should be discerned: the Grinnellian and the Eltonian niches (Devictor et al., 2010;Soberon, 2007). The Grinellian niche (Grinnell, 1917) describes the response of the species to a given set of non-interactive variables, while the Eltonian niche (Elton, 1927) describes the biotic interactions and resource-consumer dynamics through trophic variables. Despite the attempts at synthesis and unification (Chase & Leibold, 2003), the complementary concepts of the environmental niche (sensu Grinnell, 1917) and trophic niche (sensu Elton, 1927) serve as basis to assess the ecological and biogeographical dis-/similarities of species and contribute to the understanding of their distribution and diversity (Broennimann et al., 2011;Soberon, 2007). Indeed, some of the main segregation mechanisms driving species coexistence are known to occur either by means of trophic (MacArthur, 1968;Tilman, 1982) or environmental niche partitioning (Chesson, 2000). The spit up of the niche provides, therefore, an excellent opportunity to advance in the identification of segregation mechanisms, which is specially poorly understood in complex and dynamic environments such as marine ecosystems.
The study of environmental niches in marine ecosystems, for instance, faces the difficulties inherent in understanding the three dimensionality of the habitat, where physical and ecological processes occurring below the surface layer, such as subsurface thermal structure or subsurface primary production, have been found to be highly important (Kuhn, 2010;Scott et al., 2013). On the other hand, trophic niche analyses require the consideration of complex predator-prey interactions, that result from the trade-off between the energetic cost of seeking prey and the foraging profitability obtained from successful events (MacArthur & Pianka, 1966;Pyke, 1984). Traditionally, predator-prey interactions have been studied using a predominantly taxonomic approach (speciesspecies perspective), although functional characteristics related to the species' role, such as biological traits, have been proved to affect predator-prey interactions (Spitz et al., 2014). In 3D dynamic environments, addressing community ecology from functional traits perspective may be specially relevant, since prey depth (i.e. the depth at which prey is available), together with body size or energy density, can be key to understand predators' prey preferences (Boyd et al., 2017;Lambert et al., 2018;Spitz et al., 2018;Waggitt et al., 2018). Similarly, predators' physiology and morphology may also play an essential role in the foraging process by defining the metabolic cost of living or the diving capabilities of species (Spitz et al., 2012(Spitz et al., , 2014. Among marine top predators exploiting dynamic environments, seabirds are considered one of the most diverse taxa, with several evidences indicating niche differentiation among its members (Phillips et al., 2011). Most studies, however, have been conducted during the breeding season, when the competition for resources is particularly intense and segregation mechanisms are more likely to arise (Mancini et al., 2014;Navarro et al., 2013). As a result, limited knowledge exists about the potential segregation mechanisms outside the breeding season (Thiebot et al., 2012), despite being the period in which animals migrate or disperse to favourable foraging areas, aggregating in highly productive regions and coinciding with a high number of migratory species (Grecian et al., 2016).
Investigating and explaining niche segregation outside the breeding season could help complement management strategies involving all life stages, and thus, bigger efforts should be made in studying the non-breeding areas of seabirds.
The Bay of Biscay (BoB hereafter) represents such an example, as numerous seabird species stopover there during their feeding migrations attracted by a highly diverse and abundant community of small pelagic fishes (Astarloa et al., 2019). It conforms, therefore, an exceptional biogeographical area to test segregation hypothesis and provides an incomparable opportunity to understand the mechanisms that allow the coexistence of protected species. In fact, many pelagic birds visiting the BoB are protected under different international agreements, such as the Bird Directive (Council Directive 2009/147/EC) and the Convention for the Protection of the Marine Environment of the North-East Atlantic (OSPAR Convention), among others. Despite the conservation efforts, seabirds are still one of the most threatened groups, comprising rapid declining populations and critically endangered species (Croxall et al., 2012;Dias et al., 2019).
Effective conservation and management measures require the identification of variables shaping species' niches , 2020 and that is why, understanding the role that prey and environment play in species distribution turns so necessary. Within that context, the present study focused on disentangling the assembly rules of pelagic birds by addressing the following research questions: 3D environments, environmental niche, niche differentiation, pelagic birds, prey availability, trophic niche, vertical segregation 1. Does niche segregation occur among wide-ranging species during their feeding migrations?
2. Does this niche segregation occur in the environmental niche or in the trophic niche?
3. Do prey and environmental conditions on the vertical dimension influence the trophic and environmental niche of species?
To answer these questions, we modelled species environmental and trophic niches and developed a conceptual framework with the most plausible segregation hypotheses (Figure 1). This way, we aim to contribute to the understanding of protected and endangered species coexistence and provide insights on the 3D niches of species that may help advance into their management and conservation.

| Data collection: integrated monitoring schemes
Data on pelagic birds, prey species and environmental predictors were collected through the multidisciplinary oceanographic survey JUVENA, which takes place every September in the BoB by means of two different research vessels, Ramon Margalef and Enma Bardan (R/V RM and R/V EB, hereafter). Since 2013, different components of the pelagic environment (i.e. plankton, fish, megafauna, physical oceanography, marine litter) are monitored, although its main aim is to assess the population of juvenile European anchovy Engraulis encrasicolus. The sampling strategy is based on parallel transects perpendicularly arranged to the coast and spaced at 15 nautical miles, whose offshore and along-coast extension changes from year to year depending on the distribution of the European anchovy (see details in Boyra et al., 2013;Louzao et al., 2019). A schematic workflow of the entire analytical process is described in Figure 2. 2.1.1 | Pelagic birds Five phylogenetically related shearwaters were selected based on their conservation status and their high diversity and abundant records: the Cory's shearwater Calonectris borealis, the great shearwater Ardenna gravis, the near threatened sooty shearwater A. grisea, the Manx shearwater Puffinus puffinus and the critically endangered Balearic shearwater P. mauretanicus. Sightings of these species were recorded aboard R/V RM by a team of three experienced observers (2 at a time) that followed the line-transect methodology (Buckland et al., 2001). This methodology is conducted within the distance sampling framework to estimate seabirds' densities and requires the collection in the field of at least the radial distance of each observation (Heinemann, 1981), the angle of the cluster sighting with respect to the trackline (estimated with an angle meter), time of observation, species composition and group size. In addition, the behaviour of observed species (e.g. attraction) as well as environmental descriptors affecting the detectability of species (e.g. Beaufort sea-state, visibility, glare intensity or observation conditions) were collected in order to account for response bias (when animals react to the presence of the platform) and perception bias (when observer miss animals because their visibility is compromised), respectively. Sampling effort was performed during daytime, at a constant speed and under Beaufort sea-state conditions ≤6 and it was geographically located every minute with the vessel GPS ( Figure 3).

| Pelagic prey
Acoustically based biomasses were included in the analysis only for those fish species that were considered, based on bibliography, part of the diet of the seabirds (see Appendix S1: year (Appendix S1: Table S1.2). The collected acoustic data by both vessels were processed in the positive strata by layer echo integration using an ESDU (Echo integration Sampling Distance Unit) of 0.1 nmi and categorized into ten layers of varying depths according to the year (Appendix S1: Table S1.2). In parallel, mid-water trawls were performed to assign the eco-traces to species (identification purposes) and to obtain the necessary biological data (length, weight and age) to convert the acoustic back-scattered energy to fish abundance. Finally, abundance in number of individuals was multiplied by the mean weight, obtaining biomass estimates (tonnes) per age, length and depth interval (Boyra et al., 2013).

| Oceanographic and geographic environment
Oceanographic data were collected in both vessels using a CTD profiler. For each transect, a minimum of three profiles were performed (coastal, continental shelf and oceanic waters) measuring the water column from the surface (10 m, first available data) up to 200 m depth ( Figure 4). Temperature (T), salinity (SAL) and density data were directly inferred from CTD casts. Horizontal fields of these three variables were estimated every 5 db from the vertical profiles using optimal statistical interpolation scheme (Gomis et al., 2001) on a spatial grid with regular node distances of 0.15 × 0.15° latitudelongitude covering all the study area (see details in Appendix S2).
Then, geostrophic velocities (GV) were derived from the interpolated density fields following Rubio et al. (2009). Secondarily, the depth of the maximum temperature gradient (DTG, as a proxy of ocean mixed layer depth and water column stability indicator), the maximum temperature gradient (MTG, as a proxy of the strength of the water column stratification) and the sea surface temperature gradient (SSTG, as an important predictor for seabirds distributions) were estimated from temperature fields as described in Louzao et al. (2019).
In order to characterize the geographic environment, depth  (Table 1); in fact, areas of strong spatial gradient may correspond to areas where internal waves generate, which can promote an increase in primary production and small preys' availability according to Scott et al. (2010). For that, a spatial moving window of 3 × 3 cells was used in which the spatial

| Detection functions of seabirds' sightings
In line transects, it is assumed that the likelihood of detecting animals becomes smaller as the distance to the observer increases.
To account for that bias, distance sampling analyses were applied, which mainly consist on fitting a detection function to the observed distances in order to estimate the proportion of animals missed by the observer (Thomas et al., 2002). To do so, sightings of seabirds were first filtered (by removing sightings with attraction behaviour) to avoid the response bias generated when animals react to the presence of the platform; in fact, fisheries discards can attract large feeding flocks and lead to misunderstand the distribution and abundance of seabirds (Valeiras, 2003). Due to the low number of sightings per year of some species (Buckland et al., 2001), small (Balearic and Manx shearwaters) and large species (Cory's, great and sooty shearwaters) were grouped together based on their size (see Appendix S3). Once we defined the groups, the 5% of the sightings F I G U R E 4 Overview of the annual acoustic sampling accompanied by the CTD casts collected along the transects. The dash lines correspond to R/V EB, while the solid lines refer to R/V RM, which usually cover the inner and the outer section of the transects, respectively. The red circles correspond to the locations of the CTD casts. Isobaths of 200 m, 1,000 m and 2,000 m are indicated detected at the largest distances were truncated to delete outliers (i.e. by setting the truncation distance, w) (Buckland et al., 2001) and analysed using multiple covariate distance sampling (Marques & Buckland, 2004). Hazard rate and half normal functions with no adjustments were then fitted in each of the groups using the ds function from the Distance package (Miller, 2020). As covariates, only those descriptors related to the effort were considered, that is Beaufort sea-state, visibility, cloudiness, glare intensity, observation conditions and year (see details in Appendix S4: Table S4.1), which were introduced in the detection function as factor and selected by means of forward selection (Appendix S4: Table S4.2) until the lowest Akaike information criterion (AIC hereafter) was obtained (Guisan & Zimmermann, 2000;Sakamoto et al., 1986). From here, the probability of detecting an animal (Pa) was estimated, which multiplied by the truncation distance (w) provided the effective strip half width (ESW = Pa * w). ESW can be defined as the perpendicular distance in which the missing detections equal the recorded detections and serves to estimate the area effectively covered, when considering both observation sides and transect length (A = ESW * 2 * L).
In such estimations, seabirds' behaviour (on flight versus. on water) can be an important aspect to be considered, since flying individuals can lead to overestimate densities when they move faster than the observation platform (Buckland et al., 2001). In our case, most seabirds were recorded on flight (90%-95% of individuals in all species), and hence, little bias was expected between areas with flying and sitting individuals; overall overestimation, on the other hand, was considered negligible, since obtaining absolute abundances was out of our scope.

| Processing of explanatory variables
In order to understand the ecological niches of seabirds in the BoB, environmental and trophic relationships were modelled. However, the differences in the sampling coverage of seabirds, prey and environmental data did not allow for a homogenous prediction. To solve that, prey and environmental data were first processed to obtain continuous fields of explanatory variables covering the study area and then categorized by depth and size to address the multidimensionality and functionality of the environmental and trophic niches ( Figure 2; step 3).

| Prey fields
Since prey selection may be also conditioned by other factors related to prey availability such as body size or depth at which prey is available, the original biomasses (tonnes) of the seven prey species were categorized using 10 cm length classes (Lambert et al., 2018):  TA B L E 2 Prey species categorized by size and depth and above the DTG). It is known, however, that it is an important resource among procellariids (Watanuki & Thiebot, 2018), so in order to test its relevance in our seabird community, we included as predictor the biomass of the Mueller's pearlside comprised between the surface and 200 m depth (Table 2). Although seabirds are not able to dive so deep, this estimation was used as a proxy of the biomass available at dusk and dawn, that it is when the pearlside migrates close to the surface and aggregates at about 20-40 m (Kaartvedt et al., 1998), becoming available to seabirds.

| Environmental fields
As with prey fields, those environmental variables collected at different depths (i.e. T, SAL and GV) were also vertically analysed and classified into surface and deep layers

| Biologically meaningful vertical range selection
Since pelagic birds show different diving abilities, the vertical range they exploit may also differ. Without obtaining in situ diving records, we cannot test that hypothesis, but we can determine the vertical range that best explains seabirds' density patterns (biologically meaningful vertical range, hereafter). For that, two models per species were fitted: one using the conditions given by the environmental and prey variables in the surface layer, and the other using the conditions given by the same variables but in the deep layer ( Figure 2, step 4). Since this step required 3D data, the environmental variables comprised T, SAL and GV (Table 1), while the prey data included the first three axes of the PCA (explaining the 70% of the variability; Appendix S7). In all cases, individual density surface models were fitted from the previous detection function analyses using the dsm package (Miller et al., 2019). The number of individuals per unit effort was fitted by means of generalized additive models (GAM), assuming a negative binomial distribution with a probit link function (after testing with tweedy and quasipoisson families). Degrees of smoothness were limited to fit unimodal response curves and restricted to three (Bruge et al., 2016) to avoid overfitting (Burnham & Anderson, 2003). All variables were standardized to have a mean of zero and a standard deviation of one (Zuur et al., 2007) to make the comparison of effect sizes easier. The most plausible model was selected based on the lowest AIC (Guisan & Zimmermann, 2000;Sakamoto et al., 1986).
In addition, a literature survey was conducted to contrast the obtained results. For that, we focused on published biologging studies providing the average diving depth reached by the species, as it refers to the vertical range more regularly exploited (Appendix S8:

| Identification of segregation mechanisms
To identify niche differentiation mechanisms within the pelagic bird community, we modelled separately the environmental and trophic niche of species by integrating the data at their biologically meaningful vertical range (Figure 2, step 5), estimated an environmental and trophic overlap index for each pairwise species (Figure 2, step 6) and developed a conceptual framework with the most plausible segregation hypotheses (Figure 1).

| Environmental and trophic niche modelling
For the environmental niche, models combining environmental variables (both 3D and 2D, and the biomass of the Mueller's pearlside were used. In both cases, individual density surface models were fitted from the previous detection function analyses using the dsm package (Miller et al., 2019). For each species, the number of individuals per unit effort was fitted by means of generalized additive models (GAM), assuming a negative binomial distribution with a probit link function (after testing with tweedy and quasipoisson families). Degrees of smoothness were limited to fit unimodal response curves and restricted to three (Bruge et al., 2016) to avoid overfitting (Burnham & Anderson, 2003). All variables were standardized to have a mean of zero and a standard deviation of one (Zuur et al., 2007) and subsequently analysed by means of Spearman's rank correlation coefficient to identify highly correlated (|r| ≥ .6) pairwise predictors (Thiers et al., 2014). The most plausible model was selected based on the lowest AIC (Guisan & Zimmermann, 2000;Sakamoto et al., 1986).
When models differed in 2 units of AIC (ΔAIC < 2), they were considered statistically equivalent and the one with a smaller number of variables was chosen following the parsimony principle (Arnold, 2010). Once the most plausible "environmental" and "trophic" niche models were defined, seabirds' densities were predicted per year over the standard grid (latitudinal range: 43.2-48°N; longitudinal range: 1-8°W, 0.1° spatial resolution) (Figure 5a,b).

| Environmental and trophic niche overlap
In order to assess the degree of environmental and trophic niche segregation, an overlap index was calculated between pairwise species based on Ballard et al. (2012). In the case of environmental niche overlap, we first estimated the mean density for the 2013-2017 period based on the predictions obtained from the environmental niche modelling and selected only those cells containing the highest 95% of the mean density in order to avoid very low values (Figure 5c,d). After that, we assessed the degree of environmental overlap (Figure 5e) by dividing the number of cells where both species were present (i.e. cells containing both species) by the total number of cells where either species was present (i.e. cells containing one species or both) (Ballard et al., 2012). This led to a total number of 10 overlap values, derived from the pairwise combination of 5 species (C 5,2 = were posteriorly standardized so that values ranged between 0 and 1.

F I G U R E 5
Example of how environmental overlap was estimated, showing the environmental niche (standardized densities) of sooty (a) and Balearic (b) shearwaters, followed by their respective (c, d) abundant areas (after keeping only the highest 95% of mean densities) and their environmental overlap (e) To estimate the degree of trophic niche overlap, the same procedure was followed based on trophic niche modelling results.
Finally, the conceptual framework displayed in Figure 1 was developed, which describes the plausible hypotheses that may arise from the pairwise comparison of environmental and trophic overlap indexes. These main hypotheses comprise two clear segregation mechanisms, defined as environmental segregation (high trophic overlap but low environmental overlap) and trophic segregation (high environmental overlap but low trophic overlap), and two additional situations described as specific conditions (low trophic and low environmental overlap) and alternative mechanisms (high trophic and high environmental overlap). Specific conditions would refer to any situation explaining why the pairwise species found in that section do not overlap (e.g. isolation, breeding, specialization), whereas alternative mechanism hypothesis would try to find out how those species can coexist in a situation of both high trophic and environmental overlap.

| Sightings and detectability of seabirds
The most frequently observed species was the great shearwater (944 sightings), followed by the sooty (293 sightings  Note: Number of sightings refers to the final number obtained after having removed the 5% of the data detected at the largest distances.  (Table 3 and Appendix S9: Figure S9.1).

Species
In the case of small shearwaters, the function with the lowest AIC was a half normal function, with year, Beaufort and general conditions as covariates (Table 3, Appendix S9). For Cory's and great shearwaters, the best function was a hazard rate with year as covariate, while in the case of the sooty shearwater, the best function was a half normal with year as covariate (Table 3, Appendix S9).

| Biologically meaningful vertical range
According to the test conducted to identify each species' biologically meaningful vertical range, we found that the density patterns of Cory's, great and Manx shearwaters were better explained by the explanatory variables of the surface layer (10 m), while the density patterns of sooty and Balearic shearwaters were better explained by the explanatory variables integrated over the deep layer (above DTG) (Table 4). Obtained results were in agreement with the results expected from the literature review (Appendix S8: Table S8.2), in which the average diving depth records indicated a surface diving behaviour for Cory's, great and Manx shearwaters (with an average depth of 1.7, 3.2 and 5.7 m, respectively) and a subsurface diving performance for sooty shearwater (average depth 12.3 m). The only exception was the Balearic shearwater, as the average diving depth found in the literature did not agree with the obtained results from the modelling of the biologically meaningful vertical range (Appendix S8: Table S8.2).

| Environmental and trophic drivers
Trophic niche models showed a high preference for PCA1 (smallmedium fish species) in Balearic and Manx, for PCA3 (big fish species) in Cory's and for the Mueller's pearlside in great and sooty shearwaters, highlighting three main groups ( Figure 6). Environmental models, although more diverse, showed a more homogenous pattern in terms of variables' importance, suggesting a more balanced contribution of the environmental variables ( Figure 6). Nevertheless, some similarities could be found too; large shearwaters (Cory's, great and sooty), for instance, were found to rely moderately on salinity (SAL) and depth (DEP), whereas small shearwaters preferred those variables linked to land or shelf-break closeness (DCO and DHSEL).
Temperature (T), although present in most of the models, was found to be of low importance for all the species (Figure 6). In general, environmental niche models provided higher percentages of deviance explained, that ranged between 15% and 71%, while trophic niche models showed a deviance explained between 5% and 59% (see Appendix S10).

| Niche differentiation mechanisms
Overlap indexes obtained from previous environmental and trophic niche models were displayed following Figure 1 in order to assign the main four segregation hypotheses to each pairwise species ( Figure 7). In this way, the pairs composed by Balearic-great,

F I G U R E 6
Relative importance of the variables included in the environmental and prey-based models, integrated at the biologically meaningful vertical range of each species. Prey variables: the main three axes of the PCA (PCA1, PCA2 and PCA3) plus Mueller's pearlside (MAV). Environmental variables: temperature (T), salinity (SAL), geostrophic velocity (GV), depth of the maximum temperature gradient (DTG), maximum temperature gradient (MTG), distance to coast (DCO), distance to shelf break (DSHEL), depth (DEP) and depth gradient (DEPG) Manx-great and Balearic-sooty were assumed to segregate environmentally (Figures 7a and 1a), while Manx-Cory's shearwaters were presumed to segregate through trophic niche partitioning (Figures 7d and 1d). The category classified as specific conditions ( Figure 1c) was assigned to the Cory's shearwater, as it was present in all pairwise combinations showing low trophic and low environmental overlap, including Cory's-Balearic, Cory's-sooty and Cory's-great pairs (Figure 7c). On the other hand, Manx-Balearic, Manx-sooty and sooty-great pairs showed the opposite pattern, suggesting that these pairwise species coexist through an alternative mechanism in conditions of high environmental and high trophic overlap (Figures 7b and 1b).

| D ISCUSS I ON
Unlike other marine predators, air breathing seabirds are limited in prey accessibility due to their anatomy and their diving capabilities. Thus, considering the processes that concentrate prey close to the surface, prey size or depth at which seabirds can fish is essential. In this study, the incorporation of such elements has enabled us to conclude that (a) wide-ranging species coexist through environmental and trophic niche partitioning, (b) species respond differently to prey and oceanographic conditions on the vertical dimension (potential vertical segregation) and (c) phylogenetically and morphologically closer species (e.g. sooty-great or Manx-Balearic) show more similarities in their trophic and environmental niches.
These major findings were mainly extracted from the conceptual diagram, resulting from the modelling of environmental and trophic niches. Previous co-occurrence analyses conducted in the area already described some of the results found here, such as the environmental overlap between Manx-sooty or the environmental dissimilarity in Cory's-Balearic shearwaters (Astarloa et al., 2019).
However, the way in which this approach was addressed (i.e. by considering both environmental and trophic niches, prey depth and size, conditions on the vertical dimension) provided more detailed information on species assemblage and revealed four different scenarios resulting from niche overlap patterns that could not have been identified otherwise.
The two clearest scenarios were the environmental and trophic segregation (Figures 1a, 7a and Figures 1d, 7d, respectively). The former, detected in Balearic-great, Manx-great and Balearic-sooty pairs, can be understood with the spatial distribution described for these species. In fact, Balearic and Manx shearwaters are known to occupy primarily coastal waters , whereas sooty and great shearwaters show preference for shelf and oceanic waters, respectively . When they stopover in the BoB, Balearic and Manx stay closer to the coast, while great and sooty shearwaters exploit offshore areas, leading to a non-overlap pattern in their environmental niche.
On the other hand, trophic segregation was only identified in the case of Manx-Cory's pairwise species. This segregation mechanism can be explained by the results given by diet-based studies, that suggest that Cory's shearwater feeds on Atlantic mackerel and horse mackerel (Paiva et al., 2010), while the Manx shearwater mainly relies on clupeids (e.g. herring, sprat) (Thompson, 1987). This single association also indicated that overall, the community of pelagic birds in the BoB was characterized by a low trophic segregation, which can be due to the generalist behaviour of most species, known to F I G U R E 7 Comparison of trophic versus environmental overlap indexes that display pairwise species in sections of (a) high trophic overlap but low environmental overlap, (b) high trophic and high environmental overlap, (c) low trophic and low environmental overlap and (d) high environmental overlap but low trophic overlap. The most plausible segregation hypothesis explaining each section can be found in Figure 1 take advantage of available pelagic feeding resources (Bicknell et al., 2013). In fact, even the critically endangered Balearic shearwater, with a potentially more restricted foraging range compared to the remaining wide-ranging pelagic birds, feeds on the main pelagic resources of the BoB (e.g. mackerel, horse mackerel, anchovy, sardine) (Meier et al., 2017). However, it must be mentioned that, despite not having found strong evidence of trophic segregation in the conceptual diagram, the modelling of trophic niches already revealed some differences in the seabird community. Indeed, great and sooty shearwaters were associated to Mueller's pearlside, Balearic and Manx shearwaters to small-medium fish species (PCA1), and Cory's shearwaters to big fish species (PCA3), meaning that clustering prey species based on functional characteristics can help uncover subtle differences on the trophic preferences of species (Lambert et al., 2018).
The remaining two scenarios (Figure 1b,c) could not be explained by environmental and trophic requirements and instead, further information on the biology of the species was required to be untangled. Low patterns of both trophic and environmental overlap, for instance, were related to the reproductive behaviour of the Cory's shearwater (involved in all associations with low overlap). In fact, it was the only species that was breeding (they breed in the north west of the Iberian Peninsula) at the time the study was conducted (Munilla et al., 2016). During this period, seabirds act as central place foragers (Orians & Pearson, 1979), which means they have to make a balance between selecting productive areas (to obtain enough food supplies for their chicks and themselves) and performing not too long trips (to come back to the colony to feed the chicks  (Mori & Boyd, 2004;Navarro et al., 2013). In addition, our results also seem to indicate that vertical segregation mainly occurs between closely related species, such as the sooty and great shearwaters (Ardenna genus) or the Manx and Balearic shearwaters (Puffinus genus). Phylogenetically related species are expected to be ecologically similar (Losos, 2008), so we could additionally suggest that under very similar niche conditions, segregating in the vertical dimension is the only way to reduce competition. Indeed, the similarities in the niche of these species did not only arise in niche overlap analysis but also when identifying environmental and trophic drivers.
Vertical segregation hypothesis, however, cannot be confirmed without diving depth records, although it seems plausible based on the high agreement found between the biologically meaningful vertical ranges defined here and the average depths recorded by data loggers (see Appendix S8). Indeed, the only exception in which the average depth recorded in previous studies was not in agreement with our results was the case of the Balearic shearwater (see Appendix S8). This may be due to the low sample size (only diving data from one individual was obtained in Aguilar et al. (2003), but see Meier et al. (2015)) or due to the contrasting oceanographic and prey conditions between the Mediterranean Sea and the Atlantic Ocean (i.e. oligotrophic versus eutrophic conditions, respectively).
However, with the available data, no significant conclusion can be made, and we can only acknowledge that further research is needed to elucidate the diving behaviour and the vertical range of the Balearic shearwater.
Advancing in the understanding of endangered, threatened and protected (ETP) species is critical. The pelagic bird community of the BoB, characterized by a highly migratory behaviour, is protected under multiple international agreements. By combining data collected from integrated ecosystem surveys and habitat modelling, we have proved that studying seabirds outside their breeding areas can also provide useful results about their ecological niches. Indeed, we have contributed to understand the underlying environmental and trophic drivers of both environmental and trophic niches, needed to identify critical feeding grounds and high biodiversity areas in the context of marine spatial planning. Assessing the degree of overlap between such areas and anthropogenic pressures (e.g. fishing bycatch) could be, for instance, a potentially useful step to conduct in the near future that will undoubtedly help advance in the conservation of these species.

ACK N OWLED G EM ENTS
The JUVENA survey was funded by the 'Viceconsejería de