Kavala Gulf (40.88056°N, 24.41667°E) is a semi-enclosed area located on the continental shelf of the northern Aegean Sea in the eastern Mediterranean Sea (Fig. 1). It is connected to the Aegean Sea through its main mouth in the south, which is wide (20 km) and deep, and in the east, through the smaller (7.3 km wide), shallower mouth, the strait of the island of Thassos (Fig. 1). As Kavala Gulf is shallow (mean depth: 34 m, maximum depth: 50 m), with a flat seabed, the topographical variations inside it are minimal (Fig. 1). The hydrology of the gulf is influenced by the fresher and cooler Black Sea water, which oscillates seasonally in occurrence, intensity and direction (Sylaios et al. 2005), and forms a surface current flowing to the northern coastline of Greece (Poulos et al. 1997). Kavala Gulf also receives fresh water inflows from the Nestos River, though these have been significantly reduced since the construction of two hydroelectric dams in 1997 (Sylaios et al. 2005).

Mesozooplankton samples and environmental data (depth, temperature, salinity and dissolved oxygen) were collected during two sampling periods in Kavala Gulf in the summers of 2002 and 2003. The surveys were carried out across a fine scale grid of 17 stations (Fig. 1) during the period of thermal stratification, which coincides with the spawning period of European anchovy (Engraulis encrasicolus) and round sardinella (Sardinella aurita). These two small pelagic clupeoids, together with the winter spawning European pilchard (Sardina pilchardus), are the main zooplanktivorous fish in the area (Somarakis 2005, Tsikliras et al. 2005). Based on their depth and distance from the shore, the 17 stations had been previously grouped (Tsikliras and Koutrakis 2011) into shallow (maximum depth less than 30 m: stations 6, 7,8, 9, 10) and deep (maximum depth above 30 m: stations 1, 2, 3, 4, 5, 11, 12, 13, 14, 15, 16, 17).

Sea surface temperature, salinity and dissolved oxygen were measured at each station using CTD sensors (Ocean 301, Idronaut). The hydrographic sampling profiles of all these parameters are available in Tsikliras et al. (2009). For the mesozooplankton collection, a paired bongo net sampler was used, with two 0.6 m diameter frames fitted with 250 μm mesh conical nets. The cod-end consisted of a plastic container incorporating a window of 250 μm mesh net that allowed water to escape. A flowmeter was centrally mounted on the opening of each frame to the net to estimate the volume of water (m³) flowing through the net. The sampler was deployed in a double oblique tow from the surface down to 1 to 2 m from the bottom, thus forming a V-shaped dive profile. Retrieval speed was kept at the same rate (20 m min^–1) for all samplings, which were all carried out in daylight (between 09:00 h and 18:00 h). Upon recovery of the sampler, the plankton was washed into a jar and fixed with 4% buffered formalin solution.

After the removal of fish larvae for further investigation of ichthyoplankton assemblages (Tsikliras et al. 2009, Tsikliras and Koutrakis 2011, Tsikliras 2014), the major mesozooplankton taxa (e.g. Copepoda, Cladocera, Doliolidae, Appendicularia, etc) were identified using the appropriate taxonomic keys. Copepoda (only adults) and Cladocera were further identified down to species level when possible, whereas copepodites were grouped together and were not further identified. For the abundance estimation, each sample was enumerated in aliquots under the stereoscope. Subsamples were taken with a pipette, their size ranging from 1/2 to 1/30 of the total sample depending on the sample’s density. The total number of counted organisms was on average 1844.47±646.58 ind.

Mesozooplankton abundance data were expressed as number of individuals per volume unit (m³) by dividing the individuals per sample by the volume of the filtered water. The standardized numbers were used to calculate the percentage contribution of each taxon to the total abundance. In order to identify the similarity of mesozooplankton communities among stations and years, hierarchical cluster analysis based on the Bray-Curtis similarity index (Bray and Curtis 1957) was performed on the mesozooplankton taxa’s abundance matrix. Abundance data were log (x+1) transformed in order to approach normality and homogeneity of variances and to reduce bias due to highly abundant groups. CLUSTER was run using group-average linking. The similarity profile (SIMPROF) permutation test option (default settings of 999 permutations and significance level 0.05) was applied to indicate significant groups in the resulting dendrogram. The similarity analysis routines, analysis of similarity (ANOSIM) and similarity percentage analysis (SIMPER) were used to test the significance levels and sources of variance between the various zooplankton assemblages associated with the groupings identified in the hierarchical cluster analysis. The above analyses were conducted with the Plymouth Routine In Multivariate Ecological Research (PRIMER) v.5 software package (Clarke and Gorley 2001). Permutational multivariate analysis of variance (PERMANOVA), based on Euclidean distance resemblance matrix, was used to test for statistical differences in the mesozooplankton abundances between the sampled periods (July 2002 and 2003) and between water depth categories (stations at less than 30 m and more than 30 m depth). Analyses were performed using the Paleontological Statistics software package (Past v.3.15) (Hammer et al. 2001).

Direct ordination analyses were used to assess significant relationships between biological and environmental data. All environmental variables (temperature, salinity, depth and dissolved oxygen) were log (x+1) transformed. Previously, a detrended correspondence analysis (DCA) was performed and, as biological data showed a linear response with respect to environmental gradients, a redundancy analysis (RDA) was applied. All environmental parameters with an inflation factor smaller than 20 were included in the analysis as explanatory variables (Ter Braak and Verdonschot 1995); the abundance of each zooplankton taxon was included as a response variable. The statistical significance of the variation in the parameters and the overall significance of the ordination were tested with the Monte Carlo permutation test (as default settings of 499 unrestricted permutations; P<0.05). Ordination analyses were performed using the CANOCO program, version 4.5 (Ter Braak and Šmilauer 2002).

The values of sea surface environmental parameters (temperature, salinity, oxygen) and maximum depth of the sampling stations in Kavala Gulf are given in Table 2. Salinity and oxygen differed significantly between the years (ANOVA; P<0.005), whereas temperature did not (ANOVA; P>0.005).

Overall, a total of 9 main holoplanktonic and 6 main meroplanktonic groups were recorded (Table 3). Nineteen taxa were further identified, 12 of them down to species level and 7 to genus level (Table 3).

Total abundance ranged from 41069 to 164939 ind m^–3 in 2002 and from 5508 to 86747 ind m^–3 in 2003 (Fig. 2A). Total abundance values were significantly higher in 2002 than in 2003 at all sampling stations (pseudo-F=7.27, P<0.05) (Fig. 2B) except station 10, where abundance was higher in 2003 (Fig. 2A). No significant difference was observed between shallow and deeper stations (pseudo-F=0.14, P>0.05) according to the PERMANOVA procedure. Cladocera were the dominant group at all sampling stations, followed by Copepoda and Appendicularia in rank order, whose cumulative contribution exceeded 79% of mesozooplankton abundance (Table 3). Among the four identified Cladocera species, Penilia avirostris showed a marked predominance (Fig. 3A, B), followed by Evadne nordmanni Lovén, 1836, Evadne spinifera P.E. Müller, 1867 and Pseudevadne tergestina (Claus, 1877) (Table 3). Despite small differences in the species-specific distribution patterns, the bulk of Cladocera were found more or less in the eastern part of the gulf in 2002 and in the shallow stations of the southeastern part near Thassos Island in 2003 (Fig. 3A, B). Copepods, being the second dominating taxonomic group, were mainly represented by juveniles and small-sized zooplankters [e.g. the calanoida Acartia (Acartiura) clausi Giesbrecht, 1889 and Paracalanus parvus (Claus, 1863), the cyclopoida Oithona plumifera Bair 1843 and Corycaeus spp. Dana 1845, and the harpacticoida Macrosetella spp. Scott A., 1909] (Table 3), whose distribution did not follow a discrete pattern (Fig. 3C, D). Copepods’ mean contribution to mesozooplankton abundance was higher in 2003 than in 2002 (16.07% and 11.10%, respectively) (Table 3). Appendicularia were represented by the genus Fritillaria Fol, 1872, and the genus Oikopleura Mertens, 1830, which contributed 2.92% and 1.87% in 2002 and 3.01% and 2.20% in 2003, respectively (Table 3). The hierarchical cluster analysis (Fig. 4) revealed that the mesozooplankton community structure differed significantly between stations with varying maximum depth and sampling period. Four groups were generated (ANOSIM, R=0.746, P=0.001; Appendix 1), with no apparent pattern between shallow and deep stations (Fig. 4). Station 10 fell into different groups (Groups A and D in 2002 and 2003, respectively) (Fig. 4). This was mainly due to Mysidae, Ostracoda and Amphipoda abundance for Group A (SIMPER, average dissimilarity ≤27.18%; Appendix 1) and to Cirripedia, Echinoderma, Decapoda and Polychaeta abundance for Group D (SIMPER, average dissimilarity 14.48%; Appendix 1). However, according to ANOSIM, these groupings had no statistical significance (ANOSIM, Group A, R≥0.986, P>0.05 and Group D, R=0.645, P>0.05; Appendix 1). The rest of the sampling stations of 2002 and station 9 of 2003 were grouped together in the upper group of the cluster (Group B), while the stations of 2003 were placed in Group C (Fig. 4). Groups B and C differed significantly, as shown by pairwise comparison (ANOSIM, R=0.703, P =0.01; Appendix 1). Their dissimilarity was mainly due to the contribution of Ostracoda, Echinoderma, Polychaeta and Doliolidae (SIMPER, average dissimilarity 15.06%; Appendix 1), whereas the similarity among stations in the same groups was in both cases mainly due to Cladocera, Copepoda, Appendicularia and Doliolidae or Siphonophora in rank order (SIMPER, average similarity Group B 91.50% and Group C 89.72%; Appendix 1).

To assess significant relationships between surface environmental data and mesozooplankton community structure, ordination analysis of taxa assemblages, expressed in terms of biomass, was conducted (Fig. 5). In the diagnostic DCA, the highest value of the length of gradient of axis was 0.819, indicating that the relationship between mesozooplankton and environmental variables was linear (Ter Braak and Šmilauer 2002), and an RDA was performed. The significant environmental variables (P<0.05) included in the RDA were salinity and maximum depth. Although dissolved oxygen was not a significant variable (P=0.06), it was included in the analysis as an explanatory variable of zooplankton distribution. The Monte Carlo test confirmed that the selected RDA model was significant (F ratio=9.531, P=0.002). The eigenvalues of the first two axes were 0.241 and 0.072, and both of them together explained 90.4% of the variation in species-environment relation. The first axis, which accounted for a total variance of 69.6%, was negatively and strongly correlated with sea surface salinity (r=0.9768). Moreover, mesozooplankton communities in 2002, being in contrast with those of 2003, were negatively correlated with Axis 1 (Appendix 2) and consequently positively correlated with sea surface salinity. Axis 2 showed 20.8% variation and was strongly and negatively correlated with maximum depth (r=–0.9768). Dissolved oxygen was strongly and negatively correlated with the third axis (r=–0.9355), which accounted for a total of 9.6% variation. Cumulative fit indicated that the predominant mesozooplankton groups affiliated with Ostracoda and Cladocera displayed stronger correlations with Axis 1 than the rest of the groups (Appendix 3). Ostracoda, Doliolidae, Cladocera, Appendicularia and Mysidae were notable taxonomic groups with higher correlations with Axis 2 (Appendix 3).