Spatial changes in community composition and food web structure of mesozooplankton across the Adriatic basin (Mediterranean Sea)

Zooplankton are critical to the functioning of ocean food webs because of their utter abundance and vital ecosystem roles. Zooplankton communities are highly diverse and thus perform a variety of ecosystem functions; thus changes in the community or food web structure may provide evidence of ecosystem alteration. Assemblage structure and trophodynamics of mesozooplankton communities were examined across the Adriatic basin, the northernmost and most productive basin of the Mediterranean Sea. Samples were collected in June–July 2019 within the framework of the MEDIAS (MEDiterranean International Acoustic Survey) project, along coastal–offshore transects and from the surface to ca. 200 m depth, covering the whole western Adriatic side; consistently environmental variables were also recorded. Results showed a clear separation between samples from the northern-central Adriatic and the southern ones, with a further segregation, although less clear, of inshore vs. offshore stations, the latter being mostly dominated in the central and southern stations by gelatinous plankton. Such patterns were mainly driven, based on the outputs of the distance-based linear model, by fluorescence (as a proxy for primary production) for northern-central stations, i.e. closer to the Po River input, and by dissolved oxygen, together explaining 44 % of the total variance. Overall, at the basin level, the analysis of stable isotopes of nitrogen and carbon allowed for identifying a complex food web characterized by three trophic levels from filter feeders–herbivores to carnivores, passing through a general pattern of omnivory with varying preference towards herbivory or carnivory. Stable isotope signatures spatially varied between inshore vs. offshore communities and across sub-areas, with the northern Adriatic exhibiting greater δ¹⁵N and more variable δ¹³C than the other two sub-areas, likely attributable to the occurrence in the area of organic matter of both terrestrial and marine origin. Our results contribute to the knowledge of mesozooplankton community and trophic structure, at the basin scale across a coastal–offshore gradient, also providing a baseline for the future assessment of pelagic food webs within the European Council (EC) Marine Strategy Framework Directive.

1 Introduction

In an oligotrophic system, such as the Mediterranean Sea, coastal productivity largely depends on inputs from rivers, and areas of high productivity are mainly restricted to waters close to major freshwater inputs (D'Ortenzio and Ribera d'Alcalà, 2009; Ludwig et al., 2009). Here, the Adriatic basin represents an anomaly, with the northern Adriatic being one of the most productive Mediterranean areas. While the northern part is a shallow sub-basin, characterized by inputs of several rivers, with the Po representing the major buoyancy input with an annual mean discharge rate of 1500–1700 m³ s⁻¹ and accounting for about one-third of the total riverine freshwater input in the Adriatic (Raicich, 1996; Marini et al., 2008; Morello and Arneri, 2009), the southern part is characterized by highly saline and oligotrophic waters (Franco and Michelato, 1992; Boicourt et al., 1999). Thus, a trophic gradient, decreasing from northwest to southeast, is typically observed in the basin, in which the nutrient-rich waters coming from the rivers are mainly spread southward and eastward from the Italian coast (Aubry et al., 2006; Solidoro et al., 2009). Such differences may be reflected in the population dynamics of the marine biotic components (Revelante and Gilmartin, 1977; Simonini et al., 2004; Hermand et al., 2008), from zooplankton (Siokou-Frangou and Papathanassiou, 1991; Hwang et al., 2010) to fish (Wetz et al., 2011).

However, these dynamics both in terms of community composition and trophic relationships have never been investigated at the scale of the whole Adriatic basin. Zooplankton play a key role in marine ecosystems, forming the base of the marine food web because of the diversity of their functions. Zooplankton are a link between primary producers of organic matter and the higher-order consumers; they provide grazing control on phytoplankton blooms (Kiørboe, 1993) and help regulate fish stocks (Beaugrand et al., 2003), this last aspect being of crucial importance in the Adriatic basin. Because of these important zooplankton functions, a better understanding of their distribution and patterns of response to changes in the chemical and physical properties of marine waters is essential, especially under a global warming scenario, with zooplankton being a sensitive beacon of climate change (Richardson, 2008).

Moreover, trophic relationships in pelagic ecosystems are complex and complicated by the large degree of omnivory of most zooplanktonic species (Bode and Alvarez-Ossorio, 2004), which may feed on similar diets composed of a mixture of phytoplankton, detritus, and microplankton (e.g. Stoecker and Capuzzo, 1990; Irigoien et al., 1998; Batten et al., 2001). Several experimental studies allowed zooplankton (mostly copepods) to be categorized from pure carnivores to omnivores with a variety of mixtures of algae and animal prey up to strictly herbivore species (Irigoien et al., 1998; Batten et al., 2001; Halvorsen et al., 2001; see also Benedetti et al., 2016, and Hébert et al., 2016, for a review on functional traits of zooplankton). Such variety in the diet makes the quantification of flows between compartments or trophic levels difficult.

In the last decades, stable isotope analysis (SIA) has been widely used in food web studies; different studies dealt with high taxonomical groups of zooplankton (Burd et al., 2002; Blachowiak-Samolyk et al., 2008; Tamelander et al., 2008), while few investigations were focused on low taxonomical resolution (Koppelmann et al., 2003; Rumolo et al., 2017), essential to disentangle the food web structure of pelagic communities (Fanelli et al., 2011). Analysis of stable isotope composition provides indications of the origin and transformations of organic matter. Stable isotopes of carbon and nitrogen integrate short-term variations in diet and thus are less subject to temporal bias. The δ¹⁵N in tissues of consumers is typically greater by 2 ‰–3 ‰ relative to their prey and can be used as a proxy for the trophic level of organisms (Owens, 1987), while δ¹³C may act as a useful indicator of primary organic carbon sources of an animal's diet, as tissues tend to be rather weakly enriched in ¹³C at progressively higher trophic levels (1 ‰).

In this context, the main aim of this study is to analyse spatial variations in the assemblage structure and trophodynamics of mesozooplankton communities in the whole basin. Additionally, considering the complex hydrological condition of the basin, characterized by such contrasting oceanographic settings from north to south, here we explored and identified which environmental variables best explain the observed patterns.

2 Materials and methods

2.1 Study area

The Adriatic Sea is an elongated semi-enclosed basin, with its major axis in the northwest–southeast direction, located in the central Mediterranean, between the Italian Peninsula and the Balkans (Fig. 1). It is 800 km long and 150–200 km wide. It has a total volume of 35 000 km³, of which 5 % belongs to the northern basin, 15 % belongs to the middle basin, and 80 % belongs to the southern basin. The northern Adriatic is very shallow, with an average depth of 35 m with a very gradual topographic slope along its major axis, and it is characterized by strong river runoff, with the Po being the second main contributor (about 20 %) to the whole Mediterranean river runoff (Struglia et al., 2004).

Figure 1 Map of the study area with indication of WP2 net stations (black dots).

Due to this input, there is a positive water balance of 90–150 km³ that is exported to the Mediterranean. The turnover time for the whole basin is 3–4 years (Artegiani et al., 1997; Marini et al., 2017). The middle Adriatic is a transition zone between northern and southern sub-basins, with the two Jabuka (Pomo) depressions reaching 270 m depth. The southern sub-basin is characterized by a wide depression about 1200 m in depth. Water exchange with the Mediterranean takes place through the Strait of Otranto, which has an 800 m deep sill (Artegiani et al., 1997; Marini et al., 2017). The Adriatic is a temperate warm sea, with surface temperature ranging from 6 ^∘C in the northern part in winter to 29 ^∘C in summer. Even the temperatures of the deepest layers are, for the most part, above 10 ^∘C. The southern Adriatic is warmer than its central and northern parts during winter. In other seasons, the horizontal temperature distribution is more uniform (Artegiani et al., 1997; Marini et al., 2017).

Water circulation in the Adriatic is mainly driven by dominant winds (bora and scirocco) that cause a cyclonic circulation, with three closed circulation cells (one for each sub-basin). During the winter season, meteorological depressions pass over the Adriatic Sea; the first sector of the cyclone exposes the sea to warm Saharan air as the scirocco. As the cyclone passes, the winds reverse and expose the Adriatic Sea to a polar continental air mass, the so-called bora, coming from the north over central Europe and blowing the Adriatic Sea from the north and northeast. In summertime, corresponding to the time of our sampling, besides local breezes, the dominant wind, the mistral, comes from the northwest (Orlić et al., 1994). Climatological studies about the heat content of the water column (Artegiani et al., 1997) have resulted in the following definition of the Adriatic marine seasons: winter spans from January to April; spring occurs in May–June; summer goes from July to October; and autumn occurs in November–December.

Regarding temperature and salinity during the sampling period, in summer, the bathymetric effect (i.e. temperature gradients are at the same locations of topographic gradients) is evident: higher temperatures are observed in the northern part and along the western coast, and lower temperatures are in the southern part and along the eastern coast. Concerning spring, conditions are more like the summer ones (Russo and Artegiani, 1996). The distribution of salinity in the surface layer is strongly influenced (especially in the northern part and along the western coast) by river outflow, above all by the Po and other northern rivers. During summer, thermal stratification allows a wide horizontal distribution of these river waters inside the basin (vertically they are confined within the mixed layer, 10–30 m thick). The 38.0 psu (practical salinity unit) isohaline spreads southward and offshore, during spring and summer (Russo and Artegiani, 1996).

Three different water masses dominate the basin circulation: the Adriatic Surface Water (AdSW); the Levantine Intermediate Water (LIW); and the Adriatic Deep Water (AdDW), which branches out into the Northern (NAdDW), Middle (MAdDW), and Southern (SAdDW) Adriatic Deep Water. The hypersaline LIW is formed in the Levantine basin and experiences a salinity decrease on its way to the Adriatic. The AdDW is formed in the Adriatic basin, and the NAdDW is in the northern part; due to its high density, it fills up the Jabuka (Pomo) Pit and only occasionally spreads to the southern Adriatic. The MAdDW is formed in the Jabuka (Pomo) Pit area, when there is no intensive northwestward flow (i.e. during periods of low Mediterranean water inflow). The SAdDW originates in the South Adriatic Pit. During the period of the MEDIAS (MEDiterranean International Acoustic Survey) project (June), wind forcing was generally weak, and volume flux from the Po River was low, although the Po plume remained a significant feature in the northern and western Adriatic (Marini et al., 2008).

As mentioned above, the Adriatic is a very productive basin, compared to the rest of the Mediterranean. Despite being only the 5 % of the total Mediterranean surface area, the Adriatic Sea produces about 15 % of the total Mediterranean landings (and 53 %–54 % of Italian landings), with a fish production density of 1.5 t km⁻³, which is 3 times the Mediterranean density (Marini et al., 2017). This impressive feature is shaped by three main factors: river runoff, shallow depths, and oceanographic structure. River runoff is particularly strong in the northern basin and affects the circulation through buoyancy input and the ecosystem by introducing large fluxes of nutrients (Zavatarelli et al., 1998) which favour phytoplanktonic blooms and in turn cause a bottom-up effect of the whole trophic chain. Rivers can also provide suspended particulate organic matter and organic detritus that feed numerous particulate feeders and detritivores, such as bivalves (which is one of the main fisheries of the northern Adriatic Sea). The wide continental shelf favours a short trophic chain that likely improves the efficiency of energy transfer from lower trophic levels to higher ones. Moreover, the structure of the basin allows for water mixing during winter, especially in northern and middle Adriatic, transferring nutrients from sediments to the water column. From a fishery management point of view, the General Fisheries Commission for the Mediterranean (GFCM) has divided the basin in two geographical sub-areas (GSAs), the GSA 17, encompassing the northern and the middle sub-basin, and the GSA 18, including the southern part.

2.2 Zooplankton collection and analysis

Samples for this study were collected on board R/V G. Dallaporta during the acoustic survey MEDIAS 2019 GSA 17 and GSA 18 that took place in June–July 2019, in the Adriatic Sea (Leonori et al., 2021), within the framework of the MEDIAS (MEDiterranean International Acoustic Surveys) project (Leonori et al., 2021). MEDIAS coordinates the acoustic surveys performed in the Mediterranean and Black Sea to assess the biomass and spatial distribution of small pelagic fish (MEDIAS, 2019) (http://www.medias-project.eu, last access: 24 March 2022). Acoustic surveys are carried out using split-beam echo sounders working at specific frequencies which allow for discriminating target species from non-target organisms and zooplankton (see details in MEDIAS, 2019). Synoptical to acoustic sampling, fish and zooplankton are sampled, respectively, by means of a WP2 net and by a pelagic trawl.

Zooplankton samples were collected through a 200 µm mesh size WP2 net, with a circular mouth of 57 cm diameter and 2.6 m long, equipped with a MF 315 flowmeter to estimate the volume of filtered water. Vertical tows were performed with a towing speed of 1 m s⁻¹, starting from 3 m above the bottom, to the surface. Sampling stations were located along acoustic sampling transects (Fig. 1).

Zooplankton samples near the fishing hauls were subsampled and frozen at −20 ^∘C because of the requirements for SIA (see also Fanelli et al., 2009a, b, 2011, 2013; Rumolo et al., 2017, 2018). Concurrently with each vertical plankton haul, a CTD (conductivity–temperature–depth) cast was performed to acquire information on the oceanographic parameters of the chosen site. Environmental data recorded were pressure (dbar, decibar), temperature (^∘C), fluorescence (µg L⁻¹), turbidity (nephelometric turbidity unit, NTU), dissolved oxygen (expressed as mL L⁻¹ and saturation percentage), salinity, and density (kg m⁻³). For this study, the whole western Adriatic Sea has been divided in three different sub-basins or sub-areas, as described above (Artegiani et al., 1997): the northern Adriatic sub-area (NA), encompassing the stations from 4 to 38; the central Adriatic (CA), including stations 43–64; and the southern Adriatic (SA), comprising stations 5–57 (Fig. 1).

Selected zooplankton samples were analysed in the laboratory to characterize the planktonic community. First, frozen samples were defrosted and filtered with a 200 µm sieve, and the obtained mass was weighed (wet weight, WW, in g, precision 10⁻³). Then samples were quickly sorted, and larger animals were isolated first and placed in Petri dishes located on ice to preserve tissue integrity. Individuals were than identified to the lowest taxonomic level possible and stored for subsequent analysis. About 10 % of the sample was therefore weighed (WW in g, precision 10⁻⁵), and all organisms in the sub-sample were identified to the lowest taxonomic level possible (Cartes et al., 2010, 2013).

All identified taxa were then counted and weighed with an analytical weight scale to obtain abundance and biomass estimations.

2.3 Samples preparation for stable isotope analyses

The most abundant taxa in each sample were prepared for stable isotope analyses. Selected taxa were oven-dried for 24 h at 60 ^∘C. Dried samples were converted to a fine powder with a mortar and pestle. For each taxon, three replicates (when possible) were weighed (ca. 0.3–1.3 mg) and placed into tin capsules. Since it was not possible to obtain enough material of a single taxon for stable isotope analyses from stations 22 and 38, a bulk of the whole mesozooplankton community of the stations was prepared for the analyses. Acidification of samples prior to stable isotope analyses is usually regarded as a standard procedure, since inorganic carbon could lead to an increase of δ¹³C because it is isotopically heavier than most carbon of organic origin and could reflect the isotopic signature of environmental carbon (Schlacher and Connolly, 2014). However, for this study, no acidification was carried out, as this procedure generally reduces sample biomass, leading to too little matter available for isotope analyses. Moreover, some authors revealed negligible differences between acidified and unacidified samples (Rumolo et al., 2018). However, to have an indication of the possible bias, only one species was acidified, Euchaeta sp., which is a very abundant copepod in Adriatic communities. This taxon was also chosen because it has a more calcified exoskeleton, and it was abundant enough to undergo this process. Half of the sample was acidified with 1 M HCl, by adding it drop by drop to the sample until bubble cessation; then samples were oven-dried again at 60 ^∘C for 24 h. The other half, for the analysis of δ¹⁵N, was not acidified, as several studies demonstrated that the acidification procedure can alter the nitrogen isotopic signature (Kolasinski et al., 2008). Acidification of crustaceans was proved to be unnecessary, as the tested samples of Euchaeta sp. showed little and non-significant differences in δ¹³C value ($-\mathrm{21.39}\pm \mathrm{0.06}$ for untreated samples vs. $-\mathrm{21.02}\pm \mathrm{0.15}$ for acidified samples, paired t test $=-\mathrm{0.34}$, p=0.74). Then, six replicates of each sub-samples were prepared for isotope analyses. Samples were analysed through an elemental analyser (Thermo FlashEA 1112) for the determination of total carbon and nitrogen and then analysed for δ¹³C and δ¹⁵N in a continuous-flow isotope-ratio mass spectrometer (Thermo Delta Plus XP) at the Stable Isotope Ecology Laboratory of the University of Palermo (Italy). The stable isotope ratio was expressed, in relation to international standards (atmospheric N₂ and Pee Dee Belemnite for δ¹⁵N and δ¹³C, respectively), as

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\text{or}{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}:\left[\left(\left.R_{\text{sample}}/R_{\text{standard}}\right)-\mathrm{1}\right)\right]\times {\mathrm{10}}^{\mathrm{3}},\end{array}$$

where R is ${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ or ${}^{\mathrm{15}}\mathrm{N}{/}^{\mathrm{14}}\mathrm{N}$. Analytical precision based on standard deviations of internal standards (International Atomic Energy Agency IAEA-CH-6, IAEA-NO-3, and IAEA-N-2) ranged from 0.10 ‰ to 0.19 ‰ for δ¹³C and 0.02 ‰ to 0.08 ‰ for δ¹⁵N.

2.4 Community data analyses

Zooplankton abundance and biomass were standardized to a constant value. Zooplankton abundance was expressed as the number of individuals per square metre, while zooplankton biomass was expressed as milligrams of wet weight (WW) per square metre. This allows for minimizing the differences in the water column depths samples in the different stations; otherwise the use of data averaged in the water column (i.e. N or B per cubic metre) should have reduced the importance of offshore stations, as the numbers will be “diluted” in a large volume of water.

First, the Shannon–Wiener (H') diversity index of each station was calculated. Then, total biomass, total abundance, and the H' diversity index were tested by univariate PERMANOVA analyses (permutational multivariate analysis of variance, Anderson et al., 2008). Tests were run on Euclidean distance resemblance matrixes of log(x+1)-transformed data for abundance and biomass data and untransformed H^′ values (as data were normally distributed) and using a two-way design with the sub-area as a fixed factor with three levels (NA, CA, and SA, as described above) and inshore–offshore location as a fixed factor with two levels (inshore vs. offshore), crossed within each other, in order to assess the presence and significance of differences between stations. Inshore and offshore stations were selected according to Liquete et al. (2011). Univariate PERMANOVA tests were run under 9999 permutations, with the permutation of residuals under a reduced model; as the permutation method, significant p values were set at p<0.05.

To test for differences among sub-areas and inshore vs. offshore communities, a PERMANOVA test was performed on the Bray–Curtis resemblance matrix of log(x+1)-transformed abundance zooplankton data, using the same design described for univariate analyses. Data transformation is recommended for ecological data because they are often highly skewed and/or range over several orders of magnitude (as in this case) to weigh down the contributions of quantitatively dominant species to the similarities calculated between samples. This is particularly important for the most useful and commonly used resemblance measures like Bray–Curtis similarity, which do not incorporate any form of scaling of each species by its total or maximum across all samples. Here we used a severe transformation, i.e. log(x+1), which compresses large values, to take notice also of the less-abundant species (Anderson et al., 2008). A CAP analysis (canonical analysis of principal coordinates, Anderson and Willis, 2003) was then run to visualize the observed pattern on the factor found to be significant by PERMANOVA.

A SIMPER (similarity percentages) analysis was carried out according to the same sampling design to identify the most typifying taxon contributing to the average similarity/dissimilarity among sub-areas and inshore vs. offshore locations. This was conducted using Bray–Curtis similarity, with a cut-off for low contribution at 50 %.

To identify the environmental drivers of zooplanktonic communities and their structure across the sampling area, biotic data were correlated to environmental variables. Environmental data were tested for collinearity among variables by using a draftsman's plot, with fluorescence, dissolved O₂ concentration (DO, mL L⁻¹), the percentage of O₂ saturation, and turbidity data being log(x+1)-transformed to fit a linear distribution in the draftsman's plot. Finally, a DistLM (distance-based linear models, Anderson et al., 2008) was run with temperature, fluorescence, turbidity, oxygen, and salinity as environmental variables, using “step-wise” as the selection procedure and “AIC (Akaike information criterion)” as the selection criterion.

2.5 Stable isotopes data analysis

Since lipids can alter the values of δ¹³C (Post et al., 2007), samples with high lipid concentration can be defatted to avoid ¹³C depletion. However, lipid extraction can alter δ¹⁵N values; can complicate sample preparation; and can reduce sample availability, a crucial point when analysing small animals. For these reasons, the δ¹³C of samples rich in lipids was normalized according to the Post equation (Post et al., 2007):

$$\begin{array}{}\text{(2)}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\text{normalized}}={\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\text{untreated}}-\mathrm{3.32}+\mathrm{0.99}{\mathrm{C}/\mathrm{N}}_{\text{sample}},\end{array}$$

where the $\mathrm{C}/\mathrm{N}$ ratio was used as a proxy for lipid content because their values are strongly related in animals (Post et al., 2007). In particular, the normalization was applied to samples with a $\mathrm{C}/\mathrm{N}$ ratio >3, according to Post et al. (2007).

A hierarchical cluster analysis (Euclidean distance, average grouping methods) on the bivariate matrix of δ¹³C and δ¹⁵N mean values of each taxon was performed to elucidate the planktonic food web structure. Obtained clusters were also compared with literature data on the trophic guild of analysed taxa. Four main trophic groups were established a priori on the basis of literature data, where available, and by adapting the classification suggested in Hébert et al. (2016), Benedetti et al. (2016), Fanelli et al. (2011), and recent findings based on both SIA and fatty acids for some of the species here analysed (Protopapa et al., 2019). Thus, trophic groups used for the following analyses were filter feeders–herbivores (FF-HERB) considered primary consumers; omnivores with a clear tendency toward herbivory (OMN-HERB), encompassing mostly herbivore species but also those that can feed on small particles and ciliates; small carnivores (OMN-CARN), similarly to OMN-HERB but with a greater preference for small zooplankton; and carnivores (CARN), including also the parasite hyperiid Lycaea pulex. Differences among groups were tested by means of a one-way PERMANOVA test with “trophic group” (with four levels, corresponding to FF-HERB, OMN-HERB, OMN-CARN and CARN) as a fixed factor.

The trophic level of the different species was estimated according to Post (2002) as

$$\begin{array}{}\text{(3)}& \left(\left({\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_i-{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{PC}}\right)/\text{TEF}\right)+\mathit{\lambda },\end{array}$$

where δ¹⁵N_i is the δ¹⁵N value of the taxon considered; δ¹⁵N_PC is the δ¹⁵N values of a primary consumer, i.e. an herbivore or a filter feeder, used as a baseline of the food web; TEF is the trophic enrichment factor, which is considered to vary between 2.54 (Vanderklift and Ponsard, 2003) and 3.4 (e.g. Vander Zanden and Rasmussen, 2001; Post, 2002) and here is assumed to be 2.54 for low-trophic-level species, according to Fanelli et al. (2009a, b, 2011); and λ is the trophic position of the baseline, which is 2 in our case. Here, we used three different values as baselines for the food web of the three sub-areas, specifically the average values of FF-HERB taxa.

Then, differences in the isotopic composition of the overall communities by sub-area and inshore vs. offshore communities were tested by two-way PERMANOVA analysis with the same design used for assemblage analysis. The same procedure was also used to perform univariate two-way PERMANOVA and one-way PERMANOVA analyses with a pairwise test for the δ¹³C and δ¹⁵N values, separately.

Finally, maximum-likelihood standard ellipses were created for the δ¹³C and δ¹⁵N values following Jackson et al. (2011) to assess the community niche width in the different sub-areas. In addition to standard ellipse area (SEA; containing ca. 40 % of the data and representing the core isotopic niche) and standard ellipse areas corrected for small sample size (SEAc), traditional convex hulls and four Layman metrics were also estimated (Layman et al., 2007). Specifically, we calculated TA, which is the area of convex hull containing, in the case of SIBER (Stable Isotope Bayesian Ellipses in R, Jackson et al., 2011), the means of the populations that comprise the community; d15N_range, which is the distance in units between the min and max y-axis population means; d13C_range, i.e. the distance in units between the min and max x-axis population means; and CD, which is the mean distance to the centroid from the means. Ellipse sizes were compared between groups (i.e. sub-areas) using Bayesian inference techniques.

All analyses were run using the software PERMANOVA+ for PRIMER (Anderson et al., 2008; Clarke and Gorley, 2006) and within the jags and SIBER packages in R 4.1.0 (https://www.r-project.org/, last access: 10 January 2022).

3 Results

3.1 Zooplankton community and spatial changes

A total of 52 016 specimens belonging to 113 taxa were collected through the WP2 sampling (Table S1 in the Supplement). Zooplanktonic communities in the whole area were dominated by small copepods of the genus Acartia (mostly A. clausi), Oncaea, Oithona (mainly O. similis), and copepodites. Abundant large copepods were Calanoida belonging to the genera Euchaeta, Calanus, Centropages, and Temora. Since samples were frozen on board after collection for subsequent isotopic analyses, a quite considerable number of specimens (particularly amphipods and mysids and those taxa/specimens characterized by soft carapace) were damaged and therefore hard to identify at the species level. Generally, they were identified to order level or indicated as “damaged unid.” (unidentified) in Table S1. Other common crustaceans were hyperiids, such as Lestrigonus schizogeneios and Phronima atlantica, decapod larvae (mainly zoeae and megalopae), mysids, and euphausiids. Among non-crustaceans, molluscs were quite common, both as larvae of benthic organisms and adult pteropods. Chaetognatha were also locally abundant. Gelatinous zooplankton was represented mainly by thaliaceans and calycophorans, while ichthyoplankton was not very abundant, with few fish eggs and larvae found.

Zooplankton abundance and biomass varied according to geographic sub-area, decreasing from the northern to the southern Adriatic (Fig. 2a–b) and to the distance from the coasts. However, differences at the sub-area scale were significant only for abundance, while inshore–offshore differences were significant only for biomass (Table S2 in the Supplement).

Figure 2 Total abundance (number of individuals per square metre, a), total biomass (mg WW m⁻², b), and diversity (H^′, c) of mesozooplankton at each group of stations by sub-area and distance from the coast (inshore vs. offshore stations). Colours define the different sub-areas. Boxes are interquartile ranges; black lines that divide the box into two parts represent the medians; and the upper and lower whiskers represent scores outside the middle 50 %.

Diversity (in terms of H^′) increased southward (Fig. 2c), although differences were not significant for any of the investigated factors. H^′ values were on average 3.25±0.31, with only the exception of station 22, located in the NA inshore, in front of the Po delta, showing the lowest H^′ value (0.64).

PERMANOVA analysis revealed that differences in zooplanktonic communities, based on geographic sub-areas and the inshore–offshore factor were significant, while any significant differences occurred for the interaction factor (Table S3a–b in the Supplement).

The CAP plot showed a clear separation among samples from each sub-area, with the first axis separating samples from NA from those belonging to CA and SA sub-areas (Fig. 3).

Figure 3 CAP plot of the mesozooplanktonic communities of the Adriatic basin by sub-area and inshore vs. offshore location, based on abundance data. Colours indicate the sub-basins, as described in the text.

SIMPER analysis showed that Calanus-like copepods, Euchaeta sp., Euterpina acutifrons, and Evadne spinifera mainly contributed to dissimilarity between NA and CA (Table S4a in the Supplement). Bivalve and gastropod larvae, together with Acartia sp., were mainly responsible for the dissimilarity between the sub-areas CA and SA. Within NA samples, the dissimilarity between inshore and offshore zooplanktonic communities was mostly driven by Calycophorae, Calanus helgolandicus, and Chaetognatha, being more abundant at offshore stations. The cladoceran Penilia avirostris, thaliaceans, ostracods, and Calanus helgolandicus were responsible for the dissimilarity between inshore and offshore stations within CA, with P. avirostris occurring only at inshore stations and thaliaceans, ostracods, and Calanus helgolandicus being dominant at offshore ones (Table S4b). Large calanoid copepods dominated the inshore communities within the SA sub-area, while the euphausiid Meganyctiphanes norvegica was more abundant at offshore stations (Table S4b).

3.2 Environmental variables and correlation with zooplankton data

During the sampling period, temperature values were on average 18.5 ^∘C (±0.88 SD), with the lowest and the greatest values observed at inshore and offshore stations, respectively, in the central and southern sub-basins (Table S5 in the Supplement). Salinity values were on average 36 in the northern basin with the lowest value of 34.7 recorded at station 22_17 in front of the Po; salinity increased southward reaching a mean value of 38.7 in the southern basin (Table S5). Fluorescence values decreased southward from 2.45 to 0.77 µg L⁻¹, with the highest (4.9 µg L⁻¹) and the lowest (0.59 µg L⁻¹) values recorded at station 22 (in front of the Po River delta) and at station 44 (in the Otranto channel), respectively (Table S5). On the other hand, dissolved oxygen (DO) decreased southward from a mean value of 5.32 mL L⁻¹ recorded in NA stations to 4.36 mL L⁻¹ observed in SA CTD casts (Table S5). Significant variations were observed for all tested variables for sub-area and inshore vs. offshore factors and for the interaction term only for temperature and dissolved oxygen (Table S6a in the Supplement). Pairwise comparisons evidenced significant differences in salinity, fluorescence, and dissolved oxygen values between NA and CA (Table S6b). Significant differences between inshore and offshore stations occurred in the southern sub-basin for temperature, salinity, and DO; in the central sub-basin for temperature and salinity; and in the northern sub-basin only for salinity (Table S6b).

According to the results of the draftsman's plot, DO concentration (mL L⁻¹) and the percentage of oxygen saturation covaried (ρ>0.7); as well as density and pressure, therefore, only temperature, fluorescence, turbidity, DO, and salinity were used for DistLM analysis. DistLM results showed that 44 % of the variance was explained by fluorescence (33 %) and by dissolved oxygen (11 %) (Table 1, sequential test) and provided the best model solution in terms of both AIC and R² values.

Table 1 Results of the marginal and the sequential test for the DistLM model, with an indication of the best model. AIC: Akaike information criterion, Cumul.: cumulative, No. vars: number of variables, Prop.: proportion of variation explained by each variable, Res. df: residual degrees of freedom, RSS: residual sum of squares, SS(trace): portion of sum of squares relative to the analysed predictor variable.

2: fluorescence, 4: oxygen.

3.3 Stable isotope composition of zooplankton

Stable isotope analyses provided δ¹³C and δ¹⁵N values of 25 different taxa (Table 2).

Table 2 Mean values of zooplankton samples analysed for each sub-area (NA: northern Adriatic, CA: central Adriatic, SA: southern Adriatic), trophic group (TG), and trophic level (TL). Unid.: unidentified. FF-HERB: filter feeders–herbivores, OMN-HERB: omnivores with preference towards herbivory, OMN-CARN: omnivores with preference towards carnivory, CARN: carnivores. “Base” indicates the species used for the estimation of the average δ¹⁵N values of the baseline for TL calculation (see text for further details).

Cluster analysis allowed for grouping animals according to their δ¹³C and δ¹⁵N values and according to the trophic groups previously established, based on literature data (Fig. 4a). Still, the nMDS (non-metric multidimensional scaling) analysis evidenced a gradient from strictly herbivore species towards carnivore taxa (Fig. 4b). The one-way PERMANOVA test run on the factor of trophic group (TG) was significant (pseudo $F_{\mathrm{3},\mathrm{25}}=$ 13.12, p=0.0001), with significant differences between each level of pairwise comparisons across the herbivory–carnivory trophic gradient (FF-HERB vs. OMN-HERB: t=20.52, p=0.02; OMN-HERB vs. OMN-CARN: t=22.69, p=0.005; OMN-CARN vs. CARN: t=22.11, p=0.007).

Figure 4 Cluster (a) and nMDS (b) analyses on the bivariate matrix of δ¹³C and δ¹⁵N values of dominant zooplankton taxa averaged for the whole sampling area. Colours indicate trophic groups as follows. FF-HERB: filter feeders–herbivores (dark green), OMN-HERB: omnivores with preference towards herbivory (light green), OMN-CARN: omnivores with preference towards carnivory (orange), CARN: carnivores (red).

The estimates of trophic levels (TLs), considering the average δ¹⁵N value of FF-HERB for each sub-area as a baseline (from Table 2) and specifically δ¹⁵N = 3.6 for NA, of 3.4 for CA and 3.5 for SA allowed for assigning zooplanktonic taxa to three TLs from strictly herbivores located at TL 2 to carnivores at TL 4 (Table 3).

Table 3 Estimates of the convex hull (TA); standard ellipse area (SEA and SEAc, as corrected for low sample size); and Layman metrics d15N_range, d13C_range, and mean distance to centroid (CD), calculated for zooplanktonic communities from the three sub-areas. NA: northern Adriatic, CA: central Adriatic, SA: southern Adriatic. Sample size is also provided.

Overall, the δ¹⁵N of the mesozooplanktonic community was greater at NA, especially for inshore communities (Fig. 5). Conversely, the median δ¹³C value was similar among the different sub-areas; however the larger variability was observed in the inshore communities of the NA sub-area (Fig. 5).

Figure 5 Box plot of mean δ¹⁵N and δ¹³C values of zooplanktonic taxa for each sub-area at inshore vs. offshore locations. Boxes are interquartile ranges; black lines that divide the box into two parts represent the medians; and the upper and lower whiskers represent scores outside the middle 50 %.

Two-way PERMANOVA analysis on the multivariate matrix of δ¹³C and δ¹⁵N and one-way PERMANOVA analysis on δ¹⁵N values showed a significant separation according to sub-area and inshore vs. offshore factors but not for the interaction (Table S7a in the Supplement). However, the pairwise comparisons on sub-area factors did not show significant differences between contiguous sub-areas (but only between NA and SA), while the pairwise test run on the interaction factor for pairs of level of the factor “inshore vs. offshore” provided evidence for significant variations in the overall isotopic composition (δ¹³C–δ¹⁵N) and in the δ¹⁵N values between inshore and offshore communities only for taxa from CA (Table S7b). One-way PERMANOVA tests run on δ¹³C values showed significant variation for the factor inshore vs. offshore and for the interaction term (Table S7a); δ¹³C values significantly varied between CA and SA taxa and between inshore and offshore communities at CA (Table S7b).

Finally, the SIBER method for calculating ellipse-based metrics of niche width provided evidence of larger niche width for the zooplanktonic community from NA than CA and SA (Table 3 and Fig. 6). Estimated overlap by Bayesian inference evidenced low overlap among standard ellipse areas from contiguous sub-areas, being 2.71 between NA and CA and 2.07 between CA and SA. The greater d15N_range was observed for NA and SA communities, while the higher d13C_range occurred in NA communities, where also CD value was the greatest (Table 3).

Figure 6 (a) Standard ellipse areas for the three zooplanktonic communities analysed. The black circle and symbols indicate the NA community; the red ones indicate the CA community; and the green ones indicate the SA community. (b) Credible intervals for the estimated SEAc of the three communities. NA: northern Adriatic, CA: central Adriatic, SA: southern Adriatic.

4 Discussion

These are the first results on the mesozooplankton food web structure conducted at the basin scale for the Adriatic Sea. Considering that the Adriatic Sea is one of the largest areas of occurrence of demersal and small pelagic shared stocks in the Mediterranean (FAO, 2020), this study may represent an important piece to reconstruct the whole pelagic food web and spatial changes across the basin. Still, considering the increasing fishing pressure in the basin together with evidence of primary production (climate-change-related) decrease after the 1980s (Solidoro et al., 2009; Mozetič et al., 2010), this study may represent a valid baseline for future comparison on the synergic and cumulative effect of climate change and overfishing in one of the most impacted regions within the Mediterranean Sea (Coll et al., 2012; Micheli et al., 2013).

4.1 Spatial variations in zooplankton biomass, abundance, and community structure

Overall, 113 taxa and 57 species were identified during June–July 2019 in the Adriatic basin (Table S1). These values were slightly lower than those observed for the central Adriatic at 0–50 m depths, where 150 taxa were counted (Hure et al., 2018). Such differences may be only apparent and attributable to the storage method we used, as samples were kept frozen for subsequent stable isotope analyses, determining damage in many organisms, which was impossible to identify at the species or even genus level (Fanelli et al., 2011). Although this method may represent a considerable bias in species identification and biomass estimation, it allows for having an indication of both the community and food web structure (Fanelli et al., 2011, 2013; Rumolo et al., 2018). In terms of species abundance, the most representative species were Acartia clausi, Oithona similis, and Centropages typicus among copepods and the cladocerans Podon intermedius, P. polyphemoides, Penilia avirostris, Evadne tergestina, and E. spinifera, in agreement with previous studies on the mesozooplanktonic communities of the Adriatic basin (Fonda-Umani et al., 2005; Aubry et al., 2012).

Zooplankton abundances were higher, though very variable within sites, in the northern Adriatic Sea and slowly decreased moving towards the southern Adriatic, while biomass showed an increasing coastal–offshore trend, except for inshore southern stations characterized by a large within-sample variability. The abundance trend here found was also observed by Fonda Umani (1996) and can be explained by the influence of the Po River, which can determine a high nutrient input in the northern Adriatic favouring primary production and therefore zooplankton growth. Notwithstanding the general primary production reduction observed in the last years (Mozetič et al., 2010) in the northern Adriatic Sea, the area is still characterized by higher phytoplankton biomass with respect to the central and the southern basin because of the nutrients input from the Po River. Chlorophyll a concentration values from satellite data (Fig. S1, https://giovanni.gsfc.nasa.gov/giovanni, last access: 15 November 2021), analysed from 4 months before the sampling period to the survey simultaneous period (July 2019), indeed revealed a peak in primary production in May 2019, 2 months before the sampling period, in the area in front of the Po River delta, fuelling in turn zooplankton production (Aubry et al., 2012).

Although, in the northwestern Adriatic, offshore waters are less productive than inshore coastal waters and productivity of the inshore zone decreases southward away from the Po River's nutrient influx (Vollenweider et al., 1998), we did not here find significant differences in terms of abundance and biomass between inshore and offshore communities or for the interaction factors. Such differences were instead observed when we compared zooplanktonic communities' composition. Indeed, multivariate analyses evidenced a clear separation of samples as function of the sub-area and inshore vs. offshore locations, especially between the mesozooplanktonic community of the northern Adriatic from the other two. This was not surprising, as the northern Adriatic is characterized by shallower and colder waters than the rest of the basin and under the influence of riverine input, thus hosting a typical neritic community with coastal and estuarine elements. This area was dominated also by Acartia clausi, Oithona similis, cladocerans (mostly Evadne spinifera), copepodites (here comprised within the “Copepoda unid.” group), and gastropod larvae with some differences with respect to previous studies (Aubry et al., 2012) in terms of temporal shift of species maximum abundance. This could be related to the peak in primary production occurring in May 2019, although quite delayed with respect to the usual pattern of the area (Kamburska and Fonda-Umani, 2009) (see Fig. S1). Conversely, the southern Adriatic basin, except for the Gargano Promontory, being characterized by a narrow continental shelf and a steep slope, reaching high depths close to the coasts, was dominated by typical offshore species such as tunicates, chaetognaths, siphonophores, and Euchaeta spp. These results were supported by Fonda Umani (1996) that identified a clear distinction in zooplanktonic communities collected in offshore locations of the northern and central-southern Adriatic: the northern Adriatic was characterized by neritic communities, with moderate biomass, while the central and the southern Adriatic Sea were characterized by an “oceanic” community, with a higher abundance of carnivorous zooplankton, such as Euchaeta sp., a more oceanic carnivorous genus (Razouls et al., 2022), and Chaetognatha, a phylum of carnivorous animals abundant in open waters (Terazaki, 2000). Consistently, diversity was the greatest in the southern basin, with 80 taxa (out of 113) identified, likely due to the occurrence of both neritic and oceanic species in this area and comparable to other studies (Miloslavić et al., 2012) which also included deep stations.

4.2 Environmental drivers of zooplankton communities' variability

Separation among samples according to sub-areas and inshore and offshore locations were consistent with the main drivers influenced by the distance-based multivariate model, i.e. fluorescence and DO concentration, with fluorescence itself explaining 33 % of the variance. Fluorescence was strictly linked to freshwater inputs from the Po River and was likely responsible of the main separation between the northern Adriatic, with more coastal–estuarine zooplanktonic communities, and the central and southern Adriatic, with more oceanic zooplanktonic communities. Fluorescence was also found to be the main driver of zooplankton community in the northern Aegean Sea (Isari et al., 2006), another important area for small pelagic fishery. Several studies indicated that oxygen concentration could be a limiting factor for zooplankton growth and survival (Olson, 1987; Moon et al., 2006), with the inhibition of egg hatching in some copepod species (Roman et al., 1993). DO was found to be also the driving factor of zooplanktonic communities in the Strait of Sicily (Rumolo et al., 2016).

4.3 Food web structure of zooplankton communities

The trophic groups highlighted by cluster analysis fully agreed with putative trophic groups established a priori based on literature information and the previous classification on copepod functional traits (Hébert et al., 2016; Benedetti et al., 2016; Protopapa et al., 2019; Fanelli et al., 2011, and references cited therein; Rumolo et al., 2018; Conese et al., 2019). Conversely to similar works carried out on deep-sea zooplankton (Fanelli et al., 2009, 2011, 2013; Koppelmann et al., 2009), our analysis evidenced a trophic gradient from strictly herbivore species towards carnivory, with a general pattern of omnivory including taxa that may act both as primary consumers, eating phytoplankton or detritus particles or shifting to small prey, i.e. microzooplankton.

Moving from filter feeders–herbivores towards carnivores, a first group of omnivores, with phytoplankton as an important component of their diet occurred. This group contains both small-bodied calanoids that are numerically very important in the Mediterranean epipelagic (Temora stylifera; Mazzocchi et al., 2014) and also larger calanoids, some of which are strong vertical migrants, such as Calanus helgolandicus or Pleuromamma spp. (Andersen et al., 2001, 2004). These exhibit mixed feeding strategies, depending on the available food items. This is also the case of Meganichtyphanes norvegica, which can vary its diet regionally and with growth, showing a preference for phytoplankton in certain areas, in certain seasons, when juveniles (Schmidt, 2010; Fanelli et al., 2011), while preying exclusively on calanoids when adults, or depending on energy requirements (McClatchie, 1985). C. helgolandicus has been described as an herbivore species (Paffenhoffer, 1976), but some authors described density-dependent mortality through cannibalism in Calanus spp. as a form of population self-limitation (Ohman and Hirche, 2001), thus pointing out an omnivorous feeding behaviour.

Upscaling the pelagic food web, we found omnivore taxa that mostly prefer animal prey but that can shift to phytodetritus when prey was scarce or competition was high (Fanelli et al., 2011), such as Centropages typicus. C. typicus is an omnivorous copepod that feeds on a wide spectrum of prey, from small algae (3–4 µm equivalent spherical diameter) to yolk-sac fish larvae (3.2–3.6 mm length). It uses both suspensivorous and ambush-feeding strategies, depending on the characteristics of the prey (Calbet et al., 2007). Omnivorous copepods can display increased predatory behaviour in the absence of other food (Daan, 1988) and may actively target eggs even when phytoplankton is not limiting (Bonnet et al., 2004). Hyperiids (Lycaea pulex and Lestrigonus schizogeneios) also cluster with this group. Hyperiids generally use gelatinous substrate for reproduction and feeding, and some of them live in symbiosis (Gasca and Haddock, 2004) as parasites such as the genus Hyperia (now Lestrigonus). Finally, strictly carnivore species such as Euchaeta or chaetognaths clustered together with some siphonophores (Calycophorae). These species are known to prey on smaller copepods, doliolids (Takahashi et al., 2013), larvaceans (Ohtsuka and Onbé, 1989), and fish larvae (Yen, 1987).

The average enrichment between the different plankton taxa was greater than the mean value of 2.56 expected between adjacent trophic levels (e.g. Vanderderklift and Ponsard, 2003; Fanelli et al., 2011), pointing to the organization of mesozooplanktonic taxa in three trophic levels, from herbivore taxa (Nannocalanus spp., Gaetanus tenuispinus, and thaliaceans) positioned at the trophic level 2 to the highest-level species represented by large copepods and the mysid Siriella sp. located at the trophic level 4. Such results confirmed other findings (Fanelli et al., 2009, 2011) about the complexity of pelagic food webs and their lower trophic levels, calling attention to the appropriate compartmentation of zooplankton in ecosystem modelling with the final scope of small pelagic stock management (D'Alelio et al., 2016). Moreover, the predation on protozoa may have been overlooked by traditional stable isotope measurements, as phagotrophic protists do not necessarily follow the systematic ¹⁵N trophic enrichment that is well-established for metazoan consumers (Gutiérrez-Rodriguez et al., 2014). Thus, the uncertainties associated with missing one or more trophic levels using stable isotopes or other techniques significantly challenge our understanding of the pelagic food web structure.

Finally, based on our results, the isotopic composition of some species/taxa differed from the literature, as for the hyperiid Phronima atlantica. This species is reported as a carnivore, feeding on salp tissue (Madin and Harbison, 1977). However, Elder and Seibel (2015) also reported about them feeding on host mucus, which could lower their trophic position, being more similar to the basal source, i.e. the particulate organic matter (POM; Fanelli et al., 2011). Zoeae of Thalassinidea and Brachyura were also placed in this group, close to thaliaceans, that are herbivorous filter feeders (Madin, 1974).

4.4 Spatial variability in the isotopic composition of mesozooplankton from the Adriatic basin

Overall, stable isotope values of zooplankton differed significantly for both sub-areas and inshore vs. offshore factors considered, with δ¹⁵N values decreasing southward and δ¹³C showing more constant patterns across the basin but with large variability at NA. The presence of differences in isotopic signature of zooplankton between inshore and offshore locations has already been reported by other authors (Bode et al., 2003; Chouvelon et al., 2014; Espinosa-Leal et al., 2020), and it could be linked to the different contribution of terrestrial vs. marine sources of nitrogen and carbon moving from inshore to offshore waters and/or to different trophic dynamics between coastal and oceanic food webs. Here δ¹³C values were highly variable at NA (spanning from −15.9 ‰ to −22.6 ‰) in accordance with the wide array of food sources (i.e. marine and continental) available in the area due to the riverine inputs. Accordingly, the niche width of the zooplanktonic community in the area is the greatest, and SEAc decreased in CA and SA, where zooplanktonic community were likely sustained mostly by marine sources (Coll et al., 2007). Standard ellipses were mainly stretched along the x axis (δ¹³C) for NA and CA, showing a progressive decrease of the continental influence from the northern to the central Adriatic basin. The SEAc of SA was conversely mostly extended along the y axis (δ¹⁵N), likely because of the occurrence of a well-structured community with all TLs represented. The low δ¹⁵N range (and the general high δ¹⁵N values) observed for the NA community suggests a shift to omnivory in zooplanktonic communities in this area to avoid competition (Doi et al., 2010) in the high-density condition, as was generated after the phytoplankton bloom (Aubry et al., 2012) here observed in June.

5 Conclusions

This study represents the first application of the stable isotope approach to the analysis of the mesozooplanktonic food web at the Adriatic basin scale including both coastal and offshore communities. The results unveiled the presence of significant differences in zooplankton abundance, biomass, and community composition at the mid-spatial level, with the main differences observed between the northern Adriatic and the rest of the basin being due to the peculiar oceanographic conditions (i.e. cold waters) and the strong influence of the Po River. Such differences were also particularly evident in terms of isotopic composition, where a further separation between offshore and inshore communities was evident for the progressive increase of marine contribution to food sources for zooplankton in offshore communities. Such findings may represent a valuable baseline for food web studies encompassing species of a lower to high trophic level and against changes in oceanographic conditions under a climate change scenario, considering the rapid response of zooplankton communities to global warming.