Environmental Impacts on Zooplankton Functional Diversity in Brackish Semi-Enclosed Gulf

: Zooplankton as an essential component in the pelagic food web are directly linked to pelagic ecosystem functioning. Therefore, comprehension of zooplankton functional diversity (FD) and its responses to environmental changes is crucial for ecosystem-based view. To identify FD responses to environmental drivers, we analysed 25 years of summer data on the brackish mesozooplankton community (including rotifers, cladocerans, copepods, and meroplankton) from the eutrophied, shallow Gulf of Riga (Baltic Sea). We established that within the Gulf of Riga, open waters are notably different from coastal regions based on the dynamics of hydrological conditions (temperature, salinity), cyanobacterial dominance, abundance of mesozooplankton functional groups, and mesozooplankton FD indices. Competition over resources in combination with hydrodynamic features and predation by adult herring were seemingly the central structuring mechanism behind the dynamics of FD metrics (richness, evenness, divergence, and dispersion) within coastal mesozooplankton communities. Whereas predation by young herring was an important driver only for the mesozooplankton communities in the open waters. Cyanobacterial dominance, used as a proxy for food quality and availability, had no effect on summer mesozooplankton FD metrics.


Introduction
Organism survival is dependent on morphological and behavioural characteristics (traits) and associated responses to abiotic filtering and biotic interactions (predation, competition) [1,2]. Therefore, comprehension of links between trait occurrences and environmental changes is of extreme importance for ecosystem assessment and management. Functional trait approach can provide ecologically more profound information on underlying mechanisms than taxonomy-based analysis [3,4]. Functional diversity (FD) metrics consider functional traits of organisms and thus describe ecological roles and strategies of the occurring species and also characterise community assembly [5,6], ecosystem processes [7,8], and resilience [9,10].
Zooplankton are an essential component in the pelagic food web, linking primary production to higher trophic levels, thereby playing a key role in the functioning of aquatic ecosystems [11,12]. FD of zooplankton communities has been addressed in studies from different types of waterbodies worldwide [13], including the brackish Baltic Sea [14][15][16][17]. Due to the brackish environment where marine and limnic species meet their physiological limits, zooplankton of the Baltic Sea are characterised by relatively low species richness [18] and also low functional diversity [19]. However, FD of the Baltic zooplankton changes seasonally and along gradients of temperature, salinity, and depth [14,15].
The Gulf of Riga, a highly eutrophic subbasin of the Baltic Sea [20], has a strong impact from riverine discharge [21], strengthening the environmental variability. Thus, zooplankton community composition and assembly in the Gulf of Riga is strongly driven by abiotic factors [22][23][24][25]. Hydrological conditions and climate variability are also recognised as the main drivers shaping zooplankton FD in the area [15,17], albeit, the effects of biotic interactions on zooplankton FD are left underexplored.
Pecuchet et al. [17] described long-term, ecosystem-wide functional changes in three Baltic Sea subbasins, including the Gulf of Riga, and aside from abiotic factors and zooplankton, they also considered phytoplankton, zoobenthos, and fishery-related parameters. The study revealed strong causal links between multi-trophic community weighted mean, salinity, fisheries, and cod abundances, and expressed the need for further analysis of trophic interactions to assess the results explicitly. However, the study was limited to copepods, thus the results cannot be extrapolated to the whole zooplankton community. Indeed, copepods are the main food source for herring Clupea harengus [26], the main commercial pelagic fish of the region. However, rotifers and cladocerans significantly affect primary producers [27] and serve as a linkage to the microbial food web [28,29]. They provide substantial support for ecosystem production and functioning, especially in shallow, eutrophied, and coastal regions [30] such as the Gulf of Riga [31]. Hence, analysis of FD dynamics for the whole zooplankton community, rather than copepods alone, may provide additional information about the underlying processes shaping community structure.
Our study aims to identify the dynamics of mesozooplankton FD indices and their responses to abiotic and biotic drivers of summer communities (including rotifers, cladocerans, copepods, and meroplankton) in coastal and open water habitats of the Gulf of Riga. We hypothesise that FD of summer zooplankton is mainly shaped by biotic interactions due to (i) pronounced activity within and from the microbial loop [31], e.g., high cyanobacterial biomass [32,33] and (ii) intense predation pressure from herring [34,35]. Using the Gulf of Riga as a case study area, we seek answers to the questions: What are the potential structuring mechanisms of mesozooplankton FD in eutrophied shallow brackish waters and how do they change between coastal and open water habitats?

Study Area and Sampling
The Gulf of Riga is a semi-isolated basin with a north-westward salinity gradient. Salinity ranges between 0.5-2.0 PSU in its southern area, where the largest freshwater discharge occurs [21], to 7.0 PSU in northeast regions, but generally varies from 5.0 to 6.5 PSU [36]. The Gulf of Riga is located in the temperate zone, therefore, the changes in the ecosystem follow a typical temperate seasonal cycle with a productive period from March/April to October. Summer is characterised by the highest sea surface temperature, strongest stratification [37], and increased values of phytoplankton, zooplankton biomass, and diversity [15,38].
We analysed Latvian National marine monitoring data from four open water and nine coastal sites ( Figure 1). The open water sites (>40 m deep) had a stratified water column during the summer, with thermocline at approximately 15 m in depth. Coastal sites (depth < 15 m) were divided into three coastal areas based on their geographical location (western, southern, and eastern) and according to freshwater impacts (north-westward salinity gradient). Data collected in summer (25)(26)(27)(28)(29)(30)(31)(32)(33)(34)(35)(36)(37)(38) week of the year) during the period from 1993 to 2017 were used. The data frequency is shown in Table 1.
Mesozooplankton samples were collected from the whole water column following recommendations by the Baltic Marine Environment Protection Commission (HELCOM) [39] (and previous versions). The only deviation from the recommended methodology was the omission of a flow meter on the WP2 net (mesh size: 100 µm). Therefore, the sampled volume was calculated using the formula for the volume of a cylinder (diameter: 0.57 m), taking into account the depth of the sampled site and the angle of a winch wire during sampling. Mesozooplankton samples were preserved in 4% buffered formaldehyde in seawater. Zooplankton taxa ( Table 2) were identified to the lowest possible taxonomic level using atlases [40,41] and ICES leaflets [42]. The majority of holoplanktonic taxa were identified to the species level, except for Cyclopoida copepods which were identified to order. Meroplanktonic larvae were identified to class (e.g., Polychaeta, Bivalvia) and only Amphibalanus larvae were identified to genus level. Phytoplankton samples were collected from surface waters (0-10 m) using a plastic hose and processed and analysed following HELCOM recommendations [39] (and previous versions). Herring population data were obtained from an annual report [43] by the Baltic Fisheries Assessment Working Group (WGBFAS) of the International Council for the Exploration of the Sea (ICES). Sea surface temperature (SST) and sea surface salinity (SSS) were measured using a CTD water probe, calculated to obtain the weighted average of the top 10-m layer.

Mesozooplankton Traits and Functional Groups
We described FD of summer zooplankton community using functional traits on body size, feeding behaviour, and trophic role (Table 3)-a combination that is useful for determining "who eats whom" and identifying potential biotic interactions [45]. Feeding mode was defined based on classification described by Kiørboe [46], except for 'raptorial' mode that was added and combined with the C category, as it resembles 'cruising' feeding mode. 'Raptorial' feeding mode was used to describe the feeding of predatory cladocerans (Onychopoda). 'Mixed' feeding mode was assigned to organisms that show the ability to switch between two or more feeding modes. Mean individual length and prey size values were categorised. Taxa-wise values are shown in Table 3.
We applied hierarchical agglomerative clustering to define functional groups of zooplankton communities. Gower's dissimilarity matrix was used as an input for clustering. The agglomeration method and optimal cluster number were identified by evaluating cluster validation statistics (R software v3.6.1 [47]; package 'fpc' v.2.2-5 [48]). Ultimately, the 'complete linkage' method and five clusters were chosen ( Table 3). The first group (G1) compounded small-sized herbivorous filter-feeders. The second group (G2) consisted of small-sized carnivorous copepods from the order Cyclopoida but were not used for subsequent analyses due to limited occurance data. The third (G3) and fourth (G4) group were comprised of omnivorous taxa. Large omnivorous taxa feeding on larger prey were grouped in G3, whereas small omnivorous taxa feeding on small prey were assigned to G4.

Functional Diversity Indices
Four FD indices (FRic, FEve, FDis, and FDiv) were calculated by the means of the dbFD function from the 'FD' package [70,71] using R v3.6.1. Gower's dissimilarity between species based on traits (Table 3) were calculated and used as an input value. All traits were weighted equally during dbFD processing, whereas indices were weighted by the relative abundances (ind/m 3 ) of species.
Functional richness (FRic) represents the amount of niche space filled by taxa in the community. Functional evenness (FEve) describes the evenness of abundance distribution within the filled niche space. Functional divergence (FDiv) quantifies how the trait values are spread along the range of a trait space by characterising the degree to which abundance distribution in niche space maximises divergence in functional traits [72]. FDiv increases if the most abundant taxa occur towards extreme points of the filled niche space. Functional dispersion (FDis) is a multidimensional index based on multi-trait dispersion [70] and is mathematically similar to the Rao quadratic entropy [73]. The FDis metric is the average distance of individual taxa to the abundance-weighed centroid of the community trait space [5].
Indices FRic, FDis, and FDiv were adjusted with a matrix-swap null models as noted by Mason et al. [5] and Swenson [74]. Community matrix randomisation with the independent swap algorithm (1000 permutations) was performed using randomizeMatrix from the 'picante' package [75]. In further analysis, only standardised effect sizes of FRic, FDis, and FDiv indices were used, yielding SESFRic, SESFDis, and SESFDiv. FEve were not adjusted (following Mason et al. [5]) and used as calculated by the dbFD.

Data Analysis
We analysed the dynamics of mesozooplankton FD in relation to several environmental parameters: SST and SSS, the ratio of cyanobacterial biomass to other phytoplankton biomass (cyano:other ratio), total herring biomass (her-totBio), spawning stock biomass (her-totSPbio; 2+ year old herring), herring recruitment (her-Rec; one-year old herring), and abundances of mesozooplankton functional groups as defined in Table 3. The data for all environmental variables (except for herring parameters) are compiled in the Supplementary Material.
Multiple Factor Analysis (MFA) was used to identify spatial environmental gradients of the Gulf of Riga, to define differences in the study areas (open waters, eastern coastal area, western coastal area, and southern coastal area) and to describe co-varying environmental parameters. MFA input data were FD indices, SST, SSS, cyano:other ratio, and abundances of mesozooplankton functional groups, conducted using the 'FactoMineR' package [76]. Coastal sites were pooled accordingly to their location ( Figure 1) and freshwater impacts [21]. The mean value for each coastal area was used in MFA as the sites were surveyed irregularly, whereas open waters were represented site-wise. Each year was described by one value from each of the coastal areas (eastern, western, and southern; Figure 1) and four values from open waters (sites O1-O4; Figure 1). The abundance data of mesozooplankton functional groups and cyano:other ratio were normalised (boxcox transformation) prior to analysis using the 'bestNormalize' package [77].
To evaluate responses of mesozooplankton FD to variations in hydrological parameters (SST, SSS), cyanobacterial dominance (cyano:other ratio), dynamics of zooplankton functional groups, and predation pressure (her-totSPbio, her-totBio, and her-Rec), we used Generalised Additive Modelling (GAM) with the smoothing spline function constrained to five degrees of freedom. To assess the predictive performance of the relationship, a non-random 9:1 train-test split was conducted. The response was assessed on data from 1993-2014, whereas the test data were from 2015-2017. Mean values per year, per region (open waters, coastal area) were used during the analysis due to the low spatial resolution of herring data that describe the Gulf of Riga population as a whole. In addition, the MFA results supported the pooling of all three coastal areas (see Section 3.1). GAMs were calculated and visualised using functionalities of the 'INDperform' package [78].

Environmental Gradients
According to MFA results, the Gulf of Riga open waters are notably different from all of the studied coastal areas based on differences in hydrological conditions, cyanobacterial dominance (cyano:other ratio), abundance of mesozooplankton functional groups, and mesozooplankton FD indices (SESFRic, FEve, SESFDiv, and SESFDis). SST (the main driver of Dim1; Table 4) and cyano:other ratio (the main driver of Dim2; Table 4) showed evident increasing values from coastal areas to open waters ( Figure 2). Salinity gradient, although irrelevant for Dim1 and Dim2, was represented in Dim3 and Dim4 (Table 4) and it showed a negative correlation to G3 and G1 abundances ( Figure 2B).
The ensemble of mesozooplankton FD indices showed increased values in open waters compared to coastal areas ( Figure 2C,D), except SESFRic which contributed weakly to the first two MFA dimensions (<2%; Table 4). Abundances of the G1 and G3 functional groups were negatively correlated with SESFRic ( Figure 2B). Indices FEve and SESFDis contributed to Dim1 (19.4% and 15.0%, respectively; Table 4), which related to temperature gradient. SESFDiv and SESFDis were important factors for Dim2 (11.5% and 10.7%, respectively; Table 4) indicating an increase in mesozooplankton functional divergence and dispersion towards open waters.  Table 3.
Characteristics of mesozooplankton FD and other environmental factors ( Figure 2) were similar between all coastal areas. Thus, data from eastern, southern, and western coastal sites were pooled into one data set for further analyses representing the coastal region in general.

Long Term Dynamics
Herbivorous filter-feeders (G1) dominated the pelagic ecosystem of the Gulf of Riga in both coastal and open waters during the studied period ( Figure 3). Small omnivores (G4) were the second most frequent functional group but G5 (that mostly consisted of calanoid E. affinis) was found in lesser amounts, similarly as large omnivores (G3).  Table 3.
The temporal variability of mesozooplankton abundance was higher in coastal waters while the number of observed taxa was rather stable but slightly lower than in open waters ( Figure 4A). SESFRic was almost identical in coastal and open waters ( Figure 4C) and had low values in the early and mid 1990s. After a rapid increase towards the 2000s, SESFRic values stabilised and stayed approximately the same thenceforth. A minor decrease in SESFRic occurred during 2014-2017 ( Figure 4C), seemingly mirroring the decrease in the number of observed taxa ( Figure 4A) and the rapid increase in the total abundance of zooplankters ( Figure 4B). Figure 4D). SESFDiv values were high during the late 1990s and early 2000s, indicating that the mesozooplankton community consisted of functional groups that were more distinct from one another (located closer to the extreme ends of the trait space). After the year 2002, the distinction between functional groups slightly decreased implying that the abundant traits are becoming more similar ( Figure 4E). Based on SESFDis values, mesozooplankton communities in coastal areas appeared functionally less dispersed than those in open areas. However, annual dynamics were similar between area type, noting that the responses to disturbances were most likely analogous. From 2003 to 2008, SESFDis showed continuously high values indicating the most functionally diverse period (considering the traits included in the analysis) ( Figure 4F).

Mesozooplankton FD Responses to Environmental Factors
In the coastal areas, SESFRic and FEve indices responded negatively to an increase in herring biomass (both total biomass and spawning biomass) ( Figure 5A,B,E,F) and herbivorous filter-feeders ( Figure 5C) whereas high numbers of G5 zooplankters had a positive impact ( Figure 5H). Additionally, FEve decreased with increasing SST ( Figure 5G). The least responsive FD index SESFDiv exhibited only one weak relationship and that was in relation to the abundance of G4 ( Figure 5I). Conversely, SESFDis responded to cyanobacterial dominance ( Figure 5J) and to abundances of all included mesozooplankton groups (G1, G3-5; Figure 5K-N). Despite the observed wide variety of response relationships only, three of them demonstrated acceptable predictive ability (nrmse < 1.0; Figure 5A,H,I).
In open waters, SESFRic negatively correlated with salinity ( Figure 6A); although the response was almost linear (edf = 1.17) and significant (p = 0.010), its predictive ability was extremely low (nrmse = 17.541). Herring recruitment and spawning biomass demonstrated contradicting impact on the FEve index ( Figure 6B,C). Increasing herring recruitment had a positive effect on FEve, whereas increasing spawning biomass had a negative impact. Additionally, FEve was positively related to an abundance of G3 ( Figure 6D). The SESFDis index reacted to dynamics of mesozooplankton functional groups ( Figure 6E-G), but SESFDiv did not respond to any of the factors included in the analysis.  . her_totBio-total herring biomass, 1000 tonnes; her_totSPbio-herring spawning biomass, 1000 tonnes; cyano_other_ratio-cyanobacteria-to-other phytoplankton biomass ratio, box-cox transformed; G1-G5-abundance of group G1-G5, box-cox transformed (functional groups are described in Table 3); SST-sea surface temperature,°C.  Table 3); SSS-sea surface salinity, PSU.

Differences in mesozooplankton FD and their responses to environmental factors were considerable between coastal (<15 m) and open waters of the shallow brackish Gulf of Riga.
Coastal areas were lower in mesozooplankton FD but higher in total zooplankter numbers (Figure 4), whereas the mesozooplankton community in the open waters appeared to be more resilient to local disturbances as the amount of significant responses were half the number of those found in coastal areas (Figures 5 and 6).
Several studies conducted in brackish habitats have reported a decrease in diversity along the salinity gradient, linking increased plankton diversity to an intermediate salinity zone (5)(6)(7)(8) [14,79,80], i.e., the salinity of the Gulf of Riga open waters. The intermediate disturbance hypothesis [81,82] is in agreement with the observations; it states that species diversity is highest at intermediate disturbance levels due to reduced species densities that weakens interspecific competition allowing less-opportunistic and less-adapted species to survive. However, the results of the present study do not indicate salinity gradient as statistically important for the differences between open and coastal areas in the Gulf of Riga ( Figure 2B, Table 4). Besides the low salinity, the coastal waters of the Gulf of Riga are highly impacted by wind and can be defined as a profoundly fluctuating habitat, conversely to open waters. The high freshwater impact also reinforces mixing and turbulence in the area. Riverine discharge is spread along the eastern and western coastal areas almost equally during summer [83]. Most likely, these are the features responsible for the identified similarities between the studied coastal sites and their differences to open waters ( Figure 2).
Environmental filtering is particularly pronounced in dynamic and fluctuating waters, yet it has been demonstrated that they are also highly productive habitats with effective food webs and intensive biotic interactions [30]. The majority of identified relationships between coastal mesozooplankton FD and environmental drivers were responses to biotic parameters. The only exception was functional evenness (FEve) that showed a significant (p = 0.0074) and accurate (nrmse = 1.26) negative relationship to abiotic conditions, namely, SST ( Figure 5G). A decrease in the balance of the filled niche space of the mesozooplankton community (described by FEve) with increasing temperature is a direct manifestation of abiotic filtering under the fluctuation of seasonal forcing. The benefits from warmer or colder conditions differ between mesozooplankton species, consequently creating shifts in species and trait occurrences (e.g., [84,85]).
As previously mentioned, the strength of interspecific competition is also a significant aspect of FD variation. Helenius et al. [14] reported lower zooplankton FD based on feeding traits at sites where Keratella rotifers or calanoid copepods were prevailing. In the present study, coastal areas were dominated by herbivorous filter-feeders (G1 functional group), including Keratella species, and cruising small-sized omnivores (G4 functional group) that comprise Synchaeta rotifers and Copepoda nauplii ( Figure 3A), thus mirroring the patterns identified in the Gulf of Finland [14].
Both groups G1 and G4 are comparatively small-sized and body size is considered a master trait that defines the main physical abilities and constraints of an organism [86], including its power to retain horizontal and vertical position in the water mass (increasing with size) [87]. Plausibly, the aggregation of non-migrating small-size taxa in the water column is also supported by enhanced jet-like currents that are present in the Gulf of Riga within western and eastern coastal areas during the summer period [83]. Kahru et al. [88] have analysed physical-biological coupling in frontal structures in the Baltic Sea. They found indeed small-sized zooplankters, namely Bosmina and Synchaeta, as a dominating taxa in the fronts and noted the increased zooplankton abundances in the regions, which was explained by particle aggregations due to flow convergence. However, these are still a matter of speculation and only a focused study on the physical processes and their biological implications within the above-mentioned currents of the Gulf of Riga would give more profound information.
The majority of the taxa within groups G1 and G4 have short life cycles (except for Copepoda nauplii and meroplankton that are ephemeral development stages) and the ability of parthenogenetic reproduction (rotifers, cladocerans), which allows for rapid development under suitable conditions and consecutive triumph in competitive encounters [89]. Another factor affecting the abundance of coastal zooplankton can be recruitment from the dormant stages of plankters or resting eggs [90]. However, the sediment type of the coastal area, i.e., sand, does not support accumulation of resting eggs. All sedimenting material are transported to deeper accumulation zones during the autumn convection of the water column [91].
The abundance of group G4 negatively affected functional divergence (SESFDiv) of the mesozooplankton community ( Figure 5I). Functional divergence has been identified as a descriptor for niche differentiation, thus resource availability and competition within a filled trait space [72]. Therefore, the identified negative relationship builds up to the reason that the competition over resources is a significant driver for the mesozooplankton FD in the coastal Gulf of Riga, especially at times when a rapidly-developing functional group dominates the community. All other observed relationships between groups G1, G4, and FD indices showed limited predictive power (nrmse 1; Figure 5C,K,M) or were hardly interpretable (multimodal; Figure 5K) implying non-causal relationships. Although, in one instance, the G4 group had a positive effect on functional dispersion (SESFDis) ( Figure 5M).
Considering the typical summer composition of phytoplankton, prey availability and its quality was expressed as cyanobacterial dominance (cyano:other ratio) in the present study. Cyanobacterial blooms may cause a negative impact on zooplankton diversity due to inhibition [92,93] and low nutritional quality [94,95]. Moreover, feeding on colony-and filament-forming dense aggregations (e.g., Aphanizomenon flosaquae-prevailing species in the Gulf of Riga) is challenging for zooplankton and can result in mechanistic damage to the feeding apparatus [96]. Contrarily, several Copepoda species are ecologically adapted to feed on cyanobacteria [33,97]. Thus, combined effects can result in modified FD of mesozooplankton during cyanobacterial blooms as reported from freshwater ecosystems [98].
However, FD indices of open water communities did not respond to variability of cyanobacterial dominance in the present study ( Figure 6) despite the overall high eutrophication level of the Gulf of Riga. Whereas the cyano:other ratio had a positive effect on functional dispersion (SESFDis) of coastal communities when the ratio was low ( Figure 5J), i.e., non-blooming periods. At the same time, no negative impact during blooms was observed. Still, this relationship had limited predictive power (nrmse = 3.223). The functional dispersion (SESFDis) index describes niche complementarity and is a recommended metric for the detection of assembly processes [5]. Therefore, we argue that cyanobacteria had no pronounced adverse effect on mesozooplankton FD in the Gulf of Riga during the studied period. Indeed, during cyanobacterial dominance, up to 75% of the primary production can be transferred via microbial loop [31,99] meeting the energy demand for mesozooplankton needs.
Still, an important aspect of the interpretation of the present results is the long-term ecosystem view of the region. The studied data (1993-2017) cover a period after the Baltic Sea regime shift [100] and analyse a situation when planktivorous flows dominated the food web of the Gulf of Riga [101]. Beforehand, detritivorous flows were prevailing and they were supported by higher diatom biomass [38,101]. Thus, the present results describe a relatively stable period, from the perspective of the phytoplankton community, limiting the opportunities for identification of phytoplankton-induced impacts on zooplankton FD. The study by Jansson et al. [15] described mesozooplankton FD in relation to abiotic factors covering the period from 1960, including the regime shift. They identified a gradual decrease in the number (richness) of functional groups of summer mesozooplankton in the Gulf of Riga until the early 1990s, but it was followed by a rapid increase (a similar pattern was also observed in the present study; see Figure 4C). What is noteworthy is that several non-indigenous species invaded the Baltic Sea (e.g., cladocerans Cercopagis pengoi, Evadne anonyx, and polychaeta Marenzelleria viridis) in the 1990s and early 2000s [102] adding to the disturbances but, at the same time, increasing species richness ( Figure 4C) and FD [19]. Cladoceran C. pengoi, detected in the Gulf of Riga in 1992, has an evident impact on the pelagic food web [103] as it competes with planktivorous fish for larger mesozooplankton prey (groups G3 and G5; Table 3) and graze upon small-sized zooplankters (groups G1 and G4; Table 3) [104,105] whilst serving as a food source for herring [103,104].
Top-down control also appeared to be influential on mesozooplankton FD in the Gulf of Riga. Both total herring biomass and spawning stock biomass showed a negative impact on functional richness (SESFRic) and evenness (FEve) in coastal areas ( Figure 5A,B,E,F) and spawning biomass almost equally affected FEve in open waters. This pattern could be explained by the selective feeding of the herring, which tend to target larger zooplankters (G5 functional group; Table 3), especially calanoid copepods Eurytemora affinis and Limnocalanus macrurus [26]. The observed positive relationship between abundances of the G5 functional group and FD of coastal and open water mesozooplankton ( Figures 5H,N and 6G) adds to the credibility of this conclusion.
Since coastal areas serve as primary feeding grounds for larvae and young herring, we anticipated a significant impact on mesozooplankton FD in relation to herring recruitment in these areas [106]. Young herring in the Gulf of Riga prey upon various zooplankters, including cladocerans and copepods, but mainly prefer E. affinis [107]. Thus, it is plausible that larvae and young herring are not food-limited in coastal areas, which could explain the lack of a significant impact on coastal mesozooplankton FD. Alternatively, they may have already left for deeper waters by the time sampling occurred, as suggested by observed positive relationship between herring recruitment and functional evenness (FEve) of mesozooplankton communities in open waters ( Figure 6B).
The opposing impacts of spawning stock biomass or total herring biomass and herring recruitment on mesozooplankton functional evenness is unexpected. Yet, it can be explained by the broader scope of prey available for adult herring compared to younger fish resulting in more targeted feeding on the G5 mesozooplankton functional group [26,108]. In other words, adult herring will prey upon larger copepods and cladocerans only and will switch to other prey such as mysids or amphipods if the preferred zooplankton are in sub-optimal densities [108]. Conversely, young herring tend to be less capable of switching prey [108]. Additionally, the visibility is reduced in coastal areas due to the influence of opaque freshwaters reducing the abilities of visual predators [46]. Consequently, fish larvae and young fish are, most likely, forced to feed on what is available even if it is not their preferred prey, affecting the mesozooplankton community in a more balanced way than adult herring.
In conclusion, we can argue that biotic factors are important drivers shaping mesozooplankton FD in the shallow brackish Gulf of Riga during summer. Competition and resource availability in combination with hydrodynamic features were seemingly the central structuring mechanisms behind the dynamics of coastal mesozooplankton FD, although predation by adult herring was also identified as a significant driver. Predation by young herring was solely impacting mesozooplankton in open waters, implying either weaker predation on mesozooplankton or non-limiting prey availability in the coastal areas. The absence of the mesozooplankton FD response to cyanobacterial dominance was slightly surprising and our understanding of this relationship would benefit from a more comprehensive study of bottom-up effects on mesozooplankton FD in the Gulf of Riga.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/w13141881/s1, Data S1: Annual means (1993-2017) of environmental variables describing coastal (western, southern, and eastern areas) and open (stations O1-O4) habitats of the Gulf of Riga. The data describe mesozooplankton communities (number of species, functional diversity indices, total abundance and abundance of functional groups), abiotic conditions (sea surface temperature and salinity) and cyanobacteria-to-other phytoplankton taxa biomass ratio.  Data Availability Statement: The data (except for herring population parameters) presented in this study are available in the Supplementary Material Data S1.