Summer mesozooplankton assemblages in relation to environmental parameters in Kavala Gulf, northern Aegean Sea ; Comunidades zooplanctónicas estivales en relación con parámetros ambientales en el golfo de Kavala, norte del mar Egeo

Department of Zoology, School of Biology, Aristotle University of Thessaloniki, Thessaloniki, Greece. (AK) E-mail: aikatekm@bio.auth.gr. ORCID iD: https://orcid.org/0000-0003-3004-8868 (AA) E-mail: anna.artemiou@gmail.com. ORCID iD: https://orcid.org/0000-0002-7161-6879 (ACT) E-mail: atsik@bio.auth.gr. ORCID iD: https://orcid.org/0000-0002-9074-3259 (EM) (Corresponding author) E-mail: tholi@bio.auth.gr. ORCID iD: https://orcid.org/0000-0001-8194-2865


INTRODUCTION
Zooplankton play an important role in marine pelagic food webs, connecting primary producers, small pelagic fish populations and benthic communities, and thus transferring organic matter through the pelagic food webs (Fenchel 1988). Plankton and fish larvae communities are often structured in assemblages that are closely related to environmental characteristics (Cowen et al. 1993). Their response to alterations of the environmental conditions can be detected by observing changes in diversity, community structure and species distribution (e.g. Dortch et al. 1989, Richardson 2008. Moreover, zooplankton population dynamics can be indicative of environmental changes and may affect higher trophic levels (Beaugrand et al. 2003, Frederiksen et al. 2006. As an important part of the pelagic food web and fisheries ecology, zooplankton has thus been used in the evaluation and modelling of marine ecosystems (e.g. Triantafyllou et al. 2007, Petihakis et al. 2009).
The northern Aegean Sea, in the northeastern Mediterranean, is a complex and dynamic marine ecosystem (Tsagarakis et al. 2010). Its hydrodynamic complexity is closely linked to the influence of the Black Sea water masses, which outflow in the upper part of the water column through the Dardanelles Straits and are characterized by low salinity and temperature (Zodiatis and Balopoulos 1993, Zervakis and Georgopoulos 2002). Almost permanently, one branch of the Black Sea current follows a northward direction towards the Thracian Sea, where it is captured by an anticyclonic gyre formed around the island of Samothraki and is mainly restricted to the upper 0 to 20 m layer (Zodiatis and Balopoulos 1993, Zervakis and Georgopoulos 2002). The water circulation pattern in the northern Aegean has a considerable influence on the horizontal oceanographic variability of zooplankton assemblages (Isari et al. 2006, Zervoudaki et al. 2006). These hydrological conditions lead to high productivity, in contrast with the rest of the eastern Mediterranean Sea (Lykousis et al. 2002), making the northern Aegean one of the most productive areas in the region (Bosc et al. 2004).
Kavala Gulf and the surrounding area constitutes one of the most well-known fishing and nursery grounds of the northern Aegean Sea (Stergiou et al. 1997), which is particularly important for small pelagic fish populations such as European anchovy (Engraulis encrasicolus), the European pilchard (Sardina pilchardus) and the round sardinella (Sardinella aurita) (Stergiou et al. 2007, Sylaios et al. 2010, which may compete for resources (Tsikliras 2014). As zooplankton is the basic component in the diet of the early stages of fish and of the adults of small pelagic fish, zooplankton communities in Kavala Gulf are key organisms of the pelagic web (Tsikliras et al. 2005, Catalán et al. 2010, Karachle and Stergiou 2013. Although zooplankton communities have been widely studied in the northern Aegean Sea (Table 1), Kavala Gulf was not included among the study areas. Filling this gap could provide a valuable asset, adding to the theoretical knowledge of ecosystem functioning and fisheries modelling that could be used to refine or expand the ecosystem models that have already been developed for the area (Tsagarakis et al. 2010). The aim of the present study was to fill this knowledge gap regarding mesozooplankton in Kavala Gulf through the analysis of a fine grid of 17 stations during two summer cruises in 2002 and 2003. In particular, our analysis aimed to (i) demonstrate the spatial variability of major zooplankton species/ taxonomic groups and (ii) investigate the relation between sea surface abiotic factors and the structure of zooplankton communities.

Study site
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:  Siokou-Frangou and Papathanassiou 19911988-1990 Monthly 200 μm Alvanou 1999Saronikos Gulf 1984 Monthly 200 μm Siokou-Frangou 1996, Siokou-Frangou et al. 19981987 Monthly 200 μm Siokou-Frangou et al. 2009 Monthly 200 μm Christou 1998, Christou and Stergiou 1998, Siokou-Frangou et al. 2009 Monthly 200 μm Siokou-Frangou et al. 2009N Aegean 1997 Seasonal 200 μm Siokou-Frangou et al. 2002, Lykousis et al. 2002 Seasonal 200 μm Lykousis et al. 2002Lykousis et al. 1999 Snapshot 200 Ramfos et al. 2005 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).

Sampling process and laboratory analysis
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 sam-pling 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 3 ) 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.

Data analysis
Mesozooplankton abundance data were expressed as number of individuals per volume unit (m 3 ) 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 (SIM-PROF) 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 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).

Environmental parameters
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).

Community composition and structure
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 spinif-era 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 (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 (ANO-SIM, 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   (Fig. 4). Groups B and C differed significantly, as shown by pairwise comparison (ANO-SIM, 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).

Relationships between physical-chemical parameters and zooplankton community structure
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   Figure 4B. Taxa abbreviations are shown in Table 3.
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).

DISCUSSION
The present study focuses on the community structure of mesozooplankton in Kavala Gulf in the summer of 2002 and 2003. Although our dataset is relatively old, this is to the best of our knowledge the first published study on mesozooplankton in Kavala Gulf. Mesozooplankton sampling was part of a fisheries survey in the area, whose data on hydrology and fisheries have already been published (Tsikliras et al. 2005, Tsikliras 2014. Since Kavala Gulf is an important fisheries and nursery ground (Stergiou et al. 1997), data on mesozooplankton could contribute to the theoretical knowledge of ecosystem functioning in the area. Moreover, though our dataset lacks seasonality, it could be useful for comparisons with present and future data, and for further meta-analyses in fisheries and ecosystem modelling through the refining or expanding of models that have already been developed for the area (Tsagarakis et al. 2010).
The hydrological data collected during the present study were limited to the surface layer. Nevertheless, they are in accordance with other studies in the northern Aegean (e.g. Isari et al. 2008). The differences recorded between the two sampling years for salinity are attributed to the rainier days just before the 2003 survey (Tsikliras et al. 2009). The composition of the zooplankton community found in the present study is typical of coastal areas in the Aegean and Ionian Seas, such as Saronikos Gulf (Siokou-Frangou 1996), Patraikos Gulf, the northern Euboikos Gulf, Pagasitikos Gulf (Ramfos et al. 2005) and the Thracian Sea (Isari et al. 2007). Moreover, a similar zooplankton community composition has been recorded during the summer in other coastal ecosystems of the Mediterranean and Black Seas (Razouls et al. 2018). As sampling was originally planned for ichthyoplankton collection, zooplankton samples were not obtained with the typical WP-2 net but with a net of a wider mesh size (250 μm). Thus, differences in community composition are expected to be mainly in small-sized species, due to selectivity of large-sized species by the gear. Furthermore, 250 mesh size nets have also been used for studies in the northern Aegean Sea (Table 1), making our results still comparable. Seasonal fluctuations in zooplankton abundance of the coastal Mediterranean regions are characterized by two maxima, one in spring and one in autumn (Stergiou et al. 1997). In the Mediterranean coastal areas, the spring peak is common in the seasonal fluctuations of zooplankton abundance (Siokou-Frangou 1996), while the autumn peak is observed mainly in areas strongly affected by anthropogenic activities, often during the period June-September (Siokou-Frangou 1996, Alvanou 1999. We acknowledge that our dataset is limited to the summer period and is not able to reflect the annual fluctuations of mesozooplankton communities. Nevertheless, the mesozooplankton abundance values recorded, after comparison with summer data in other coastal areas of the Mediterranean Sea, seemed to correspond to the autumn peak (Kovalev et al. 2003). In order to broaden our knowledge on Kavala Gulf a more intensive annual sampling scheme should be designed in the future.
Total mesozooplankton abundance values in both sampling years were relatively higher than those of samplings in the surrounding area of the Thracian Sea (Isari et al. 2006) and in other coastal systems in the Mediterranean [e.g. Gulf of Naples (Mazzocchi and Ribera d'Alcalà 1995); Gulf of Trieste (Lipej et al. 1997); Ionian and central Aegean Sea (Ramfos et al. 2005)] and the Black Sea (e.g. Sazhina 1964). This may be linked to the special conditions occurring in this semi-enclosed gulf, namely the effects of the Black Sea water, the nutrient rich intake from the nearby river outflows and the urban and industrial activity in the area (Friligos 1985, Kardaras 2005, Sylaios et al. 2005, Isari et al. 2006. A water circulation study that was synchronous with our sampling in Kavala Gulf states that northern Aegean water enters westwards in the Gulf from the Strait of Thassos (near station 10, Fig. 1), forming a cyclonic sea surface pattern (Tsikliras et al. 2009). Moreover, in Kavala Gulf cold water masses with low salinity are provided by the outflow of the Nestos River during winter and spring (Kardaras 1998), whereas in the summer period water circulation is influenced by the low-salinity water masses entering from the Dardanelles (Kardaras 1998, Tsikliras et al. 2009). The influence of the Black Sea water on zooplankton communities has been previously studied in the northern Aegean Sea, showing that it favours the dominance of Cladocera and small-sized Copepoda due to its lower salinity (Isari et al. 2006(Isari et al. , 2007(Isari et al. , 2011, as is the case with the mesozooplankton communities of Kavala Gulf. Moreover, the abundance recorded in 2002 was higher than in 2003, while the reverse pattern is observed in the ichthyoplankton data (Tsikliras et al. 2009). This reverse relationship may be indicative of fish larvae predation on the mesozooplankton communities (Cushing 1990). However, inter-annual variability of the mesozooplankton-ichthyoplankton trophic link and longer time-series would be needed to support this hypothesis with more evidence.
The dominant zooplankton group in both sampling years were Cladocera, followed by Copepoda and Tunicata (Fig. 4). Cladocera dominate coastal areas in the summer period according to similar studies, contributing to the summer peak [e.g. Christou and Stergiou 1998;Strymonikos and Ierissos Gulfs (Michaloudi 1999); Thracian Sea (Isari et al. 2007)]. They are favoured by high temperatures, low salinity (Moraitou-Apostolopoulou andKiortsis 1973, Christou andStergiou 1998) and low depth (Siokou-Frangou 1996), which are conditions encountered in Kavala Gulf during the summer. This trend was reflected in the RDA (Fig.  5), with stations being grouped by depth and surface salinity. In 2002, mesozooplankton communities were positively correlated with the sea surface salinity values, which were significantly higher, whereas in 2003 they were negatively correlated. Although salinity was significantly differentiated between the years studied, it was not an actual forming factor of the communities' distinctness because of the narrow range of the measured values (1.89±0.95). According to Siokou-Frangou et al. (1998), temperature is the main driving factor of mesozooplankton community formation.
In our study, the Cladocera community consisted almost entirely of P. avirostris, which is a highly efficient filter feeder feeding mainly on phytoplankton, such as diatoms, nanoflagellates (HNF) and bacterioplankton (Turner et al. 1988, Kim et al. 1989, Atienza et al. 2006. P. avirostis peaks may be favoured by the presence of Cymodocea and Posidonia meadows in Kavala Gulf (Orfanidis et al. 2010), taking advantage of the dissolved organic carbon enrichment in the pelagic food web and the microbial loop that is favoured in such habitats (Barrón et al. 2004, Barrón andDuarte 2009). Moreover, massive populations of P. avirostris may be a result of the combination of its ability to exploit the available food resources and of its life cycle characteristics (a short life cycle and parthenogenetic reproduction) (Colton 1985, Valentin andMarazzo 2003). In conclusion, the zooplankton community of a semi-enclosed gulf in the northern Aegean Sea during the summer was mainly dominated by Cladocera, small pelagic Copepoda and Tunicata. The high abundance values might have been favoured by the Black Sea water and river input in the upper layer of the water column. Our results on the mesozooplankton communities of Kavala Gulf may provide valuable information concerning the function of the pelagic food web in the gulf and contribute to the construction and/or improvement of models used in fisheries management of economically important fish stocks of the area. APPENDICES Appendix 1. -Results of the similarity percentage analysis and the pairwise analysis of similarity of variance among the four groups indicated by the hierarchical cluster analysis (Fig. 4B). Taxa with a cumulative contribution of more than 50% are indicated in bold. Taxa abbreviations are shown in Table 3.