Materials and methods

Sampling

Organisms were collected by epi- to bathypelagic trawling in canyons of the Bay of Biscay continental slope (North-East Atlantic) during EVHOE scientific cruises (“Evaluation Halieutique de l’Ouest de l’Europe ”; https://doi.org/10.18142/8) that took place between 2002 and 2021. Trawls were conducted at night between 25 m to 2000 m depth at 25 stations (Figure 1). The trawl net was 192 m long with a headline of 76 m and a foot rope of 70 m. The average vertical mean mouth opening was about 24 m and the horizontal opening of about 58 m. The mesh size gradually decreased from a very large 8 m (stretched mesh) at the mouth to 20 mm (stretched mesh) in the cod-end. To allow the capture of very small specimens, the trawl was also equipped with a 7.5 m long sock with a 12 mm mesh size. Each haul was made at a specifically chosen immersion depth. Once the trawl reached the selected depth it was towed horizontally (i.e., constant immersion depth) for 1 hour at 4 kn.

Datasets

Two different datasets were used to study ontogenetic changes in deep-sea pelagic fish species from the Bay of Biscay. The trawling (first) dataset included all data collected by trawls (i.e. number of individuals per species per sampling depth and total body length of each individual; n = 4165). To study the trophic aspects of ontogeny (see methodology below), muscle sampling was performed on 12 species of the trawling dataset to access the δ15N values of individuals (n = 682). This constituted the isotopic (second) dataset. The size measured for the individuals sampled for the isotopic dataset was the standard length. The size distribution of the individuals composing the species included in the isotopic dataset was representative of the size distribution observed in the trawling dataset (see appendix 1).

Nitrogen stable isotope analysis

A total of 682 muscle samples belonging to 12 of the most abundant species (seven migratory and five non-migratory (Loutrage et al., 2023) were collected (table II). For each individual, the standard length (cm) was measured on board and a small piece of muscle was collected and frozen at -20°C. To have sufficient material for stable nitrogen isotope analysis, the muscles of the smallest individuals were pooled. Within each of these pools, the individuals were of equivalent size and were sampled at the same depths. At the laboratory, muscle samples were freeze-dried (72h). To reduce the samples to a fine powder, samples containing a single individual were manually homogenised, while samples containing a pool of individuals were homogenised using a ball mill (MM400 Retsch®) with zirconium oxide-coated bowls and balls. A fraction of this powder (0.50 ± 0.05mg dry mass) was weighed in tin cups. Analyses were then performed with an isotope ratio mass spectrometer (Delta V Advantage with a Conflo IV interface, Thermo Scientific) coupled to an elemental analyser (Flash EA, 2000; Thermo Scientific). Results are presented in the usual δ notation relative to the deviation from an international standard (atmospheric nitrogen, forδ 15N values), in parts per thousand (‰). Based on repeated measurements of USGS-61 and USGS-62 used as laboratory internal standards, the experimental analytical precision was< 0.15‰.
The isotopic dataset included individuals sampled from different years (i.e. between 2007 and 2021), which could have affected our data. However, more than 90% of the muscle samples were collected between 2019 and 2021 and nearly 75% in 2021, which reduces the potential inter-annual effect. Furthermore, if a significant temporal change inδ 15N values had occurred, we expect that a common pattern would have been observed for the different species considered, which was not the case. Thus, we hypothesize that changes in size and depth have a much greater influence on theδ 15N values than does the temporal aspect. Details on the sampling by year for each species are presented in the appendix 2.

Relationships between size distribution and depth

The different depth layers were defined as follows: the epipelagic zone between 25 and 175m, the upper mesopelagic zone between 175 and 700m, the lower mesopelagic zone between 700 and 1000m, and the bathypelagic zone below 1000m. This division corresponds to the one used in the literature (Sutton, 2013) and is congruent with the depth structuration observed in the canyons of the Bay of Biscay (Loutrage et al., 2023). To study the changes in size distribution with depth, the trawling dataset was used. At both community and specific levels, a linear model was performed, the sampling depth corresponding to a continuous explanatory variable. Results were considered significant when the p-valuewas ≤0.05.

Relationships betweenδ 15N values and size

The relationship between δ 15N values and individual size was explored with the isotopic dataset at both community and species levels. In the case of pooled samples for nitrogen isotope analysis, the size data is the mean size of all pooled individuals. A linear model was used in each case (species and community levels). For non-significant relationships between fish size andδ 15N values (i.e. p-value ≥0.05) at the species level, coefficients of variation were calculated to explore the dispersion of values.

Variance partitioning

Variance partitioning was used to calculate the variance explained by the different variables included in a model (Borcard et al., 1992; Legendre and Legendre, 2012). This is done by developing a set of partial models (in a multivariate or univariate framework) created using a subset of predictor variables. Here, the objective was to test to what extent the individual size and the sampling depth influence theδ 15N values at the specific level. Due to the restricted depth range at which Aphanopus carbo and Stomias boa were captured (≤ 100m), the variance partitioning was not performed on these two species. The model results are composed of the proportion of δ 15N values influenced by size and depth separately, and a third fraction representing the shared fraction of variation explained when both variables are included in the model. An ANOVA-type permutation test was performed for each model to test the significance of the influence of each variable (depth and size) on δ15N values. Since the third fraction is deduced from the sum of variances, it cannot be tested statistically. The R packagevegan was used to perform the tests (Oksanen et al., 2022). All the graphics were performed with the ggplot2 R package and all statistical analyses were performed in the R environment version 4.3 (Wickham et al., 2016; R Core Team, 2023).