The Gulf of Naples (GoN) is a semi-enclosed embayment in the central-southern Tyrrhenian Sea along the Campania Region (southern Italy) ( Figure 1 ). The area has a typically temperate Mediterranean climate with marked seasonality. A complex geomorphology and wind-driven water movements result in a regular alternation between coastal and offshore waters within the gulf (Cianelli et al., 2015). From October to May, intense meteorological events in the area are triggered mainly by westerly and southern air flows, driven over by cyclonic circulations resulting from the oscillation of the polar front. In summer (June-September), however, the area is characterized by almost stable high pressure (Capozzi et al., 2023; Saviano et al., 2023; Saviano et al., 2022). During spring, the strong northwestward Tyrrhenian current, which is the prevailing driving force of the circulation patterns within the GoN during winter, becomes unstable due to the weakening of the cyclonic wind-stress curl. As a result, the local forcing becomes the main driver of the water masses within the basin, with the formation of an anticyclonic circulation within the GoN (de Ruggiero et al., 2020; Gifuni et al., 2022). Within this highly dynamic system, planktonic communities appear to have a stable composition. Phytoplankton form assemblages that occur consistently in time and are formed mainly by typical coastal species of temperate areas (Bosso et al., 2025; Bellardini et al., 2024; Zingone et al., 2023). Similarly, heterotrophic and mixotrophic protists as well as zooplankton communities include mainly coastal species, although the number of offshore species seems to have increased in recent years (Mazzocchi et al., 2023; Del Gaizo et al., 2021).

Pelagia noctiluca were collected within the GoN at the sea surface using a dip net during the 4 outbreaks recorded in March, June, July and November 2019 ( Figure 1 ). Each specimen was placed in a plastic bag with seawater to allow gut evacuation. Plankton samples, which encompassed all potential prey of P. noctiluca within the GoN, were collected weekly throughout 2019 within the monitoring at the Long-Term Ecological Research site MareChiara (LTER-MC) (Kokoszka et al., 2023) by vertical tows (from about 70 m depth to surface) using 20 and 200 µm plankton nets and transferred into 0.5L plastic jars. To isolate the smallest planktonic fractions, water was sampled by deploying a 5L Niskin bottle at the subsurface (1 m depth) and then transferred into plastic bins. All samples were kept fresh and transported to the laboratory on ice. Because the turnover rate of scyphomedusae is about 20 days (D’Ambra et al., 2014), only the samples collected 3 weeks before each outbreak of P. noctiluca were considered in the data analysis within the present study.

In the laboratory, each P. noctiluca was sized (bell diameter as the distance between two opposite rhopalia, cm), rinsed with sterilized filtered seawater to remove detritus and/or other planktonic organisms and frozen at – 30°C. Water and plankton samples were processed as described in detail in Merquiol et al. (2023), in order to isolate 4 distinct size fractions: 0-20 µm, which included pico- and nanoplankton; 20-200 µm (microplankton); 200-2,000 µm (mesozooplankton); > 2,000 µm (macrozooplankton). Seawater samples were pre-filtered through a 20 µm mesh nand the fraction <20 µm was concentrated onto pre-combusted Whatmann GF/F 0.7 µm filters to isolate pico- and nanoplankton (<20 µm), Net samples were concentrated onto a 20 and a 200 µm meshes and transferred into vials. Macrozooplankton were found only in March within the 200 µm sample and were sorted at the lowest taxon possible. Individuals were then pooled into two separate groups, Siphonophora and Salpida, to ensure sufficient organic matter to determine SIs. All samples were freeze-dried and the powder was homogenized using a mortar and a pestle.

Carbon and nitrogen SI ratios were determined for each scyphomedusa and for each planktonic size fraction and taxonomic group. According to the diverse organic content of each group, 4.0 ± 0.1 mg of scyphomedusae and macrozooplankton, and 2.0 ± 0.1mg of the other groups were packed into tin capsules and sent to the Stable Isotope Facility at the University of California in Davis (USA) for determination of δ¹³C and δ¹⁵N (for details about the instrument and standards see Merquiol et al. (2023). Duplicate aliquots of 26 samples were determined in order to ensure reproducibility, which was ± 0.04‰ for δ¹³C and ± 0.06‰ for δ¹⁵N.

The C:N of P. noctiluca and Siphonophora were lower than 3.5, while ratios were higher in all other groups ( Supplementary Table S1 ). Therefore, the δ¹³C of pico- and nanoplankton and Salpida was corrected according to Post et al. (2007) due to the lack of a specific equation, while the δ¹³C of mesozooplankton was normalized following Smyntek et al. (2007).

Bell diameters of scyphomedusae were compared among the 4 outbreaks using an analysis of variance (ANOVA). The relationship between the δ¹³C and δ¹⁵N and the bell diameter of scyphomedusae was defined using Pearson’s correlations.

Numerical simulations were performed using the Regional Ocean Modeling System (ROMS), a 3-D free-surface hydrostatic primitive-equation and finite-difference model widely used by the oceanographic scientific community for a wide range of applications (Haidvogel et al., 2008). A one-way nesting approach was adopted to obtain high resolution for the Campania coastal area. Firstly, simulations were performed on a grid covering the whole Tyrrhenian Sea with a 2 km resolution and 60 vertical levels, using the reanalysis data of the Mediterranean Forecasting System (https://medforecast.bo.ingv.it/ as initial and boundary conditions with a resolution of ~ 6–7 km and available with a daily frequency (Dobricic et al., 2005). Results of this first simulation were then used as initial and boundary conditions for a finer-grid computation, covering the Campania coast with a 500 m resolution and 50 vertical levels [see Kokoszka et al. (2022)]. Lateral conditions at the open boundaries were the same as those proposed by Marchesiello et al. (2001), with an adaptive algorithm where inward and outward fluxes were treated separately. The turbulence model used for the vertical mixing in the vertical direction was the K-Profile Parameterization (KPP) introduced by Large et al. (1994). In the horizontal field a constant value (equal to 10 m² s^-1) was attributed to the diffusion coefficient for velocity. Surface forcing was calculated from atmospheric values using the bulk parameterization of Fairall et al. (1996) based on the Coupled Ocean-Atmosphere Response Experiment (COARE) algorithm. The river inflows were modelled as point sources, with daily flow rate data obtained from a variety of sources, in particular from the Italian State agencies in charge of monitoring the environment (personal communications). The numerical configuration adopted in this work is an improved version of the GoN high resolution model (GNAM) validated in (Kokoszka et al., 2022). The major improvement regards atmospheric data, such as winds and air temperature, that were previously obtained from ERA-Interim reanalysis data (Dee et al., 2011) at horizontal resolution of about 16 km, while, in the present version, they were derived from dynamically downscaled ERA5 reanalysis data at a 2.2 km resolution. Lagrangian trajectories were simulated using the ROMS floats module, which allows the release and tracking of numerical particles, using a fourth-order Milne predictor and fourth-order Hamming corrector. Floats were all treated as geopotential, i.e. they were kept at the same depth during all simulations. The time step used in all simulations was 60 s, while particle locations were stored every 6 hours. While P. noctiluca may exhibit vertical migrations, our simulations assume passive surface transport due to limited data on vertical behavior in the GoN, thus the approximation of a purely horizontal displacement of the particles was selected for all the numerical simulations, in line with other studies using a Lagrangian approach for P. noctiluca in the Mediterranean Sea (Malačič et al., 2007; Bergamasco et al., 2022).

The following set-up was used for all Lagrangian simulations: 20 particles were released in random positions within each grid cell contained in the study area [lon = (13.6, 14.7), lat = (40.25, 41.05)], for a total of 892,960 particles. The release area was chosen based on initial tests (not shown here) to cover a sufficiently large area around the target zone. Particles were all treated as geopotential; thus, they were kept at the same depth during all the simulations. The time step used in all simulations was 60 s, while particle locations were stored every 6 hours. . In each of the 4 simulation scenarios, all particles were released at the same time, 20 days before the outbreak date and after 6 weeks from the start of the simulation, to allow for an initial spin-up of the velocity field from initial conditions. The choice of starting simulations 20 days before the outbreaks considers the surface dynamics of the study area (Cianelli et al., 2015) and the turnover rate of isotopes in scyphomedusae (D’Ambra et al., 2014). Post-processing of the simulations was made by individuating the trajectories that cross a target area in the 24 hours of the outbreak date (from 0:00 h to 24:00 h). The choice of the size of the target area, centered on the position of jellyfish outbreaks, was made by balancing the need to identify trajectories passing very close to the outbreak’ locations and the resolution of the numerical model, in order to gain a realistic accuracy. After a series of tests for all the 4 outbreaks presented here, a square target area of 1 km size was chosen in all simulations, corresponding to two grid cells at the 500 m resolution of the model.

The accuracy of the model trajectories was tested by a wind analysis which was carried out for 20 days before each outbreak to align with models and SI determinations, using in situ wind data (speed (m/s) and direction (°)).”Wind data were collected by two weather stations. The wind data used for the outbreaks in June, July and November were collected by the weather station operated by ISPRA that is located in the harbor of Naples (Molo del Carmine: lon:14.2699220345E; lat: 40.84149055N; anemometric sensor height: 10 m a.m.s.l.; data freely downloadable at http://www.mareografico.it/; date of access: 19 May 2024) ( Figure 1 ). Due to a lack of data in the ISPRA dataset, we used for the event in March the data retrieved by the weather station located at the University of Naples “Parthenope” (lon: 14.2452361E; lat: 40.8321083N) owned by Fondazione Osservatorio Meteorologico Milano Duomo (FOMD Network), a private and professional network of urban meteorological stations in Italy developed since 2010 (Curci et al., 2017). The weather stations are equipped with ultrasonic anemometer sensors located approximately 3 km away. Hourly averaged data were used in the analysis. Stick plots of the wind data for each outbreak were used to directly indicate the magnitude and direction of the wind during the 20 days before each outbreak ( Figure 2 ).

Overall, the mean bell diameter of scyphomedusae collected during the 4 outbreaks in 2019 was 7.6 ± 1.7 cm and varied across collections ( Supplementary Table S2 ). While the δ¹³C of scyphomedusae was -18.9 ± 0.7 ‰ and it was not correlated with the bell diameter (R² = -0.069, p = 0.056), the δ¹⁵N averaged 4.7 ± 1.6 ‰ and decreased with increasing bell diameter (R² = -0.856, p < 0.01). This finding is unusual because the correlation is positive as larger medusa can capture larger prey with higher δ¹⁵N (Schaub et al., 2023).

When looking at the relative position of P. noctiluca and plankton collected within the GoN, the δ¹³C and δ¹⁵N of scyphomedusae did not reflect the SI composition of any size fractions and/or taxonomic group ( Figure 3 ). During the outbreaks in March, June and July ( Figures 3A–C ), the SI composition of scyphomedusae fell within the range of SI values of their potential prey. However, if the diet-tissue fractionation factor (DTDF), which accounts for the trophic enrichment between the predator and the prey, is applied, P. noctiluca will remain outside the polygon of the potential prey. This suggests that prey with a very low δ¹⁵N are required to close the polygon around scyphomedusae. Conversely, in November ( Figure 3D ), the δ¹⁵N of scyphomedusae was already lower than the δ¹⁵N of all planktonic groups. By accounting for the DTDF, the δ¹⁵N would be similar to values determined in plankton collected in pelagic areas (Montoya et al., 2002).

The simulations were conducted for each outbreak observed in 2019, and the trajectories of particles were analyzed up to 20 days before the event.