**Species Richness and Taxonomic Distinctness of Zooplankton in Ponds and Small Lakes from Albania and North Macedonia: The Role of Bioclimatic Factors**

#### **Giorgio Mancinelli 1,2,3, Sotir Mali <sup>4</sup> and Genuario Belmonte 1,5,\***


Received: 13 October 2019; Accepted: 11 November 2019; Published: 14 November 2019

**Abstract:** Resolving the contribution to biodiversity patterns of regional-scale environmental drivers is, to date, essential in the implementation of effective conservation strategies. Here, we assessed the species richness S and taxonomic distinctness Δ+ (used a proxy of phylogenetic diversity) of crustacean zooplankton assemblages from 40 ponds and small lakes located in Albania and North Macedonia and tested whether they could be predicted by waterbodies' landscape characteristics (area, perimeter, and altitude), together with local bioclimatic conditions that were derived from Wordclim and MODIS databases. The results showed that a minimum adequate model, including the positive effects of non-arboreal vegetation cover and temperature seasonality, together with the negative influence of the mean temperature of the wettest quarter, effectively predicted assemblages' variation in species richness. In contrast, taxonomic distinctness did not predictably respond to landscape or bioclimatic factors. Noticeably, waterbodies' area showed a generally low prediction power for both S and Δ+. Additionally, an in-depth analysis of assemblages' species composition indicated the occurrence of two distinct groups of waterbodies characterized by different species and different precipitation and temperature regimes. Our findings indicated that the classical species-area relationship hypothesis is inadequate in explaining the diversity of crustacean zooplankton assemblages characterizing the waterbodies under analysis. In contrast, local bioclimatic factors might affect the species richness and composition, but not their phylogenetic diversity, the latter likely to be influenced by long-term adaptation mechanisms.

**Keywords:** crustacean zooplankton; species richness; phylogenetic diversity; bioclimate; freshwater ponds

#### **1. Introduction**

Lentic freshwaters are acknowledged to play a crucial role in regulating the global ecosystem functions e.g., carbon cycle [1] and they are among the Earth's most threatened habitats in terms of intensity of anthropogenic pressures, biodiversity loss, and non-indigenous species introduction [2–4].

They include an extreme variety of habitats differing in ecological characteristics and fragility [5,6]. Surface area represents one of the most apparent differentiating properties: indeed, lentic environments

(304 million water bodies; 4.2 million km<sup>2</sup> in total area [7]) include the lake Superior (82,000 km2), together with small ponds, i.e., waterbodies less than 0.05 km<sup>2</sup> in area [8].

Ponds and small lakes (hereafter PSL) have significant ecological functions [9,10]: among others, they provide a considerable contribution to inland water CO2 and CH4 emissions [11]. In addition, they are, to date, recognized as important biodiversity hotspots, especially in mountainous regions, supporting a high species richness and contributing a high degree of rare species to regional pools ([12–15]; see also [16] for an example on planktonic Calanoida). Noticeably, PSL are threatened by a number of anthropogenic pressures, including nutrient loading, contamination, acid rain, and invasion of exotic species [17]. In addition, infilling (both natural and caused by direct habitat destruction), land drainage, decline in many of their traditional uses, and changes of function determine at a regional scale the drastic reduction in PSL number and connectivity [12].

In the last decade, several investigations have focused on the diversity of benthic invertebrates, as they are excellent bio-indicators of PSL ecological integrity [18,19]. Local factors that are related with e.g., hydroperiod, environmental harshness, water chemistry, spatial connectivity, habitat heterogeneity, and presence of predators, have been recognized to influence the biodiversity of macroinvertebrate assemblages ([20] and literature cited). At a regional scale, attention has been primarily given to the influence of waterbodies area [21–23], while assuming, within the general theoretical background provided by the species-area relationship (SAR) hypothesis [24], that basins' size correlates with the number of microhabitats within the basin itself and with populations' abundance, and thence inversely correlated with the likelihood of random extinctions. However, resolving the contribution to biodiversity patterns of environmental factors acting at a regional scale is, to date, essential to the implementation of effective conservation strategies in the face of e.g. deforestation and climate change ([25,26] and literature cited). Accordingly, several attempts have been made to model biodiversity of freshwater environments by means of regional bioclimatic factors [27,28].

In the present study, we focused on the diversity of crustacean zooplankton assemblages in 40 ponds and small lakes differing remarkably in terms of origin, extension, and altitude from a relatively wide region comprising part of Albania and North Macedonia. A recent faunal inventory focusing on ponds and lakes in the area [29] provided the starting reference information on the taxonomic characteristics of the assemblages.

A number of studies have generally indicated a positive relationship between the surface area of lacustrine environments and zooplankton diversity (e.g., [30–32]; but see [13]); accordingly, crustacean zooplankton has been shown to have higher species richness in small ponds as compared with lakes [16,33,34]. This notwithstanding, we hypothesized that area alone may not be an adequate predictor, and that local bioclimatic conditions may ultimately contribute in explaining diversity variations across waterbodies by influencing their physical-chemical characteristics, as observed in recent investigations on freshwater macroinvertebrates and macrophytes [27,28,35]. This could be particularly true for waterbodies in mountainous habitats, where temperature and precipitation regimes intensely reflect the chemical-physical characteristics and hydroperiod of the waterbodies themselves [36], regulating the harshness and stability of the aquatic environments and, in turn, the diversity of the biota living in them ([22] and literature cited).

To verify the hypothesis and test whether bioclimatic factors can predict assemblages' diversity, we identified a minimum adequate model (MAM) predicting assemblages' diversity across the different waterbodies by means of a heuristic multiple regression approach and Bayesian Information Criterion model selection method while using satellite-derived bioclimatic variables as predictors. Multiple indices are, to date, available to quantify different aspects of biodiversity [37]. Here, we identified predictive MAMs estimating the diversity of crustacean zooplankton assemblages in terms of species richness and average taxonomic distinctness.

Species richness is the most classical measure of biodiversity across ecosystems that has been extensively used in studies on lentic habitats (see references cited above). This index provides an incomplete understanding of biological variability, because it neglects information on the identity and taxonomic relationship among species, and it is hampered by a number of critical limitations [38,39]. Accordingly, we used the average taxonomic distinctness Δ+ [40] to compare the taxonomic relatedness of species in the crustacean assemblages of every water body. In addition, we tested the influence of bioclimatic factors on crustacean assemblages in terms of species composition. To this end, multivariate approaches that are based on a canonical analysis of principal coordinates were used to model the changes in the structure of the assemblages as affected by bioclimatic variables, and identify relationships between the latter and specific groups of zooplankton taxa.

#### **2. Material and Methods**

#### *2.1. Sampling Sites and Collection Procedures*

A total of 40 sites were selected among those (53) that were surveyed between 2005 and 2017 by Belmonte and colleagues [29] in an area comprised between 39◦55 22–42◦04 30 N, and 19◦24 30–20◦47 36 E. (Figure 1, Table 1).

**Figure 1.** Location of the 40 waterbodies in Albania and North Macedonia.



#### *Water* **2019** , *11*, 2384



Sampling procedures are described in detail elsewhere [29]. As the water bodies varied in terms of area, altitude, origin (natural, artificial), as well as in hydrology (permanent, temporary), banks morphology, and degree of aquatic vegetation cover, it was not possible to apply a standardized sampling protocol across all of the sites.

In the selected sites, sampling operations were carried out in spring–early summer (i.e., between April and July) always by the same operators, while using a hand-held plankton net (200 μm mesh-size, mouth diameter, 30 cm). The collection (by plankton net towing from two opposite edges of the pond) covered the whole water body when its diameter was smaller than 100 m. For larger water bodies, sample collection was carried while using a canoe. In the case of small ponds, the collection procedure was repeated three times (each sample derived from the execution of three collections). In the case of larger water bodies, a sample collection was carried out in three different stations; the three different samples were ultimately cumulated.

After collection, the samples were fixed in situ in 90–96% ethanol. In the laboratory, taxa were identified to the species level while using a compound microscope (×30–×300 magnifications) that was equipped with a *camera lucida*.

The quantification of the abundance of each taxon was not performed, and only presence/absence data were considered due to the huge variability of water volumes in each pond, which made impossible the comparison among the concentrations of plankton of different sites.

#### *2.2. Landscape-Climate Variables*

Together with altitude, we used area and perimeter of the water bodies together with bioclimatic factors (i.e., temperature and precipitation) to represent the landscape-climate variables. Water bodies were geo-referenced in Google Earth Pro version 7.3.2., where their surface (in km2) and perimeter (in km) were measured while using the software tools. Measurements were performed by preferentially choosing images that were taken in spring or summer between 2008 and 2017, assuming that negligible variations in the water bodies size and morphology occurred during this period.

Nineteen climatic layers with a 30-second spatial resolution (0.93 <sup>×</sup> 0.93 = 0.86 km<sup>2</sup> at the equator; approximately 0.92 <sup>×</sup> 0.70 = 0.64 km2 within the study area), including temperature and precipitation variables, were extracted from the WorldClim v2 data set [41] (Table A1 in online material). Besides climate, there is a growing recognition of the importance of vegetation cover in characterizing the spatial environmental heterogeneity in a given area at the meso- and topo-scales, in turn affecting climate, soil composition, hydrology and geomorphology, and, ultimately, biological processes that were related with species richness and community complexity [42,43]. Accordingly, two vegetation variables (i.e., percent tree cover and percent non-tree cover) were extracted from the Terra MODIS Vegetation Continuous Field (VCF) product (available as MOD44B v006 https://lpdaac.usgs.gov/products/mod44bv006/). Percent tree cover included all forest types and age classes, while percent non-tree cover included meadows, regeneration areas, and clear-cut areas. MODIS tiles of the study area were re-projected and re-sampled to meet the coordinate system and resolution of WorldClim layers; percent cover data were subsequently obtained for the years from 2006 to 2017, and averaged. In addition, the third VCF component of ground cover, i.e., percent bare soil (including bare soils and rocks) was extracted according to the aforementioned procedures, and was used together with tree and non-tree vegetation percent cover data to estimate the Shannon's diversity index (H) as a proxy of habitat heterogeneity.

#### *2.3. Data Analysis*

The values in the text are expressed as averages ± 1SE; for parametric statistical analysis, data were tested for conformity to assumptions of variance homogeneity (Cochran's C test) and normality (Shapiro–Wilks test) and transformed when required.

The taxonomic diversity of crustacean assemblages in each water body was estimated in terms of species richness S and average taxonomic distinctness Δ+. The index can be used as a proxy for phylogenetic diversity and it measures the mean path length through a taxonomic tree connecting every species [40]. Here, mean taxonomic distinctness values were calculated assigning equal weighting to branch lengths from a linear Linnaean classification while using eight taxonomic levels (i.e., class, subclass, order, suborder, infraorder, family, genus, and species). The taxonomic classification tree was built according with the World Register of Marine Species (WoRMS, available at https://www.marinespecies.org) and the Integrated Taxonomic Information System (ITIS, available at https://www.itis.gov).

For the sake of completeness, other taxonomic diversity indices were calculated, including the total taxonomic distinctness sΔ+ and the variance in taxonomic distinctness Λ+ [44,45], the average phylogenetic diversity Φ+, and the total phylogenetic diversity sΦ+[46]. S resulted in being significantly related with sΔ+, Φ+, and sΦ+ (*r* = 0.98, −0.89, and 0.95, respectively; *P* always < 0.01, 38 degrees of freedom), while Δ+ scaled negatively with Λ+ (*r* = −0.44, *P* = 0.004); conversely, S and Δ+ were characterized by a non significant negative correlation (*r* = −0.26, *P* = 0.12, 38 d.f.); thus, the two indices were chosen for further analyses.

We verified the influence on both indices of potential artefacts, due to (i) possible differences in sampling procedures and (ii) the different number of total collected samples per water body. To this end, we performed a one way permutational analysis of variance (PERMANOVA; [47]) based on Euclidean distances and 999 permutations with "sampling procedure" as a fixed factor (two levels, "hand", or "canoe") and the number of collected samples as the covariate (P and N hereafter). Both factors exerted negligible influences on the S and Δ+ estimations (S: Pseudo–FP = 4.42, P(perm)P = 0.08, Pseudo–FN = 0.48, P(perm)N = 0.47, Pseudo-FP×<sup>N</sup> = 0.12, P(perm) <sup>P</sup>×<sup>N</sup> = 0.72; Δ+: Pseudo–FP = 1.67, P(perm)P = 0.21, Pseudo–FN = 0.01, P(perm)N = 0.94, Pseudo–FP×<sup>N</sup> = 1.07, P(perm) <sup>P</sup>×<sup>N</sup> = 0.31). Consequently, the S and Δ+ values were assumed to provide robust estimation of planktonic Crustacea diversity across the studied water bodies. PERMANOVA was further used to test the effects of the factors "origin" (two levels, "natural" and "artificial") and "typology" (two levels, "pond" and "lake") on the diversity indices. As water bodies varied greatly in elevation (Table 1), the latter was included in the analyses after log-transformation as a continuous covariate.

The georeferenced locations of the sampling sites were used to extract climatic and vegetation data from environmental layers. The final data set included 19 climatic, two vegetation (% tree cover, % non-tree cover), four geomorphological (elevation, perimeter, surface, and perimeter/surface ratio), and one habitat heterogeneity variable (Table A1). They were log-transformed and z-scaled; subsequently, their original number (25) was reduced while using an iterative variance inflation factor (VIF) analysis [48,49]. In brief, if a strong linear relationship links a variable *x* with at least another variable *y*, the correlation coefficient would be close to 1, and the VIF for *x* would be large. Here, diversity measures with VIF factors that were larger than 10 were excluded. Variables with VIF factors larger than 10 were discarded. The identification of a minimum adequate model (MAM hereafter; [50]) linking diversity measures with environmental variables was based on the heuristic generation of alternative regression models (see [51] for complete details on the procedure). Model selection was performed while adopting an Information Theoretic criterion [52]; the second-order Akaike Information Criterion *AICc* [53,54] was calculated for each combination of *n* explanatory variables and used to identify the best MAM among the alternative regression models that were generated by the procedure. For model comparison, *AICc* values were used to estimate a set of positive Akaike weights *wi* summing 1:

$$w\_i(AIC) = \frac{\exp\left[-1/2(AIC\_i - minAIC)\right]}{\sum\_{1}^{K} \exp\left[-1/2(AIC\_i - minAIC)\right]} \tag{1}$$

With *K* = number of models. The model showing the highest *wi* was accepted as the best candidate; other candidate models were accepted if characterized by *wi* values within 12.5% of the highest [55–57]. The model building and MAM identification procedures were performed while adopting Fox and Weisberg [58] as a general reference.

Given the non-conclusive outcomes of the analyses that were performed on taxonomic distinctness (see Results section), we verified whether bioclimatic factors influenced planktonic Crustacea assemblages in terms of species composition. Species incidence data were used to calculate a Jaccard similarity matrix across the different waterbodies. In addition, a similarity matrix that was based on Euclidean distances was constructed for bioclimatic variables, and the consistency of the two matrices was tested while using the Kendall coefficient of concordance (W). Subsequently, a canonical analysis of principal coordinates (CAP) [59] was performed to model changes in assemblages composition, as affected by bioclimatic factors. The appropriate number of principal coordinates *m* was chosen as to minimise the *P* value from the permutation test based upon the trace statistic and maximizing the leave-one-out allocation success [60]. Post-hoc PERMANOVA tests were performed to confirm the results of the ordination for both bioclimatic factors and species; SIMPER analyses were performed on the Euclidean distance matrix of the former to assess the percentage contribution of each factor to the dissimilarity between the groups of waterbodies that were identified by the CAP procedure. Furthermore, the Spearman's rank correlations were estimated to identify the bioclimatic variables and the crustacean species that most effectively described the groups of waterbodies that were identified by the CAP procedure. Only the variables with a Spearman rank correlation coefficient *r* > 0.55 were considered.

All of the analyses were implemented in the R statistical environment v3.6.1 [61] while using a suite of packages including *taxize* (for taxonomic information retrieval from online databases) *vegan* (for diversity measures and multivariate analyses), *raster*, *rgdal*, and *maptools* (for environmental layers manipulation), *usdm* (for VIF analysis), *car*, *leaps*, and *HH* (for MAM identification).

#### **3. Results**

#### *3.1. General Features*

The 40 water bodies analysed, varied remarkably in terms of altitude, area, and perimeter (Figure 2). They showed an average altitude of 1168.3 m a.s.l. (± 141.9 m SE), ranging from 2 m a.s.l. (Narte Zvernec pond, ZVE in Figure 1) to 2,435 m a.s.l. (SHR 5). The average surface extension was 0.09 km2 <sup>±</sup> 0.03 SE, ranging from 9.1×10−<sup>4</sup> to 0.86 km2. The average perimeter was 1.01 <sup>±</sup> 0.21 km, ranging between 0.035 and 6.1 km. For both of the variables, the minimum and maximum values corresponded with a small artificial pond in the karst highlands of Progonat (PRO1) and the Lake Seferan in the Dumre region (SEF).

In the 40 water bodies, 79 Crustacea species were identified in total, being almost equally distributed between the classes Branchiopoda and Hexanauplia (41 and 38 species respectively).

Among Branchiopoda, Anomopoda outnumbered the other two orders Ctenopoda and Haplopoda (38 vs. 2 and 1 species). *Daphnia*, *Ceriodaphnia*, and *Moina* were the genera that were characterized by the highest number of species (nine *Daphnia* species, four *Moina* species, and three *Ceriodaphnia* species) together representing the majority of all the sampled Anomopoda species. Ctenopoda were represented by the congeneric *Diaphanosoma brachyurum* and *D. lacustris* while Haplopoda by the single species *Leptodora kindtii*. The class Hexanauplia (*alias* Copepoda) was dominated by species belonging to the order Cyclopoida (31) and to a minor extent Calanoida (7). The genus *Cyclops* (seven species) together with *Acanthocyclops*, *Paracyclops* (four species each), and *Mesocyclops* (three species) constituted to the majority of the species in the order; Calanoida were represented by the genera *Eudiaptomus* (three species), *Mixodiaptomus* (two species each), and by *Arctodiaptomus salinus* and *Neodiaptomus schmackeri*.

**Figure 2.** Frequency distribution of waterbodies area, perimeter, and altitude. The histograms for first variables adopt geometric scale increments (×2), while for altitude (elevation) a linear scale is used.

The species showing the widest distributions were the Anomopoda *Bosmina longirostris* (16 sites), *Chydorus sphaericus* (15 sites), and *Daphnia longispina* (14 sites) (Table 2); in addition, the Cyclopoida *Mesocyclops leuckarti*, the Calanoida *Mixodiaptomus tatricus*, and the Ctenopoda *Diaphanosoma brachyurum* occurred in 12 sampled sites.

**Table 2.** Summary of PERMANOVAs on species richness S and taxonomic distinctness + of crustacean zooplankton assemblages testing for the effects of waterbodies' origin and typology including elevation as a continuous covariate \*\*: *p* < 0.01.


As regarding high level taxa, only Anomopoda (Cladocera) were present in every site. Cyclopoida (Hexanauplia) were present in 38 sites (95% of the total), Calanoida (Hexanauplia) in 30 sites (75%), Ctenopoda (Cladocera) in four sites (10%), and Haplopoda (Cladocera) in three sites (7.5%).

#### *3.2. Diversity Patterns and Bioclimatic Correlates*

On average, 6.7 ± 0.4 species per water body were found, ranging between three (JAH, PIL, RRE) and 14 species (KOR 1). The taxonomic distinctness Δ+ of the different planktonic assemblages was on average 72.1 ± 0.78, ranging between 60 (SHB 2) and 88.9 (GRA).

The factor "origin" was the only exerting significant effects of the species richness of waterbodies (Table 2), with the six artificial water bodies being included in the study characterized by lower S values as compared with natural basins (3.5 ± 0.22 vs. 7.29 ± 0.38, respectively). Conversely, negligible effects were generally observed for the taxonomic distinctness Δ+ (Table 2).

The 25 predictor variables (Table A1) were reduced to a set of nine characterized by negligible collinearity (TableA2). They included five climatic variables (i.e., Isothermality, Temperature Seasonality, Mean Temperature of Wettest Quarter, Annual Precipitation, and Precipitation of Coldest Quarter), % tree and % non-tree vegetation cover, habitat heterogeneity, and water body surface. Besides water body surface, the variable Mean Temperature of Wettest Quarter showed the greatest among-water bodies variability (Table A2), ranging from a minimum of −5.32 ◦C (SHB 2) to a maximum of 11 ◦C (ZVE). It was followed by % tree cover (varying between 1 and 70%, BEL and DEG, respectively) and % non-tree cover (ranging between 20 and approx. 81%, BEL-DEG and SHR 2, respectively).

The heuristic search procedure identified a Minimum Adequate Model (MAM) predicting the variation of species richness S across the different water bodies relying on the three explanatory variables % non tree cover, Temperature Seasonality, and Mean Temperature of Wettest Quarter (Figure 3; multiple *r* = 0.58, *P* = 0.002, d.f. = 3, 36). The MAM was characterized by the lowest *AICc* value, and by an Akaike weight *wi* approximately eight times larger than the second-best candidate, based on the variables Temperature Seasonality, Mean Temperature of Wettest Quarter, and Habitat heterogeneity (Table 3). All of the predictors provided significant contributions to S variation across water bodies (minimum absolute *t* value = 2.34, *P* = 0.02 for the variable Mean Temperature of Wettest Quarter); the contributions of both % non-tree cover and Temperature Seasonality were positive (b = 0.27 ± 0.12 and 0.64 ± 0.15, respectively), while the Mean Temperature of the Wettest Quarter provided a negative contribution (b = −0.26 ± 0.14). Noticeably, none of the first ten best-performing models included water body area (Table 3).

**Figure 3.** Species richness S of crustacean zooplankton assemblages in the 40 waterbodies under analysis plotted against the values predicted by the best Minimum Adequate Model (MAM; see Table 3 for predictor variables) identified adopting a multiple regression procedure (top). Dashed curves are 95% CI. For the sake of completeness the predicted taxonomic distinctness values Δ+ by the best MAM are also reported (bottom).

**Table 3.** Summary of the heuristic multiple regression analysis followed by a parsimonious selection procedure of the Minimum Adequate Model (MAM) predicting species richness (S) and of crustacean zooplankton assemblages by means of bioclimatic variables; only the first 10 best MAMs are reported. For the sake of completeness, results of MAM analysis are reported also for taxonomic distinctness (Δ+), even though the statistical power of the models was negligible (see Results). K: number of predictors included in the model; AICc: second-order Akaike Information Criterion; *wi*: Akaike weight. For predictor abbreviations see Table A1.


In contrast with species richness, the heuristic search procedure was unable to identify a MAM with a significant predictive power for taxonomic distinctness. The single variable % tree cover, resulted the best predictor of Δ+ (Table 3); however, the correlation resulted in being non-significant (*r* = 0.25, *P* = 0.11, d.f. 1,38; see also Figure 3). Other models showed an even worst performance, independently from the number of variables involved (Table 3), indicating, in turn, that the taxonomic distinctness of the crustacean assemblages cannot be predicted by bioclimatic, landscape-scale factors.

Canonical analysis of principal coordinates (CAP), followed by a confirmatory PERMANOVA test identified two main groups of waterbodies significantly different in terms of species composition (Figure 4; Pseudo–F = 7.2, P(perm) = 0.001). The variables Isothermality, Mean Temperature of Wettest Quarter, Annual Precipitation, and Precipitation of Coldest Quarter showed a correlation (*r* > 0.65) with the canonical axis 1 (Figure 4) and significantly differed between the two groups of waterbodies (PERMANOVA, Pseudo–F = 27.4, P(perm) = 0.001). A Simper procedure indicated that the variable Mean Temperature of Wettest Quarter contributed by 32.2% to inter-group differences, followed by Isothermality and Annual Precipitation (28.6 and 24.2%, respectively).

**Figure 4.** Canonical analysis of principal components (CAP) analysis (*m* = 6, misclassification = 13%, *P* = 0.0001) testing the differences in crustacean zooplankton assemblages across the 40 waterbodies included in the study as affected by bioclimatic factors. *Markers* represent waterbodies labelled with the respective code (Table 1); their color categorizes the two groups (i.e., red for group1 and blue for group2) showing significant differences in species composition (post hoc PERMANOVA, P(perm) < 0.001). Vector overlay are Spearman correlations of bioclimatic factors and species with canonical axes with *r* > 0.55.

In group1, the Mean Temperature of Wettest Quarter was remarkably lower than that determined for group2 (0.06 ± 0.69 vs. 8.44 ± 0.09 ◦C; *t*-test for separate variances: *t* = 9.41, *P* < 0.0001, 25.53 d.f.). Similarly, Isothermality showed lower values in group1 (33.91 ± 0.23 vs. 38.42 ± 0.24, *t* = 9.19, *P* < 0.0001, 29.52 d.f.), while the Annual Precipitation showed an inverse pattern (1098 ± 10.42 vs. 1012 ± 1.71 mm, *t* = −6.43, *P* < 0.001, 26.09 d.f.).

The analysis of the relationships between species occurrences and the canonical axis 1 (see Table A2 for Spearman correlations for the complete list of species) indicated that *Eucyclops serrulatus*, *Chydorus sphaericus*, *Mixodiaptomus tatricus*, and *Daphnia rosea* in the first group were correlated mainly with Precipitation. The occurrences of *Neodiaptomus schmackeri*, *Diapahnosoma brachyurum*, *Cyclops scutifer*, *Ergasilus* sp, *Bosmina longirostris*, and, to a lesser extent, *Mesocyclops leuckarti* in the group2 were

related with the variables Isothermality, Mean Temperature of Wettest Quarter, and Precipitation of Coldest Quarter.

#### **4. Discussion**

The earth is undergoing an accelerated rate of native ecosystem conversion and degradation and there is increased interest in measuring and modelling biodiversity while using landscape-scale, remotely-sensed predictors. Here, we made an attempt towards this direction while using crustacean zooplankton. This group of organisms, ubiquitous in lentic habitats, has been recently subjected to renewed interest as an effective bio-indicator of the environmental status of ponds and small lakes [62–65]. The heuristic procedure was used here to identify minimum adequate models predicting the diversity the crustacean zooplankton assemblages across the 40 waterbodies under analysis provided non-univocal results. On one hand, they showed that a subset of bioclimatic variables could effectively predict the variation in species richness and composition across the different waterbodies. On the other hand, they also indicated that the assemblages' taxonomic distinctness Δ+ is unrelated with landscape-scale environmental drivers.

Before discussing these results, it must be considered that the emphasis we put on landscape and bioclimatic drivers of zooplankton diversity by no means imply that other chemical, physical, and biotic characteristics of the waterbodies, such as nutrient concentration, pH, depth, predators, aquatic vegetation, etc. are to be considered of secondary importance. A number of studies have unequivocally indicated that they can directly affect zooplankton species richness ([66,67] and references cited in the introduction; but see also further in this section). In addition, lake age [68,69], connectivity, and, in general, the spatial arrangement of the habitat have been acknowledged to play an important structuring role (e.g., [70]; see also [71] for a marine example). However, in the present study, a hierarchy of effects at different spatial and environmental scales is implicitly assumed, with landscape and bioclimatic drivers indirectly affecting the characteristics of the biota (including crustacean zooplankton) by affecting the chemical and physical conditions of the waterbodies. Indeed, the limited extension of the basins that were included in our study (Figure 1) actually implies for them a low thermal inertia, and thus the ability to rapidly respond to the external climatic conditions [72]. An indirect support to this view is also provided by the negligible predictive power of waterbodies' area for assemblages' S, Δ+, and species composition, confirming the results of investigations performed on the benthic fauna of high-altitude ponds [22].

The best MAM included as predictors the degree of non-arboreal vegetation cover of the land areas neighboring the waterbodies, temperature seasonality, and mean temperature of the wettest quarter. The positive influence of the non-arboreal vegetation cover on species richness is consistent with the results of studies that were focused on macrobenthos in lotic habitats (e.g., [73] and literature cited). This could be ascribed to a positive, indirect effect of lateral trophic enrichment on aquatic primary producers, increasing zooplankton diversity by a phytoplankton-mediated bottom up effect or, alternatively, promoting habitat heterogeneity through an increase in aquatic vegetation [66,74,75].

Noticeably, the two temperature variables had contrasting effects on species richness: the lowest number of species was predicted to occur in waterbodies that were subjected to minimum temperature variability during the year and to maximum temperatures during the wettest months, i.e., in winter. The positive influence of temperature variability on patterns of species diversity has been acknowledged for zooplankton and other freshwater invertebrates [76,77], and can generally be ascribed to an effect of habitats environmental heterogeneity on empty niches availability, and, in turn, on species richness [78]. The negative influence of the mean temperature of the wettest quarter (MeanTwet) on species richness can be explained while considering that the former scales negatively with the altitude of the water bodies (*r* = −0.95, *P* < 0.001, 38 d.f.). Thus, MeanTwet maximum values were observed for basins that were located at low altitudes, generally in highly anthropized areas. Accordingly, the lower species richness characterizing these environments might be actually determined by the interplay of a spectrum of anthropogenic perturbations, such as pollution, water caption, or introduction of

predatory fish. Indeed, until 1990, the vast majority of low-altitudes waterbodies in Albania have been stocked with both native and non-indigenous fish species [79], and fish predation has been repeatedly recognized to influence the species richness of zooplankton in lentic habitats [80,81].

The "assemble first, predict later" modelling strategy that was used in the present study was successful in predicting species richness and is generally acknowledged to have several advantages, among others an enhanced capacity to synthesize complex data into a form more readily interpretable by scientists and decision-makers [82]. However, it is apparent that it presents important limitations, as testified by the failure in modeling taxonomic distinctness. Δ+ resulted in being not predictable, confirming the results of several investigations that have found weak or negligible relationships of the taxonomic distinctness of macrobenthic communities with environmental factors [27,83,84]. Species richness S and taxonomic distinctness Δ+ are not conceptually (or mechanistically) related, and they behave differently [44]. The lack of congruence between S and Δ+ and the negligible predictability of the latter is because S is likely to respond to short-term environmental changes in the waterbodies under analysis, while Δ+ is a proxy for phylogenetic diversity, reflecting a complex set of intrinsic and extrinsic traits and expressing evolutionary long-term adaptations to local environmental conditions [85]. The canonical analysis of principal coordinates allowed for us to partially overcome these limitations, applying an "assemble and predict together" strategy [82] in order to model changes in the species composition of planktonic assemblages and provide an advanced resolution of species-specific relationships with bioclimatic factors. The CAP analysis (Figure 4) distinguished two distinct group of waterbodies, showing different climatic characteristics in terms of isothermality, mean temperature of the wettest quarter, and annual precipitation. The first (red circles in Figure 4) was mainly constituted by high-altitude ponds and lakes (elevation 1725.3 ± 123.7 m, mean ± 1SE) distributed throughout the study area characterized by Crustacea species (e.g., *Eucyclops serrulatus*, *Mixodiaptomus tatricus*) that are peculiar of pristine alpine environments [86]. The second group (blue circles in Figure 4) comprised low-altitude karst waterbodies that were located in the Dumre area (Figure 1; elevation 239.9 ± 86.2 m, mean ± 1SE), where they are subjected to several anthropogenic pressures, including agricultural and urban pollution, eutrophication, and the introduction of non-indigenous fish species [87]. Accordingly, the group is characterized by the occurrence of *Neodiaptomus schmackeri*, an Australasian species of Chinese origin that was recently introduced in Albanian lentic habitats through fish stocking [88] and the copepod *Ergasilus* sp., whose adult females are ectoparasites of fish [89].

Regarding the three isolated waterbodies in Figure 4 (i.e., JAH, PIL, and RRE), they are artificial reservoirs located at 1330, 701, and 280 m a.s.l., respectively, being heavily affected by cattle frequentation (Belmonte, personal observation). Their isolation in the CAP diagram is due to the low species richness (three species in all the waterbodies), that might be ascribed to cattle-induced eutrophication conditions [90]. However, it is worth noting that copepods vary their body size (at the community level and even for single species) inversely with the eutrophication level [91,92]. Thus, by using a plankton net with a mesh size of 200 μm, we may have underestimated the smaller component of the planktonic assemblages thus biasing the species count.

A final consideration deserves a brief mention. In this study, the spatial resolution of the bioclimatic layers (approximately 0.64 km2) was lower than the area characterizing most of the ponds and lakes included in the analysis (Figure 2). In other words, the bioclimate spatial grid "matched" the dimensions of the waterbodies, the latter being completely included (and described in terms of climate and vegetation cover) within the same grid cell. Additional studies including a wider size range of waterbodies as well as bioclimatic layers resolved at different spatial resolutions are necessary to provide a more complete picture of the actual relationships linking bioclimatic factors and the diversity of lentic zooplankton at multiple regional and environmental scales.

**Author Contributions:** Conceptualization, G.B. and G.M.; methodology, G.B., S.M. and G.M.; software, G.M.; validation, G.B.; formal analysis, G.B. and G.M.; investigation, G.B. and S.M.; resources, G.B. and S.M.; data curation, G.M.; writing—original draft preparation, G.B. S.M. and G.M.; visualization, G.M.; supervision, G.B.

**Funding:** This research received no external funding.

**Acknowledgments:** Many collaborators allowed, during 12 years, the collection of samples. Among them we are particularly grateful to Bilal Shkurtaj (Univ. of Vlore, Albania), and Giuseppe Alfonso (Univ. of Salento, Italy). The authors finally thank two anonymous reviewers whose comments greatly improved the original manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Table A1.** List of the 25 bioclimatic variables used in the study. Layers coded MeanT- PrecColdQ were obtained from the Wordclim v2 dataset [41] available at http://www.worldclim.org/, and refer to average monthly climate data relative to the period 1970-2000 with a 30 arc-second spatial resolution (approx. 0.92×0.70 km within the study area). The vegetation layers %tr and %nontr (% tree cover and % non-tree cover) were extracted from the MODIS Vegetation Continuous Field (VCF) product (https://lpdaac.usgs.gov/products/mod44bv006/; [93]). In the table, a third VCF variable—% bare soil (%bare, in italics)—is included, as it was used together with variables %tr and %nontr only for the estimation of habitat heterogeneity (see text for further details). The R packages *raster*, *rgdal*, and *maptools* were used for the manipulation of bioclimatic layers [94–96].



**Table A2.** List of the 79 crustacean species sampled in the 40 waterbodies under analysis. Linnaean classification of species using eight taxonomic levels (i.e., class, subclass, order, suborder, infraorder, family, genus and species) and the total number of occurrences are included.

**Table A2.** *Cont.*







**Table A4.** *Cont.*

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Dynamics of Mesozooplankton Assemblage in Relation to Environmental Factors in the Maryland Coastal Bays**

#### **Efeturi U. Oghenekaro and Paulinus Chigbu \***

NSF CREST-Center for the Integrated Study of Coastal Ecosystem Processes and Dynamics in the Mid-Atlantic Region, and NOAA Living Marine Resources Cooperative Science Center, Department of Natural Sciences, University of Maryland Eastern Shore, Princess Anne, MD 21853, USA; euoghenekaro@umes.edu **\*** Correspondence: pchigbu@umes.edu; Tel.: +14-10-621-3034

Received: 16 August 2019; Accepted: 4 October 2019; Published: 14 October 2019

**Abstract:** The mesozooplankton composition and dynamics in coastal lagoons of Maryland, mid-Atlantic region, USA have received little scientific attention despite the fact that the lagoons have undergone changes in water quality in the past two decades. We compared mesozooplankton abundance and community structure among sites and seasons, and between 2012, a year of higher than average salinity (33.4), and 2013 with lower than average salinity (26.6). It was observed that the composition, diversity, and abundance of mesozooplankton in 2012 differed from those of 2013. Barnacle nauplii were abundant in 2012 contributing 31% of the non-copepod mesozooplankton abundance, whereas hydromedusae were more dominant in 2013 and contributed up to 83% of non-copepod zooplankton abundance. Gastropod veliger larvae were more abundant in 2013 than in 2012 while larvae of bivalves, polychaetes, and decapods, in addition to cladocerans and ostracods had higher abundances in 2012. The abundance and diversity of mesozooplankton were explained by variations in environmental factors particularly salinity, and by the abundance of predators such as bay anchovy (*Anchoa mitchelli*). Diversity was higher in spring and summer 2012 (dry year) than in 2013 (wet year). The reduction of salinity in fall 2012, due to high freshwater discharge associated with Hurricane Sandy, was accompanied by a decrease in mesozooplankton diversity. Spatially, diversity was higher at sites with high salinity near the Ocean City Inlet than at sites near the mouth of tributaries with lower salinity, higher nutrient levels and higher phytoplankton biomass. Perhaps, the relatively low salinity and high temperature in 2013 resulted in an increase in the abundance of hydromedusae, which through predation contributed to the reduction in the abundance of bivalve larvae and other taxa.

**Keywords:** mesozooplankton; salinity; abundance; distribution; diversity; Maryland Coastal Bays

#### **1. Introduction**

Planktonic organisms ranging in size 0.2–20 mm are referred to as mesozooplankton [1]; examples are copepods, decapod zoea, rotifers, gastropod veliger, barnacle nauplii, larvaceans, cladocerans and bivalve larvae. Mesozooplankton are major grazers of primary producers, and important predators of microzooplankton, including ciliates, thereby acting as a pathway for transfer of materials and energy from microbial food web to higher trophic levels [2,3]. Copepods are the most abundant of, and contribute the most biomass to, the mesozooplankton group in marine and estuarine zooplankton [4].

Although non-copepods do not constitute the bulk of mesozooplankton in United States mid-Atlantic estuaries, they occasionally become abundant [5–7] and play an important role in the diet of zooplanktivores [8,9]. Mesozooplankton composition, abundance and diversity in estuaries and coastal lagoons vary seasonally, spatially and between years in response to abiotic and biotic factors [10–16], and have been well documented in large river-dominated estuaries of the mid-Atlantic such as the Chesapeake Bay [12,17] and Delaware Bay [18,19]. Nevertheless, we know little about mesozooplankton composition and dynamics in Maryland Coastal Bays (MCBs), which are poorly flushed, eutrophic, largely polyhaline lagoonal systems that receive relatively small amount of freshwater directly from tributary rivers and creeks. Because of differences in the hydrography, the composition, abundance and distribution of mesozooplankton in MCBs may be different than those of river-dominated estuaries, and are expected to vary between years due to variations in climatic factors.

The coastal bays support commercially and recreationally important species of shellfish such as blue crabs (*Callinectes sapidus*) and hard clams (*Mercenaria mercenaria*) which rely on the successful recruitment of their meroplanktonic larvae to sustain the adult populations. The blue crab is the most abundant shellfish in the MCBs and requires high salinity waters for spawning and larval development [20,21]. In the Chesapeake Bay, *C. sapidus* spawn near the mouth of the estuary, or in the nearby coastal ocean [6,22]. The larvae (zoeae) are then transported into adjacent coastal waters on the continental shelf [23–25] where they develop through seven or eight zoeal stages before subsequently returning as post-larvae and juveniles to the estuary [25], while the adults remain resident in the system. Due to the high salinity (>26) levels in a large portion of the MCBs, it has been hypothesized that resident adult populations of the *C. sapidus* spawn in the high saline waters of the coastal bays and the larvae complete their development within the Bays [26]. In the 2014 updated Coastal Bays Blue Crab Fishery Management Plan Implementation document, the need for studies to determine the "level of localized reproduction and entrapment of larvae" of blue crabs in the MCBs was emphasized.

Furthermore, salinity in the MCBs increased after the opening of the Ocean City Inlet in 1933, which favored the development of hard clams (*M. mercenaria*) and, by the 1960s, the fishery superseded the commercial landings and value of oysters in the bay [27]. The population of hard clams is currently dominated by large and older adults; the recruitment is low and intermittent except in areas close to the Ocean City Inlet, and densities are lower than historical levels [27] although increasing. Efforts to restore the population including stoppage of the hydraulic clam dredging activities have not resulted in significant increases in clam density, the cause of which is unknown, although intense predation by species such as blue crab, and the frequent occurrence of brown tide, *Aureococcus anophage*ff*erens*, blooms that negatively impact bivalve larvae have been speculated to be the cause [27]. Information on the distribution and abundance of bivalve larvae, and the environmental factors influencing their abundance would be useful for designing effective management strategies to enhance bivalve populations in the MCBs.

Within the embayments in the MCBs, nutrient levels, phytoplankton biomass and fish abundance vary spatially, seasonally and inter-annually [28–30], which in addition to variations in the abundance of gelatinous zooplankton predators can influence the distribution and abundance of mesozooplankton. For example, mortality of bivalve larvae due to predation by gelatinous zooplankton can potentially slow recovery of hard clams despite the implementation of effective management plans [31,32].

The objectives of the study were:


The study addressed the following questions: Are mesozooplankton abundance and diversity higher at sites closer to the Ocean City Inlet with higher salinity, than at sites farther from the Inlet with lower salinity? Are mesozooplankton abundance and diversity higher in spring and summer 2012, a year with relatively higher salinity and lower temperature, than in 2013 when salinity was lower, and temperature higher? Do blue crabs spawn within the MCBs?

#### **2. Materials and Methods**

#### *2.1. Study Area*

The MCBs located in the Mid-Atlantic region on the east coast of the United States is separated from the Atlantic Ocean by barrier islands. This lagoonal system consists of five embayments (Figure 1) that are well mixed mostly by wave and tidal actions that facilitate the exchange of water and materials between the MCBs and the nearby coastal ocean [29,33]. The system, however, is poorly flushed, has low freshwater input from groundwater, rivers and tributaries, and therefore high salinity that is close to that of the ocean except in upstream areas of rivers and creeks [28,29]. The MCBs are relatively shallow (average depth about 1.2 m) with a maximum depth of 9.3 m at the Ocean City Inlet [29]. Some characteristics of the northern bays watershed include low forest covers, non-point sources of pollution, dead-end canals, high stream nitrate, and nutrients from agricultural activities and development [34]. Increased nutrient input has resulted in eutrophic conditions in the northern bays and tributaries, although a similar condition exists in the Newport Bay located in the southern MCBs.

**Figure 1.** Map of the study area showing sampling locations within the Maryland Coastal Bays.

#### *2.2. Sample Collection*

Water quality parameters such as temperature, salinity, dissolved oxygen (DO) concentration, and pH were measured in situ using YSI 6600 sonde. Zooplankton samples were collected monthly from February 2012 to December 2013 from 13 sites (Figure 1). Sites 1–5 are located in Chincoteague Bay (CB), Site 6 in Newport Bay (NB), Sites 7 and 8 in Sinepuxent Bay (SB), Site 9 in Isle of Wight Bay (IWB), Site 10 at the mouth of St. Martin River (SMR) and Sites 11–13 in Assawoman Bay (AB). Zooplankton sampling was conducted using a plankton net, with 60-cm diameter mouth opening, 240-cm-long conical section, a mesh size of 200 μm, and a cod-end made of PVC with a volume of 1000 mL and diameter of 11 cm. Because of the shallow nature of the lagoons (average depth < 2 m), the net was towed horizontally in the upper meter of the water column similar to how zooplankton samples were collected in the shallow Barnegat Bay, NJ [35]. Sampling was for a duration of 2 min. The mouth opening of the net was equipped with a General Oceanic flow meter (Model 2030R) attached to the frame to measure the volume of water filtered by the net. Zooplankton samples were preserved in 5% buffered formalin solution immediately after collection. Biovolumes of gelatinous zooplankton (*Mnemiopsis leidyi*) were determined with a measuring cylinder.

Preserved plankton samples were processed in the lab using a dissecting microscope. A Folsom splitter was used to split samples when necessary. Samples were diluted to a known volume, and a minimum of 5 mL subsample was extracted with a Stemple pipette and 200–500 zooplankton individuals were counted. Cladocerans and crab zoeae were identified to species level when possible using a compound microscope and by reference to various publications [1,36,37]. Data on bay anchovy (*Anchoa mitchelli*) catch per unit effort (CPUE) were obtained from the Maryland Department of Natural Resources (MD DNR). Water samples for analyses of chlorophyll *a* and 19'-Butanoyloxy-Fucoxanthin (But-fuco), a diagnostic pigment for pelagophytes such as the Brown tide, *Aureococcus anophage*ff*erens,* were collected from 0.5 m below the surface using a peristaltic pump, and placed in 2-L acid-rinsed polyethylene amber bottles. The samples were kept in ice until arrival at the lab. Water samples were then filtered using a Whatman GF/F 47 mm filter (GE Healthcare, Fisher Scientific, Hanover Park, IL, USA) and a vacuum pump for 1–4 min. A minimum of 1 L by volume of water was filtered except in situations when the filter was clogged. Chlorophyll *a* and 19'-Butanoyloxy-Fucoxanthin (But-fuco) concentrations were determined, as part of the integrated monitoring of the MCBs, using a High Performance Liquid Chromatography [38]. Freshwater discharge data for Birch branch at Showell, Maryland, USA were downloaded from the United States Geological Survey (USGS) website (http://waterdata.usgs.gov/nwis/monthly).

#### *2.3. Data Analyses*

The Mann–Whitney U test was used to compare mean abundance of mesozooplankton between years for each season. One-way ANOVA was used to compare differences in mean abundance of mesozooplankton among sites and embayments after log(x+1) transformation of data. The relationships between non-copepod taxa distribution and environmental factors were examined by canonical correspondence analysis (CCA), using CANOCO version 4.5 (Microcomputer Power, Ithaca, New York, USA.) [39]. Non-copepod mesozooplankton abundance was log(x+1) transformed; rare species were down-weighted when necessary and environmental data were standardized. Non-copepod mesozooplankton abundance matrices included only taxa that appeared in >10% of the sampling sites, representing >1% in 2012 and >0.1% in 2013 of the non-copepod community. Shannon diversity index [40] was computed using PRIMER 6 (PRIMER-e, Auckland, North Island, New Zealand. [41]. Spearman's rank correlation was applied to identify the environmental factors (water quality variables and biota) that were associated with diversity.

#### **3. Results**

#### *3.1. Environmental Factors*

Mean monthly water temperature in 2012 varied from 5◦C in February to 25◦C in July and then decreased to about 5 ◦C in December (Figure 2a). In 2013, mean temperature was 30.7 ◦C in July and 7 ◦C in December. Water temperature was on average lower in 2012 than in 2013. Lower water temperatures were observed at sites close to the inlets, while St. Martin River had the highest mean water temperature (Table 1).

**Figure 2.** Seasonal patterns of environmental variables expressed as mean values ± SE: (**a**) Temperature ( ◦C); (**b**) salinity; (**c**) Freshwater discharge (m<sup>3</sup> s<sup>−</sup>1); (**d**) dissolved oxygen (mg L<sup>−</sup>1); (**e**) Chlorophyll-α (μg L<sup>−</sup>1); and (**f**) pH.

Mean monthly salinity was higher in 2012 (33.4) than in 2013 (26.6) during the period of March to August, but the reverse occurred from September to November (Figure 2b). This salinity pattern reflected the pattern of freshwater discharge by Birch Branch at Showell, MD (Figure 2c). From January to July, Birch Branch mean discharge was higher in 2013 (0.45 m<sup>3</sup> s<sup>−</sup>1) than in 2012 (0.11 m3 s<sup>−</sup>1), and the reverse was the case from August to November (Figure 2c). Discharge was highest in October 2012, the month that Hurricane Sandy occurred in the area. Mean monthly salinity was between 23.6 and

34.7. Newport Bay had the lowest mean salinity level and Sinepuxent Bay had the highest salinity (Table 1). Salinity was also lower at sites in St. Martin River and Assawoman Bay relative to sites in Chincoteague and Isle of Wight Bays.

Monthly dissolved oxygen (DO) levels on the average were high in fall and winter, but relatively lower in spring and summer, although dissolved oxygen levels were > 4.6 mg L−<sup>1</sup> at all sites throughout the duration of this study (Figure 2d). The lowest DO level was recorded in July (4.6 mg L<sup>−</sup>1), and the highest in March (12.7 mg L<sup>−</sup>1).

The seasonal patterns of chlorophyll *a* concentrations are similar in 2012 and 2013. Monthly levels of chlorophyll *a* averaged between 0.2 and 3.4 μg L−<sup>1</sup> (Figure 2e). Chlorophyll *a* levels were relatively high in winter but decreased in spring, particularly in April. Chlorophyll *a* peaked in summer (July) in both years. Site 6 (Newport Bay) had the highest mean level of chlorophyll *a* in the MCBs during this study (Table 1). Monthly mean pH levels in the MCBs ranged from 7.8 to 8.2 and did not vary much within the system (Figure 2f).

**Table 1.** Mean physico-chemical and biotic factors measured at sampling stations within embayments of the Maryland Coastal Bays. Embayment codes are: CB, Chincoteague Bay; NB, Newport Bay; SB, Sinepuxent Bay; IWB, Isle of Wight Bay; SMR, St. Martin River; AB, Assawoman Bay. Temp, temperature; Sal, salinity; Chl *a*, chlorophyll *a*; But-Fuco, 19'-Butanoyloxy-Fucoxanthin; Cteno, Ctenophores; CPUE, Catch Per Unit Effort.


#### *3.2. Temporal Pattern of Mesozooplankton Abundance*

Mesozooplankton abundance expressed as individuals m−<sup>3</sup> (ind. m−3) varied from 5 ind. m−<sup>3</sup> (Site 12) in May 2012 to 410,509 ind. m−<sup>3</sup> (Site 8) in September 2013. The average monthly abundance was high (24,543 ind. m<sup>−</sup>3) in February, but peaked in early March (36,878 ind. m−3) during winter of 2012; it declined to a relatively low level in spring and then gradually increased from summer through fall (Figure 3). Relatively high abundance (30,032 ind. m−3) of mesozooplankton was observed in December 2012. A gradual decline was also observed from winter 2013 through early spring. Densities remained low (378 ind. m−3, May 2012) in spring and early summer, peaked (34,828 ind. m−3) in September 2013 and then declined in fall.

Mesozooplankton abundance differed significantly between years (Table 2) in summer and fall (p < 0.001), but not in winter and spring (p > 0.05). In fall, abundance was significantly higher in 2012 than 2013 (Mann–Whitney U = 261, T = 2040.0, p ≤ 0.001), but in summer, it was lower in 2012 than in 2013 (U = 392.0, T = 1909.0, p ≤ 0.001).

Seasonal and spatial patterns of abundance of mesozooplankton are presented in Figure 4. For each of the seasons examined, there were no significant differences among sites in mesozooplankton abundance (ANOVA, p > 0.05) except in summer 2012 (ANOVA, p = 0.05) when the highest density was recorded (Figure 4 and Table 3) at Site 7 (Tukey test, p < 0.05). When mesozooplankton abundance was compared among bays, higher abundance was recorded in Sinepuxent Bay than in any other bay in winter and summer 2012 (ANOVA, p < 0.05), as shown in Table 3.

**Figure 3.** Temporal pattern in the abundance (±SE) of mesozooplankton in the Maryland Coastal Bays.

**Table 2.** Seasonal average densities (ind. m<sup>−</sup>3) of mesozooplankton in the Maryland Coastal Bays. Winter, January–March; spring, April–June; summer, July–September; fall, October–December.


**Table 3.** One-way ANOVA result of analyses comparing mesozooplankton density among sampling sites and embayments of the Maryland Coastal Bays in the different seasons in 2012 and 2013. d.f, degree of freedom; Fs, F-value; P, probability value; \* 0.05 > p > 0.01; n.s, p> 0.05. There were insufficient data to perform ANOVA in winter 2012.


**Figure 4.** Spatial patterns in the densities (+SE) of mesozooplankton in the Maryland Coastal Bays in 2012 (year of below average freshwater discharge) and 2013 (year of above average freshwater discharge): (**a**) winter; (**b**) spring; (**c**) summer; and (**d**) fall.

*3.3. Community Composition, Relative Abundance and Distribution of Non-Copepod Mesozooplankton in MCBs*

Mesozooplankton community in the MCBs was represented by 14 major taxa dominated by copepods, but only groups that are relatively abundant are presented in Figure 5. Non-copepods

contributed about 5% of the total mesozooplankton abundance. A detailed description of copepod species composition and abundance in the MCBs has been made [42].

**Figure 5.** Mesozooplankton community composition in the Maryland Coastal Bays based on samples collected in 2012 and 2013.

The non-copepod community on the average was dominated by hydromedusae (3%), followed by gastropod veliger (0.5%), and barnacle nauplii (0.5%) (Figure 5). Cladocerans (0.1%) and larval stages of decapods (0.2%), polychaetes (0.2%), bivalves (0.1%), and fish (0.1%) were also relatively abundant within the non-copepod community. Larvaceans, nematodes, amphipods, and isopods made up <1% of the zooplankton community.

The percentage contribution of non-copepods to total zooplankton abundance was as high as 98% at Site 12 (AWB) in the summer when copepod relative abundance was low. Monthly averaged contribution varied between 0–31.9% in 2012 and 0.4–76.5% in 2013 (Figure 6). Non-copepods were most represented in May 2012 and June 2013.

**Figure 6.** Percentage contribution by number of non-copepod taxa to total mesozooplankton abundance (+SE) in the Maryland Coastal Bays.

In 2012, barnacle nauplii were the most abundant non-copepods, with mean density of 810 ind. m−<sup>3</sup> in Newport Bay (Figure 7a,b). Annual peaks were observed in April 2012 (368 ind. m−3) and

September 2013 (102 ind. m<sup>−</sup>3). The maximum barnacle nauplii density (4645 ind. m−3) was attained in April 2012 in Newport Bay and mean abundance of barnacle nauplii was significantly higher in Newport Bay than in other embayments (ANOVA, p = 0.02). Although, barnacle nauplii were absent in fall 2012, they were present in all seasons in 2013 (Figure 8a–d). Maximum density in 2013 was 1119 ind. m−<sup>3</sup> in September at Site 3 (Chincoteague Bay). Barnacle nauplii were more abundant in spring and summer than in the fall and winter (ANOVA, p < 0.001).

**Figure 7.** Spatial distribution, abundance (**a**,**c**) and percentage composition (**b**,**d**) of non-copepod mesozooplankton in the MCBs during 2012 (below average) and 2013 (above average) freshwater discharge from winter to summer. CB, Chincoteague Bay; NB, Newport Bay; SB, Sinepuxent Bay; IWB, Isle of Wight Bay; SMR, St. Martin River; AB, Assawoman Bay.

**Figure 8.** Temporal patterns in abundance (**a**,**c**) and percent contribution (**b**,**d**) of non-copepod mesozooplankton in the Maryland Coastal Bays during 2012 (year of below average freshwater discharge from winter to summer) and 2013 (year of above average freshwater discharge from winter to summer).

In 2013, hydromedusae were the most abundant non-copepod mesozooplankton (Figure 7c,d). Mean densities as high as 1367 ind. m−<sup>3</sup> and 5156 ind. m−<sup>3</sup> were recorded in Assawoman Bay and St. Martin River, respectively, in 2013. Hydromedusae were not observed at St. Martin River and Newport Bay sites in 2012 and also were not observed in Sinepuxent Bay in 2013 (Figure 7b,d). They were present from March to August in 2012 and from March to September in 2013, but were not found in May 2013 (Figure 8a–d). Average hydromedusae abundance was less than 20 ind. m−<sup>3</sup> each month

and at all sampling sites in 2012. Monthly mean densities peaked in April, 2012 (16 ind. m<sup>−</sup>3) and in March 2013 (5252 ind. m−3). The highest density recorded in 2013 (50,399 ind. m−3, March, Site 10) was 298 times greater than the highest density (181 ind. m<sup>−</sup>3, March, Site 5) in 2012.

Gastropod larvae were intermittently present and their monthly mean density reached a peak in May 2012 (158 no m<sup>−</sup>3, Figure 8a). A secondary peak (103 ind. m<sup>−</sup>3) occurred in July of the same year. Mean density in 2013 peaked in June (721 ind. m<sup>−</sup>3). Gastropod larvae were abundant (107 ind. m−3) in Sinepuxent Bay, but were not observed in Newport Bay in 2012 (Figure 7a–b). A maximum density (1299 ind. m−3) of the larvae was observed in May 2012 at Site 7 (Sinepuxent). In 2013, maximum density (1852 ind. m−3) occurred at Site 4 (Chincoteague Bay) in June. Gastropod larvae were least abundant in Assawoman Bay (2 ind. m<sup>−</sup>3) in 2013.

Polychaete larvae occurred from February through August in 2012 and were most abundant (266 ind. m<sup>−</sup>3) in March (Figure 8a). Mean polychaete larval densities never exceeded 10 ind. m−<sup>3</sup> in 2013, occurred sporadically, and were most abundant in March (9 ind. m<sup>−</sup>3). Polychaete larvae were more abundant in 2012 than 2013 (Mann–Whitney U-test, p = 0.0002). Peak numbers were 2384 ind. m−<sup>3</sup> in March 2012 at Site 10 located in St. Martin River and 60 ind. m−<sup>3</sup> in March 2013 at Site 5 in Chincoteague Bay (Figure 7a–d).

Bivalve larvae were present from February through August in 2012 (Figure 8a), but their occurrence was irregular in 2013 (Figure 8c) with mean monthly densities < 5 ind. m<sup>−</sup>3. Bivalve larvae were better represented in plankton samples from 2012 than 2013 (Mann–Whitney U-test, p = 0.001) and relatively abundant in March during both years (Figure 7). Average density was 161 ind. m−<sup>3</sup> in March 2012, and 4 ind. m−<sup>3</sup> in March 2013. Maximum densities of 1490 ind. m−<sup>3</sup> at Site 10 (St. Martin River) in 2012 and 58 ind. m−<sup>3</sup> at Site 8 (Sinepuxent Bay) in 2013 were recorded in March. Bivalve larvae were only found in Sinepuxent Bay and Isle of Wight Bay (Site 9) in 2013, but their distribution in 2012 was widespread except at Sites 3 and 6 where they were not observed. Densities did not differ significantly amongst embayments (ANOVA, p = 0.12).

Cladocerans (*Evadne nordmanni*, *E. spinifera*, *Pseudoevadne tergestina, Pleopis polyphemoides*, *Podon intermedius*) were most abundant in Sinepuxent and Isle of Wight Bays (ANOVA, p = 0.002) during summer and were absent from Newport Bay (Figure 7b).

Decapod larvae (shrimps and crabs) were prevalent in Sinepuxent Bay where they occurred from March through November in 2012 and 2013 (Figure 7; Figure 8). In July, they comprised 42% and 67% of non-copepod zooplankton abundance, and 10% and 11% of total mesozooplankton abundance in 2012 and 2013, respectively. Larval stages of crabs *Panopeus herbstii*, *Neopanope texana*, *Uca* spp. and *Pinnixia* spp. were observed in zooplankton samples collected in 2012. Zoeae of *Callinectes sapidus*, *Rhithropanopeus harrisii*, *Ovalipes ocellatus*, *Cancer irroratus*, *Emerita* spp., *Hemigrapus* sp., *Libinia dubia*, *Ocydopode spp., Petrolisthes armatus*, and *Lepidopa websteri* were also present.

Fish larvae were first observed in April 2012 (mean: 16 ind. m<sup>−</sup>3) when they were most abundant (Figure 8). They were observed in December 2013 (0.1 ind. m−3), but mean densities were very low. The highest density of fish larvae (202 ind. m<sup>−</sup>3) was observed at Site 6 in Newport Bay (April 2012).

Larvaceans occurred from April to June 2012 and in August 2013. Average densities ranged from 1 to 5 ind. m−3. The maximum density of larvaceans was 23 ind. m−3, at Site 8 (Sinepuxent Bay) in August, 2013. Other non-copepods irregularly found in samples included amphipods, isopods, nematodes, and ostracods.

#### *3.4. Relationships Between Non-Copepod Mesozooplankton and Environmental Factors*

In 2012 (Figure 9a), the first two CCA axes explained 47.1% of the cumulative percentage variance of taxa-environment relationship and the correlations were 0.89 and 0.89 for the first and second axes, respectively. The biplot in 2012 showed that the main variability (29.6% of variance) was due to the positive association of larval shrimp and bivalves with salinity (p = 0.046). In contrast, barnacle nauplii and ostracods were negatively related to salinity (p < 0.05). Gastropod and polychaete larvae showed

a positive relationship with pH (p = 0.022), while crab zoeae, *Evadne* spp. and *Pleopis* sp. related inversely to pH. Hydromedusae neither associated closely with salinity nor pH in 2012.

**Figure 9.** Results of CCA analysis showing environmental variables and taxa relationships in: 2012 (**a**); and 2013 (**b**) in MCBs. BA, bay anchovy; DO, Dissolved oxygen; Sal, Salinity; Chl *a*, Chlorophyll *a.*

In 2013 (Figure 9b), the first two axes of the CCA explained 70% of the cumulative variation in taxa-environment relationship, and the correlation was 0.98 for the first axis and 0.91 for the second axis. The taxa-environment biplot (Figure 9b.) showed that the variability (63.7% of variance) was accounted for by the positive association of gastropod veliger with salinity (p = 0.026), and inverse relationship of hydromedusae to salinity (p < 0.05). Barnacle larvae were negatively aligned to Chl. *a* (p = 0.034); bay anchovy larvae were positively correlated to DO (p = 0.04), which was also plotted inversely with cladocerans, decapod zoae and polychaete larvae. Temperature, pH, and ctenophores correlated strongly with DO and were not selected by the model.

#### *3.5. Mesozooplankton Diversity*

The diversity of mesozooplankton varied temporally and spatially and was affected by salinity such that diversity was higher at higher salinity than at lower salinity (Figure 10). The Shannon–Wiener diversity index was high (≥2.5) in 2012 at all sites except at Site 6 (H = 2.2) in Newport Bay with the lowest diversity (Figure 10). Diversity index value in Newport Bay (Site 6) in 2012 (H = 2.2) was slightly lower than in 2013 (H = 2.3). On the average, diversity was significantly higher (t = 6.2; p ≤ 0.001) in 2012 than in 2013. Sinepuxent Bay supported the highest diversity (H = 2.8 in 2012; H = 2.6 in 2013) of mesozooplankton in both years. Diversity decreased from Sinepuxent Bay (Sites 7 and 8) close to the Ocean City Inlet to Assawoman Bay (Site 13) in the northern MCBs in both 2012 and 2013. Likewise, diversity in Chincoteague Bay decreased from Site 5 to Sites 1 and 2 in both years.

**Figure 10.** Spatial (**top**) and seasonal (**bottom**) patterns in Shannon–Wiener diversity index of mesozooplankton in the MCBs in 2012 and 2013 from winter to summer.

Diversity was higher in spring and summer (H > 2.4) than in fall (H < 1.9) of 2012. In 2013, diversity in March and from September to December was higher (t = −2.3; p ≤ 0.06) than diversity in 2012. Periods of maximum and minimum diversity index values were different in 2012 and 2013. Mesozooplankton was most diverse in June (H = 2.9) and March (H = 2.5), but least diverse in October (H = 1.4) and July (H = 1.8) in 2012 and 2013, respectively.

When examined spatially, diversity indices values in 2012 were significantly correlated negatively with Chl. *a* (r = −0.7, p = 0.01) and positively with salinity (r = 0.7, p = 0.02), but not with But-fuco (r = −0.5, p = 0.08). In 2013, diversity also correlated positively with salinity (r = 0.3, p < 0.050). All other environmental factors were weakly correlated with mesozooplankton spatial diversity in 2013, and no significant relationship was observed (p > 0.05). With respect to monthly diversity, H was significantly correlated with ctenophores (r = 0.7, p = 0.02) in 2012, but not in 2013 (r = − 0.1, p = 0.9). Other variables such as temperature (r = − 0.4, p = 0.2) and salinity (r = − 0.4, p = 0.3) showed weak negative correlations with diversity in 2013.

#### **4. Discussion**

This is the first study to describe mesozooplankton assemblage dynamics in the MCBs. The mesozooplankton composition, abundance and diversity varied spatially and temporally in the MCBs, driven in part by fluctuations in environmental factors, particularly salinity. Lower salinity due to increased freshwater input coupled with higher ambient temperature in 2013 (Figure 2) relative to that of 2012, perhaps promoted high densities of hydromedusae in winter and early spring. Gelatinous zooplankton can exert high predatory pressure on crustacean zooplankton species as well as on gastropod veligers, barnacle nauplii, polychaete larvae, and ichthyoplankton (fish larvae) [32,43–47]. The feeding rate of ctenophores, in particular, increases with water temperature [48], possibly contributing to the decline in annual peak abundance observed for some non-copepod taxa (polychaete larvae, ostracods, bivalve and fish larvae). Nevertheless, the high freshwater inflow associated with Hurricane Sandy in fall 2012 accompanied by the relatively low salinity that occurred through spring 2013 might have contributed to the observed low abundances of the taxa in the MCBs.

#### *4.1. Mesozooplankton Community Composition and Taxa-Environment Relationships*

The percent contribution of non-copepods (e.g., decapod larvae, cladocerans) to mesozooplankton abundance in this study was highest in late spring and early summer when copepod density was relatively low [42]. Similar observations were made in the lower Narragansett Bay [16] and in other mid-Atlantic estuaries [6,19,35,49]. The increase in non-copepod taxa relative abundance that occurred when copepod density was low is important since they serve as alternate prey for zooplanktivorous fishes. Bay anchovy is abundant in summer in MCBs [50], and the juveniles usually make up the majority of estuarine bay anchovy population biomass in summer and fall [51]. Although copepods contribute the bulk of bay anchovy's diet, non-copepod taxa, especially larvae of barnacles, bivalves, and polychaetes, as well as amphipods, crab zoeae, mysids, ostracods, cladocerans, and ichthyoplankton seasonally dominate the diet in mid-Atlantic estuaries such as Long Island Sound, Cheseapeake Bay and MCBs [8,9,52–55].

Temporal variations in freshwater discharge and salinity influenced the relative abundance of mesozooplankton taxa. For example, Barnacle nauplii dominated the non-copepod community, especially in Newport Bay from January to August 2012. During this period, freshwater discharge into MCBs was relatively low and salinity was high. Naupliar stages of barnacles were absent in fall of 2012, which might have been due to the dramatic reduction in salinity caused by excessive rainfall that was associated with Hurricane Sandy that increased freshwater inputs into the system possibly increasing mortality of the larvae and perhaps the adult stages. Mortality of barnacle larvae is higher at low salinity (<6) than at higher (>14) salinity [56], and high rainfall with associated high freshwater discharge has been reported to depress the abundance of some species of barnacles. Higher densities of barnacles were collected during dry season than wet season in coastal lagoons in Panama [56]

and in Lagos Harbor, Nigeria [57]. The salinity observed in the MCBs in spring of 2013 (wet year) was much higher (23–26) than the levels observed during the rainy season in these other lagoons, though the salinity in MCBs during the period Hurricane Sandy occurred might have been much lower than the values we recorded after the event. Spring and summer peaks of barnacle naupliar observed in the MCBs were also reported for the Delaware River estuary [19]. In 2013, a year of above average freshwater discharge into MCBs and lower salinity, hydromedusae dominated non-copepod community particularly in Assawoman Bay and at the mouth of St. Martin River. Larvae of polychaetes and bivalves, in addition to ostracods, and larvaceans, were better represented in the plankton in 2012 when salinity was higher, than in 2013 when salinity was lower. In estuaries that experience much horizontal gradient in salinity, seawater penetration can bring in zooplankton species from the coastal ocean during periods of low freshwater inputs [11].

There were seasonal changes in the relative abundance of mesozooplankton taxa in the MCBs. Hydromedusae were present from winter through summer, but scarce in fall (October–December). Hydromedusae were abundant in Delaware Bay River estuary during winter [18] and in Long Island estuary in April [49]. Polychaete and bivalve larvae were also abundant in winter to spring, with a shift in seasonal dominance from polychaetes to barnacles, gastropod larvae, crab zoea, barnacles, shrimps, and finally polychaete larvae in 2012. This seasonal pattern was different in 2013 when a shift in dominance of the major taxa occurred from hydromedusae to gastropod larvae, crab zoea, barnacles, and finally polychaetes.

The spatial distributions of mesozooplankton taxa were related to variations in salinity. Larvaceans were positively correlated to salinity, and found only in Sinepuxent Bay close to the Ocean City inlet in 2013. Bivalve larvae also showed a similar preference for high salinity bay areas in 2013, and gastropod larval densities decreased at sites near the mouths of MCBs tributaries. The relatively high abundance of barnacle nauplii in Newport Bay and St. Martin River in both years as well as the occurrence of polychaete and bivalve larvae in high abundance in St. Martin River in the dry year (2012) and within the open waters of Sinepuxent Bay in the wet year (2013), suggests that salinity is an important factor that influenced their distributions and abundance. As salinity increases towards the ocean, clams, crustaceans and polychaete worms dominate the benthos of the coastal bays [27], thus accounting for the high abundance of meroplanktonic larvae of these organisms around Sinepuxent Bay. The occurrence of gastropod and polychaete larvae in high abundance in embayments closest to the coastal ocean (downstream) is contrary to the findings in Pagan River, a sub-estuary of Chesapeake Bay [17] and in Mondego estuary, Western Portugal [58] where the larvae were abundant upstream, although this might be due to differences in the species reported in the systems. In the Delaware River estuary, polychaete larvae occurred from the entrance of the river to the middle of the bay and peaked in abundance from winter to spring [19]. The lowest salinity at which they were observed was 24.1. In this study, annual peaks in the abundance of polychaete larvae were recorded at salinity levels >30.

#### *4.2. Occurrence of Blue Crab Larvae in the MCBs*

Blue crab, *C. sapidus* zoeae, were common in July in this study which is similar to the month during which they were observed by other investigators in nearby estuaries [18,37]. The study confirms what was previously suspected based on the capture of gravid females by the MD DNR that blue crabs spawn in the MCBs. Blue crab zoeae were most prevalent in Sinepuxent Bay and the northern bays and rarely occurred in southern bay areas with adequate salinity for larval development, although no sampling occurred near the Chincoteague Inlet in Virginia. This supports results of previous studies that indicated that blue crabs hatch their eggs near the mouth of estuaries [6,22]. The occurrence of zoeae at sites located in Assawoman Bay close to tributary creeks, however, suggests that spawning occurs not only close to the Ocean City Inlet, but also in areas that are some distance away from the inlet. It is also possible that crab larvae that hatched near the inlets were dispersed into Assawoman Bay by waves and tidal action. The proximity of spawning locations to the inlet enhances easy transportation of zoeae from the MCBs into the nearby continental shelf waters where further development takes

place [22,25]. Additional studies are needed, however, to determine if blue crab larvae complete their development from zoeae to megalopae within the MCBs.

#### *4.3. Diversity of Mesozooplankton*

Mesozooplankton assemblage was most diverse in spring and summer of 2012 when salinity was relatively high, especially in areas close to the Ocean City Inlet, connecting the bays to the Atlantic Ocean. Intrusion of seawater in estuaries can act as a major driver for zooplankton diversity [59]. In Pearl River estuary China, zooplankton diversity and species richness were enhanced in the middle and lower portions with high diversity at salinity level >25 [11]. The higher zooplankton diversity in spring and summer months in 2012 might have been contributed by the increase in the abundance of cladocerans, decapod larvae, and other meroplanktonic forms in the system during that period, some of which were likely transported into the bays from coastal ocean by tidal action.

The most diverse assemblage was observed in areas with high salinity, but low chlorophyll *a*. The observed negative correlation between chlorophyll *a* and diversity was because chlorophyll *a* levels were higher in areas close to MCBs tributaries (Newport Bay and St. Martin River) with lower salinity and high nutrient levels due to freshwater inflow [38,60], which also had lower mesozooplankton diversity than in areas close to the inlets. Howson et al. [35] observed higher zooplankton diversity in the southern bay of the Barnegat Bay, NJ, with more oceanic influence and less anthropogenic impact, than in the more nutrient enriched northern bay.

#### *4.4. Implications of Temporal Variability in Environmental Factors on Mesozooplankton*

Variability in freshwater discharge altered salinity and influenced the abundance and distribution of mesozooplankton in the MCBs. High freshwater discharge that began in fall 2012 and continued through spring 2013 lowered salinity that likely favored hydromedusae whose abundance was highest in March and April, but perhaps disfavored bivalve and barnacle larvae. During this period, salinity was about 23 (March) and 26 (April) compared with about 30.5 (March) and 34 (April) in 2012. Salinity has been reported to influence the abundance and distribution of hydromedusae in other estuaries [61,62]. In Mississippi Sound, the incidence of hydromedusae was highest at 25.1 to 30 below and above which the abundance decreased [61]. With the occurrence and high abundance of hydromedusae, ctenophores [42] and scyphozoan jellyfish in the MCBs, and their documented negative impacts on the standing stock of planktonic communities [31,45,47,48], the slow rate of recovery of the hard clam population abundance in the MCBs is not surprising. This is because, peak abundance of bivalve larvae (March) coincided with the increase in the abundance of hydromedusae (March and April), which was followed by ctenophores (April–June) especially in 2013 when salinity was relatively low [42]. Hydromedusae are known predators of zooplankton [44] and when numerically abundant, are able to exert their potential maximum clearance impact [63]. Ctenophores as predators of bivalve larvae [31,43,64] cause high mortality, with a potential of consuming up to 94.1% of the larvae during spawning periods [32].

Non-copepod zooplankton are important prey items for carnivorous zooplankton, larval fish, and adult planktivorous and forage fish such as bay anchovy. These organisms are relatively abundant when copepod abundance is low in the MCBs. During years with high freshwater input, the MCBs salinity levels in some areas may not be favorable for the development of some decapod crab zoeae and bivalve larvae (*Mercenaria mercenaria*) that prefer high salinity. In addition, the temperature range (~17–30 ◦C) tolerated by *M. mercenaria* larvae is reduced when salinity decreases from the optimum range of 26–27 [65]. Thus, increased temperature, reduced salinity as well as high densities of gelatinous zooplankton (hydromedusae and ctenophores) can affect the breeding success of finfish and shellfish species in the MCBs. Gelatinous zooplankton can change food web dynamics and decrease recruitment of fish and shellfish populations by preying on the larval forms or via competition for zooplankton prey [45,46].

#### **5. Conclusions**

Results of this study indicate that variations in freshwater discharge affect the abundance, distribution and diversity of mesozooplankton in MCBs. These changes in mesozooplankton structure likely influence trophic relationships and recruitment of finfish and shellfish larvae to the adult populations. The percent contribution of non-copepods to mesozooplankton abundance in MCBs was highest in late spring and early summer when copepod density was relatively low. The composition and abundance of mesozooplankton differed between 2012 and 2013 such that barnacle larvae were the most abundant of non-copepod mesozooplankton, whereas hydromedusae were the most abundant group in 2013. Mesozooplankton diversity varied temporally and spatially, and was higher at higher salinity than at lower salinity. Seasonally, diversity was higher in spring and summer 2012 than in fall. Additionally, in 2012 when salinity was comparatively high because of low precipitation, diversity was higher than in 2013 when salinity was relatively low. Blue crab zoeae were common in the MCBs in the summer, especially in Sinepuxent Bay and the northern bays confirming that the adults spawn in the MCBs, although further studies are needed to determine if they complete their life cycle in the bays. Additional studies are also needed to investigate the predation rates of ctenophores, jellyfish and zooplanktivores such as bay anchovy on various zooplankton taxa. This will be useful for understanding the role of non-copepod community as alternate prey for fish in the MCBs during spring and summer minimum of copepods, and will also provide more information on the impact of predation on the survival of bivalve larvae and their recruitment to the adults.

**Author Contributions:** This study was conducted as a joint effort of all authors. Conceptualization, P.C. and E.U.O.; Methodology, E.U.O. and P.C.; Software, E.U.O.; Validation, E.U.O. and P.C.; Formal Analysis, E.U.O.; Investigation, E.U.O. and P.C; Resources, P.C.; Data Curation, E.U.O.; Writing—Original Draft Preparation, E.U.O.; Writing—Review and Editing, E.U.O. and P.C.; Visualization, E.U.O.; Supervision, P.C.; Project Administration, P.C.; and Funding Acquisition, P.C.

**Funding:** This publication was made possible by the National Science Foundation Center for Research Excellence in Science and Technology (NSF-CREST) award number (1036586) to the Center for the Integrated Study of Ecosystem Processes and Dynamics, University of Maryland Eastern Shore, and in part by the National Oceanic and Atmospheric Administration, Office of Education, Educational Partnership Program award number (NA16SEC4810007) to the Living Marine Resources Cooperative Science Center at the University of Maryland Eastern Shore. Its contents are solely the responsibility of the award recipient and do not necessarily represent the official views of the U.S. Department of Commerce, National Oceanic and Atmospheric Administration.

**Acknowledgments:** We would like to thank Kam Tang and Jamie Pierson for their valuable comments on the earlier draft of this manuscript, and Capt. Christopher Daniels (Boat Captain), Ozuem Oseji, Andres Morales-Nunez, Kingsley Nkeng, Chinwe Otuya, Abena Okyere Acheampong and Jessica Mazile for help with collection and enumeration of plankton samples. Data for bay anchovy were provided by Steve Doctor, Maryland Department of Natural Resources (MD DNR). Young Tsuei (DOEE) and Tracy Bishop (UMES) assisted with creation of map of the sampling sites.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript or in the decision to publish the results.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
