Plankton Taxonomic and Size Diversity of Mediterranean Brackish Ponds in Spring: Influence of Abiotic and Biotic Factors

Abstract

In this study, performed in Mediterranean brackish ponds during spring season, we assessed the effects of biotic interactions and abiotic factors on the size and taxonomic structure of the phytoplankton and zooplankton. We used a taxonomic and a size diversity index as a descriptor of the community structure. We predicted that the size diversity of each trophic level would be mainly related to biotic interactions, such as size-based fish predation (in the case of zooplankton) and food resource availability (in the case of phytoplankton), whereas taxonomic diversity would be more affected by abiotic variables (e.g., conductivity, pond morphology). Our results showed a negative relationship between phytoplankton size diversity and food resource availability leading to low size diversities under food scarcity due to dominance of small species. Conductivity also negatively affected the phytoplankton size diversity, although slightly. Regarding zooplankton size diversity, none of predictors tested seemed to influence this index. Similar fish size diversities among ponds may prevent a significant effect of fish predation on size diversity of zooplankton. As expected, taxonomic diversity of phytoplankton and zooplankton was related to abiotic variables (specifically pond morphometry) rather than biotic interactions, which are usually body size dependent, especially in these species-poor brackish environments.

1. Introduction

Mediterranean ponds are ecologically very important ecosystems that support relevant hydrological, chemical and biological processes and are biodiversity hotspots in terms of both species composition and biological traits [1]. They are also very vulnerable habitats, because they are threatened by several anthropogenic pressures [2,3]. Nevertheless, these habitats have received less scientific attention than other water bodies such as lakes or rivers [4,5].

Most ecological studies of Mediterranean ponds have been focused on their taxonomic diversity [2], the dynamics of a certain species and communities [6,7,8,9] and its relationship with nutrients dynamics and hydrological patterns [10,11,12,13], as well as with anthropogenic pressures [14]. However, there are not many studies considering size-based interactions among adjacent trophic levels (predators and preys) of the food web [15,16,17] despite such interactions play a key role in the trophic structure and functioning of aquatic ecosystems [15,16,17,18]. It is worth to mention that size-based interactions are particularly relevant in the species-poor communities like those of Mediterranean brackish ponds where trophic interactions are mainly body-size dependent [19,20]. A way to analyze the size structure of a community is through the size diversity index, which is analogous to Shannon diversity index but adapted for continuous variables, such as body size [21,22]. In Mediterranean brackish ponds, it has been shown that size diversity provides complementary information about the community structure to classical taxonomical approaches. While taxonomic diversity of zooplankton and macroinvertebrates communities was more sensitive to abiotic factors such as nutrient availability, size diversity was mainly related to biotic interactions (e.g., predation or inter-and/or intraspecific competition) [19,23]. However, abiotic factors can also affect the size structure of a community. For example, in other European brackish ecosystems, an increase in conductivity has been related to a decrease in zooplankton mean size and size diversity due to the replacement of large cladoceran species by copepods, small cladoceran species and rotifers, which are overall more salt-tolerant [24,25,26].

Size structure within a trophic group of the planktonic food web can be determined by the size structure of the adjacent trophic levels since both predation and food selection are size-dependent [16]. However, each trophic group may be affected in different ways by top-down (e.g., size-based predation) or bottom-up (e.g., size diversity of resources) controls. Several studies in aquatic ecosystems have shown that fish predation is the most important driver of zooplankton size structure, as an increase in fish density is related to a decrease in density, mean body size and size diversity of zooplankton [27,28,29]. However, the effects of size-structured predation (i.e., predation by individuals of different sizes) on the individual size structure of prey in natural food webs are less well understood, although this information brings insight into the strength of the interactions between adjacent trophic levels and into the biomass transfer through the food web [15,16,30,31]. A recent study [30], found a negative relationship between predator (mesozooplankton) and prey (nanomicroplankton) size diversities which was explained due to the enhanced strength of top-down control at increasing predator size diversity. Nevertheless, the opposite situation has also been found—positive relationship between the size diversity of predators (fish) and prey (zooplankton)—suggesting that a higher diversity of sizes in consumers may promote diversification of resources by size [16]. On the other hand, size diversity of a trophic group may also be determined by bottom-up control since low resource availability often results in higher size diversity of consumers, as has been found in zooplankton and fish communities [21,28]. Thus, when resources are scarce, the competitive pressure for resources between small and large individuals can be reduced by size-based selection of food, where large and small predators tent to prefer large and small prey respectively, allowing the coexistence of a wide range of organism sizes [21,28,32,33]. Concerning phytoplankton, recent studies have shown that its size diversity is mainly determined by changes in resource availability, rather than by predation, as due to the prevailing influence of abiotic factors in their nutrient uptake [16,28,31,34].

Whereas size diversity is usually more sensitive to biotic interactions, taxonomic diversity had been found to be more related to abiotic factors, such as nutrient concentration [19,23,35], conductivity [23,35] and water body size [36]. In shallow ecosystems, zooplankton and benthic macroinvertebrates taxonomic diversity was found to be negatively related to increasing concentrations of phosphorus [19,37] and total organic carbon [23] respectively, whereas nitrogen and phosphorous were found to decrease phytoplankton taxonomic diversity [38]. Regarding conductivity effects, in Mediterranean brackish communities a decrease in taxonomic diversity of zooplankton [23,35], macroinvertebrates [8], and phytoplankton [39] was observed at highest conductivities. In these communities, conductivity plays an important role for shaping species composition and food web interactions [24,25,27,40,41]. Increased conductivity is usually related to reduced richness due to low osmoregulatory ability of several species [9,24,26,42,43]. Significant positive regression between area of the water body and the number of species had been found for ponds gastropods [44], for macroinvertebrates in streams and ponds [36,45], and for phytoplankton and zooplankton in lakes [37,46,47,48]. In a study on Danish shallow lakes, lake depth was also positively related to phytoplankton species richness [37].

Although both taxonomic and size-based approaches can be used to complement each other, studies dealing with their combination when determining plankton community-structuring factors are still scarce [19]. The combination of both approaches becomes especially relevant in species-poor communities, like those in Mediterranean brackish ecosystems, where harsh environmental conditions make size-based interactions more relevant [23]. In this study, performed in Mediterranean brackish ponds, we aimed to determine the effects of biotic interactions (i.e., predation, food resource availability) and abiotic factors (environmental variables) on the phyto- and zooplankton taxonomic and size diversity indexes, and explore if they are affected by the same factors. To analyze the effects of biotic interactions we searched for relationships among diversity metrics of adjacent trophic levels (phytoplankton, zooplankton, and fish) as previous studies did [16,30]. The study was carried out during the spring season since the plankton peak is broader [49], and fish activity is higher [8,50]. We expected that taxonomic diversity for both phyto- and zooplankton would be more related to abiotic factors (e.g., conductivity, depth, area, etc.) [19,23] whereas size diversity would be mainly related to biotic interactions (top-down and bottom-up processes). Regarding zooplankton size diversity, we expected that it would mainly be driven by top-down effects (i.e., fish predation) due to size-based predation pressure by fish. At increasing predator size diversity, top-down control would be enhanced [30,51] and, therefore, we predicted a negative relationship between the zooplankton and fish size diversity. In the case of phytoplankton, we expected that the size diversity would be more affected by bottom-up effects (i.e., changes in resource availability) rather than by top-down effects (i.e., zooplankton predation), as it has been observed in previous studies [16,31,52], and therefore we predicted a negative relationship between the size diversity of phytoplankton and its resource, due to competitive interactions for resources in less productive systems [16].

2. Materials and Methods

2.1. Study Area

The study was performed in brackish ponds located in two protected areas in the Empordà coastal wetlands (Figure 1) that are situated between 42°01′42″ N–3°11′18″ E and 42°15′58″ N–3°08′17″ E of Catalonia region (NE Spain) and are characterized by Mediterranean climate with hot, dry summers and mild, wet winters. The hydrological pattern of these coastal wetlands, free from tidal influence, is mainly determined by sudden and irregular flooding events (rainfalls or marine intrusion during sea storms, the latter not more often than twice a year) followed by prolonged periods of confinement without water inputs (winter anticyclone and summer droughts) [53]. Therefore, hydrological connection among ponds and/or to the sea or rivers takes place only during the strong flooding events. During the rest of the hydroperiod, ponds remain confined because of the presence of dunes, sand bars, and/or artificial levees [12,53]. In contrast to many temperate lakes, in these confined Mediterranean ponds nitrogen rather than phosphorus is usually the limiting nutrient for primary production (e.g., [19,54,55]) due to strong denitrification processes, low water turnover and the high internal load of phosphorous [56,57].

2.2. Sampling and Analysis

A total of 13 permanent ponds (Table S1) were sampled once during the spring season (May to early June 2016). We chose to sample the ponds in spring since the plankton peak is broader [49], and fish activity is higher [8,50]. Mean water column depth (cm) was calculated from repeated measures of water column depth obtained in situ using a two-meter rule. Total area (m²) was estimated in each pond by using ‘Google Maps Area Calculator Tool’ [58]. Physicochemical variables such as temperature (°C), pH, conductivity (mS·cm⁻¹), total dissolved solids (mg·L⁻¹) and dissolved oxygen (mg·L⁻¹) were measured in situ using a multiparameter probe (Hanna Instruments, Woonsocket, RI, USA). Water transparency was estimated as Secchi depth (cm) out of maximum water column depth (cm) as it has been previously used in shallow waterbodies [16]. Water samples were analyzed for total phosphorus (mg·L⁻¹), soluble reactive phosphates (mg·L⁻¹), total nitrogen (mg·L⁻¹) and nitrates (mg·L⁻¹), according to Koroleff 1973 [59], adapted by Seal Analytical to an integrated system of CFA digester. Chlorophyll-a was measured using spectrophotometry after 90% acetone extraction following Parsons and Strickland, 1963 [60].

In each pond, water samples for planktonic organisms were taken through the water column by means of a 6 L container. Several subsamples were obtained in each pond from different sites and subsequently mixed to overcome the expected patchy distribution of plankton. To obtain a phytoplankton sample, 250 mL of unfiltered water were stored in 4% acid Lugol’s solution, and for the zooplankton sample 5 L were filtered through a 50 µm mesh-size net and preserved in 4% acid Lugol’s solution. Zooplankton individuals (including rotifers, copepods and cladocerans) were counted, identified to species level (whenever possible) and measured (total length in µm) using a stereoscope and an inverted microscope (Utermöhl method). For the individual counting and identification, we analyzed the whole sample, whereas for the measuring, we measured the first 100 individuals (when possible) assuming all individuals were equitably distributed in the observed sample. Phytoplankton was counted and identified to specie level under inverted microscope using Utermöhl chambers following the protocol for phytoplankton identification described in the EU project ‘WISER’ (Water bodies in Europe: Integrative Systems to assess Ecological status and Recovery) [61] When identification at species level was not possible, the different species of the same genus were numbered differently in order to respect the diversification. For both zooplankton and phytoplankton, at least 100 individuals of the most abundant taxa were counted. For phytoplankton, biovolume was estimated through the calculation of similar geometric models [62].

Fish Sampling

Different methods of fish sampling were selected to cover efficiently the different microhabitats of ponds. Shallow littoral areas were sampled using plastic minnow traps (PMT) and fyke nets (FN). Multi-mesh gillnets (GN) were set on central and deepest areas of ponds that were at least 1.5 m deep since it is the height of the gillnet. PMT were made with 2 L soda plastic bottles as described in Fouilland & Fossati, 1996 [63] and Clavero et al., 2006 [64]. The upper piece of each bottle was cut and inverted, acting as a funnel (main characteristics: 21.5 mm of funnel diameter, 72 cm² of interception area, 22.5 cm of length and 9.6 cm of height). FN consisted of a semicircular entrance ring followed by three smaller circular rings surrounded by a net (3.5 mm mesh) and had two consecutive funnels (120 mm of funnel diameter, 1050 cm² of interception area, 98 cm of length, 30 cm of height and 95 cm of wing length). GN followed the European standard [65] composed of 12 mesh-size panels of length 2.5 m each ranging between 5 mm to 55 mm (length 30 m, height 1.5 m, knot to knot dimensions following a geometric series: 5, 6.25, 8, 10, 12.5, 15.5, 19.5, 24, 29, 35, 43, and 55 mm). PMT and FN were set during 24 h. GN were set for approximately 12 h overnight to include both evening and morning phases of high fish activity [65]. The total number of fish traps set in each of the studied ponds varied according to their area and depth (Table S1). Captured fish were identified, measured for total length (mm) and released. In each studied pond, captures per effort unit (CPUE) were calculated for each fish species and each sampling method by dividing the captures for the number of traps, fyke nets or gillnets respectively.

2.3. Calculation of Diversity Metrics

For each phytoplankton, zooplankton and fish samples we obtained two diversity measures: the Shannondiversity index (taxonomic) and the size diversity index (non-taxonomic). The taxonomic index (H′) was calculated using the numerical abundance of each identified taxon as following [66]:

$$H\prime =-{\displaystyle \sum }_{i=1}^{S}pi\ln pi$$

(1)

where pi is the proportion of individuals belonging to the ith taxon, and S is the total number of identified taxon. Shannon diversity index (H′) was calculated with the ‘vegan’ R-package [67].

The size diversity index (µ) was calculated using individual size measurements as proposed by Quintana et al., 2008 [22]. This index is computed based on the Shannon diversity expression adapted for a continuous variable, such as body size. This measure takes the form of an integral involving the probability density function of the size of the individuals described by the equation

$$\mu =-\underset{0}{\overset{\infty }{{\displaystyle \int }}}{\rho }_x\left(x\right){\log }_{2 }{\rho }_x\left(x\right)dx$$

(2)

where ${\rho }_x$(x) is the probability density function of size x. The nonparametric Kernel estimation was used as a probability density function, which is applicable to any type of size distribution. Before computing size diversity, data were automatically standardized by division of each size value by the geometric mean of the size distribution. The size diversity index is the continuous analogue of the taxonomic Shannon diversity index, and it produces values in a similar range to those of the Shannon index. For each trophic level, we randomly measured at least 100 individuals in each sample thatrepresents a size diversity error estimation lower than 10%. Size diversity index was computed following an R code provided by Quintana et al., 2008 [22].

2.4. Data Analysis

We used general linear models (GLMs) to test the effects of biotic interactions (e.g., top-down and bottom-up effects) and abiotic factors on the phyto- and zooplankton taxonomic and size diversity indexes. Thus, at each trophic level (phytoplankton and zooplankton) two GLM models were carried out considering taxonomic and size diversity indexes as response variables respectively. Since the small sample size (N = 13) limited the number of independent variables to be included in each GLM model, we used the smallest number of predictors which allowed us to test for our hypothesis: two variables reflecting biotic interactions (top-down and bottom-up effects, respectively) and three abiotic environmental variables.

As predictor variables reflecting biotic interactions, we used the two diversity indexes of the adjacent trophic levels, taxonomic or size-based depending on the response variable (i.e., size-based metrics when size diversity was the response variable and vice versa). Therefore, predictor variables used to test top-down effects (i.e., predation from upper trophic level) in phytoplankton GLMs were the zooplankton diversity indexes, and in zooplankton GLMs were the fish diversity indexes. For what concern the fish community, we only considered those planktivorous fish species since we are testing predation effects on zooplankton. We only used captures from PMT since they are more effective than FN in capturing small juveniles [64]. To test for bottom-up effects (i.e., food resource availability), predictor variables used in zooplankton GLMs were the phytoplankton diversity metrics whereas in phytoplankton GLMs we used the “phytoplankton biomass: total nitrogen” ratio. Similar ratios (e.g., Chlorophyll-a: TP, Chlorophyll-a: TN) have been previously used in studies dealing with phytoplankton yield, as a measure of the phytoplankton resources limitation (ratio values are high under resource limitation) [16,26,68].

As predictor abiotic variables for both phyto- and zooplankton GLMs, we considered conductivity, pond area, and depth because of the absence of correlation among them and with the diversity metrics (Pearson’s r < |0.31| and p-value > 0.30 for all cases. See Supplementary Table S2 for correlations details). These three abiotic variables have already been shown to play an important role in structuring Mediterranean pond communities [15,16,31,32]. We dismissed using other abiotic variables that, although not correlated, did not show enough variation among ponds (e.g., temperature, pH). The rest of the measured environmental variables were significantly correlated among them (Pearson’s r > |0.56| and p-value < 0.05 for all cases).

We looked for the most parsimonious of the full models by performing an automatic backward selection of one predictor variable at a time by minimizing the Akaike information criteria (AIC). The most parsimonious model was the combination of variables having the strongest impact on outcomes. To compare the relative strength of the significant predictors, we calculated their standardized (beta) coefficients, and adjusted R² was used as a measure of the variability explained by the model. Variables were log-transformed when necessary to improve linearity and reduce heteroscedasticity, and residual plots were inspected to detect violations of regression assumptions. The variance inflation factor (VIF) was used to check collinearity, and Cook’s distances were analyzed to check for outliers and influential cases. We are aware that a low number of observations (N = 13) prevent making reliable predictions when performing multiple regression models taking the risk to incur into Type I error. Nevertheless, to make an estimation of the statistical power of the GLM models obtained, we computed the statistical power analysis of each model using the G*Power software [69]

All data have been analyzed using R version 3.4.2 (R core Team, 2017, Boston, MA, USA). All plots have been created with ‘ggplot2’ R-package [70].

3. Results

Physico-chemical variables measured in the studied ponds are shown in

Table 1. Mean, minimum, maximum, and standard deviations (SD) values of the physico-chemical and morphometric variables measured in the study ponds (N = 13). The “phytoplankton biomass: total nitrogen” ratio used in phytoplankton GLMs models to test the bottom-up effects (i.e., food resource availability, see text Section 2.4) is also included.

	Mean	Minimum	Maximum	SD
Mean water column depth (cm)	55.08	16.00	150.00	37.99
Pond area (m2)	9457.90	147.90	68,150.00	18,162.80
Water transparency	0.89	0.64	1.00	0.12
Temperature (°C)	25.56	21.14	32.20	2.97
pH	8.47	7.39	9.97	0.72
Conductivity (mS·cm−1)	42.61	10.07	69.10	18.74
Total dissolved solids (g·L−1)	0.03	0.01	0.05	0.00
Dissolved oxygen (%)	115.04	41.10	262.50	56.66
Chlorophyll-a (µg·L−1)	11.92	0.31	43.41	14.54
Total nitrogen (μmol·L−1)	97.00	58.55	234.40	47.46
Total phosphorus (μmol·L−1)	3.19	0.44	7.57	2.15
Nitrate (μmol·L−1)	0.21	0.07	0.66	0.16
Soluble reactive phosphate (μmol·L−1)	1.27	0.09	6.94	2.01
Phytoplankton biomass: total nitrogen	0.35	0.09	1.58	0.42

. Since we sampled once during the spring season, water temperatures (°C) were relatively similar between ponds (with a mean value of 25.56 °C and a standard deviation of 2.97), as well as pH. Studied ponds are brackish, with conductivity values ranging from 10.7 mS·cm⁻¹ to 69.10 mS·cm⁻¹, and shallow, with mean depths ranging from 16 cm up to 150 cm. In addition, they differed quite a lot also in their areas, with values ranging from 147.90 m² to 68,150 m². The measured dissolved oxygen content in water showed a situation very far from anoxia, with values from 41.10 mg·L⁻¹ to 262.50 mg·L⁻¹. Chlorophyll-a values ranged from a minimum value of 0.31 μg·L⁻¹ to a maximum value of 43.41 μg·L⁻¹, and the mean water transparency, measured as Secchi depth: maximum depth, was 0.89 cm. With regard to nutrient concentrations, inorganic forms of nitrogen and phosphorus (nitrate and soluble reactive phosphate) were always lower than total nitrogen and phosphorus, which included organic and inorganic forms (Table 1). Nitrates were especially lower in the studied ponds with an average concentration of 0.21 μmol·L⁻¹. Values of the ratio “phytoplankton biomass: total nitrogen”, as a measure of resource limitation for the phytoplankton, ranged from 0.09 (high resource availability) to 1.58 (low resource availability).

With regards the phytoplankton community in the studied ponds, a total of 39 species were identified belonging to seven classes (Figure 2A). Prasinophyceans (mostly marine species) dominated the community (>60%) under high conductivity (ponds PI: ‘Bassa del Pi’, LLU: ‘Bassa de la llúdriga’ and AN: ‘Bassa de l’anguila’), whereas bacillariophyceans (diatoms) were more abundant in the 3 ponds that showed lower conductivities, with the specie Navicula sp. (Bory de Saint-Vincent, 1822), except for FAN (Fangassos) pond where diatoms were probably composed by marine species. In the rest of the ponds, phytoplankton was dominated by cryptophyceans and dinophyceans (dinoflagellates). Phytoplankton individual length ranged from 2.55 µm (cyanophiceans) to 112.81 µm (prasinophyceans).

Regarding the zooplankton community, a total of 17 species were found in the studied ponds being rotifers the dominant group in 8 of the ponds (>50% of total zooplankton abundance, Figure 2B). In the other five ponds, the zooplankton was dominated (>70% to 100%) by the group of copepods, where planktonic larval stages (copepodites + nauplii) were the dominant followed by adult individuals of calanoids and harpacticoids. Cladocerans were the less abundant and frequent (present only in two ponds) and did not dominate the community in any pond. Zooplankton individual length ranged from 0.39 mm (rotifers) to a maximum of 5.54 mm (cladocerans).

For what concern the fish community, we found 8 species including both continental and marine species, being the latter the less abundant. The endemic Aphanius iberus (Valenciennes, 1846), an omnivorous an euryhaline species, was well distributed along the conductivity gradient and was the most abundant species reaching a maximum of 957.6 CPUE in all the ponds together. It was present in 11 out of 13 ponds and dominated the community (>70%) in 8 ponds (Figure 2C). The second most abundant fish was the freshwater, invasive planktivorous Gambusia holbrooki (Girard, 1859) that reached a maximum total abundance of 705.9 CPUE in all the ponds together. It was only present in the less salty ponds (Figure 2C) except in the pond BF (Bassa Fartet) where the species never arrived. In order of abundance, the following species were the planktivorous Atherina boyeri (Risso, 1810), with a maximum CPUE of 102.9, the benthivorous Pomatoschistus sp. (Gill, 1863) with a maximum CPUE of 62.2 and Mugil cephalus (Linnaeus, 1758) with a CPUE of 41.8, all them present in five ponds at most. The less abundant (≤4.8 CPUE) and frequent (<2 ponds) species were, in order, the marine fish Solea solea (Queusel, 1806), Anguila anguila (Linnaeus, 1758) and the invasive Lepomis gibbosus (Linnaeus, 1758). Fish individual length ranged from 6 mm (A. iberus juveniles) to a maximum of 60 cm (A. anguila).

Diversity metrics calculated for the three trophic levels (phytoplankton, zooplankton, and fish) are shown in

Table 2. Diversity metrics of phytoplankton, zooplankton and fish assemblages in the studied ponds (N = 13). The descriptive statistics are the mean, minimum, maximum and standard deviation (SD). Fish diversity metrics were computed using only those planktivorous species (A. iberus, G. holbrooki, A. boyeri).

	Mean	Minimum	Maximum	SD
Phytoplankton
Size diversity	1.79	0.31	2.95	0.61
Taxonomic diversity	0.94	0.02	2.14	0.53
Zooplankton
Size diversity	1.95	0.70	2.84	0.67
Taxonomic diversity	0.38	0.00	1.13	0.34
Fish
Size diversity	1.42	0.72	2.08	0.37
Taxonomic diversity	0.08	0.00	0.40	0.14

. Size diversity values ranged from a minimum of 0.31 up to a maximum of 2.95, both found in phytoplankton assemblages, while taxonomic diversity ranged from 0.00 (in those zooplankton and fish assemblages composed by one species), to a maximum of 2.14 (in phytoplankton assemblages).

Variables significantly affecting size and taxonomic diversity in each of the GLM models performed are shown in

Table 3. Results of the GLMs (N = 13) showing the predictor variables that affect size diversity and taxonomic diversity of phytoplankton and zooplankton assemblages. Only the most parsimonious significant models are presented. For each one, intercept (estimate and standard error, S.E.), Beta coefficients (standardized), t-value, significance (p-value) and degrees of freedom (df) are shown. For each model, the global p-value, R² and the statistical power are shown. Phytoplankton biomass: TN is the ratio of “phytoplankton biomass: total nitrogen” used as a measure of food resource availability for phytoplankton.

Response Variable	Predictor	Estimate	S.E.	Beta Coefficients	t-Value	p-Value	df	Global p-Value	R2	Statistical Power
Phytoplankton
Size diversity	Conductivity	0.02	0.01	0.48	2.04	0.06	11	0.04	0.47	0.55
	Log Phytoplankton biomass:TN	−3.71	1.42	−0.62	−2.62	0.02
Species diversity	Log Pond Area	0.32	0.11	0.46	2.96	0.01	11	<0.01	0.76	0.99
	log Mean Depth	−1.44	0.31	−0.73	−4.73	<0.01
Zooplankton
Species diversity	Log Pond area	0.29	0.10	0.65	2.86	0.02	11	0.02	0.42	0.8

, whereas full models are shown in Supplementary Table S3. According to GLMs results (Table 3, Figure 3 and Figure 4), in the phytoplankton assemblages, size diversity was slightly positively related with conductivity (marginally significant p-value = 0.06), and negatively related with the ratio “phytoplankton biomass: total nitrogen” (p-value = 0.02) (Figure 3). Beta coefficients showed values of 0.48 for conductivity and −0.62 for the ratio meaning that the latter had a stronger effect on the dependent variable. This model explained 47% of the variation in the size diversity (p-value = 0.04, statistical power = 0.55). Taxonomic diversity in the phytoplankton was positively related with area and negatively related with depth of the pond (Figure 4). In this second model, beta coefficients showed values of 0.46 for pond area and −0.73 for the pond depth showing a stronger effect of pond depth. This model explained 76% of the variance in the phytoplankton taxonomic diversity (p-value < 0.01, statistical power = 0.99). Regarding zooplankton assemblages, none of predictor variables tested were selected in the GLMs to be significant for size diversity. Instead, taxonomic diversity was positively related with the area of the pond (Figure 4). A Beta coefficient of 0.65 showed a high strength of the effect of pond area on the dependent variable. The model explained 42% of variance (p-value = 0.02, statistical power = 0.80).

4. Discussion

In the studied ponds, the physico-chemical characteristics of the water, as well as the composition of planktonic communities, were the expected for Mediterranean coastal brackish ponds during spring season. Thus, they show a variable conductivity (ranges within expected for brackish waters) and low concentrations of inorganic nutrients (nitrates and phosphates) due to the absence of water inputs during the sampling period. Inorganic nutrients are quickly up taken by phytoplankton as they enter into the pond, and then transformed into organic form [55,57]. Therefore, high concentrations of total nutrients (both nitrogen and phosphorus) are typical in these confined ecosystems, which even tend to increase in spring and summer due to evaporation phenomena [53]. Regarding biological communities, the phytoplankton species found were also characteristic of the brackish coastal ponds during spring season, with a high presence of marine species (dinophyceans) that can stand high salinities and entered the ponds during sporadic flooding events (e.g., sea storms) [39,53]. Zooplankton composition was also similar to the composition found in spring in previous studies in brackish ecosystems [12,57] with the dominance of salt tolerance rotifer species and planktonic larval stages of copepods. Regarding fish community, marine species were found in low abundance in ponds with higher salinity, due to sporadic flooding events, whereas small planktivorous fish (A. iberus and G. holbrooki) were the dominant and more abundant species in all the studied ponds. Ranges of taxonomic and size diversity values found in each trophic level were also characteristic of the spring season in Mediterranean brackish ponds [71,72].

Regarding the effects of abiotic and biotic factors on plankton diversity indexes, our results agreed with our hypothesis regarding taxonomic diversity since it was significantly affected by an abiotic factor (pond morphometry) for both phytoplankton and zooplankton assemblages. Thus, in the studied Mediterranean brackish ponds, our results showed that when increasing pond area both the phytoplankton and zooplankton taxonomic diversity also increased. Phytoplankton taxonomic diversity increased also when decreasing pond depth. These results are supported by the biogeographical principle that a larger area supports more species, [73], and that area is one of the best single predictors of species number [74]. Actually, one of the most commonly observed patterns in aquatic ecology, and for a wide variate of taxa (planktonic and benthonic organisms), is the positive relationship between the diversity metrics and the area of the water body [44,48,73,75,76]. According to those studies, in a bigger pond there is likely more habitat heterogeneity, favoring the coexistence of more species. In the present study, a bigger area was translated to a more diverse littoral and aquatic habitats (e.g., presence of Phragmites australis, Arthrocnemum, Juncus maritimus, Ruppia sp., algal mats, unvegetated open waters, etc.) In a study on biodiversity in Danish lakes [37], lake depth was also related to phyto- and zooplankton richness.

Even though conductivity varied among the studied ponds during the spring sampling (from 10.07 to 69.10 mS·cm⁻¹), probably due to different evaporation rates (different surface, wind influence) [53,57] or water inputs during previous flooding events (fresh- or seawater inputs) [53,77], we did not find a significant effect of conductivity on phyto- and zooplankton taxonomic diversity. This may suggest that conductivity values did not differ enough among ponds to have an appreciable effect on the taxonomic diversity of plankton. Conductivity has been widely recorded as a main environmental factor shaping taxonomic structure and biodiversity in Mediterranean planktonic communities [6,20] having negative consequences in species composition, such as a reduction in species richness and diversity, and also affecting the food web interactions [26,27].

Concerning the effects of the analyzed environmental factors on the plankton size diversity in the spring season, results obtained in the present study partially agreed with our hypothesis that size diversity (in contrast to taxonomic diversity) was more related to biotic interactions (predation and/or food resource availability). However, this was only true for phytoplankton assemblages since for zooplankton none of the abiotic nor biotic predictor variables tested significantly affected the size diversity. As we predicted, our results showed that phytoplankton size diversity was mainly affected by bottom-up effects (i.e., food resource availability) rather than top-down effects (i.e., predation), as it has been found in previous studies in lakes [16,28,31] and marine ecosystems [78,79,80]. In the present study, resource limitation (i.e., high values of the “phytoplankton biomass: total nitrogen” ratio) would lead to a dominance of small phytoplankton cells and, therefore, to a decrease in size diversity under low nutrient availability. Small phytoplankton is more efficient than the large one due to their low resource requirements and high cell’s surface-to-volume ratio [81,82]. Our results agree with previous studies in lake and marine ecosystems where changes in phytoplankton resource supplies, like total phosphorous [16] and a resource supply index [34], were the main drivers of phytoplankton size structure since smaller sizes are favored at low nutrient availability due to the higher surface: volume ratio or lower resource requirement. Nevertheless, there are some studies which found both positive [83] and negative [30] relationships between zooplankton and phytoplankton size diversities thus indicating a top-down control on the phytoplankton size structure. In the case of positive relationships, prey size diversity was found to increase predator size diversity, promoting diversity of consumers [83] whereas in the case of negative relationships, an increase in predator size diversity was found to enhance the strength of top-down control, reducing prey size diversity [30]. It is worth mentioning that, in the present study, a slightly positive effect of conductivity on phytoplankton size diversity was also observed. The additional presence of phytoplankton marine species in high-conductivity ponds (i.e., prasinophyceans and some bacillariophyceans) likely increased the phytoplankton size diversity. The entrance of marine phytoplankton in salt marshes was also observed in López-Flores et al., 2006 [84] after a sea storm, leading to a change in the specie composition. In the present study, the pond with the less size diversity was dominated by small-sized species (e.g., cryptophyceans) suggesting that the prevalence of small-sized individuals leads to a low size diversity.

Against our predictions, zooplankton size diversity in these Mediterranean ponds during the spring season was not affected by fish size diversity (i.e., top-down effects) although previous studies have found significant relationships between predators and prey size diversity (e.g., zooplankton and planktivorous fish [16], planktivorous and piscivorous fish [85]). It is important to highlight that in these previous studies a wide geographic range was considered (including different ecoregions and altitudes), and predators size diversity ranges were wider than the ones found in the present study. For example, in Brucet et al., 2017 [16] fish size diversity ranged from a minimum of −0.81 to a maximum of 2.42, whereas in our study, where geographical variation was negligible, size diversity of planktivorous fish ranged from 0.72 to 2.08. This was probably due to the similar sizes of the planktivorous fish present in the studied ponds (A. iberus, G. holbrooki, and A. boyeri, size ranges between 6 mm and 5.8 cm) that lead to a small range of size diversity values. This, together with the fact of not including in our analysis the presence of other predators, that although they were observed in situ (e.g., jellyfish Odessia maeotica, the amphipod Gammarus aequicauda) they could not be sampled properly, may prevent finding significant top-down effects on zooplankton size structure. We cannot ignore the fact that in the Mediterranean brackish ponds both the jellyfish and the amphipod have been found to act as planktonic top predators [86,87] causing changes in lower trophic levels through cascading effects [88,89].

In conclusion, our results suggested that in Mediterranean brackish ponds during spring season both the taxonomic diversity of phyto- and zooplankton were mainly influenced by abiotic factors (pond morphometry variables), whereas size diversity was mainly affected by biotic interactions but just in the case of phytoplankton. Nutrient limitation (i.e., bottom-up effects) rather than zooplankton predation (i.e., top-down effects) affected the phytoplankton size diversity leading to small phytoplankton size diversity when resource availability was low. Regarding zooplankton size diversity we did not find the expected top-down effects (i.e., fish predation) probably due to the similar sizes of planktivorous fish lead to similar fish size diversities among ponds. We are aware that our findings may have been limited by the low number of ponds studied, but also by the fact of considering the classical three-level food chain (phytoplankton, zooplankton and fish) neglecting the effects of omnivory and intraguild predation [28,90,91], and the presence of the microbial loop [92,93]. It is worth mentioning that our results are limited to spring plankton communities in Mediterranean brackish ponds and we cannot rule out the effect of seasonality in our results since, in the Mediterranean region, it leads to greater environmental variability which is also reflected in planktonic community dynamics and composition [26,53,57].