The relative importance of biological and environmental factors on the trophodynamics of a pelagic marine predator, the blue shark

Marine ecosystems have been significantly altered by the cumulative impacts of human activities. Pelagic sharks have become vulnerable to increases in mortality rates caused by fishing. The decrease in number of these top predators could have substantial cascading effects on wider marine communities. Concerns about these potential impacts


Introduction
Human-induced environmental changes are notably affecting marine systems by homogenizing species assemblages, with negative effects on biodiversity and ecosystem functioning (Worm et al., 2006). Pelagic shark populations are considered to be particularly vulnerable to large-scale pelagic fisheries (Baum et al., 2003;Juan-Jordá et al., 2022), because of their slow growth, late age of maturity, and low reproductive rate (Myers et al., 2007;Myers and Worm, 2005). In the Atlantic Ocean, several large shark species have declined by more than 75% between 1986 and 2001 (Baum et al., 2003). Similarly, in the Mediterranean Sea, there have been declines in the abundance of several species of sharks, as well as reductions in shark diversity and body size (Ferretti et al., 2008(Ferretti et al., , 2010. This is of special concern because these species play an ecologically important role, exerting a top-down influence on the ecosystem structure (e.g. Stevens 2000;Bascompte et al., 2005;Dulvy et al., 2014).
Concerns about the viability of shark populations and the critical need for the effective management of these populations has increased the interest in understanding pelagic shark trophic ecology, in order to predict the sensitivity of these species and the potential consequences of environmental changes (Bird et al., 2018;Carrier et al., 2012;Coll et al., 2013;Cortes, 1999). Despite considerable progress in this ecological knowledge, anticipating such outcomes remains challenging due to the complexity of extrinsic and intrinsic factors that shape shark trophic ecology , and due to our fragmented understanding induced by spatially and temporally restricted studies (Bird et al., 2018). Furthermore, although the diet of marine predators has been traditionally analysed using stomach content analysis (SCA; Hyslop, 1980), sharks often have empty stomachs, and the prey that are recovered are often skewed towards those that are difficult to digest (Hyslop, 1980). SCA requires large sample sizes to quantify the trophic ecology of organisms, which are difficult to obtain for endangered sharks, and especially at a large scale (Hussey et al., 2012;Navarro et al., 2014). As an alternative to SCA, the use of ecological tracers such as stable isotopes analysis (SIA) of Nitrogen (δ 15 N) and Carbon (δ 13 C) have been used to investigate the trophic ecology of these pelagic predators (Bird et al., 2018;Hussey et al., 2012;Madigan et al., 2021;Shiffman et al., 2012). This approach is based on the fact that δ 15 N and δ 13 C values are transformed from dietary sources to consumers in a predictable manner (Layman et al., 2012). δ 15 N values are a proxy of trophic position, while δ 13 C values show little change with trophic levels, but instead are a useful indicator of the primary source of carbon (Layman et al., 2012).
The aim of this study was to investigate the relative effect of biological factors, spatial location, and environmental features on the spatial trophodynamics of blue shark (Prionace glauca), assessed by analysing the δ 15 N, δ 13 C and trophic level values from individuals sampled along a wide geographic area between south of the Canary Islands in the Atlantic Ocean, and the north-western Mediterranean Sea. The blue shark is a clear example of a globally distributed pelagic shark that is highly impacted by fisheries, both directly and as bycatch (Clua, 2020;Druon et al., 2022;Queiroz et al., 2016). In fact, due to overfishing, the blue shark is considered Near Threatened at the global scale and Critically Endangered in the Mediterranean Sea by the International Union for Conservation of Nature (IUCN) . Previous trophic studies with the blue shark indicated spatial variation of its trophic ecology with diet changes, such as the relative importance of pelagic fish, cephalopods and crustaceans Navarro et al., 2017;Rosas-Luis et al., 2017). However, the extent to which this spatial variation is related to environmental or intrinsic factors has been largely overlooked for this pelagic marine predator. We investigated the relative importance of biological (sex and body size), environmental (sea-surface height anomaly, chlorophyll-a, pelagic productivity, sea-surface temperature, bathymetry and human impact), and geographic factors in determining the isotopic and trophic position variability in the study area using deviance partitioning analyses (Bocard et al., 1992;Navarro et al., 2016). These factors were chosen as they have previously been described as important determinants of the spatial and trophic ecology of the blue shark (Druon et al., 2022;Estupiñán-Montaño et al., 2019;Loor-Andrade et al., 2017).

Sampling and stable isotope analysis procedures
Between the years 2017 and 2018, a total of 180 blue sharks were collected between south of the Canary Islands in the Atlantic Ocean, and the north-west of the Mediterranean Sea (Fig. 1). All blue sharks were captured as by-catch by Spanish commercial longline vessels targeting swordfish. From each captured individual, the sex (visually) and the fork length (FL, error of ±0.1 cm) were recorded, and a small portion of muscle (without skin) was collected from the ventral region and stored at − 20 • C until stable isotopic determination.
Urea content from the muscle samples was removed prior to the stable isotope determination to prevent potential bias in the estimation of the stable isotopic values (Kim and Koch, 2012;Shiffman et al., 2012).
To extract urea, muscle samples were placed in 2 ml vials and rinsed twice in 1.5 mL of deionized water. A rinse entailed sonication for 15 min and then decanting the supernatant. After this process, all muscle samples were freeze-dried and powdered. Then 0.3-0.4 mg of each sample was packed into tin capsules, and the stable isotopic determination was performed at the Laboratory of Stable Isotopes of the Estación Biológica de Doñana (LIE-EBD, Sevilla, Spain). Encapsulated samples were combusted at 1020 • C using a continuous flow isotope-ratio mass spectrometry system by means of a Flash HT Plus elemental analyser coupled to a Delta-V Advantage isotope ratio mass spectrometer via a CONFLO IV interface (Thermo Fisher Scientific). The isotopic composition was reported in the conventional delta (δ) per mil notation (‰), relative to Vienna Pee Dee Belemnite (δ 13 C) and atmospheric N 2 (δ 15 N). Replicate assays of standards routinely inserted within the sampling sequence indicated analytical measurement errors of ±0.1‰ and ±0.2‰ for δ 13 C and δ 15 N, respectively. The standards used were EBD-23 (cow horn, internal standard), LIE-BB (whale baleen, internal standard) and LIE-PA (razorbill feathers, internal standard). These laboratory standards were previously calibrated with international standards supplied by the International Atomic Energy Agency (IAEA, Vienna). The C:N ratio of all tissues was always lower than 3.5‰, and hence, no correction of the δ 13 C values was required to account for the presence of lipids in muscle samples (Logan et al., 2008).

Trophic position of blue sharks
The estimated trophic position (TP) of blue sharks was calculated based on their δ 15 N values according to the algorithm proposed by Vander Zanden and Rasmussen (2001): where δ 15 N consumer and δ 15 N basal were, respectively, the δ 15 N values of blue shark individuals and the δ 15 N values of the phytoplankton at the location where each blue shark was collected. δ 15 N basal values were obtained from a global ocean model of phytoplankton Somes et al., 2017). Δδ 15 N was the trophic discrimination factor (3.4%, following Pethybridge et al., 2018) and the TP basal was 1 for phytoplankton .

Environmental and human impact variables
We selected 6 variables as descriptors of environmental variability (Table 1): (1) surface chlorophyll-a concentration (Chl-a, mg⋅m − 3 ), (2) sea-surface temperature (SST, • C), (3) sea surface height anomaly (SSHa, Fig. 1. The study area and capture locations (white dots) of blue sharks sampled for this study. Blue shark drawing by À lex Mascarell. m), (4) the Ocean Productivity available to Fish (OPFish) based on the tracking of productivity fronts (Chl-a horizontal gradient), (5) bathymetry (BAT, m) and (6) human impact (IMP). Chl-a values were acquired from MODIS-AQUA satellite sensor as monthly values at 1/24 • resolution (https://oceancolor.gsfc.nasa.gov). Chl-a and SST are oceanographic variables related to ocean productivity and therefore have been described as a good proxy for prey availability for pelagic sharks. SSHa provides information about the presence of eddies, oceanographic structures that could also play a role in upwellings and nutrient supply to the ecosystem and have been related to the distribution of pelagic predators (Godø et al., 2012;Vandeperre et al., 2014). SST and SSHa levels were extracted from the EU-Copernicus Marine Environment Monitoring Service global model (https://marine.copernicus.eu/acc ess-data). The SST and SSHa data were linearly interpolated from the original grid at 1/12 • resolution, to the habitat grid at 1/24 • resolution. PEL was taken from OPFish as monthly values (https://fishreg.jrc.ec.eur opa.eu/fish-habitat). PEL is a plankton-to-fish index that represents the distribution of Ocean Productivity available to Fish (OPFish, or potential fish production) (Druon et al., 2021), which provides insights about the prey availability for blue sharks (Druon et al., 2022). BAT was obtained from ETOPO1, a global relief model provided by NOAA's National Geophysical Data Center (NGDC, http://www.ngdc.noaa.gov/mgg/g lobal/). Different bathymetric regimes are generally characterized by different levels of productivity (e.g., open ocean vs. continental shelf). Moreover, marine predators that feed on deep-water prey may appear to have elevated trophic positions relative to species foraging predominantly on shallow dwelling prey, due to isotope fractionation via heterotrophic degradation of sinking or suspended particulate organic matter. Human impact data was a composite measure of different anthropogenic factors such as fishing, pollution, and climate change (it ranges from 0 to 1, with 0 representing no impact and 1 maximum impact; Halpern et al., 2008). Human impacts are expected to lead to changes in the structure of food webs and therefore in TPs of fish species that occur in open fishing areas, such as blue sharks (Dell et al., 2015;Smith et al., 2011).
For consistency, the spatial resolution of all environmental variables was 1/24 • . Preliminary analyses considered increasingly broader resolutions (1 • , 2 • , and 3 • ) in order to account for the large areas that blue sharks can cover (Campana et al., 2011;Queiroz et al., 2012). However, these resolution changes did not affect model fit (see Supplementary Material S1). As stable isotope values in the muscle of large sharks potentially reflect their diet during the 6 months prior to the sampling (Logan and Lutcavage, 2010), for each individual and oceanographic variable we extracted the average of the six monthly values previous to sampling. To avoid problems related to multicolinearity, we performed Pearson correlation matrix analysis using all explanatory variables. Two pairs of variables had correlation coefficients higher than 0.70 (Chl-a and SST, r = − 0.74, and Chl-a and PEL, r = 0.77) and were not included together in GLM models (see Statistical analyses).

Statistical analyses
To assess the factors influencing variation in δ 15 N and δ 13 C and TP of blue sharks, we fitted generalized linear models (GLM, normal error distribution and identity-link functions) with different combinations of predictor sets. We started by fitting separate models for specific predictor sets (hereafter 'partial models'), including geography (sampling area -Canary Islands, Southern Portugal and western Mediterranean Sea-, latitude and longitude and their interaction 'latitude × longitude'), biological influences (FL and sex) and environmental variables (SST, SSHa, CHL-a, PEL, BAT and IMP). Both the linear and quadratic form of the environmental variables and FL were considered in models to test for potential parabolic trends (i.e. increases/decreases in δ 15 N, δ 13 C or TP values at intermediate values). We subsequently constructed models considering twofold combinations of predictors sets (namely, geography + biology, geography + environment, biology + environment) and a final 'combined model' (geography + biology + environment). We included year as a control fixed factor in all models, to account for potential non-independence of the data across years. All GLMs were conducted using a multimodel inference approach based on Akaike's information criterion adjusted for small sample sizes (AICc) (Burnham, 2004). Final models were calculated as averaged values of sub-models receiving high support (ΔAICc ≤2) (Burnham, 2004). Model-averaged parameter estimates and their standard errors were calculated, along with variable weights, which indicate the relative importance of each variable in the average model. Multimodel inference was implemented in the R software with the functions 'dredge' and 'model.avg' from the 'MuMIm' library. We checked for spatial autocorrelation in model residuals using Moran's I tests with 100 permutations using the function 'moran.mc' from package 'spdep'' (version 1.1-8) in R. We used variation partitioning analyses based on R 2 to quantify the independent and joint contribution of different target predictor sets. Only records with complete data (N = 173) were used in all analyses.

Models for δ 15 N and TP values
Predictions from the models combining all the variables fit observed δ 15 N (R 2 = 0.64) and TP (R 2 = 0.63) well, and better than any of the partial models investigated (Table 2). Both biological (FL and sex) and specific environmental variables were relevant in explaining the variation of δ 15 N values (Table 2; Fig. 2) and TP (Table 3; Fig. 2). The independent contributions of biological and environmental predictor sets accounted for 8.2% and 8.9% of variation in δ 15 N, and 6.2% and 12.3% of variation in TP, respectively (Fig. 2). Among biological variables, FL had high support in both partial and combined models of δ 15 N and TP. δ 15 N and TP increased in a curvilinear fashion with FL, with higher δ 15 N and TP at intermediate FL values (Fig. 3). δ 15 N and TP values were higher for females according to combined and partial models of δ 15 N and TP, respectively (Tables 2-3 and Fig. 3). For the environmental variables, SST and PEL had high support as explanatory variables in combined models of δ 15 N (Tables 2-3). In contrast, Chl-a received higher support in combined models of TP (Table 3). Note that SST and PEL were highly correlated with Chl-a and were not allowed in the same models  (Fig. 3). TP increased in a curvilinear fashion with Chl-a (Table 3 and Fig. 3). Geographic variables were significantly related to δ 15 N and TP, but most of the geographic variables contribution to variation explained was due to joint effects with biological and environmental variables (Fig. 2). No spatial autocorrelation was observed in model residuals (δ 15 N: Moran I test statistic = 0.016, p-value = 0.18; TP: Moran I test statistic = 0.050, p-value = 0.10).

Models for δ 13 C values
Neither geographic, biological, nor environmental variables alone were relevant to explain variation of δ 13 C values. The performance of partial models considering those variable sets separately was poor ( Table 4). The combined model fit the data better than any partial model investigated, but model fit was moderate (R 2 = 0.14). The independent contribution of environmental and geographical variables to the combined models was almost null (Fig. 2). Among the biological variables, only sex received high support in the combined model (Table 4). δ 13 C values were higher in females (Fig. 3). No spatial autocorrelation was observed in model residuals (δ 13 C: Moran I test statistic = − 0.05, pvalue = 0.94).

Discussion
This study examined the relative impact of several biological, environmental, and geographic predictors on the stable isotopic values and trophic position of the blue shark, over a large area between south of the Canary Islands in the Atlantic Ocean, and the north-western Mediterranean Sea. Statistical models provided new information on the main factors that influence the trophodynamics of this pelagic predator, showing that biological, environmental, and geographic variables all play a relevant role. According to models that considered stand-alone predictor sets, environmental variables followed by geographic zones, appeared to better fit δ 15 N and TP variation than biological variables. However, when the independent contribution of the three predictor sets was partitioned, only environment and biological predictors appeared to be relevant. This indicates that most of the geographic variation observed was accounted for these two groups of factors. No clear pattern Table 2 Model-averaged results for δ 15 N of the blue shark according to environmental, intrinsic variables and geographic variables. The table indicates the relative importance (weights, i.e. selection probability in the top-ranked models, ΔAIC≤2) of each variable, as well as full average parameter estimates for predictors receiving support (weight>0). Significant average parameter estimates (P < 0.05) are in bold. Model fit (R 2 ) and the number (N) of the top-ranked models are shown. N = 173.

Variables
Partial models b Combined model b  was observed for δ 13 C values. Among the biological variables, body size (body length measured with the fork length, FL) appeared to play a relevant role in δ 15 N variation and trophic position. However, contrary to previous studies that indicate a positive relationship between the body size and the trophic position of other pelagic predators (Cohen et al., 1993;Costa, 2009;Dhurmeea et al., 2020), we found that the higher δ 15 N or trophic position values were found at medium body sizes (see Martínez-Rincón and Acosta-Pachón 2022 for tuna species). Specifically, blue sharks with a medium body size (90-160 cm) had more enriched δ 15 N values, Table 3 Model-averaged results for trophic position of the blue shark according to environmental, intrinsic variables and geographic variables. The table indicates the relative importance (weights, i.e. selection probability in the top-ranked models, ΔAIC≤2) of each variable, as well as full average parameter estimates for predictors receiving support (weight>0). Significant parameter estimates (P < 0.05) are in bold. Model fit (R 2 ) and the number (N) of top-ranked models are shown. N = 173.

Variables
Partial models b Combined model b  occupying subtle higher trophic position values than larger individuals (>180 cm). In addition to body size, the sex also affected trophodynamics of the blue shark. In particular, females showed higher TP than males. These results suggest that blue sharks change their feeding behaviour ontogenetically as body size is a proxy of their age (Estupiñán-Montaño et al., 2019), and there is also a difference in feeding behaviour between females and males . Overall, blue sharks have been shown to have a plastic foraging behaviour, being able to alternate between feeding on fish in the surface waters or diving deep to feed on cephalopods depending on prey availability and the abilities of the individuals (Estupiñán-Montaño et al., 2019;Queiroz et al., 2012). The high δ 15 N values and trophic position found in blue sharks with medium body size and in females in the present study may reflect that they consumed prey occupying high trophic positions, or they are present in deep-water habitats, as deep-water species have enriched δ 15 N values in comparison to shallow organisms (Hannides et al., 2013;Pethybridge et al., 2018). Alternatively, these body size and sexual isotopic differences could also be associated with a different spatial distribution, as suggested by the joint contribution of these variables with geography and environment in models (Fig. 2). This last potential explanation might be consistent with the differences in δ 13 C values between male and female blue sharks. In marine environments, δ 13 C values present natural spatial gradients related to the latitude or the distance from the coast (inshore vs offshore) (Cherel and Hobson, 2007;Navarro et al., 2013). In addition to the biological variables, our results also showed that SST, Chl-a and pelagic ocean productivity PEL index (OPFish) predicted the δ 15 N values and trophic position of blue sharks. In marine ecosystems, Chl-a and gradients have a direct effect on marine productivity and, thus, on prey availability for blue sharks (Druon et al., 2022). In fact, these two environmental factors broadly explain the stable isotope spatial variability of other pelagic predators (e.g. Pethybridge et al., 2015;Dhurmeea et al., 2020;Martínez-Rincón and Acosta-Pachón 2022). However, in our case, each environmental factor showed contrasting curvilinear relationships with the stable isotope and trophic position values of blue sharks. In particular, at intermediate SST levels (18-20 • C) and intermediate Chl-a values, blue sharks showed the highest δ 15 N values. These patterns are likely related to the distribution of the main prey consumed by blue sharks, with larger sized or more δ 15 N enriched prey, such as cephalopods, found at depth in the intermediate SST values (Navarro et al., 2015).

Conclusions
This study provides new information on the factors that influence the trophodynamics of the blue shark, showing that biological, environmental, and geographic variables all play a significant role. The independent contributions of environmental variables and biological factors seemed to be more important than geographic location for δ 15 N and TP. δ 15 N and TP increased in a curvilinear fashion with body size, and TP was higher for females. Among the biological variables, sex was a better predictor than body size for δ 13 C models, suggesting that male and female blue sharks segregate their main feeding areas, confirming the habitat segregation result in Druon et al. (2022). Among the environmental variables, Chl-a, PEL and SST proved to be reliable predictors, particularly for δ 15 N and TP, most likely due to their relationship with productivity and prey availability. These results increase our understanding of the variables that influence the isotopic signature of the blue shark in the Iberian region and thus, on food-web dynamics along the study area. They point to differential feeding behaviour between size, sex, and spatial areas, and highlight the importance of large scale and comprehensive studies of trophic ecology for marine predators.  Model-averaged results for δ 13 C of the blue shark according to environmental, intrinsic variables and geographic variables. The table indicates the relative importance (weights, i.e. selection probability in the top-ranked models, ΔAIC≤2) of each variable, as well as full average parameter estimates for predictors receiving support (weight>0). Significant average parameter estimates (P < 0.05) are in bold. Model fit (R 2 ) and the number (N) of top-ranked models are shown. N = 173.

Variables
Partial models b Combined model b a Canary Island taken as reference level. b All models were conducted including the factor year as control.

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

Data availability
Data will be made available on request.