Spatio-temporal connectivity and dispersal seasonal patterns in the Adriatic Sea using a retention clock approach

1 Introduction

In marine systems, the persistence and recovery of a population are influenced by the dispersal and transport processes of individuals within it, which are controlled by the hydrodynamic circulation of a given study area (Williams and Hastings, 2013). When complemented by favorable environmental conditions that support the survival of the dispersed organisms, this dynamic process gives rise to the concept of population connectivity, a phenomenon that encompasses the exchange of individuals within and between geographical subregions through the transport of water bodies (Cowen and Sponaugle, 2009; García-Lafuente et al., 2021). For most marine species, exchanges occur through the pelagic dispersal of early life stages, when eggs and developing larvae can be treated, at least to some degree, as plankton (Largier, 2003). However, circulation can still play a key role in the success of organisms at advanced developmental stages, by conditioning their settlement and dispersal in particular locations, or by influencing the reproductive stock, e. g., altering the spawning location and time (Ciannelli et al., 2014). Estimating these transport pathways is a necessary step towards understanding the regional functioning of pelagic ecosystems, as well as for marine ecosystem management, the design of marine protected areas (Lester et al., 2009; D'Agostini et al., 2015) and the optimization of vulnerable fisheries resources (Fogarty and Botsford, 2007; Gaines et al., 2010; Nadal et al., 2022), among others.

A key component to estimate hydrodynamic connectivity, this understood as the potential of a hydrodynamic field to connect different sub-areas through the exchange of its individuals (García-Lafuente et al., 2021), is the spatio-temporal context in which sub-populations are connected (Treml et al., 2012). To this aim, hydrodynamic numerical models have become robust tools capable of representing circulation features and, when coupled with appropriate mathematical approaches, effectively analyzing dispersal pathways by simulating the lagrangian transport of virtual particles (Van Sebille et al., 2018). This is particularly relevant for coastal environments, which normally are highly dynamic regions dominated by tides, upwelling or small-scale turbulence processes, harsh to study via observational methodologies (Pineda, 1991). Ideally, a numerical model of a continental shelf should capture these small-scale coastal features on a high-resolution grid (at least, 0.5-1 km horizontally), with accurate near-shore tidal and meteorological forcing implementation and satisfactory calibration and validation.

An example of such complex application with a strong linkage to its coastal system is found in the Adriatic Sea (AS hereafter), a marginal water body characterized by a peculiar topography, having a very shallow northern area gradually deepening towards the southeast, and by a large number of freshwater sources (Russo and Artegiani, 1996). This paper explores the role of AS hydrodynamic circulation in governing population connectivity and entails the understanding of key physical processes and temporal-scales controlling the drifting of individuals, including spawning, dispersal and final settlement. The relevance of each species-specific behavior within the hydrodynamic circulation of the area makes every specific case study challenging. Thus, this study did not aim to replicate the behavior of any specific species, but rather focused on resolving the physical mechanisms controlling dispersal of passive elements, representing propagules of generic abundant species within the specified region. A high-resolution circulation model, already used in previous studies of this particular region, coupled to a lagrangian tracking module, is used as a numerical tool for the estimation of population connectivity and dispersal. A sophisticated concept introduced by Defne et al. (2016), consisting of a “retention clock” that highlights the evolution of source-sink connections through time, is employed to compare simultaneous spatio-temporal analysis of connectivity. The proposed methodological approach is applied on a tested year as a case of study.

The paper is organized as follows: Section 2 presents the hydrodynamic characteristics of the study area and the biodiversity of the region, including a list of ichthyofauna inhabiting the basin to which the results of our research may apply. The model framework used to represent the main circulation features of the AS and the lagrangian approach used to study the transport of virtual tracers, as well as the methods used to estimate connectivity from model results, are described in Section 3. The results of the lagrangian experiments on the tested year are discussed in Section 4, which includes subsections examining the influence of seasonality on virtual particle dispersal pathways and factors related to developmental aspects of abundant fish species in the AS, with a particular focus on pelagic larval duration. Finally, Section 5 contains the concluding remarks of this work.

2 Study area

The AS is a semi-enclosed basin in the northernmost Mediterranean (Figure 1). It is narrow and elongated in a NW–SE direction, with a length of ∼800 km and a maximum width of ∼200 km. It is bounded by an extremely complex rocky coast along the Balkans peninsula and by a sandy and regularly shaped shore along the Italian coastline (Ružić et al., 2019). It is divided into three main sub-basins, explained by the diverse bathymetric profiles (Russo and Artegiani, 1996). The Northern Adriatic (“NA” in Figure 1), extending from the Gulf of Venice and Trieste to the 100 m isobath, is characterized by an extremely shallow mean depth (∼30 m) with a very weak bathymetric gradient, and occupies the flooded seaward extension of the Po Delta. The Middle Adriatic (“MA” in Figure 1), still has relatively shallow waters, and constitutes the transition between the northern and southern sub-basins. MA is characterized by two depressions, the Pelagosa Sill and the Pomo (or Jabuka) Pit, with depths of ∼170 m and ∼270 m, respectively. The southern part (“SA” in Figure 1) extends to the Strait of Otranto, the channel connecting the Adriatic with the Ionian Sea, and is characterized by the depression of the Otranto Basin, with a maximum depth of ∼1270 m. The Otranto Channel, which is ∼75 km wide and ∼800 m deep, acts as an exchange area for water masses, and export of nutrients from the AS (Civitarese et al., 1998; Fanelli et al., 2023).

Figure 1 Map of the AS showing bathymetry from EMODnet, and sketching its general surface circulation: the Western and Eastern Adriatic Currents (WAC and EAC, respectively) and the North, Central and South Cyclonic Gyres. Orange dashed lines show the virtual boundaries of the three sub-basins, Northern, Middle and Southern Adriatic (NA, MA, SA, respectively), while semitransparent gray bands show the limits of the two Geographical Subareas (GSA17 in the north and GSA18 in the south; see Section 2.2). Maps were created using ArcGIS^® 10.8 software from ESRI.

2.1 Relevant hydrologic features

The surface flow dynamics of the AS are characterized, with a certain seasonal variability, by a general cyclonic circulation (Russo and Artegiani, 1996). Levantine intermediate waters and Ionian surface waters flow northward from the Otranto Channel along the Balkans coast within the Eastern Adriatic Current (EAC), and return southward with the Western Adriatic Current (WAC) along the Italian peninsula (Orlić et al., 1992). Coastal currents exhibit seasonal variability, with the WAC generally stronger in summer, and the EAC stronger in winter (Zore, 1956; Poulain, 2001). Prevalently during winter, part of the EAC recirculates, shaping the nearly-permanent Southern Adriatic Gyre in the southern sub-basin (López-Márquez et al., 2019), and less frequently, the Northern and Central Adriatic gyres (Martin et al., 2009) in the upper and middle sub-basins. These eddies partially favor the cross-basin (east-to-west) transport via their northern rims (López-Márquez et al., 2019), and also result in the pumping of highly productive waters that entail rich areas favorable for mesopelagic spawners (Specchiulli et al., 2016).

Tides in the AS are amplified northwesterly (from the Strait of Otranto to the northernmost Adriatic), reaching amplitudes up to 1 m in the Venice shoreline, an exceptionally large value in the Mediterranean basin (Medvedev et al., 2020). The influence of tides and local winds, namely Bora and Sirocco, induce an important variability in the main circulation pattern (Orlić et al., 1994). The dry, cold and strong northeasterly Bora, more frequent in winter storms, provokes a double gyre response in the North Adriatic circulation, consisting of a cyclonic loop drifting northward and an anticyclonic drifting southward (Kuzmić et al., 2006). Otherwise, the wetter and warmer southeasterly Sirocco, typical from late autumn to early spring, tends to accumulate water near the northernmost coasts (Pasarić et al., 2007; Molinaroli et al., 2023), usually strengthening the EAC (Book et al., 2007) and weakening or even reversing the WAC (Bignami et al., 2007). Summer winds (e.g. the northwesterly Mistral) are weaker and more stable, which favors a slower and steadier circulation during this period (Pasarić et al., 2009).

2.2 Ichthyo-biodiversity in the AS

The AS is a recognized biodiversity hotspot in the Mediterranean, driven by the nutrient transport and runoff by rivers and the consequent phytoplankton primary production, and marked by the strong seasonality of the Po River discharges (Cozzi and Giani, 2011). It is a highly productive fishing ground (Cavraro et al., 2022), grouped into two Geographical Subareas (GSA17 in the north, and GSA18 in the south, as shown in Figure 1) according to the General Fisheries Commission for the Mediterranean (GFCM, 2009). GSA17 has a wide variety of seabed habitats, from shallow and muddy bottoms in the east to steep and rocky bottoms in the west, whereas GSA18 is predominantly deep, steep, and rocky, supporting sensitive marine habitats that are under high fishing pressure (Grati et al., 2018). Despite the differences, both GSAs share important fish stocks (UNEP/MAP-RAC/SPA, 2015), which has led to the proposal of a large transboundary marine protected area (Bastari et al., 2016). Particularly, the AS supports regionally important fisheries of small pelagic stocks such as sardine (Sardina pilchardus) and anchovy (Engraulis encrasicolus) (Morello and Arneri, 2009), which are widely distributed throughout the whole AS. Over the last two decades, sardine specimens in the AS have tended to spawn from early autumn until late spring (Zorica et al., 2020). The peak of reproductive activity occurs between November and February, primarily depending on environmental parameters such as temperature (Zorica et al., 2019). In contrast to the AS, the spawning period for S. pilchardus appears to end earlier in other Mediterranean areas (Tsikliras et al., 2010; Basilone et al., 2021). The Adriatic anchovy spawns from the end of March (winter) to October (autumn), with its peak spawning season occurring in July (Zorica et al., 2020). According to Somarakis et al. (2004) and Basilone et al. (2006), the species’ reproductive peak is in the summer in the central Mediterranean. In the Bay of Biscay, anchovy spawn earlier and for a shorter period, peaking in May and June (Motos, 1996).

As for the bottom trawl fishery, the European hake (Merluccius merluccius) is the most important target species in terms of both landed weight and value for the fleets involved (Arneri and Morales, 2000). Unlike the other demersal species, it has two recognized spawning phases in the AS: the first occurring in winter in deeper waters, and then, a second occurring in summer after an adult migration to shallower waters (Ungaro et al., 1993; Vrgoč et al., 2004). Based on a recent study conducted on specimens sampled from the GSA17 (Candelma et al., 2021), the hake reproductive season peaks from April to July. However, spawning females can be found throughout the year, indicating a protracted spawning period (Zorica et al., 2019). It should be noted that the reproductive peak of M. merluccius varies in different geographical areas. Along the Tunisian coasts, the main peak is recorded from June to October, with minor peaks in January and April (Khoufi et al., 2014). In the Gulf of Lion, spawning is highlighted in autumn (Morales-Nin and Moranta, 2004). Along the Portuguese coast, spawning peaks are observed in March, May, and August (Costa, 2013). In North Atlantic waters, the spawning peak lasts from January to April (Alvarez et al., 2004; Murua and Motos, 2006).

Medium-sized pelagic fish species such as the Atlantic mackerel (Scomber scombrus), horse mackerel (Trachurus trachurus) and Mediterranean horse mackerel (Trachurus mediterraneus), despite their relatively lower commercial value compared to other pelagic fish species, play also a key trophic role as mesopredators in the basin (Da Ros et al., 2023), with their reproductive activity peaking in winter, late winter and summer, respectively (Jardas et al., 2004; Zardoya et al., 2004). High variability in range and peak spawning season is associated with latitude, with T. trachurus spawning season extending up to 8 months, with a peak in spring, in both the Atlantic and the Mediterranean (Abaunza et al., 2003). Other relevant fish stocks targeted by the bottom fleet are the red mullet (Mullus barbatus) and the common pandora (Pagellus erythrinus), with a spawning peak during summer (Carbonara et al., 2015; Muntoni, 2015; Zorica et al., 2020), as well as flatfish species such as the common sole (Solea solea), with a spawning peak in winter (Fanelli et al., 2022).

Some of the authors previously cited in their works indicate that some of the most important spawning and nursery grounds on the AS for these species are the Gulf of Trieste, Po Delta, Gargano Cape, Manfredonia Gulf, Kvarner Gulf, and the northern and the central Italian coastlines (see locations in Figure 1). Management measures intended to protect marine resources in the AS were additionally set out by the (GFCM, 2021) in the Pomo Pit (Figure 1), an area recognized as a critical habitat for demersal species and identified as an Ecologically or Biologically Significant Area (EBSA) under the 1992 Convention on Biological Diversity.

3 Experimental procedure

3.1 Hydrodynamic model

Modeling the hydrodynamic processes of the Adriatic is challenging due to the large number of spatial scales involved and the extremely intricate coastal morphology (McKiver et al., 2016). Such complex application was carried out with a high-resolution numerical model, whose code is based on the Shallow Water Hydrodynamic Finite Element Model (hereafter SHYFEM; Umgiesser et al., 2004). SHYFEM solves the shallow water equations of motion for complex morphology and wide bathymetric gradients on unstructured grids. The effectiveness of SHYFEM is well-documented in various applications in Europe (Bajo et al., 2014; Umgiesser et al., 2014 among others) and in the Adriatic basin (Bellafiore and Umgiesser, 2010; Ghezzo et al., 2015, Ghezzo et al., 2018; Umgiesser et al., 2022 and citations therein). The version used in this study (SHYFEM-Tiresias, Ferrarin et al., 2019) covers the entire AS, from 12.05° E to 19.92° E and from 39.95° N to 45.80° N, including the small lagoons of Grado-Marano, Lesina and Venice and the Po Delta (Figure 1).

The full numerical grid consists of 96,392 triangular elements with horizontal resolution varying from 5 km in the open-sea to a few dozens of meters along the coastlines (Figure 2A), and 72 vertical Z-layers of uneven discretization (Ferrarin et al., 2019). It includes the contribution of five boundary conditions, namely, the sea level, current velocity, temperature and salinity fluxes at the Strait of Otranto, and the freshwater discharges from 17 tributaries. Temperature and salinity fluxes were provided by the oceanographic fields of Mediterranean Forecast System (Tonani et al., 2008), and the freshwater discharges were obtained from monthly and annual climatological values at the river boundaries (for the location of tributaries, see Ferrarin et al., 2019 and references therein). Atmospheric forcing, i.e., pressure and wind, is given by the ECMWF ERA-5 atmospheric reanalysis (Hersbach et al., 2023).

Figure 2 (A) Numerical grid and bathymetry of SHYFEM-Tiresias (i.e., the SHYFEM application in the AS; Ferrarin et al., 2019). (B) Modeled surface velocity (color scale) and direction (arrows), in cm/s, in winter (January-March). (C) Same as panel (B) for summer (July-September).

The model outputs in this study are hourly surface values of zonal (u) and meridional (v) components of drift velocity, sea level, temperature, and salinity for the year 2018. The selected year is as a representative scenario for the general hydrographic conditions of the AS and is accessible in NetCDF through the Institute of Marine Sciences of the National Research Council of Italy (CNR-ISMAR) OpenDAP catalog (https://iws.ismar.cnr.it/thredds/catalog/emerge/catalog.html). It serves as a case study, examining spawning areas and intra-annual connectivity patterns, rather than assessing inter-annual variability, which is currently being investigated but not shown for simplicity in this first application. Supplementary Figure S1 displays the monthly averaged surface velocity data from SHYFEM and the ECMWF ERA-5 wind fields for the simulated months (January to December 2018). Winter and summer averaged surface circulation is presented in Figures 2B and C, respectively.

The synoptic fields in Figure 2 reveal similar circulation patterns, with the velocity field exhibiting greater variability in winter than in summer. The winter mean (Figure 2B) shows a more prevalent WAC in the northern sub-basin, and a stronger EAC along the north-flowing eastern rim of the southern gyre. The circulation is weaker in summer (Figure 2C), although the WAC is stronger along the south-flowing rim of the southern gyre in this season. Despite not being the most typical configuration reported in the AS, which is characterized by the presence of the EAC and the WAC and the three reported cyclonic gyres, the slight weakening of the currents during summer is not unusual. The prevailing northwesterly Mistral winds during summer intensify the east-to-west cross-basin transport, leading to stronger southeastward coastal currents and increased export to the northern Ionian Sea. During winter, the effect of both easterly and southeasterly winds disrupts the unidirectional shore transport, resulting in a more dispersed and irregular pattern within the basin.

3.2 Particle tracking module

Dispersal and connectivity were assessed by implementing a lagrangian particle tracking algorithm based on the open-source software OpenDrift (version 1.10.6, Dagestad et al., 2018). OpenDrift consists of several particle-based sub-models that can be used to predict the transport and fate of different types of elements. This study applies the main sub-model, OceanDrift, that uses neutral buoyant particles as tracers. Virtual particles represent propagules (i.e., eggs and larvae) with the typical biological traits of generic marine organisms (Table 1). The analysis focuses on basin-scale larval connectivity from an ecosystem perspective, rather than on specific target species.

Horizontal trajectories were simulated by integrating the zonal and meridional velocity field using a 4^th order Runge-Kutta advection scheme, where particle positions are bilinearly interpolated using the model output data in the hydrodynamic grid. Vertical velocities, mainly associated with wind-driven upwelling/downwelling in certain AS areas, are of small magnitude (less than mm/s to cm/s) compared to horizontal velocities (cm/s to dm/s). Vertical motions are influenced by other phenomena such as diel cycles, feeding patterns, egg buoyancy or settlement due to changes in fat content, diurnal temperature variations affecting buoyancy, etc. As a result, the vertical velocity is poorly determined, with an uncertainty in its value that could even change its sign, for which it has not been considered in the advection module. Another simplification adopted in our method is that the tracers are fully-passive, meaning that complex larval traits, such as swimming capability, migrations, and natality and mortality rates are not considered in the approach. Under these assumptions, transport of eggs and developing larvae is determined by the PLD (i.e., the time propagules spend drifting with currents), the duration of spawning, and the time-varying horizontal velocity. To prevent particle stranding near the release site, a strategy of relocating particles to the open sea was implemented upon encountering the continental coastline or islands when oceanographic conditions allowed. The coastline was represented using the Global Self-Consistent Hierarchical High-Resolution Geography in its full resolution version (GSHHG_F, version 2.3.7). A 30-minute time step with non-additional diffusion nor wind field were added to the trajectories, as they were already included in the hydrodynamic simulation run of SHYFEM-Tiresias.

To investigate the spatial connectivity patterns, the AS basin was divided into 40 sub-zones, 21 along the coastlines and 19 distributed by the open sea (Figure 3). These areas were partially selected based on the experience gained from studies on the circulation and ecological characteristics of the region of interest (Coll et al., 2010; Lipizer et al., 2014; UNEP/MAP-RAC/SPA, 2015; Bastari et al., 2016). Each area played, simultaneously, the role of source and destination of particles. Two hundred lagrangian particles as a representative number of propagules are seeded in each box randomly distributed in the surface layer, so each release allocates 8,000 new particles in the AS, and their position is computed every 30 min. To avoid any potential bias that could be introduced in trajectories if particles were only released at a single time of the day, particles are released four times per day for 14 days (56 releases, 448,000 particles in total). All the particle tracks were stored in netCDF format and the post-processing and visualization of the simulations were performed in Matlab (version 2023a).

Figure 3 Map of the study area with the subareas of source and destination of particles and bathymetric features. A summarizing table that illustrates the code numbers, associated names and surface areas in km² of each box is included.

In addition to the spatio-temporal analysis of the trajectories of eggs and larvae, several factors related to the biology and development of common fish stocks in the AS with influence on the dispersal pathways of particles were tested. Table 1 provides the biological parameters considered in the lagrangian experiments. The first one is the spawning time, understood as the time of particle release. Simulations were run for winter and summer, representing the most active spawning periods of the abundant species inhabiting the region. The second parameter being tested is the pelagic larval duration (PLD), which is the time between when a pelagic propagule leaves the spawning site and when it finally settles. In practice, the duration of egg and larval drift of marine organisms varies from weeks to months, depending on the ontogeny of fish larvae. The value of this parameter per species has been inferred from the literature cited in Table 1. Overall, given the uncertainties in determining the transition from passive to active behavior, we considered a wide PLD of 60 days as a good compromise between a stable value for connectivity in the AS and an affordable computational time. To test the sensitivity to PLD, connectivity was recalculated for time windows ranging from 15 to 90 days (Section 4.3).

Table 1 Literature review of biological traits for some emblematic species relevant for modeling in relation to the AS oceanographic conditions.

3.3 Construction and analysis of connectivity

Two different methods were followed to address hydrodynamic connectivity (Figure 4). First approach consisted in the so-called “connectivity matrix” (Figure 4A), which represents the probability of larval exchange among geographically separated sites (Cowen and Sponaugle, 2009). In this matrix, each cell is the number of particles (p) released from a certain source i (along the vertical axis) and collected in a certain destination j (along the horizontal axis), so C_i,j is the probability that an individual or group of individuals of population i will move to population j after a certain tracking time ($C_{(i,j)}= \frac{p(j)}{p(i)}$). Diagonal cells of this matrix (where i = j) stand for the self-recruitment, i.e., the number of individuals that remain in the same region from which they were originally released, including particles that eventually return to the sourcing area (Cowen et al., 2006). Connectivity ratios (C, with values ranging from 0 to 1), therefore, quantify the strength of connectivity between different sites, with higher values indicating more favorable connections. Source-sink dynamics of larval dispersal can be elucidated with this approach, revealing the direction and intensity of connections. It further enables to estimate the potential and effectiveness of a given area as a nursery ground for developing larvae. Despite this, connectivity matrices only offer a partial depiction of larval connectivity complexity, neglecting crucial information such as the identification of areas that are rapidly dispersing, both dispersing and recruiting, or only recruiting for individuals through time.

Figure 4 (A) An example of a connectivity matrix estimated in ten artificial subareas with fictitious ratios of connectivity varying from 0 (blank cells, null connectivity) to 1 (brownest colored cells, maximum connectivity). Rows correspond to the release zones of virtual larvae, whilst columns correspond to the destination areas of the released larvae. Diagonal stands for self-recruitment, which is the number of virtual larvae retained within the area from which they were released. (B) RCM applied to the same example, with an inset that zooms an example of a box connection with itself. In the inset, time progresses clockwise from 0 to a time scale T (e.g. 60 days) as concentration of particles decreases from 1 to 0, with a given temporal resolution of Δt (e.g. 12 slices of 5 days each).

Prompted by this lack of temporal information in the network, the second method used a more-recent approach developed by Defne et al. (2016). The method, called a “retention clock matrix” (RCM, hereafter), uses a clock in each cell to track the temporal changes in source-sink connections, evaluating the time-dependent connectivity (see Figure 4B). Specifically, each retention clock describes the release event as a circular clock that tracks the number of particles (p) over the entire time scale of interest (T) (see inset of Figure 4B). The time scale is discretized into slices with a temporal resolution of Δt_n ($C_{(i,j,t_n)}=\frac{p_{ }(j,t_n)}{p(i,t_0)}$). As in the conventional connectivity matrix, in all RCMs, the horizontal axis represents the boxes from which particles are released, and the vertical one is where the particles end up. The strength of connectivity between each possible combination of source-destination is depicted with a color scale, in which slices with the darkest color intensity indicate the larger fraction of particles moving from their origin to the destination. Naturally, domains with different retention characteristics have unique retention clocks, which reflect the rate at which particles are retained or lost within a domain through time. In slowly dissipating environments, the particle concentration, hence the clock, would gradually approach zero, indicating that individuals are being spawned at a slower rate than they are being dispersed. In rapidly dispersing domains, the clock would approach zero faster and a larger portion of it would remain at null values, suggesting that tracers are subjected to strong currents that result in direct connections to other domains. The opposite behavior is depicted by retentive domains, identified by clocks with connectivity values close to 1 at all-time scales, which indicate that spawned individuals are effectively retained due to physical or biological features that prevent their advection.

To account for the discrepancy in sizes and shapes of source and destination polygons (see Figure 3), the estimated connectivity in both cases was normalized by the areas (a) of each cell involved in the connection $C = C\times [\frac{a(i)}{a(j)}]$, so that the maximum value of the matrix is 1. In both methods, connectivity was only considered significant if the ratio C satisfied a minimum threshold of 0.01, meaning that at least 1% of the particles released from the source reached the destination polygon. Thus, the presence of either colored cells (in method 1) or clocks (in method 2) indicate a connection between the source and destination polygon, while a blank cell indicates that the connection is zero or less than 1%, and thus, negligible. The connectivity ratios C are hereinafter interpreted in the text as S# - D#, indicating particle transport from source (row number) to destination (column number), respectively. A numerical code made available to the scientific community by Defne et al. (2016) was adapted to calculate the RCMs in this study.

4 Results and discussion

4.1 Connectivity matrix

The mean connectivity matrix estimated at winter and summer for the year 2018 is illustrated in Figures 5A and B respectively, and the result of its difference (summer minus winter) is provided separately in Supplementary Figure S2.

Figure 5 Average connectivity matrices representing the mean exchange of particles between sub-zones in winter (A) and summer (B). The color scale represents the particle concentration at the selected time scale. A small map of the study area showing the sub-areas is displayed in the summer scenario to facilitate interpretation.

A characteristic feature of both winter and summer scenarios is the prevalence of high connectivity rates along the diagonal of the matrix. This indicates that these regions are characterized by a certain level of auto-retention of propagules. Self-retention probability differs greatly between the distinct boxes, ranging from 60% in the Kvarner Gulf (S13-D13) to 5% in Bosnia (S33-D33) for winter (Figure 5A), and from 55% in Manfredonia (S26-D26) to 4% in Brindisi (S34-D34) for summer (Figure 5B). This spatial difference is foreseeable, as the complex system of islands in the Kvarner Gulf, and the existence of a small cyclonic eddy in the Gulf of Manfredonia (see Figure 1) act as trap of propagules (Specchiulli et al., 2016; Sciascia et al., 2018), whilst Brindisi and Bosnia coasts, nearby the Strait of Otranto, are subjected to stronger advection by the EAC-WAC circulation system. Overall, self-retention in the northern and middle Adriatic sub-basins (areas 1-25) is more likely than in southern Adriatic sub-basin (areas 26-40), denoted by the darker color intensity of the former diagonal cells compared with the latter elements (~28% in areas 1-25, compared to ~13% in areas 26-40, in both winter and summer scenarios).

Considering also the connectivity values out of the diagonal, Gargano (D22) and Manfredonia (D26) during summer (Figure 5B), display a large fraction of particles arriving from S14-20 and S18-22, with mean connectivity values of 10% and 20%, respectively, and maximum reaching 50% (S18-D22) and 70% (S22-D26). This spatial pattern is to be expected, since the WAC acts as the main advection process along the Italian shelf, although it is also subject to seasonal fluctuations (Supplementary Figure S1). Contrastingly, Istria (D9), Kvarner Gulf (D13), and the Kornati (D17), Korcula (D25), and Elaphiti (D29) Islands, exhibit external particle reception rates less than 5%. Extreme case is provided by D17 (Figure 5B), which shows zero reception of individuals from external boxes. These areas harbor relatively closed populations, receiving particles from a limited number of grounds, but suggest their potential as self-sustaining areas. During winter (Figure 5A), isolated self-sustaining areas are less obvious, although the shielding capacity of the Kornati Islands (D17), northern Croatia (D21) and the Korcula Islands (D25) is noteworthy, as they receive nearly null abundance of particles from the northern part of the basin. Areas D2 (Venice), D3 (Trieste), D5 (OpenSea2) and D9 (Istria) only receive particles from a small part of the southern sub-basin, while source areas S12 (OpenSea6), S13 (Kvarner Gulf), S17 (Kornati Islands), S21 (North Croatia), S24 (OpenSea12) and S25 (Korcula Islands) do not send any particles to the southern sub-basin. A striking example is the Kvarner Gulf (D13), which paradoxically shows minimal external particle receptions in summer, but high external particle exchange in winter, while remaining self-retaining (elevate connection with itself). This suggests that this area may play a hatchery and nursery role for marine species, as suggested by the previously cited authors (Zorica et al., 2020), and suggest the existence of physical processes determining the temporal evolution of the particles in the basin.

Based on the averaged result, we can derive macro-regions as combinations of boxes of our network. For instance, more evident in summer than in winter, areas 1-9 and 26-40 may be identified as two isolated regions with a slight or null inter-connection. A third area can be identified in 10-25 (central Adriatic), where the particle reception is mainly limited to the upper part of the connectivity matrix (northern and central Adriatic), with limited or no exchange with the southern sub-basin during summer. During winter, the prevailing currents over the AS inject more energy in the basin, leading to a wider dispersion of particles and a less distinct delineation of macro-regions. The seasonality is also evident by the significantly higher number of valid elements, i.e., cells with connectivity probabilities higher than 1%, in winter (~900 out of 1600, Figure 5A) compared to summer (~600 out of 1600, Figure 5B). The likely reason is the atmospheric forcing, which is less variable and more stable in summer, resulting in diminished velocities and consequently shorter drifter paths in this period. The much more dispersed transport because of the more variable atmospheric forcing during winter may create otherwise less structured connectivity patterns within the basin.

4.2 Time-dependent connectivity

Temporally averaged connectivity matrices (Figure 5; Supplementary Figure S2) highlight the role of circulation structures on particle dispersion and provide predictable connectivity patterns that show seasonal variability. However, it also shows that the selection of a temporal snapshot is subjective and can limit the analysis due to loss of information. It is more adequate to preserve the temporal information of connectivity, while benefiting from the simplicity of a connectivity matrix (Defne et al., 2016). To this end, the RCM for all source-destination pairs in the AS basin has been estimated for the same temporal configurations (winter and summer). A maximum time scale of 60 days is used for all model scenarios. Each slice on all retention clocks represents 5 days. The color scale in each clock represents the particle concentration at the selected time slice, with darker colors indicating larger fraction of particles moving from source to destination.

4.2.1 Auto-connectivity

Auto-connectivity, the ability of a system to receive and retain particles from and within itself, is a pivotal determinant of population persistence. It represents the self-sustaining capacity of a population structure, a factor intricately influenced by the physical processes that mediate the propagules transport (Cowen et al., 2006). The quantification of self-recruitment probability, a direct metric of auto-connectivity, is made by evaluating the diagonal cells of the connectivity matrix (where i = j), with higher values signifying a more robust tendency for particles to remain within the area from which they were initially released. The cells corresponding to these diagonal matrices are extracted from the estimated RCMs and depicted in Figure 6, illustrating the time-dependent self-recruitment probabilities for both winter and summer seasons.

Figure 6 Clocks of auto-connectivity extracted from the diagonal of RCMs estimated at winter (A) and summer (B) and displayed spatially in the centered positions of the 40 defined subareas.

In the winter experimental configuration (Figure 6A), a low self-sustained clock pattern is evident along the Italian coast, indicating rapid southward transport within the western Adriatic shelf, with particles leaving the source area within a maximum of 10 days. This is consistent with the winter WAC pattern shown in Figure 2B. Regions located above the northern Adriatic gyre show retentive characteristics, as indicated by elevated clock values in both intensity and duration, particularly observed in Po (1), Venice (2), OpenSea5 (11), and OpenSea6 (12). The summer scenario (Figure 6B) presents a stark contrast, with the exception of the Croatian coast, which shows similar auto-retention characteristics. Specifically, the Kornati Islands (17) and northern Croatia (21) retain their winter characteristics for an extended period of up to 50 days, and the Gulf of Kvarner (13), along with the islands of Korcula (25) and Elaphiti (29), exhibit clock patterns almost identical to the winter configuration throughout the 60-day tracking period. This suggests that the coastal topography of the eastern Adriatic coast, rather than the seasonal regime, plays the primary role in influencing its auto-retention behavior. Po (1) and Venice (2), which contain complex lagoon systems, also maintain their behavior in both seasons due to their topographic conformation. Conversely, the coastal areas of Gargano (22) and Manfredonia (26) along the Italian shelf, exhibit higher auto-retention values in summer than in winter, in both intensity and duration, which persist for almost the entire 60-day period. Similarly, the onshore regions along the western Italian coast south of the Po estuary (6, 7, 10, 11) show high retention rates, in agreement with previous findings in this region (Revelante and Glimartin, 1992; Bray et al., 2017).

4.2.2 Winter connectivity

When the RCM estimated in winter (Figure 7) is analyzed, more granular temporal information becomes available with respect to the time-averaged analysis. For instance, among the previously identified macro-regions, the one corresponding to the northern Adriatic Sea (boxes 1-9), not only shows low ratios but also exhibits poor retention capacity with decreased persistence of connectivity over time. This is especially noticeable in the OpenSea3 (S7), which exhibits high average connectivity ratios (~40%, Figure 5A) but, when evaluated temporally, reveals consistently short and variable durations, with peak rates occurring at 40 days of tracking for particles originating from Trieste (S3-D7) and OpenSea2 (S5-D7).

Figure 7 RCM representing the time-dependent exchange of particles between sub-zones in winter, with the diagonal cells highlighted in blue color. A small map of the study area showing the 40 sub-areas and their corresponding codes is displayed to facilitate the interpretation of results.

The central macro-region (10-25) shows a more extensive distribution of connections, encompassing both the northern and southern sub-basins, with particle transfers to the southern Adriatic being significantly stronger and more persistent throughout the tracking period. Retention areas in this part of the basin occur in the Kvarner Gulf (D13), OpenSea10 (D20), OpenSea12 (D24) and the Korcula Islands (D25). It is particularly noteworthy in the case of the Kvarner Gulf (D13), which, despite the extremely complex rocky coastline, shows the highest connectivity receptions within the RCM (S16-D13, S15-D13, S12-D13 with a ratio of >60% sustained over longer periods). This pattern is likely due to enhanced atmospheric variability during winter, which induces rapid and direct transport of particles into these coastal areas. A strong and persistent rate of particle reception is further in the OpenSea6 (D12), which maintain a strong connectivity throughout the 60-day tracking period with the OpenSea8 (S16-D12), the Kornati Islands (S17-D12), the OpenSea10 (S20-D12) and North Croatia (S21-D12). These four areas shape a continuously well-connected water parcel compatible with the effect of the central Adriatic gyre, and likely identify one of the most productive areas for fish and important spawning and nursery grounds for commercially valuable fish species. Other important connections with the southern subbasin occur in the OpenSea9 (D19) and in the OpenSea10 (D20), which receives particles from Gargano (S22), OpenSea11 (S23), OpenSea12 (S24), and Korcula Islands (S25). These areas build a corridor of connection south of the Pelagosa Sill, with connectivity percentages of 45% throughout the 60-day clock period. The area D20 receives particles secondarily from S26-40, with special intensity in a period from 25 to 45 days from areas S29, S33, S37 which represent the pattern of the EAC. The connection of the central Adriatic with the northern is relatively weak and variable, as evidenced by the low (<30%) and delayed (>40 days) ratios in the upper portion of the matrix, attributed to the northward transport of currents towards the end of the tracking period. A relevant connection in this part of the basin occur in the OpenSea7 (D15) for particles arriving from the northernmost Adriatic (S1-D15, S2-D15, S4-D15, S5-D15, S7-D15, S10-D15 and S11-D15), indicating a particularly retentive area south of the Po River, corroborating previous findings in the region (Bray et al., 2017).

Among the three identified macro-regions, the southern one (26-40) emerges as the most retentive. It maintains a generally persistent connectivity, being particularly significant in specific regions of the central macro-region (S28-D24, S33-D24), with a high rate of connections (>70%) persisting 35 days. Additionally, it exhibits a high degree of inter-connectivity, with the most extensive and enduring connection occurring between areas 31, 32 and 35, 36. This pattern is likely attributable to the presence of the south gyre, which is particularly prominent during winter (Figure 2B; Supplementary Figure S1).

4.2.3 Summer connectivity

The rather scattered and variable connections seen in winter contrast with the much more unidirectional connectivity pathways in summer (Figure 8). This becomes apparent in the less heightened variability observed in both direction and duration of connections compared to the winter RCM. The northern macro-region (boxes 1-9), which showed a reduced ratio and persistence over time in winter, shows not only the highest but also the most persistent connectivity percentages in summer, denoted by the clocks with probabilities above 20% on all time scales (S2-D1, S3-D2, S4-D2, S5-D2, S4-D3, S2-D3, S2-D4, S3-D4, S3-D5). The darker color intensity in the matrix indicates a significant influx of individuals into destination 7 (OpenSea3), with a concentration of elements ranging from 20% to 60% between S1-D7 (Po-OpenSea3) and from 30% to 90% between S8-D7 (OpenSea4-OpenSea3) during the tracking time. Within the same box, particles originating from source boxes 3 (Trieste), 4 (OpenSea1), and 5 (OpenSea2), once arriving, maintain connectivity percentages above 70% for the rest of the tracking period. The particles reaching these areas exhibit rapid travel times (5 to 15 days), which is within the range of the PLD of the relevant species (Table 1) and indicates their potential for retention. Overall, the entire northern macro-region requires a longer time (approximately 30-60 days) to establish a connection with the central part of the basin (D10-25), with slight or null transport towards the southern macro-region (D26-40). This pattern is consistent with the slower flow conditions seen during summer (Figure 2C), but also provides important information regarding species strategy, as individuals arrive in that region at advanced developmental stages.

Figure 8 Same as Figure 7 for the summer experimental configuration.

Regarding temporarily sustained connectivity, the southern macro-region 26-40 follows, with the highest and longest connectivity observed between the boxes 35 and 36 (OpenSea17 and OpenSea18). These two sites also exhibit strong connectivity for the entire 60-day tracking period to boxes 31 and 32 (OpenSea15 and OpenSea16), forming a well-connected water parcel consistent with the isolation effect of the southern Adriatic gyre. Paradoxically, a notable connectivity is observed between Bosnia (S33) and Montenegro (D37), with 67% of particles flowing in the opposite direction of the EAC. This unexpected behavior is attributed to the reversal of currents very close to the shore due to instability created by the weakening of the coastal current during summer. In the southern macro-region, a mirrored relationship emerges in comparison to the northern counterpart along the south shore, characterized by a limited number of connections with the northern areas at longer times of arrival (>45 days). This pattern confirms that the northern and southern Adriatic, while internally connected, act as isolated water bodies within the region, which maximizes their potential as unique nursery environments.

The opposite situation is seen in the areas of the central macro-region 10-25, which reveal connections with a wider portion of the AS and demonstrates both retention and dispersal traits. Destination areas D10-12,14-16,18-19 consistently exhibit longer arrival times (>45 days) and poor retaining capacity for particles originating from the northern sub-basin compared to other boxes. In contrast, Gargano (D22) and Manfredonia (D26) are identified as top retaining areas, with connectivity consistently exceeding that of other destinations throughout all time slices. The majority of particles observed in the Manfredonia Gulf (D26), previously identified as a recruiting area in the averaged analysis (Figure 5B), are primarily and nearly continuously transported from the adjacent Gargano (S22) during the simulation period, with secondary and slightly delayed contributions from S23 and S24 (>15 days). This pattern arises from east-to-west shore transport, particularly noticeable during summer due to the displacement of the Adriatic central gyre (see Figure 2C). In the same area, particles require more than 30 days to arrive both from North Croatia (S21) and the Kornati Islands (S17), and more than 40 days from Ancona (S10) and the OpenSea areas of the northern Adriatic (S11, S12, S15, S16). This temporal pattern is to be expected: as the distance increases, the time of arrival of particles also does. But the point to note here is that the high rate of particles reaching Manfredonia stays and recirculates nearly continuously within the same box for the rest of the simulation period. This observation confirms the Gulf’s role as a particle receptor during the summer and subsequently as a nursery ground for larvae, a finding that could not be obtained in the time-averaged analysis (Figure 5B).

The spatial pattern observed between Gargano and Manfredonia (S22-D26) is consistent with patterns observed between San Marino and Ancona (S6-D10), and similarly with Pescara (S10-D14), and Vasto (S14-D18). This is likely attributable to the WAC, the primary advection mechanism transporting particles southward along the Italian shelf. Particularly during this season, the wind regime indicated a Mistral event blowing northeasterly (Supplementary Figure S1, July-August), which further supports the local displacement of particles onshore, despite the usual tendency of the coastal current to be detached from the shore during summer (Vilibić et al., 2012).

4.3 Effects of PLD on connectivity

The temporal component of connectivity between specific areas must be directly related to the species-specific PLD to ensure the practical application of the research. Determining the time that propagules spend drifting with currents is a key issue in shaping the dispersal potential and population connectivity of a given area. As mentioned in Section 3.2, the range of plausible durations is broad and uncertain (Table 1). While continuing to use the numerical model as a tool, it is worth exploring the dependence of connectivity on PLDs. To address this, we investigated the rate of particle reception in each area across a range of PLD windows from 15 days to 90 days, with intervals of 15 days (Figure 9).

Figure 9 Pie charts showing the proportions of connectivity at each area at the distinct temporal bands of PLD (see legend for interpretation) for the winter (A) and summer (B) experimental configurations.

In the winter configuration (Figure 9A), subareas show particle receptions evenly distributed across the PLD windows, with the most significant connections forming after at least 30 days of simulation. This is particularly noticeable along the Italian shelf (areas 6, 10, 14, 18, 22, 26, 30, 34), where the probability of receiving particles is less than 5% within the initial 15 days and remains below 20% within the first 30 days. This hinders the possibility of early life stages persisting in these areas, but suggests the eventual occurrence of individuals at more-advanced developmental stages, a phenomenon partially supported by works reporting the Gulf of Manfredonia’s role as a nursery ground (Sciascia et al., 2018). An exception to this trend is evident in the northernmost Adriatic, specifically in Venice (1) and Trieste (2), where the highest rate of particle reception occurs within the initial 15 days (35% and 53%, respectively), indicating a higher probability of occurrence of early life stages in these regions.

During the summer (Figure 9B), the arrival times are remarkably short, characterized by an average probability of reception of 30% and a maximum of 70% within the first PLD window (0-15 days). The most favorable scenario is observed in Istria (9), where the highest probability of receipt is concentrated in the first 15 days, with no further particle arrivals after 30 days of tracking. This observation indicates a prevalence of individuals in early life stages in this subarea, aligning with prior research highlighting Istria’s role as a spawning ground (Zorica et al., 2020). Other potential spawning grounds are discerned in the southeastern Adriatic, specifically in OpenSea17 (35), Montenegro (37), and OpenSea19 (39). Favorable scenarios are also evident along the northern edges of the Adriatic gyres, favoring cross-shore connections, with peak reception occurring between 15 and 30 days of individual tracking.

From a biological perspective, these findings offer valuable insights into spawning strategies. According to our results, species with shorter PLD, such as the red mullet or anchovy (<40 days), would likely thrive along the Croatian coast in the summer, benefiting from the shorter arrival times and the more favorable oceanographic conditions. Conversely, species with longer PLD, such as sardine or hake (>40 days), would find the northern Adriatic and Italian shelf during winter advantageous.

4.4 Implications for the species inhabiting the AS

Identifying the unique retentive characteristics of different areas in the AS is a crucial step in comprehending the distribution of marine species and designing effective conservation strategies in the basin. Our results indicate that particles exhibit a strong preferential direction along the coasts. The behavior in the Italian shelf was notably diverse, with some areas maintaining high self-sustaining ability during the summer and values shifting offshore in the winter. The Istrian coast (Figure 1) is identified as a crucial habitat for species spawning success, particularly in the summer, where arrival times fall within the ranges of PLD estimated for abundant species in the basin, such as anchovy, red mullet, and common pandora (Table 1). Areas distributed along the Italian shelf can be identified as essential nursery grounds, particularly when the region’s productivity is influenced by wind-induced upwelling and river discharges, creating favorable environmental conditions for larvae. During the summer, when the velocity is reduced, retention is promoted, coinciding with the spawning period of anchovy, Mediterranean horse mackerel, and hake (Table 1). The variability observed during winter is attributed to environmental fluctuations, specifically the prevalence of strong winds that significantly affect the displacement of drifters. This phenomenon provides advantageous opportunities for marine species with longer PLD and spawning periods concentrated in winter, such as sardines and common soles (refer to Table 1). Many particles arrive after undergoing recirculation, a process that requires longer times and results in extended arrival durations.

5 Conclusions

In this study, a hydrodynamic model integrated into the AS, coupled with a lagrangian module, was employed to systematically identify and evaluate sub-areas with distinct particle retention characteristics. While traditional connectivity matrices offer a valuable visualization tool for assessing the probability of individual exchange between zones, the challenge lies in extending this analysis across multiple time scales, often overlooked in existing approaches (Cowen et al., 2006; Crochelet et al., 2016; Gamoyo et al., 2019). To overcome this limitation, our study introduces an innovative application of the “retention clock” methodology, originally developed by Defne et al. (2016), and previously applied only in a few back-barrier estuaries (Barnegat Bay, in New Jersey, USA; Defne et al., 2016; Goodwin et al., 2019 and Chincoteague Bay, in Maryland/Virginia, USA; Beudin et al., 2017). Therefore, our research represents a pioneering effort and the first large-scale application of the temporal approach in the AS and, to our knowledge, in Europe.

The proposed methodology was tested in 2018 as a case study with representative hydrographic patterns of the AS, and in areas with known ecological characteristics. This allowed the verification of the distribution patterns of specific areas in the AS with the existing scientific literature (Morello and Arneri, 2009; Coll et al., 2010; UNEP/MAP-RAC/SPA, 2015; Fanelli et al., 2022). It also established a correlation between the findings and the life cycle traits of different organisms. Future research should aim to improve the coupling of connectivity and biological knowledge at larger spatiotemporal scales. The spatial discretization, as shown in Figure 3, followed a logical subdivision of the domain into a limited number of releasing sub-regions. A more detailed subdivision of space may be necessary in future applications, especially in areas where prior knowledge is lacking. This will require significant computational effort, which can be reduced by using computer clusters. Similarly, the analysis of connectivity over multiple years requires the development of an algorithm-based method to synthesize information and highlight the most relevant aspects. Such a study is currently in progress. All in all, the tool will improve the accuracy of the AS connectivity estimates which could, for example, inform the design of networks of marine protected areas (MPAs), fisheries assessment tools, and management policies.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://iws.ismar.cnr.it/thredds/catalog/emerge/catalog.html.

Author contributions

IN: Conceptualization, Data curation, Formal analysis, Visualization, Writing – original draft. MP: Data curation, Formal analysis, Writing – review & editing. FF: Conceptualization, Data curation, Methodology, Supervision, Writing – review & editing. JG-L: Conceptualization, Formal analysis, Supervision, Writing – review & editing. SS: Investigation, Methodology, Visualization, Writing – review & editing. MG: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Supervision, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. The publication of this work was financially and scientifically supported by the research activities of CNR ISMAR.

Acknowledgments

Authors are grateful to the Institute of Marine Sciences of the National Research Council of Italy for providing the model data necessary to conduct the numerical experiments (accessed at https://iws.ismar.cnr.it/thredds/catalog/emerge/catalog.html), with special thanks to Dr. Christian Ferrarin for his availability. IN acknowledges a predoctoral fellowship from the Spanish Ministry of Science and Innovation under the project BLUEMARO (PID2020-116136RB-100). The constructive suggestions and advices of the reviewers are gratefully acknowledged.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2024.1360077/full#supplementary-material

References