Mountain ponds of glacial origin are one of the most remote and undisturbed aquatic environments in Europe (Čiamporová-Zaťovičová et al. 2010). These high alpine ponds are particularly sensitive to global change, especially to climate warming, and constitute ideal systems for long-term monitoring of global changes (Gurung 2005; Oertli et al. 2008). Among the biological communities present in these lakes, benthic invertebrates are excellent indicators of local as well as global temperature changes (Fjellheim et al. 2000) and they are one of the most common groups of organisms used to assess the health of aquatic ecosystems (Sharma, and Rawat 2009).
On the other hand, these small aquatic systems are rich and diverse habitats and play a key role in safeguarding aquatic biodiversity (Biggs et al. 2005; Shieh, and Chi 2010). These systems are considered to support a high richness of organisms, particularly macroinvertebrates (Oertli et al. 2002; Williams et al. 2004), both on a local and regional scale (Toro et al. 2006). Understanding patterns of biodiversity distribution is essential to conservation strategies (Gaston 2000), and resolving the relative contributions of local and regional processes might provide a key to understanding global patterns of biodiversity (Gaston, and Spicer 2004). Therefore, it is crucial to achieve a good understanding of the richness patterns in these systems, especially the environmental factors influencing their diversity.
Several studies document clear relationships between composition or richness of macroinvertebrate communities and a variety of ecologically relevant gradients in ponds, such as hydroperiod (Collinson et al. 1995), surface area (Céréghino et al. 2008), water chemistry (Friday 1987; Boix et al. 2008), pond connectivity (Gascón et al. 2008; Oertli et al. 2008), habitat heterogeneity (Della Bella et al. 2005), presence of large predators (Zimmer et al. 2001) and altitude (Lin et al. 2006). However, only a few of them have specifically addressed the issue of which of a large set of factors really determine richness of ponds (see Rundle et al. 2002), and only rarely in alpine systems (Hinden et al. 2005; Füreder et al. 2006; Collado, and De Mendoza 2009). Ecological researches on littoral macroinvertebrates from Spanish mountain ponds have only recently been undertaken. Good examples are the contributions by Catalan et al. (2009) and De Mendoza, and Catalan (2010) in Pyrenean lakes, as well as those by Granados, and Toro (2000) and Toro et al. (2006) in the Central mountain system.
The effects of ecological factors on macroinvertebrate communities (and on their richness in particular) are complex and difficult to disentangle, partially because richness–environment relationships are typically masked by other physical or geometric constraints (De Mendoza, and Catalan 2010). We might expect richness to be depressed in harsh environments, such as high-altitude (Hoffman et al. 1996), ultraoligotrophic (Wetzel 2001) or acid ecosystems (Schell, and Kereks 1989), as well as in heavily impacted systems, as those subjected to salmonid introductions (Carlisle, and Hawkins 1998; Martínez-Sanz et al. 2010). In contrast, it should be high in heterogeneous ponds supporting a high habitat diversity. The characteristics of the habitat (substrate type, vegetation cover, etc.) are also likely to influence richness, although in unknown ways. All these interactive factors give a complex scene making richness prediction a difficult task.
A good understanding of how the combined action of factors determines richness of mountain ponds would aid in management of ponds and biodiversity preservation. It would be helpful to know which environmental factors may be relevant and which secondary for a particular system. In particular, we must gain knowledge on the effect of human activities compared to natural factors. In this study, we aim to generate, by means of a prediction model, a gradient of factor relevance on macroinvertebrate richness in mountain ponds from a relatively wide region. We have also individually analysed the hypothetical influence on richness of several factors traditionally considered relevant for pond macroinvertebrates.
Materials and Methods
Fifty-one mountain ponds of quaternary glacial origin were selected for this study (Fig. 1). All of them are located in Castilla y León, a wide Spanish region consisting of a vast central plateau (from 600 to 1000 m a.s.l.) surrounded by mountains (altitudes up to 2600 m a.s.l.). The pond selection was intended to record a wide gradient of environmental conditions of altitude (from 1360 to 2195 m a.s.l.), area (between 3000 m2 and 0.12 km2), depth (maximum depths between 0.3 and 15 m), littoral slope, habitat type (substrate and type and coverage of aquatic vegetation), and water permanence (temporary and permanent systems).
All the ponds were sampled once in late June or early July between 2004 and 2008. Benthic macroinvertebrate samples were collected from the littoral zone (at depths below 1 m) with a pond net (FBA standard, mesh size 500 μm) following a multihabitat time-limited sampling (Collinson et al. 1995; Biggs et al. 2000; Briers, and Biggs 2005). Each pond was sampled during three (ponds below 0.01 km2), four (0.01-0.05 km2) or five minutes (>0.05 km2). Total sampling time was proportionally distributed among the main habitats according to their surface in the pond. Macroinvertebrates were separated from the plant material, counted under a binocular microscope (10×) and identified to genus except for Oligochaeta (identified to class); Dolichopodidae, Empididae, Ephidridae, Sciomyzidae, Simuliidae, Tabanidae, Tipulidae (to family); Limnephilidae, Ceratopogonidae, Chironomidae and Limoniidae (to sub-family).
Twenty-four environmental variables were measured in this study. Several of these environmental variables were directly measured in each pond (maximum depth, surface, conductivity, pH and Secchi disk). The granulometric nature of the substrate was estimated as percentage of stones, gravel or fine substrate in the sampled area. Land use was accounted for by measuring percentage of the basin covered by bushes or wood. Other variables considered were altitude, longitude, percentage of macrophyte cover, concentration of chlorophyll, alkalinity and nutrients in the water column (nitrate, total nitrogen, soluble reactive phosphorus and total phosphorus), and slope. Water persistence was expressed as one of two classes (temporary or permanent) and littoral slope was qualitatively estimated by assigning a value between 1 (nearly plain bottom) and 5 (steep shores).
In addition, we generated a discrete variable to asses the impact of introduced fishes (from 0, no fish, to 4, maximum impact). Each pond was assigned to a fish impact class on the basis of previous information on fish stocking in some ponds (Granados, and Toro 2000a; Martínez-Solano et al. 2003), expert knowledge of local anglers, personal observations during the survey and the results of a fyke-net sampling carried out during summer 2009 to obtain rough estimates of species abundances. Fish density as inferred from fyke-net sampling and personal observations during the survey and number of fish species present (also used as a continuous variable in the analysis) have been key criteria to assign a pond to an impact class. The highest impact value (4) was attributed to ponds with high intensity of salmonid stocking or extremely high densities of other species. These population levels are always kept by repeated introductions. Ponds with no fish were assigned to class 0. The remaining ponds were ranged between these two extremes according to the data provided by the fyke-net sampling.
There is one more variable which could prove critical in richness prediction: habitat heterogeneity. We addressed this issue by identifying five possible meso-habitat types in the study area study (hydrophytes, helophytes, stones without vegetation, sand/gravel without vegetation and fine substrate without vegetation) and counting how many of them were present in a given pond.
This variable set accounts for natural among-pond differences as well as for the main impacts in the study area (fish stocking and changes in land use). Eutrophication, which can be consequence of these and other environmental pressures, is also considered. Of all these variables, only a subset was used for further analysis. Many variables were ruled out because they were highly intercorrelated (significant Pearson or Spearman’s correlation, p=0.05). Several variables were surrogates for others. For example, pH was correlated with conductivity. Other variables as surface and depth were correlated with littoral slope, water persistence and percentage of macrophyte. Some evident relationships among total nitrogen and total phosphorus with soluble reactive phosphorus, chlorophyll and Secchi depth were also detected. Twelve variables were finally included in the model: longitude, altitude, pond area, depth, alkalinity, pH, nitrate, total nitrogen, total phosphorus, fish richness, fish impact, and number of mesohabitats (Tab. 1). Altitude was a surrogate for substrate type and land use, as also shown by de Mendoza, and Catalan (2010), while pond area and depth were significantly correlated with littoral slope and macrophyte cover. Maximum depth, on the other hand, is a good measure of the hydroperiod length with the additional advantage on water permanence of being measured as a continuous variable (and not just as one of two values, temporary or permanent).
General prediction model. We performed multiple linear regression (MLR) analysis generated by a step-wise selection technique to quantify the magnitude, direction and significance of the relations (Robinson et al. 2004). Thirty-nine ponds were randomly selected for MLR analysis while the remaining ones were kept as an independent data set to test the model.
The linear regression model assumes a linear relationship between the dependent variable and each predictor. We examined a scatterplot of richness by each independent variable to determine whether a linear model is reasonable for these variables. However, in practice, this assumption can virtually never be confirmed; fortunately, multiple regression procedures are not greatly affected by minor deviations from this assumption. Linear regression model involves several additional assumptions. Firstly, the error term has a normal distribution with a mean of 0. Again, even though most tests (specifically the F-test) are quite robust with regard to violations of this assumption, it is always necessary, before drawing final conclusions, to review the distributions of the major variables of interest. We tested this assumption by using the Kolmogorov-Smirnov test over the error term. Secondly, the variance of the error term is constant across cases and independent of the variables in the model. This homoscedasticity of the error term was tested by a residual scatterplot in the analysis. Thirdly, the value of the error term for a given case is independent of the values of the variables in the model and of the values of the error term for other cases. We tested autocorrelation among error terms by Durbin-Watson test. Besides, we removed from the regression analysis intercorrelated variables (Spearman’s test, p=0.05). The statistical processing was performed with the SPSS Statistic 17.0 package.
We tested the MLR model by applying it to an independent data set (12 ponds of the study area) and comparing predicted with observed richness values by means of a one way analysis of variance (ANOVA) of repeated measures and Scheffé test.
Finally, we used a generalized additive model (GAM) to explore the graphical response of richness to the variables that were significant in the linear regression model (Schimek 2000; Guisan et al. 2002). Through a cubic spline scatterplot smoother (along with the 95% confidence bands) generated by the model, we can evaluate the nature of the relationship between the predictor and the residualized (adjusted) dependent variable values, and hence the nature of the influence of the respective predictor in the overall model. Computational details regarding the cubic spline smoother, and comparisons to other smoothing algorithms, can be found in Hastie, and Tibshirani (1990) and Schimek (2000). The statistical processing was performed with the STATISTICA 8.0 package.
Particular prediction models. We generated an additional multiple regression model on ponds of the study area without human impact (fish impact: 0 or 1; 37 ponds) to explore the possible relation between observed richness and other natural predictors undetected in the general model. Multiple linear regression was also applied to additional data subsets in order to remove the effect of both habitat heterogeneity and fish impact. Two models were generated, one on temporary (fishless) ponds (n=11) and another on ponds with one single habitat (stones) and low (0-1) fish impact (n=8)
Influence of selected variables. On the same pond set, we explored by means of box-plots and one–way analysis of variance (ANOVA) hypothetical relation-ships between three variables (altitude, habitat type and water permanence) and richness.
General prediction model
The MLR model included two variables as the best predictors of richness: number of habitats (NH) and impact by fish introductions (FI) (R2=0.579, F2.29=19.913, p<0.001). The relationship between richness and the predictors was described by the following equation: S=4.922(NH) – 2.359(FI)+11.227, being number of habitats the first variable selected by the model. The significance value of the F statistic is less than 0.01, which means that the variation explained by the model is not due to chance. The regression and residual sums of squares were approximately equal (863 and 628, respectively), which indicates that about half of the variation in observed richness is explained by the model. The large value of the multiple correlation coefficient (R=0.761) indicates a strong relationship between the observed and modelpredicted values of the dependent variable and error of estimate is considerably low, about 4.656. The tolerance collinearity statistic (T), percentage of the variance in a given predictor that cannot be explained by the other predictors, was very near to one (0.982), and variance inflation factor was below two (1.018). Therefore, there were not multicollinearity effects between predictors and the standard error of the regression coefficients will not be inflated. On the other hand, about assumptions of MLR, the error terms of this regression had a normal distribution (Kolmogorov-Smirnov test); they were constant across cases and independent of the variables in the model (there were not heteroscedastic) and Durbin-Watson showed not autocorrelation between them (1.719). Therefore, all the assumptions of the multiple linear regression model were fulfilled.
The model performance was tested on a set of 12 new ponds. ANOVA showed no significant difference between observed richness and richness predicted by the MLR model (p=0.525) in this new data set.
The graphs generated by GAM give an insight into the response of richness to the two variables in the MLR. Richness and number of mesohabitats showed a positive relation, but not all along the gradient. In ponds with only one and two habitats, apparently the richness was the same (Fig. 2). In contrast, richness was inversely related to fish impact (Fig. 3), although this relationship was only evident for high impact values (3 or 4).
Particular prediction models
In order to analyse only the influence on richness of natural variables, we carried out another multiple linear regression on a subset of the available samples: those with no or low fish impact (0 or 1). In the model generated, only the number of habitats significantly explained the response of richness (S=4.533(NH) + 10.699). This model was separately applied to two groups of samples: ponds with no or little fish impact (impact value 0 or 1) and ponds with high fish impact (3-4). Differences between observed and predicted richness values were significant for ponds with high fish impact (p=0.006) but not for low fish impact (p=0.99). This result helps to show that negative effect of fish is only evident at high impact level (Fig. 3).
Multiple linear regression was also applied to additional data subsets in order to remove the effect of both habitat and fish impact. The model applied to permanent ponds with stony substrate and low fish impact proved that no variable significantly explained richness values. When applied to temporary, vegetated ponds with low fish impact, a significant contribution of depth and, secondly, total nitrogen to the model was observed, with F significant value of 0.02. This result seems to indicate a direct relationship between hydroperiod length (as indicated by depth) and richness.
Influence of selected variables
Finally, we explored by means of box-plots and ANOVA hypothetical relationships between three variables (altitude, habitat type and water permanence) and richness. Only ponds with low fish impact (impact level 0 or 1) were considered. The results are shown in Figs 4 to 6. No significant differences were detected except for water permanence (p=0.031): richness proved to be slightly lower in temporary ponds.
These results might be conditioned to some extent by the effect of habitat heterogeneity in the case of permanent ponds. Therefore, we tried to cast some light on the effect of altitude on richness by selecting a small but homogeneous data subset: nine permanent ponds with just one mesohabitat. The selected ponds were otherwise heterogeneous in altitude (1480 to 2195 m a.s.l.) and richness (from 13 to 22 taxa). The Pearson correlation coefficient showed no significant relationship between altitude and richness (R2=0.1589).
The only relevant variables influencing macroinvertebrate richness in the study area were habitat heterogeneity (measured as number of mesohabitats in a pond) and intensity of fish impact. The relationship between these variables and richness is well known and its ecological basis well understood. The habitat heterogeneity classic hypothesis proposes that an increase in the number of habitats leads to an increase in species diversity in a landscape (MacArthur, and MacArthur 1961) because of an expansion in the number of partitionable niche dimensions (Willig 2000). On the other hand, several researches in Europe, North America and Oceania have reported negative consequences of fish introductions on macroinvertebrate communities, whether on their composition or on richness (Crowder, and Cooper 1982; Carlisle, and Hawkins 1998; Oertli et al. 2001; Knapp et al. 2001; Zimmer et al. 2001; Tate, and Hershey 2003). Similar results have been reported in Spanish lakes and ponds after the introduction of lake trout (Salvelinus fontinalis) in Peñalara Lake (Toro et al. 2006) or of rainbow trout (Onchorhynchus mykiss) in mountain ponds of the Iberian System (Martínez-Sanz et al. 2010).
The result is relevant because it stresses the predominant influence of heterogeneity and man-induced impacts over that of variables traditionally considered to influence richness in mountain ponds and lakes. Some of these factors, apart from heterogeneity and fish stocking, are altitude (Hinden et al. 2005; Čiamporová-Zaťovičová et al. 2010; Mendoza, and Catalan 2010), pH (Schell, and Kerekes 1989; Dumnicka, and Galas 2002; Mousavi 2002), and water permanence (Wissinger et al. 2009). Different habitat types (stone, gravel, cobbles, sand, silt, aquatic plants, or woody debris) are also known to support macroinvertebrate communities differing in size structure and species composition (Oertli 1995; Gardner et al. 2001; Tolonen et al. 2001; Weatherhead, and James 2001). Of these variables, pH was not expected to influence richness in our study. Only pH values below 5.0 seem to have an appreciable effect on lake macroinvertebrate richness or community composition (Mossberg, and Nyberg 1979; Harve, and McArdle 1986; Schell, and Kereks 1989). Sometimes, the acidification threshold causing richness depletion is even much lower than 5.0 (Crisman et al. 1980). Therefore, it is hardly surprising that pH did not proved to be relevant in our study area, where pH values were always above 4.9 and usually above 6.0.
Neither was pond size expected to play a significant role. The assumption of increasing richness with increasing habitat size dates back to the island biogeograpical theory proposed by MacArthur, and Wilson (see Rosenzweig 1995 for a review) and has often been invoked to explain richness patterns. However, the relation between macroinvertebrate richness and pond size is far from evident. Several studies have failed to find such relationship or found it only for selected taxa (Gee et al. 1997; Oertli 2002). Of course, this does not fully refute the richness-size hypothesis. We should state, instead, that size effect is often overruled by other factors, such as vegetation and fish presence (as proposed by Scheffer et al. 2005) or heterogeneity (as suggested by Heino 2000, Céréghino et al. 2008, or by the results shown in this study). In a sense, our results do not completely oppose the theory of island biogeography, which predicts an effect of size on richness because larger islands have more habitats, as well as larger populations (see Rosenzweig 1995). Nevertheless, the possibility of a detectable effect of pond size itself cannot be ruled out provided that the range of areas is wide enough, as supported by the data provided by (Aho 1978) whose study comprised ponds and lakes between 0.03 and 147.3 km2.
The reason why some other variables (altitude and water permanence) were not significant in our MLR model was not so easily explained and required an additional data processing to check whether they might play a secondary role (unravelled in the general model) or not. Besides, habitat type could not be included in the MLR model but deserved consideration too.
Water permanence has sometimes considered to shape macroinvertebrate richness in a wide range of pond types (Collinson et al. 1995; Della Bella et al. 2005; Waterkeyn et al. 2008), but not always (Boix et al. 2008). We did find temporary ponds to harbour less taxa than permanent ones, although the relationship was weak (ANOVA), and, in any case, superseded by heterogeneity and fish impact (as shown by the global regression model). Nevertheless, no much importance must be attached to this result. It might perfectly be a simple consequence of the low habitat heterogeneity in temporary ponds. In fact, when they are compared with permanent ponds having similar number of environments (one or two) and fish impact (0-1), differences disappear (richness values between 13 and 21 in temporary ponds, and between 14 and 22 in permanent ones, except for one pond with 27 taxa).
There is no doubt that different habitats usually harbour different macroinvertebrate taxa (Weatherhead, and James 2001; Della Bella et al. 2005). However, the effect of habitat type on richness is less known. This is a relevant issue because of the implications when interpreting variable-richness relationships. If there were substantial among-habitat differences in richness, trends attributed to environmental gradients (altitude, for example) might partly be a consequence of varying environments or pond type (vegetation cover, type of substrate, etc.) making it difficult to draw sound conclusions. Vegetated environments have sometimes been reported to have higher richness values than unvegetated ones (Wetzel 2001) and, on the other hand, sediments environments frequently have low taxonomic richness (Oertli et al. 2005). No significant among-habitat differences in richness were found in this study in spite of the wide range of habitat types involved (stones, gravel, silt, helophytes, hydrophytes). Within-habitat variability is very high but average values are similar across types. We must assume, therefore, that factors other than habitat type are responsible for differences in richness.
Altitude has often been considered the most important driver force of macroinvertebrate biodiversity in mountain ponds (e.g., Oertli et al. 2000). Why was not altitude relevant in our study? Certainly, the altitudinal range available (1360-2195 m a.s.l.) is not as large as those in the Alps or the Pyrenees; however, it should be enough to detect changes if richness decreased gradually along the altitude gradient. One possible explanation is that richness drop is only apparent above an altitude threshold. The conclusions obtained by Mendoza, and Catalan (2010) and Čiamporová-Zaťovičová et al. (2010) are in this line. Mendoza, and Catalan (2010) showed that the inverse relationship between the number of macroinvertebrate groups relative to altitude was expressed best by a simple linear regression model (p<0.001) when considering only lakes above 2500 m. Čiamporová-Zaťovičová et al. (2010) found that the strongest response of macroinvertebrate fauna will likely occur in shallow alpine lakes above 2000 m a.s.l., when the duration of ice cover shortens below the 190-day threshold. Similarly, Lin et al. (2006) found only 4.6 taxa per pond in eight alpine ponds above 3000 m a.s.l. in central Taiwan, indicating that taxon richness in Taiwan ponds was related to altitude, with increasing altitude being associated with a reduction in species richness In addition, Hinden et al. (2005) showed an altitude-related decline in taxon richness in ponds of Switzerland (1860-2757 m a.s.l.). However, an analysis of their results show no richness decline when only ponds below 2500 m a.s.l. are selected. In the light of these data, it is quite possible that effects on richness are only evident above 2000 or 2500 m a.s.l., near or above the upper altitude limit in our study. Most authors explain this relationship as a consequence of harsh conditions in high altitude ponds and the specific adaptations required to survive in such extreme conditions (Čiamporová-Zaťovičová et al. 2010; Mendoza, and Catalan 2010). The results obtained in our study might suggest another possible explanation. High altitude ponds might support a lower number of habitats (see Hoffman et al. 1996). It is likely that many ponds from high altitude zones are dominated by a single meso-habitat (probably stones or sand) while lower lakes may incorporate other mesohabitats (silty sedimentation areas, vegetation of different types, etc.). As shown in our study, littoral macroinvertebrate richness may be closely related to number of habitats (but not to habitat type). In our study area there is not a trend towards habitat impoverishment at high altitudes. This fact might be responsible for the lack of altitude-richness relationship if heterogeneity (but not altitude) was the real determinant factor. Habitat homogeneity and harsh conditions might be complementary explanations of low richness in high-altitude ponds.
For a correct pond management, we must gain knowledge on the simultaneous effects of factors and their relative importance for richness. Furthermore, in areas where there are few (if any) pristine systems left, such as most of Europe, it is necessary to define to which extent the effect of man-induced impacts can overcome that of natural factors. As deduced from literature, macroinvertebrate richness in mountain ponds can be expected to depend on a number of interacting variables, such as altitude, heterogeneity, habitat type (substrate, vegetation), water permanence, pond size and, when extreme values happen, pH. However, heterogeneity is likely to be an essential factor relegating the remaining ones to a background position. Apart from these natural variables, richness can greatly depend on man-induced impacts (Wilcove et al. 1998; Steinman et al. 2003; Martínez-Sanz et al. 2010). The role played by fish in determining macroinvertebrate richness is particularly outstanding for mountain pond management. As deduced from the present study, several relevant evidences arise. Firstly, fish introductions may become a severe impact for littoral macrobenthic communities, even more than any other disturbance like land use changes. Secondly, low density fish communities, even though artificially introduced, seem to have no evident effect on macroinvertebrate richness (lower than other factors, in any case, and, thus, undetectable). In contrast, high densities caused by fish stocking can turn fish into an essential factor overcoming most other environmental variables. Nevertheless, sound conclusions on this issue require further, detailed analyses taking into account fish community composition and density.