Seascape connectivity of European anchovy in the Central Mediterranean Sea revealed by weighted Lagrangian backtracking and bio-energetic modelling

Ichthyoplanktonic dataset

Spatial distribution of anchovy early-life stage was estimated from ichthyoplanktonic surveys in the GSA 16 from 2009 to 2012. During each survey, the same systematic sampling was carried out following a regular grid of stations (1/10° × 1/10° on the continental shelf and 1/5° × 1/5° further offshore). Mesozooplanktonic samples were collected in the study area using a bongo net (40-cm opening) towed obliquely from the surface to a 100-m depth, equipped with a 200-µm mesh size net. In each tow, the volume of filtered water was estimated using mechanical flowmeters (General Oceanics Inc., FL, USA). Over the period included in the analyses for anchovy eggs data, the average number of sampling stations was 154.

Samples were immediately fixed after collection and preserved in a 70% alcohol solution for further analysis in laboratory by stereomicroscopy. Anchovy larvae were identified in a land-based laboratory and otoliths were extracted for the age determination. Both sagittae were extracted using insect pins and fixed on a slide using a mounting medium for microscopy. Daily rings were counted at 1000× magnification (oil immersion) starting from the spawning mark³¹; each otolith was read blind (i.e., with no size or collection location information available) by a single reader⁴³. When all otoliths were completed, the process was repeated. If the difference in age between the 2 readings was < 5%, one was chosen randomly as the final read. If there was a discrepancy in age of ≥ 5% between the 2 readings, the otolith was read a third time. If after a third read a high discrepancy in the readings remained, the otolith was discarded.

No use of live vertebrates has been required for this study and no specific permissions were needed for the sampling activities since Engraulis encrasicolus is commercially harvested (neither endangered nor protected) and it was caught in areas where fishing is allowed. Sampling and analysis protocols are in accordance with the European Union regulations (EU Dir 2010/63).

Lagrangian modeling

Lagrangian approaches use numerical particle trajectories, as obtained from velocity field of marine currents. What matters most to the problem we are dealing with is the role of the spatial and temporal resolution of the velocity fields. In our case the Mediterranean Forecasting System (hereafter MFS)^44,45,46. The MFS model output consists in daily fields defined on a regular horizontal spatial grid, with 1/16 × 1/16 degree step, extended over 72 unevenly spaced vertical levels in depth. Considering that the Rossby first radius of deformation is O(10) km for the Mediterranean Sea, the MFS spatial resolution is not accurate enough to reconstruct mesoscale and sub-mesoscale structures in the study area.

MFS spatial resolution (~ 6.5 km × 6.5 km) does not allow an accurate reconstruction of mesoscale and sub-mesoscale structures (see Methods). Numerical experiments on predictability issues, related to Lagrangian simulations with the MFS model, demonstrate that the relative dispersion between tracer trajectories in the mesoscale and sub-mesoscale range is, indeed, significantly underestimated when compared to Mediterranean drifter data²⁴. To compensate this drawback, we hire a kinematic modelling technique (see Methods), which has been recently set up and applied to numerical studies of transport and dispersion in the Mediterranean Sea^17,23,24,47.

Therefore, to the MFS Eulerian fields for horizontal marine currents, we add two different versions of the Kinematic Lagrangian Model (hereafter KLM). The two-dimensional version, called KLM2, is a superposition of N_m independent modes, each associated to a given wavenumber k_n, which form the following velocity components u_klm2 (zonal) and v_klm2 (meridional):

where: x and y are zonal and meridional coordinates, respectively, of a Lagrangian particle; A_n is the velocity scale associated to the mode k_n; δ_n and ω_n are amplitude and pulsation, respectively, of the periodic oscillations of the velocity modes, and φ is phase difference (common to all modes) between x and y assigned to these oscillations²⁴. Here, for the KLM2, the range [k₁, k_Nm] corresponds to the inertial range of mesoscale turbulence of the Mediterranean Sea, as observed from experimental data. The velocity scales of the modes refer to the corresponding spatial scales via a Kolmogorov ‘− 5/3’ spectrum (two-dimensional inverse cascade): A_n ∝ (εk_n)^−1/3 where ε indicates the mesoscale 2D turbulent dissipation rate, of the order of 10^–9 m² s⁻³. This value is obtained from the horizontal dispersion that was estimated from experimental Lagrangian buoys²⁴, released in the Mediterranean Sea (https://marine.copernicus.eu/, product: INSITU_MED_NRT_OBSERVATIONS_013_035).

The oscillation parameters, δ_n and ω_n, control the ‘level of chaos’ of the Lagrangian trajectories: if they are both (or even only one of them) set to zero, the trajectories are regular and confined to periodic orbits; if they are assigned suitable values, Lagrangian Chaos can play its role at full regime⁴⁸.

For the mixed layer Lagrangian modelling, we use a three-dimensional version, called KLM3, which has three velocity components, i.e., u_klm3 (zonal), v_klm3 (meridional), and w_klm3 (vertical), defined as superposition of N′_m modes:

Here, x, y, and z are zonal, meridional, and vertical coordinates, respectively, of a Lagrangian particle; φ and φ’ are phase differences of the periodic oscillations between the three directions; [k₁,k_N’m] is the inertial range of the model; A_n and k_n are related by some scaling property, e.g., a Kolmogorov-type of spectrum, in this case, describing a 3D direct cascade^23,24. We observe that both KLM2 and KLM3 are conservative dynamical systems, i.e., the divergence of the velocity field is everywhere null. The number of modes in KLM3 does not have to be large since, in actual facts, what concerns most to our study is a good efficiency of vertical mixing, and this is accomplished just with a few velocity components.

By indicating with U_mfs, V_mfs, and W_mfs the three components of the ocean model velocity field, the full equations, which determine the evolution of the Lagrangian coordinates can be written as:

To mimic daily vertical migration (DVM) we also use the additional term²³.

where transfer rate γ is ≃ 0.1 h⁻¹ and z₀ is set to 3 (day time) and 100 m (night time), depending on daylight.

Otolith georeferenced data, which ranged between 1 and 17 days post hatching, are used to generate a large number (i.e., 100 individuals) of virtual larvae/eggs with a given age distribution. The backtracking during the larval period assumes that the growth rate is a normal variable with 0.49 mm/day (mean), and 0.05 mm/day standard deviation. The temporal extension of backward simulation is defined according to the planktonic stage duration, i.e. the sum of the egg and larval period³¹. Regarding the egg stage, we considered representative a period of 48 h, which is the incubation time of this species at the temperature of ~ 22 °C⁴⁹.

Biological constraint for Lagrangian trajectories

In nonlinear dynamics, the chaotic behavior of trajectories enhances dispersion over large scales. The same occurs for backward motion if the system is with good approximation conservative (i.e., Hamiltonian-like), as the ocean is. Hence, the identification of the origin of a tracer is largely affected by this problem and, consequently, concentration of traced particles may spread over unreliable large areas when integrated backward-in-time.

We therefore evaluate the amount of Chlorophyll encountered by each i trajectory during the simulation time by integrating the Chlorophyll concentration, i.e., ({rho }_{i}^{Chl}left(overrightarrow{x}right)), along trajectories (overrightarrow{{x}_{i}}left(tright)):

We then calculate the pdf of “zero-age positions” that satisfies the condition ({overline{{rho }^{Chl}}}_{i}>{rho }_{min}^{Chl}), where ({rho }_{min}^{Chl}=frac{1}{2}K) and (K) = 0.09 mg m⁻³, i.e., half of the saturation value for anchovy larvae feeding⁵⁰.

This idea is justified by the relationship between the pelagic growth rate and the Chlorophyll concentration: only those individuals that travel within Chlorophyll-rich marine waters have a significant survival probability. In this way, we have a criterion to rule out all the final positions of the simulated backward trajectories, which have poor or null Chlorophyll content.

Bioenergetic modeling

DEB represents a reliable and powerful tool to mechanistically describe the whole life cycle of an organism individual performance and to make predictions of life-history traits^{27, 40,50,51,52}. In particular, DEB is a complex mechanistic model that relies on several differential equations that are solved to obtain the final amount of energy allocated to vital functions, such as metabolic maintenance and to growth and reproduction⁵³. Here we adapted this model to the European anchovy and we adopted the standard version of the DEB model, which considers one reserve, one structure compartment and isomorphic growth⁵⁴.

The energy gathered through feeding processes is stored in a reserve pool, from which it is allocated according to the κ-rule: part of the energy (κ) sustains the somatic tissues and the growth of structures, while the rest (1 − κ) maintains the maturity level and maturation or gamete production in adults. Temperature controls the rates of all energetic flows and it follows the Arrhenius rules within the thermal-tolerance range^26,54,55. The Type II functional response⁵⁶ instead models the relationship between food density and ingestion rate⁵⁴. DEB theory therefore allows, through the explicit modeling of energy and mass fluxes through organisms, to derivate individual performance in terms of the most important life-history traits of a species such as, for instance, the total reproductive output (TRO) and maximum length.

In the present study, we followed a well-tested spatially explicit contextualized approach already successfully adopted in several companion studies^27,40,57. This approach consists in running DEB models in each spatial pixel of the study area using organismal body temperature (which is assumed to be similar to the Sea Surface Temperature, as extracted through satellite imagery for every single pixel) and environmental food availability. Food availability was expressed as density (wet mass mg m⁻³), which for anchovy primarily comprise zooplankton, and obtained as a spatially continuous dataset on the distribution of food throughout the study area⁴⁰. Nevertheless, to run DEB model in each pixel with local temperature and food density is a long computational process. Thus, due to the large number of pixels of the study area, we moved the standard DEB model from the original Kooijman’s Matlab code (https://www.bio.vu.nl/thb/)—once adapted it using anchovy DEB parameters⁴⁰—to an R code. In doing so, we were able to automatize and speed up the process since R coding improve the computational effort and allows running DEB models at larger spatial and temporal scales^27,40,58. The coding re-arrangement was performed by the Ecology Lab of the University of Palermo, which is one among the DEB node of the DEB world net (https://www.debtheory.org/wiki/index.php?title=DEBnet).

Our simulations were restricted to the continental shelf, based on depth (from 0 to 200 m below sea level). A vector polygon grid feature class of 346 square cells (having a size of 0.11° × 0.11° [~ 150 km²]) covering the study area was used. Food availability is an important forcing of the model and, for anchovy, primarily comprises zooplankton⁵². Since locally collected data for zooplankton were spatially and temporally fragmented due to sampling effort, we followed the recent approach proposed by Strömberg et al.⁵⁹ and applied by Mangano et al.^28,40 to obtain a continuous (space and time wise) dataset on the distribution of food throughout the study area and across time (i.e. weekly Net Primary Productivity was transformed into wet mass of zooplankton). Due to the short life span of the anchovy (~ 4 years), we extracted daily sea surface temperatures (SST; 1 km resolution) from JPL MUR SST data (2010) (https://podaac.jpl.nasa.gov/Multi-scale_Ultra-high_Resolution_MUR-SST) over a time range of 4 years (2009–2012) for each cell (see the recent proof of concept by Mangano et al.⁴⁰ for a detailed presentation of the model and its validation).

TRO values are validated by using in situ data, collected during ad-hoc oceanographic surveys⁴⁰. This approach assumes stationarity in biological parameters (i.e. DEB parameter values estimated for populations in one location/time are valid for populations elsewhere). We adopted the DEB parameters designed for the Mediterranean anchovy by Pethybridge et al.⁵⁰, the degree of uncertainty of our simulations was low (see Fig. 2 in Mangano et al.⁴⁰) and sufficiently robust to allow reliable predictions of anchovy life-history traits. We are aware that phenotypic plasticity and/or local adaptation have the potential to increase the degree of uncertainty of modeling outcomes and we suggest the use of DEB parameters values that, to the extent possible, realistically match those of local populations rather than global (species specific) parameters.

Source: Ecology - nature.com