Zooplankton Community Dynamics in Temporary Mediterranean Wetlands: Which Drivers Are Controlling the Seasonal Species Replacement?

Abstract

Temporary Mediterranean wetlands are characterized by both intra and interannual variations in their environmental conditions. These inherent fluctuations in limnological features affect the seasonal variation in the structure and dynamics of the aquatic communities. In this study, we hypothesized that zooplankton community is coupled to seasonal changes of the environmental variables along the hydroperiod. To get this purpose, the study was focused in monitoring, by collecting monthly samples during an annual period, seven temporary Mediterranean ponds located in the south-eastern region of the Iberian Peninsula (Alto Guadalquivir region, Andalusia). The relationships between zooplankton community and the different limnological variables were analyzed based on two approaches: a Spearman correlation analysis and a correspondence canonical analysis (CCA). The results have shown that chlorophyll-a concentration, Secchi depth, total nitrogen concentration, wetland area and depth were the variables with a greater influence on the zooplankton community, explaining the zooplankton species replacement. Moreover, optima and tolerance of the zooplankton species were obtained from the position of species within CCA diagram, allowing the separation of different groups of zooplankton along the hydroperiod. We finally highlight that the monitoring of zooplankton community and environmental conditions are necessary to evaluate how theses singular and endangered aquatic ecosystems will be affected by anthropogenic activities in the future.

1. Introduction

Mediterranean ecosystems drastically differ from the colder-temperate aquatic ecosystems of Northern European in both structure and function [1]. Mediterranean wetlands are characterized by uneven hydroperiods with high intra and interannual water level variations, which lead to spatial and temporal fluctuations. Consequently, one of the most definitory aspects of these ecosystems is the temporality of many of them, which normally begin to fill in the autumn and dry up in the summer. This also conditions in many cases intra-annual differences in salinity values, so that they contain waters with lower salinity at the beginning of the hydrological cycle and more saline waters at the end of the hydroperiod [1]. This situation is responsible for changes in their communities and ecological processes, making them ecologically rich ecosystems supporting a high diversity [2]. This peculiar feature has increased the interest in Mediterranean wetlands, and a considerable number of studies have been focused on processes and the structure of their communities [3,4,5]. Despite this, more studies are still required to know the structure and dynamics of aquatic communities that temporary Mediterranean wetlands support, which actually are threatened by habitat degradation [6,7,8]. This concern has induced a change in the direction of research in the European framework trying to ensure the conservation of these ecosystems, with the implementation of diverse management and protection directives (Habitat Directive Natura 2000 network, Water Framework Directive 2000/60/EC).

Within aquatic ecosystems, plankton communities are an essential element, playing an important role in nutrient recycling and energy flow [9]. Among plankton, zooplankton is recognized as a fundamental component in Mediterranean wetlands [10]. Zooplankton are excellent components for monitoring aquatic environments [11,12,13] due to their small size, short generation time, high sensitivity to environmental changes and low-cost sample collection.

Multiple studies, in a variety of shallow Mediterranean wetlands, have recognized that temporal and spatial variability of these systems can significantly influence the structure and composition of these aquatic invertebrates [2,3,14]. In fact, zooplankton composition and abundance are controlled by many factors such as water hydrochemistry, season, lake morphology, presence of macrophytes, predators and lake production [15]. However, the seasonal relationships between environmental factors and zooplankton assemblages are neglected or only beginning to be explored.

In this study, and considering the seasonal environmental variability that characterize Mediterranean wetlands, we hypothesized that zooplankton community is coupled to seasonal changes of the environmental variables along the hydroperiod. To test the hypothesis, we identify the most important variables in the dynamics of these Mediterranean aquatic ecosystems and the influence in the structure of zooplankton communities (cladocerans and copepods) during a period of flooding and desiccation.

2. Materials and Methods

2.1. Study Area and Sample Collection

The study was carried out in seven endorheic ponds located in the South-East of the Iberian Peninsula (Figure 1). These ecosystems are all of them small size, shallow, fishless and temporary, with drastic water level fluctuations in relation to rainfall. Studied ponds present some differences such as the geology and the percentage of olive groves in their watersheds. Three of the ponds (Quinta, Casillas and Hituelo) are located in the Guadalquivir valley and with a high percentage of olives tree (>80%) in their catchment area. Three other ponds (Ardal, Santisteban and Castillo) are located on siliceous materials, while the other pond (Orcera) is located on calcareous materials. Of these four ponds, Castillo is located on an intermediate agricultural land (between 40% and 60% of olive groves) and finally the other three ponds are located in lands with no or low presence of olive groves (less than 10%).

Samples were collected during an annual period, from December 2009 to December 2010. All studied ponds were dry in late spring or summer 2009, filling in autumn-winter 2009–2010. Indeed, all ponds (except Castillo) were dry at the beginning of the sample period. Physicochemical variables and biological samples were taken monthly in each pond. Pond surface area at each sampling date was calculated using measurements of the major and minor axes fitted to simple shapes (circle or ellipse). The hydroperiod length (days of pond inundation—DWI) was measured from the first time when the pond basin was waterlogged.

Water depth at the deepest point was measured in situ in each sampling date with a ruler. Temperature, pH, oxygen and conductivity were determined in situ with a multiparameter probe (YSI 556 MPS). Turbidity was measured in situ as nephelometric turbidity units (NTU) with a Hanna HI93703 probe. Water transparency was determined using a Secchi disk. Many times, and due to the shallowness of ponds, the bottom could be seen from the surface and then Secchi disc depth could not be measured. To tackle this problem, we used the index of the transparency [16], which permit us to obtain the percentage of Secchi disc depth with respect to pond depth. Thus, if Secchi depth is equal to the depth (which means that the bottom could be observed from the surface), the value of the index is 100%; lower values of percentage mean lower transparency and difficulties for light to reach the bottom. Chlorophyll-a concentration (Chl-a) was obtained from fluorescence measurements taken in situ with a field fluorometer Aquafluor Turner Design Handheld, previously calibrated by fluorometry. The percentage covered by submersed macrophytes (SM—charophytes and vascular plants) was estimated as a percentage of pond surface based on visual site mapping.

In the laboratory, unfiltered integrated water samples were digested with potassium persulfate for analysis of total P (TP) and with sulfuric acid for analysis of total N (TN). TP and TN were spectrophotometrically quantified by using the molybdenum-blue method [17] and the ultraviolet spectrophotometric screening [18], respectively. Detection limits were 1 μg P L⁻¹ (TP) and 40 μg N L⁻¹ (TN). Organic matter (OM) was estimated by the loss on ignition at 450 °C during 3 h.

Zooplankton samples were taken along two longitudinal transects, one from the shoreline to the open waters and another along the shoreline, with a plankton net with a mesh size of 63 µm and a mouth diameter of 27.5 cm. Both samples were combined to obtain an integrated sample and fixed in situ with 4% formalin (final concentration). In each sampling time two zooplankton samples were collected, one for estimating the abundance of copepods (calanoids and cyclopoids) and branchiopods, and another one for measuring total biomass. Once in the laboratory, zooplankton taxa were firstly identified to species level using taxonomic keys [19,20,21], and later abundance data was also obtained. As a result, for each zooplankton sample, the following parameters were calculated: (i) community parameters: species richness (R), Shannon–Wiener diversity index (H’), Pielou evenness index (J) and Simpson dominance index (D); (ii) ratio calanoid copepods abundance (Cal) versus cyclopoid copepods (Cycl) plus cladocerans abundance (Cla—[22]); (iii) ratio of abundance of large cladocerans (Daphnia and Simocephalus) versus the total abundance of cladocerans [10]; (iv) total zooplankton biomass; obtained from the dry weight values at 105 °C. Community parameters values were obtained using the software PAST 2.8b.

2.2. Statistical Analysis

Zooplankton assemblages were analyzed based on three different approaches: (i) a permutational analysis of variance (PERMANOVA), performed with 9999 permutations on the Bray–Curtis measure to test differences in the environmental variables and zooplankton communities across ponds and date, (ii) a Spearman correlation analysis between the four arrays above mentioned (community parameters, different ratios between zooplankton groups, zooplankton biomass and environmental variables) and (iii) a canonical correspondence analysis (CCA), a multivariate direct gradient analysis commonly used in ecology to obtain a linear combination of environmental variables assuming an unimodal species–environment relationship [23,24]. CCA was selected because a preliminary detrended correspondence analysis (DCA) indicated that the data set exhibited a length of gradient >3 [25].

Later, the results of CCA were plotted using the obtained linear combination scores (LC scores [26]). In the ordination graph only the first two axes were used because they showed the highest eigenvalues and amount of variation by CCA ordination. A Monte Carlo permutation test (999 permutations) was used to assess the statistical significance of the ordination axes and of the species–environment correlations (Intraset correlations; [27]).

In turn, the position of species within the CCA diagram indicates approximately their environmental optima [27]. In this sense, the CCA is a method based on weighted averaging (WA) where the optimum for an environmental factor of a species is estimated as:

$$u_k=\frac{\sum y_{ik}{x}_i}{\sum y_{ki}}$$

(1)

and the tolerance of a species as:

$$t_k={(\frac{\sum y_{ik}{(x_i-u_k)}^2}{\sum y_{ki}})}^{1/2}$$

(2)

where u_k is the optimum of species k; y_ki is the abundance of species k at site i; x_i is the value of the environmental factor at site i.

Prior to analyses, the data (species abundances and environmental variables) were tested to normality with a Kolmogorov–Smirnov test. When the assumption of normality was not fulfilled, a log10 (× + 1) transformation was used to reduce the influence of the most dominant species and to remove the influence of different scales of measured variables. Rare species (<10% of abundance in the samples) were removed from the analysis. Spearman correlation analysis was performed by using the software Statistica and CCA with PC-ORD.

3. Results

Although in general, strong differences in their environmental conditions among the seven studied ponds emerged (see

Table 1. Mean and standard deviation of limnological variables measured in the studied ponds.

	Castillo	Orcera	Ardal	Santisteban	Hituelo	Quinta	Casillas
Conductivity (µS cm−1)	144 ± 55	444 ± 126	116 ± 52	454 ± 128	343 ± 74	561 ± 199	236 ± 44
pH	8.3 ± 0.6	8.1 ± 0.5	8.1 ± 0.8	7.8 ± 0.2	8.1 ± 0.4	8.5 ± 0.7	8.5 ± 0.6
Turbidity (NTU)	28 ± 67	6 ± 4	31 ± 25	71 ± 96	20 ± 24	11 ± 6	32 ± 44
Chl-a (µg L−1)	8.91 ± 7.46	7.64 ± 4.57	7.23 ± 3.55	8.55 ± 2.97	7.76 ± 10.08	8.82 ± 8.56	6.98 ± 5.06
SM (%)	80 ± 19	88 ± 14	65 ± 17	56 ± 27	54 ± 14	66 ± 23	64 ± 28
Area (ha)	1.94 ± 1.24	1.68 ± 1.24	0.14 ± 0.09	5.10 ± 5.20	16.82 ± 4.82	22.81 ± 3.45	3.62 ± 1.14
Depth (cm)	148 ± 61	172 ± 78	29 ± 9	27 ± 17	179 ± 52	20 0 ± 51	242 ± 87
Secchi depth (cm)	94 ± 54	172 ± 78	29 ± 9	15 ± 6	111 ± 65	146 ± 69	113 ± 76
Temperature (°C)	17.1 ± 8.2	14.4 ± 7.6	16.3 ± 6.5	14.7 ± 5.9	20.5 ± 7.8	20.4 ± 7.0	17.8 ± 7.1
O2 (%)	87.8 ± 17.9	77.6 ± 17.1	116.4 ± 32.8	74.5 ± 14.4	92.0 ± 10.0	82.7 ± 15.6	82.6 ± 21.7
OM (mg L−1)	3.9 ± 2.8	3.5 ± 2.5	8.2 ± 5.4	8.5 ± 3.7	4.2 ± 2.0	6.1 ± 4.9	3.8 ± 3.9
TP (µg L−1)	57.8 ± 44.1	35.4 ± 23.5	150.9 ± 87.0	122.3 ± 80.4	55.5 ± 33.6	49.1 ± 22.1	206.7 ± 155.3
TN (mg L−1)	1.1 ± 0.3	1.1 ± 0.8	2.0 ± 1.0	1.6 ± 0.7	0.7 ± 0.3	1.0 ± 0.6	1.0 ± 0.3
DWI (days)	414	353	192	171	343	340	343
Coordinates (UTM)	30SWH2357	30SWH3542	30SVH4721	30SVH8134	30SVG0679	30SUG8679	30SVG1084
Altitude (m)	780	1270	400	637	476	289	442

), some general patterns can be described. Conductivity values indicate that all ponds are fresh-hyposaline during the study period (sensu [28]), ranging from 57.7 in Ardal to 955 μS cm⁻¹ in Quinta, both in November 2010. For all studied ponds, Chl-a concentration was extremely low (<10 μg L⁻¹), which additionally to the high-water transparency, lastly promotes that macrophytes may cover a high percentage of pond area. Contrary to it was expected based on the low algal biomass, TP concentrations were generally high, ranging from 35.4 in Orcera to 206.7 μg L⁻¹ in Casillas. Finally, flooding period widely differ among the studied ponds, experiencing the shortest flooding time in ponds, which seasonally dried up (Ardal and Santisteban), and the longest in the first to be flooded (Castillo).

Attending to the zooplankton community parameters, the studied ponds also shown obvious differences (see

Table 2. Annual mean and standard deviation of zooplankton community parameters in the studied ponds.

	Total Cumulative Richness	Copepod Abundance (ind L−1)	Branchiopod’s Abundance (ind L−1)	Shannon-Wiener (H’)	Equitability (J)	Simpson (D)	Cal:Cycl+Cla	GranCla:TotCla	Total Biomass (mg L−1)
Castillo	9	47.06 ± 127.21	1.51 ± 4.16	0.71 ± 0.29	0.49 ± 0.20	0.61 ± 0.16	0.0023 ± 0.0070	0.411 ± 0.420	0.104 ± 0.149
Orcera	8	65.73 ± 148.89	22.07 ± 67.51	0.60 ± 0.37	0.44 ± 0.27	0.67 ± 0.21	0.2299 ± 0.5309	0.151 ± 0.205	0.123 ± 0.201
Ardal	6	419.93 ± 678.49	45.62 ± 86.90	0.58 ± 0.22	0.50 ± 0.28	0.64 ± 0.18	---	0.002	0.001 ± 0.001
Santisteban	9	809.82 ± 956.02	18.36 ± 34.63	0.30 ± 0.33	0.20 ± 0.19	0.85 ± 0.20	---	0.0347 ± 0.0482	0.303 ± 0.195
Hituelo	13	39.72 ± 100.37	0.97 ± 2.67	0.46 ± 0.43	0.29 ± 0.27	0.74 ± 0.27	52.9825 ± 118.8370	0.059 ± 0.092	0.074 ± 0.119
Quinta	12	289.70 ± 938.60	11.56 ± 32.24	0.56 ± 0.43	0.40 ± 0.34	0.70 ± 0.25	0.001 ± 0.001	0.234 ± 0.259	0.015 ± 0.014
Casillas	12	21.36 ± 51.69	1.96 ± 4.00	0.47 ± 0.35	0.32 ± 0.26	0.75 ± 0.19	11.9280 ± 23.0303	0.078 ± 0.143	0.039 ± 0.048

). The total cumulative species richness ranges from a minimum of six species in Ardal to a maximum of 13 species in Hituelo pond, with minimum and maximum values in Casillas pond (one species in summer and ten species in winter). As a general rule, the richness of branchiopods species is higher than the copepods. The low values of diversity and equitability are also noteworthy.

Table 3. Seasonal mean abundance and standard deviation of copepods and branchiopods (ind L⁻¹) in all studied ponds. Additionally, shown Cal:Cycl+Cla, GranCla:TotCla ratios and total mean and standard deviation for biomass values (mg L⁻¹).

Pond	Season	Copepod Abundance (ind. L−1)	Branchiopoda Abundance (ind. L−1)	Cal:Cycl+Cla	GranCla:TotCla	Total Biomass
Castillo	Winter	13.86 ± 28.91	1.35 ± 3.07	0.0003 ± 0.0002	0.77 ± 0.37	0.08 ± 0.13
	Spring	5.74 ± 7.91	0.57 ± 0.78	0.009 ± 0.01	0.55 ± 0.38	0.01 ± 0.01
	Summer	115.12 ± 238.43	0.26 ± 0.32	0.0021	0.45 ± 0.47	0.10 ± 0.11
	Autumn	47.75 ± 49.63	5.92 ± 10.25	---	0.014	0.19 ± 0.23
Orcera	Winter	1.32 ± 1.20	0.08 ± 0.08	0.82 ± 0.92	0.071	0.02 ± 0.02
	Spring	8.05 ± 7.13	1.75 ± 3.94	---	0.041 ± 0.030	0.004 ± 0.003
	Summer	312.95 ± 302.73	79.39 ± 126.66	---	0.12 ± 0.01	0.23 ± 0.30
	Autumn	58.11 ± 74.21	6.51 ± 8.76	0.103 ± 0.080	0.47 ± 0.10	0.24 ± 0.24
Ardal	Winter	135.1 ± 156.12	90.74 ± 134.91	---	---	0.0004 ± 0.0001
	Spring	1171.47 ± 1656.33	71.73 ± 49.78	---	0.002	0.001
	Summer	---	---	---	---	---
	Autumn	359.31 ± 392.70	5.92 ± 7.19	---	---	0.001 ± 0.001
Santisteban	Winter	709.65 ± 881.71	11.97 ± 28.41	---	0.07 ± 0.07	0.28 ± 0.15
	Spring	1580.43 ± 833.22	28.17 ± 40.01	---	0.03 ± 0.01	0.20 ± 0.04
	Summer	---	---	---	---	---
	Autumn	603.03 ± 1040.81	25.06 ± 47.76	---	---	0.40 ± 0.54
Hituelo	Winter	136.87 ± 192.97	1.28 ± 1.48	1.18 ± 0.03	0.23 ± 0.06	0.21 ± 0.19
	Spring	7.98 ± 13.59	1.41 ± 3.87	38.31 ± 33.30	0.010 ± 0.005	0.07 ± 0.06
	Summer	4.08 ± 6.18	0.01 ± 0.003	172.21 ± 216.15	---	0.009 ± 0.01
	Autumn	28.37 ± 40.41	0.87 ± 2.81	0.24 ± 0.34	0.08 ± 0.08	0.009 ± 0.009
Quinta	Winter	989.41 ± 1790.86	27.23 ± 52.55	0.001 ± 0.001	0.60 ± 0.05	0.017 ± 0.02
	Spring	4.67 ± 9.92	0.24 ± 0.25	0.003	0.37 ± 0.10	0.02 ± 0.02
	Summer	93.66 ± 132.21	1.57 ± 1.78	---	---	0.006 ± 0.008
	Autumn	57.50 ± 69.46	9.90 ± 20.42	0.0003	0.02 ± 0.02	0.02 ± 0.008
Casillas	Winter	45.52 ± 87.95	3.08 ± 5.93	0.53 ± 0.57	0.04 ± 0.06	0.08 ± 0.05
	Spring	21.83 ± 26.58	1.67 ± 3.21	5.40 ± 0.36	0.24 ± 0.32	0.08 ± 0.04
	Summer	0.46 ± 0.59	1.81 ± 3.29	0.32 ± 0.22	---	0.005 ± 0.006
	Autumn	6.56 ± 12.03	0.85 ± 1.94	41.75 ± 33.35	0.12 ± 0.11	0.0003 ± 0.0002

shows the seasonal abundance data for copepods and branchiopods in all studied ponds, and the ratios Cal:Cycl+Cla, GranCla:TotCla and total biomass values. The results obtained show a high variability between stations and ponds in all community parameters, being remarkable the greater abundance of copepods on branchiopods.

The results of the PERMANOVA have shown that considering the seven studied ecosystems, samples from different ponds are statistically different attending both, environmental variables (F_(1,6) = 7.771; p < 0.001) and zooplankton communities (F_(1,6) = 10.44; p < 0.001). When the seasonal differences between all the ponds are taken into account, the results show the same results, both for environmental variables (F_(1,3) = 7.801; p < 0.0001) and for the zooplankton community (F_(1,3) = 1.611; p < 0.01), reflecting a different seasonal dynamic. Furthermore, when the seasonal changes are analyzed for every pond individually, the results shown that for environmental variables exists significant statistical differences in Castillo (F_(1,3) = 2.514; p < 0.05), Orcera (F_(1,3) = 3.768; p < 0.0008), Hituelo (F_(1,3) = 4.777; p < 0.05), Quinta (F_(1,3) = 7.307; p < 0.0001) and Casillas pond (F_(1,3) = 6.045; p < 0.05), while no differences have been shown in Ardal (F_(1,1) = 3.501; p = 0.2011) and Santisteban (F_(1,2) = 1.53; p = 0.2595). For zooplankton communities the results are quite different, with statistical differences in three of the same ponds (Orcera (F_(1,3) = 4.015; p < 0.005), Hituelo (F_(1,3) = 1.768; p < 0.05) and Casillas (F_(1,3) = 2.241; p < 0.05)), and without differences in the other four ponds (Castillo (F_(1,3) = 1.356; p > 0.05), Ardal (F_(1,1) = 0.6502; p > 0.05) and Santisteban (F_(1,2) = 0.749; p > 0.05) and Quinta (F_(1,3) = 1.479; p > 0.05)).

Spearman correlation between zooplankton community parameters and environmental variables shows that among all variables, Chl-a, Secchi depth, pond area, depth and TN were those with a greater influence on the parameters of the community (

Table 4. Coefficients of the Spearman correlation among zooplankton parameters and environmental variables. * p < 0.05; ** p < 0.01; *** p < 0.001.

	Conductivity	pH	Turbidity	Chl-a	SM	Area	Depth	Secchi	Temperature	O2	MO	TP	TN	DWI
Total richness				−0.290 *	−0.318 **	0.256 *	0.284 *	−0.350 **					−0.250 *
Copepod richness	−0.403 ***						0.438 ***	−0.270 *	−0.272 *				−0.387 **
Branchiopod’s richness				−0.316 **	−0.318 **	0.304 **	0.279 *	−0.321 **
Shannon-Wiener (H’)				−0.246 *
Equitability (J)						−0.250 *
Simpson (D)				−0.231 *
Cal:Cycl+cla				−0.252 *		0.350 *	0.535 ***	−0.270 *			−0.301 *		−0.403 ***
GranCla:TotCla							0.270 *
Total biomass

). Chl-a was negatively related to total zooplankton and branchiopods richness, to zooplankton Shannon–Wiener diversity, to zooplankton dominance, and to the ratio Cal:Cycl+Cla. TN concentration was negatively correlated to total zooplankton and copepod richness and to the ratio Cal:Cycl+Cla. Positive relationships between Secchi depth and richness were also observed, and also with the ratio GranCla:TotCla. Finally, morphometric parameters such as pond area and depth were, in general, positively and significantly related to zooplankton richness (i.e., total and branchiopods richness), and to the ratio Cal:Cycl+Cla.

The relationship between zooplankton assemblages and environmental variables in CCA (

Table 5. Summary of CCA analysis among zooplankton assemblages and environmental variables. Bold letters reflect the significant related variables for each axis.

	Axis1	Axis 2	Axis 3
Axis summary statistics
Eigenvalue	0.677	0.594	0.401
p-value	0.001	0.001	0.001
% of variance explained	10.3	9	6.1
Cumulative % explained	10.3	19.3	25.4
Pearson Correlation, Spp-Evnt*	0.920	0.884	0.842
p-value	0.001	0.001	0.001
Correlations of environmental variables and canonical axes (Intraset correlations of terBraak, 1986)
Conductivity	0.056	0.232	−0.490
pH	−0.052	0.331	0.140
Turbidity	0.470	−0.114	0.195
Chl-a	−0.135	−0.173	−0.384
SM	−0.531	−0.051	0.207
Area	0.633	0.711	−0.038
Depth	−0.117	0.898	0.252
Secchi depth	−0.350	−0.293	−0.431
Temperature	0.209	0.007	−0.091
O2	0.235	0.007	0.319
OM	0.199	−0.383	−0.430
TP	0.361	−0.298	−0.039
TN	−0.044	−0.481	−0.287
DWI	−0.565	0.318	−0.476

) was significant (Monte Carlo test, p < 0.001). The first three CCA axes explain 25.4% of the total variance of the species–environment relationship. The first axis (10.3%) was positively and significantly correlated with turbidity and negatively with SM and DWI. The second axis (9%) was positively correlated with area, depth and Secchi depth and negatively with TN. Finally, the third axis 3 (6.1%) was only negatively correlated with conductivity. The CCA ordination (Figure 2) shows that on the first axis, samples with positive coordinates that represent the initial stages after the flooding event are dominated by the copepods Tropocyclops prasinus and Arctodiaptomus wierzejskii. These species appear in sites with high values of turbidity and obviously with low percentages of macrophytes. On the other hand, samples located in negative coordinates (with an increase of the flood period, macrophytes coverage and reduction of turbidity) are represented by the copepods Diaptomus cyaneus, Macrocyclops albidus and by the cladocerans Chydorus sphaericus and Simocephalus vetulus. The second axis distributes the samples according principally to pond morphometry and nutrient concentration (TN). Samples located in the upper part of this axis are present in ponds with higher depth and the species associated are the cladocerans Daphnia magna, Ceriodaphnia dubia, Pleuroxus aduncus and the copepod Acanthocyclops venustus. In the lower part, appear samples that have greater nutrient concentrations, and the species related are the anostraca Chirocephalus diapahanus, the cladocerans, Pleuroxus letourneuxi, Macrothrix hirsuticornis and the copepod Cyclops sp.

Finally, Figure 3 shows optima and tolerance of the most important variables that explain the zooplankton distribution in the CCA analysis (turbidity, SM, DWI, area, depth, Secchi depth, TN and conductivity). The position of the different species in each of these figures provides a very interesting information about their ecological niche. Thus, for example, the predominant position in the lower parts of the figures of the anostraca species Chirocephalus diaphanus indicates, and thus agrees with the bibliographic data [21], that it is a species that lives in small and shallow ponds, with short hydroperiods and with transparent and fresh or subsaline waters. In the same way, these results also show that calanoid diaptomids appear in ponds with a medium-high surface area and depths greater than 1 m, being also a group that does not tolerate very high concentrations of nutrients (TN). Finally, it is also interesting to note the predominance of brachiopods in ponds with a high percentage of submerged macrophytes.

4. Discussion

The Mediterranean region includes high-diversified wetlands with unusual limnological characteristic that makes them different from the temperate aquatic systems [1]. These wetlands are valuable in terms of salinity, which range from freshwater to hypersaline and, also, in terms of the hydroperiod comprising a gradient from permanent to temporary systems [1,29]. In this context, it is essential to recognize that zooplankton represent a model community for monitoring changes in these ecosystems [30].

The results obtained in this study shows that local conditions, like geology and the land uses in the watersheds between others, are manifested in differences in the environmental variables and also in the different biotic communities (zooplankton community) that ponds present. These results coincide with those obtained previously in the same study area [31,32]. Moreover, these differences are also globally pointed out between seasons. In a more detailed study and considering the environmental variables, these differences are not patent in ponds with a lower hydroperiod (Ardal and Santisteban), but are clearly manifest in ponds located in the Guadalquivir valley (Hituelo, Quinta and Casillas) and in the other two ponds (Castillo and Orcera). Regarding the zooplankton community, the situation is unlike, since only in three ponds significant statistical differences have been found (Orcera, Hituelo and Casillas). In these ponds we could indicate that seasonal differences in environmental conditions are reflected in zooplankton communities. However, the absence of differences in two of the other ponds would be clearly justified by the similarity in the environmental variables (Ardal and Santisteban), while in the other two ponds (Castillo and Quinta) it would be more difficult to find a justification.

Our results clearly show the high intra and interannual variability of Mediterranean wetlands [1]. Thus, if this study were carried out on a multiyear scale, the differences that have not been manifested in some studied ponds would surely become evident. Among the most obvious differences we would have the changes in salinity, which would clearly manifest themselves in different zooplankton communities [33], a community that shows adaptations to unpredictable aquatic habitats [34].

The results found in this study are similar to those shown in the bibliography, which exhibit high seasonal variability in populations of crustaceans in Mediterranean temporary ponds, with marked differences between the filling, intermediate and desiccation phases [35,36,37]. However, there are some differences, such as that normally anostracans and diaptomid copepods are the most abundant organisms in the filling phase, cladocerans in the intermediate stage—coinciding with a clear water stage—and some cyclopoids and small cladocerans in the drying stage [35]. In our case, this pattern is not so evident because: (i) except in two ponds, there is no drying phase; (ii) the late and short duration of the flood period could affect the composition of the zooplankton assemblages [36,37] and so, could surely also limit individual and population growth, and restrict the presence of some species [38]; (iii) during the study period, copepods were the most abundant group in all seasons and ponds (except in the Casillas pond in summer), without great abundance of branchiopods in any case. This situation could be due to the high trophic level of the ponds studied [31], a situation that is supported by the low values of the Cal:Cycl+Cla ratio, which has been associated with eutrophic systems [22].

Our results obtained in the correlation analysis reflect the suitability of use total zooplankton and cladocerans species richness and the ratio Cal:Cycl+Cla as biological monitoring indicators. Our results confirm the results obtained in previous studies carried out in coastal Mediterranean wetlands [39]. Attending to the factors that regulates this previous result, several studies [40,41] had reported a negative relationship between zooplankton richness and eutrophic conditions, a result that is also observed in this study (Table 3). Similar results have been previously found indicating the same relationships between TN concentration and zooplankton richness in other Spanish Mediterranean wetlands [39].

On the other hand, and alike to [42], we found a positive and significant influence of pond area and depth on richness of total zooplankton and cladocerans richness. This positive relationship might be explained by a greater number of available niches and larger habitat heterogeneity [43]. However, and contrary to the results observed by many other authors [44], in our case we have not found a relationship with the macrophytes coverage, which support an increment of the habitat heterogeneity, and zooplankton richness. At this point, it is worth to note that this correlation between zooplankton richness and macrophytes coverture has been mainly obtained in lakes characterized by the presence of fish where macrophytes act as a refuge zone [45]. However, a great majority of Mediterranean wetlands, and all of the studied ponds are fishless ecosystems, and thus, an increase in macrophytes biomass could cause a reduction in the phytoplankton biomass, which lastly promote a decrease in the biomass and richness of the zooplankton. Furthermore, macrophytes may also host a diverse and abundant macroinvertebrate assemblage, and it is likely that zooplankton will also face predatory macroinvertebrates once they enter the littoral zone [46]. Therefore, if this phase of high abundance of macrophytes is prolonged in time, as has happened in the study ponds, it translates into a decrease in the total richness of zooplankton.

The results of the CCA have evidenced that eight variables play an important role for determining the zooplankton composition. The extraction of these variables and plotting in a graph have led to a well-defined separation of three groups, which represent a transition of distinctive species assembly along the hydroperiod. In the initial stages, an increase in nutrients concentration occurs as a result of nutrient release from sediments linked, among others, to desorption processes and also due to the organic matter mineralization [47,48,49]. Moreover, an increase in turbidity values also occurs as a consequence of the resuspension of sediment particles during the filling process and by the growth of phytoplankton promoted by an increase in nutrient availability. The species related with these conditions appears in the lower-right quadrant of the CCA, represented among others by Chirocephalus diaphanus, an anostraca, typical of initial stages after pond rewetting [50]. As we indicated previously, Figure 3 shows the same results, that is, a species that lives in small and shallow ponds, with short hydroperiods.

In the intermediate and successive stages of the hydroperiod, nutrient concentrations decreased mainly due to the biological uptake (e.g., macrophytes). It is well-known that the establishment of macrophytes lead to maintain a clear water state because they are effective competitors for nutrients and light with phytoplankton [51]. During this stage, two groups of species can be distinguished. The first is placed in the upper-right quadrant and it is related to intermediate stages characterized by large cladocerans as genus Daphnia. It plays an important role in establishing water clarity by top-down control (grazing) of phytoplankton [52,53,54]. The second group appears in the left part being related to maturate stages after flooding events, which in the case of the study wetlands is related to a longer flood period and therefore a greater coverage of macrophytes with flowering and fruiting stages. This period can be lengthened more or less depending on the hydroperiod of each year, and it is characterized by Dyaptomus cyaneus, Macrocyclops albidus and Chydorus sphaericus. Figure 3 confirms the presence of these species in ponds with high macrophyte coverage.

Finally, it is necessary to indicate that changes in the physical and chemical variables are important drivers of the dynamics of the aquatic community of temporary wetlands, including the zooplankton community. Therefore, an adequate assessment of these wetlands requires a previous step of understanding the organization, structure and functioning of the ecosystem and the composition of these aquatic communities [14,50]. The results obtained under a regional approach facilitate the future application of measures to address possible threats in these temporary ecosystems. Among all of them, the most significant threat is climate change, being its most notable outcome the unpredictability in rainfall patterns that ultimately affects the hydroperiod of temporary aquatic systems. Accordingly, in a scenario of global change, drastic changes are expected in physical, chemical and biological variables, which together with the inherent unpredictable functioning of Mediterranean wetlands will have a greater effect on these ecosystems. In this sense, the results obtained have clearly revealed the role of the zooplankton community as a useful biological indicator of the functioning of the aquatic ecosystem. All in all, our results support the need to include the zooplankton community in monitoring studies as they allow a rapid identification of ecosystem changes, helping to propose management tools for these ecosystems.