Integrity of fluvial fish communities is subject to environmental gradients in mountain streams , Sierra de Aroa , north Caribbean coast , Venezuela

We examined physical habitat and fish assemblages in rivers of the Aroa Mountains (Venezuela) with different levels of environmental protection due to the creation of Yurubí National Park within the drainage. We developed an Index of Biotic Integrity (IBI) and evaluated it using principal components analysis (PCA) and canonical correspondence analysis (CCA). Tributary rivers were divided into classes according to their origin (protected by the park) and physical characteristics of each, including substrate. Fishes were captured using standardized electrofishing. Fish communities showed greater species richness in heterogeneous habitat and protected rivers but overall abundance was higher in unprotected and impacted rivers. The IBI was sensitive to these differences and the scores were higher in protected rivers. The IBI detected degree of disturbance of fish communities without direct consideration of habitat parameters measured. The PCA revealed a gradient in substrate heterogeneity. Similarly, CCA revealed differences in fish assemblage composition along the environmental gradient and that varied with protection status of the river. The relationship between PCA and IBI scores was highly significant (r = 0.61, P < 0.0001). The PCA and CCA analysis moderately validated the structure and predictability of IBI; but it is still necessary to refine the model and to extend its application for more time and over a wider area.


Introduction
The conservation status of a watershed has direct and measurable effects on the fluvial ecosystem and fish communities.Colonization and residence of organisms are determined by a physical and biological gradient in the hydrosystem continuum and by local habitat conditions.The assemblage of species that make up fish communities is attributed mostly to substrate composition, water type and velocity, but also, the degree of conservation and land use in the basin (Roth et al., 1996;Naiman & Decamps, 1997;Johnson et al., 2005).In contrast to basins with significant human activities, protected areas (such as na-tional parks or riparian forest) are known to act as buffer regions that reduce the entry of sediment into the channels, provide allochthonous organic material and permit the persistence of diverse communities (Berkman & Rabeli, 1987;Poff et al., 1997;Lammert & Allan, 1999;Growns et al., 2003).Furthermore, point disturbances (e.g.dredging, sand or gravel extraction, etc.), even in pristine environments, can impede migration corridors thus reducing populations of sensitive species.These regional and local effects are considered the principal forces influencing the biotic integrity of aquatic ecosystems.
The changes in species richness and composition of fish communities in rivers with different levels of environmental impact can be determined in multiple ways.Attributes of communities (e.g.species richness and diversity) provide basic information about qualitative changes (Magurran, 2003), but the recognition of emergent attributes (structure, composition and organization) permit better assessments of environmental effects on communities (Karr, 1981).The combination of these attributes is used to predict the biotic integrity of aquatic ecosystems.
A pioneer model to estimate the health or biotic integrity is the IBI (index of biotic integrity) proposed by Karr et al., (1986).The IBI considers variation in attributes of fish communities in rivers with different degrees of perturbation relative to a regional baseline or reference system with minimal exposure to anthropogenic perturbations.This model was originally developed in the midwestern United States and has been tested successfully in different regions and ecological situations, for example in Africa (Toham & Teugels, 1999), Mexico (Lyons et al., 2000), Belgium (Belpaire et al., 2000) and Argentina (Hued & Bistoni, 2005).Other methods, such as multivariate or multidimensional analysis, explain how variation in species assemblages are associated with environmental gradients (Bond & Lake, 2003;Gerhard et al., 2004), and how these relate to the conservation status or integrity of hydrosystems.
The effective understanding of the attributes in fish communities and their consequent biotic integrity require data with spatial and temporal amplitude (e.g.complete annual hydrologic cycle) in order to recognize the changes and the relations in the communities respect to disturbance gradients, but these still are limited.However, it is possible that the IBI or multivariate analyses applied on a smaller scale (e.g.dry season) can reflect major changes in fish communities associated with specific gradients of disturbance.Such approaches may not explain all effects on fish communities, but can contribute useful information for management and conservation of natural resources.Most neotropical regions still suffer the disadvantage of only having short term historical databases to document variation in the freshwater ecosystems and their organisms, which limits research and the application of conservation programs.This situation is well known in protected areas and their buffer zones which are usually threatened by human impacts.For this reason there is an urgent need for short term evaluation techniques that identify the main effects of major disturbances on freshwater biological communities.This paper presents a comparative analysis of the biotic integrity of fish communities using an index developed for piedmont rivers and multivariate analysis of parameters that describe environmental gradients and apparent conservation status in rivers of the Aroa Mountains, Caribbean versant of Venezuela.
Study area.The Aroa river basin (2450 km 2 ) drains the western Caribbean coast of Venezuela (Fig. 1).Most tributaries originate in the Aroa Mountains and are protected by the Yurubí national park (YNP, 23,670 h;Gaceta Oficial de Venezuela 26210, March 1960).In the piedmont, climate is seasonal with rainfall (800 -1500 mm) concentrated in July-August (Lentino y Bruni, 1994), and the driest period from January-April.In the mountains, humid and cloud forest predominate, and in the lowland plains there are still remnants of tropical deciduous forest (Huber 1997).Vegetation in the areas bordering the park has been extensively modified by colonial and current human intervention, principally in the piedmont transition zone between mountains and the plains, where deforestation, extensive farming and urbanized areas have dramatically modified both terrestrial and aquatic ecosystems (Rodríguez-Olarte et al., 2006).
Agricultural expansion is extensive even in areas bordering the national park and generates social pressure for the use and administration of natural recourses.The different degrees of conservation (park-protected vs. unprotected watersheds) is manifest in the low percentage of remaining riparian forest, the large amount of water extracted for irrigation, extensive dredged areas of riverbed, and urban pollution.Some tributaries, for example the Tupe river, have intermittent flow in the lower course due to water extraction, whereas others (e.g.Crucito) are dredged to remove excessive sediments to prevent flooding.

Methods and Materials
The sampling stations located in the tributaries originating in the Yurubí National Park were classified as group 1 (G1, n = 14): which includes the following rivers Crucito (Cr: stations Gy and Pt), Tesorero (Te), Guarataro (Gu: stations Gu and Tt), Zamuro (Za) and Carabobo (Ca).Rivers lacking protection: Tupe (stations Ta, Tb and Tu), Cumaraguas (Cu) and Oro (Or) were included in group 2 (G2, n = 10).All stations of G1 were located no more than 2 km from the limits of the park, except for the Crucito river (~ 3 km).We sampled in two reaches in each river, except for the Guarataro (n = 4), Crucito (n = 4) and Tupe rivers (n = 6).Samples were separated by > 200 m and located within the piedmont region of geomorphological transition between mountains and lowland plains.The degree of disturbance associated with each river was estimated as low: (Zamuro and Carabobo rivers), medium (Guarataro, Tesorero and Oro), high (Cumaragua and Tupe) and very high (Crucito) based on visual evaluation of the intervention on the riparian vegetation, channel and human activities.
Since access to the rivers and the use of electrofishing during the rainy season is very limited, fishes were sampled between December and April of 2002-2003.We used electrofishing gear consisting of a portable electrical generator and hand nets with fine mesh (< 5 mm).Captures were made in river sections (~ 50 m) at depths < 1.5 m and from all available mesohabitats (riffle, pool).The temporal effort for collecting each sample (35 min) was determined based on successive captures and accumulative species curve.Most fishes were identified at the time of capture, recorded and returned alive, except for a few individuals that were preserved in formalin [10%] for taxonomic confirmation in the laboratory.
Channel depth (cm) and width (m) of each sample site was measured.For the substrate characterization, we estimated the percent coverage of components inside 1 m 2 squares arranged along transects in the beginning, middle and end of each survey station.We used the classification system proposed by Bain & Stevenson (1999) for the granulometry [sand (0.1-1 mm), fine gravel (0.1-5 cm), coarse gravel (5-25 cm), stones (25-100 cm) and rocks (> 100 cm)].In addition, we estimated the cover of leaves, branches and logs.Also, the concentration of total dissolved solids (TDS, mg/l) was measured for each river, according to APHA (1985).

Analysis
To normalize the environmental data, the percent cover of substrate was expressed as a ratio before being transformed using arcsine square-root (p 0.5 ), and the width, depth and total dissolved solids variables were log(X) transformed (Sokal & Rohlf, 1995).Because some species had differences in their abundance, this variable was square-root transformed to normalize the data.This reduced the impact of the extremes in species abundances for a better interpretation of the data (Pegg & McClelland, 2004;Gratwicke & Speight, 2005).We used oneway ANOVA tests to detect differences in means for habitat and fish variables between river groups.The Pearson coefficient was used to correlate the habitat variables and fish community descriptors (Helms & Feminella, 2005).The abundance variable was not transformed for the IBI procedures.
We classified the fishes according to their known life history attributes or, for poorly known taxa, those of related species.Trophic guilds (omnivores, herbivores and carnivores) were assigned following Winemiller & Taphorn (1989), Taphorn (1992), Rodríguez-Olarte & Taphorn (1995), Arrington et al. (2002) and Keith (2003).The assignments of habitat of occurrence were based on capture data and field observations.Baseline information for each species is presented in appendices.
Index of biotic integrity.The general IBI model is expressed using categories based on metrics (e.g.species richness and composition) and other numerical values obtained by measuring different parameters of the fish community present in each river classified.To permit comparison, the maximum possible value for each parameter evaluated is applied to the optimal condition of that metric (e.g. the highest possible species richness) and a minimum value is assigned the poorest situation encountered in rivers with the most deteriorated condition.The final result is a sum of all the metrics considered and these are then related to a specific biotic integrity class (Karr et al. 1986, Fausch et al. 1984).The rivers of the Aroa Mountains lack detailed information on taxonomy and ecology of several fish species, and because of this we modified the IBI model.We also eliminated the metrics about morphological anomalies and hybridization, the first because they are extremely uncommon, the latter for lack of information.
The baseline conditions considered optimum for species and others metrics were obtained from rivers with the best conservation values (group G1).The Crucito river was excluded from this analysis because it has high levels of impact even though its headwaters are protected by the park.The tolerance class (tolerant, medium tolerant and intolerant) was assigned in concordance with Rodríguez-Olarte y Taphorn (1995) and refined by the reference values determined for a comparison of abundances according to the coincidence of interquartile ranges in box plot diagrams (Silveira et al., 2005).The coincidence, or lack of coincidence, for these ranges was tested for statistical significance to develop a base list of species and to determine the confidence intervals for the IBI metrics (Table 1 and appendices):

Species richness and fish composition category:
1. Sample fish richness: the total number of species in sample; this included the amphidromous species.2. Diversity: (complementary to richness).Although usually the richness shows positive association and redundancy with the diversity index, the disturbance gradient produces changes in relative abundance when comparing diverse communities to simple assemblages (Fausch et al., 1990;Harrison & Whitfield, 2004).We used the Simpson index (D) according to Magurran (2003).This index emphasizes dominance by heavily weighting the most abundant species in the sample, but is less sensitive to species richness.3. Richness of intolerant species: considered taxa with specialized diets (e.g.invertivorous) and with requirements for a clean substrate for feeding or reproduction.A species could be considered intolerant if its abundance or distribution has been drastically reduced or if it is restricted to high quality habitat sites (Hued & Bistoni, 2005).4. Richness of long lived species: included species with annual reproduction and/or long life (e.g.Hoplias malabaricus).
Trophic composition and habitat use category: 5. Proportion of Loricariidae fishes: These catfish are benthic and their diets include autochthonous (algae, roots) and allochthonous (wood) materials.Chaetostoma anomalum was excluded from this metric because it is tolerant to adverse conditions based on the only slight variation in its abundance (see appendices).6. Proportion of invertivores: several fishes feed on benthic invertebrates, but others select drifting aquatic invertebrates and/or those of riparian (terrestrial) origin.7. Proportion of piscivores: included fish predators (e.g.Crenicichla geayi, Roeboides dientonito) that feed on fishes as adults.8. Proportion of drift-feeding omnivores: species associated with currents and that feed on allochthonous materials (Astyanax spp, Hemibrycon jabonero and Pimelodella odynea).9. Proportion of tolerant species: dominant, generalists and/ or omnivorous species (e.g.Creagrutus lassoi) with increase or maintenance of its abundance in degraded rivers could be considered tolerant.
Fish abundance category: 10.Density of fishes per sample (individuals /m 2 ): was estimated for the sample area and compared with the mean width of river in each locality.The fish with lowest abundance (< 5%) were excluded for the IBI (see appendices).
The sum of these ten metrics multiplied by two gives the final result for the IBI, a value that falls within an interval 0-100, according with Ganasan & Hughes (1998).The numerical value thus obtained for the IBI was then assigned to a biotic integrity class (Karr et al. 1986), but modifying their intervals: poor (< 20 points), bad (20-40), regular (>40-60), good (>60-80) and excellent (>80).The IBI values for groups G1 and G2 were compared using the non-parametric Kruskal-Wallis test.
We applied principal component analysis (PCA) to detect environmental gradients among river stations.The PCA techniques reduce the dimensionality of the habitat variables and examine initially the influence of habitat characteristics (Legendre & Gallagher, 2001).Cross-products matrix in PCA were centered and standardized by standard deviation.We retained the first three by the Brocken-Stick criteria (eigenvalue > 1) and the variables with their highest loadings in this axis.Since the ordination of samples and species is constrained by their relationships to environmental variables, we applied canonical correspondence analysis (CCA) to contrast the relationships between the fish community composition and fluvial habitat.In the CCA the total amount of variation in species and sampling stations that can be explained by the environmental variables is calculated by dividing the sum of canonical eigenvalues by the sum of all unconstrained eigenvalues (Ter Braak, 1994).We tested the statistical significance in CCA analysis using the Monte Carlo test with 999 permutations.For the joint plots only variables with r values > 0.500 (cut-off level = 0.5) were shown.The PCA and DCA axes were correlated with primary attributes of fish communities and IBI values to determine and to validate the relationships between the biotic integrity and environmental gradients.Statistical tests were done using STATISTICA version 6 (StatSoft, 2004).For the multivariate analysis we used PC-ORD, Version 4.25 (McCune & Mefford, 1999).
The IBI model identified a gradient in values related to conservation and protection status of the rivers (Fig. 2).IBI values were lower in unprotected rivers with obvious impacts.Biotic integrity class for all protected rivers except the Crucito was good.Mean IBI was higher in group G1 and it was significantly different from the mean for group G2 (Kruskall-Wallis = 16.19,P < 0.0001).The IBI showed high positive correlation only with river width (r 2 = 0.60, P < 0.0001) and negative correlation with total dissolved solids (r 2 = 0.71, P < 0.0001) and fine gravel substrate (r 2 = 0.38, P = 0.0024).
Principal component analysis for habitat variables produced three significant axes and explained 69.29% percent of accumulated variance (Table 2).The PCA 1 axis represented 32.56% of the explained variance and described a gradient determined positively by fine gravel and total dissolved solids, and negatively by width and branches and logs.For the second axis (22.71%) the gradient was dominated by sand, coarse gravel and stones.The third axis was determined by rocks and leaves.The recognized gradients (Fig. 3) demonstrate that in protected and better conserved rivers the presence of total dissolved solids, sand and fine gravel were smaller.These rivers were wider and deeper.The disturbance gradient was dominated mainly by the increase in the sand cover, fine gravel and total dissolved solids, better recognized in Crucito and Tupe stations.
CCA analysis showed a functional relationship between fish species, habitat variables and localities.The correlation species-habitat was significant (P = 0.0010) and elevated for all axes (Table 3).The three axes in CCA accounted for 40.3% of the variation in species abundances and habitat variables.The CCA 1 resumed 20% of the variance and was determined by the substrate particle size expressed as the sand and fine gravel, as well as by depth and total dissolved solids.In the rest of the axes the substrate particle size also predominated (sand and fine gravel) but their explained variances were low (CCA 2: 11.9%; CCA 3: 8.5%).The CCA plots (Fig. 4) can be associated with the geographic location of many rivers and localities.The diagram indicated a probable separation between rivers according to protection status (G1 and G2 groups) and the degree of watershed disturbance estimated initially.In the CCA the width, depth, fine gravel and total dissolved solids defined the fish assemblages.A gradient in the complexity of the aquatic habitat determined the composition of fish communities; that is, from communities with greater richness and environments dominated by more components (e.g.Zamuro river), to less complex communities with low richness and homogenous sub-strates (e.g.Tupe river).Nevertheless, intermediate assemblages also were detected; for example, the rivers Crucito and Guarataro showed communities with extremes in richness but with altered substrates.These results show that the habitat variables contribute moderately to the interpretation of the species richness and distribution.
The relationships between the environmental gradient with fish attributes were evident but only the axis PCA1 could partially explain species richness (r 2 = 0.34, P = 0.0030, Fig. 5A).The geographic location and protection status (G1 and G2 groups) for the rivers were associated with the PCA and CCA axes (Fig. 5B    degree of conservation and protection status for the rivers.Thus, the greater positive values of the PCA 1 and CCA 1 were recognized in rivers without protection by the national park. Variations in IBI scores were significantly attributable (r 2 = 0.61, P < 0.0001, Fig 6A ) to the variations in the environmental gradient shown by PCA analysis, suggesting that the IBI model could discriminate the changes in the fish communities and that it was independent of the environmental physical characteristics.The graphical interpretation suggests greater scores of the IBI associated with substrates with little cover of fine gravel and low concentration of total dissolved solids.On the other hand, smaller values of IBI were associated with a reduction in the river width and the presence of abundant wood (branches and logs).Even though the comparison of the IBI with the CCA (Fig. 6B) could be redundant, the correlation was only moderate.Our data indicate that the biotic integrity index and it the differences detected among studied fish communities are attributable to substrate heterogeneity and composition in fluvial environments.

Discussion
Differences observed among habitat variables and primary fish attributes from rivers in the Aroa Mountains were detected by simple comparisons, demonstrating differences in groups respect to habitat heterogeneity, fish assemblages, protection status and the degree of watershed disturbance associated with each locality.Also, the differences among fish communities were recognized with better refinement using both the IBI and multivariate analysis.Some fish community attributes (e.g.richness) were most sensitive to habitat variables, but the IBI were most responsive to overall gradients.The effect of environmental gradients and the degree of disturbance of the aquatic habitat on the IBI was evident in the rivers of the Aroa Mountains.The relationship between the IBI and habitat disturbance gradients have been demonstrated by many authors; for example, Toham & Teugels (1999) showed that variation in the IBI in relation to a gradient of human impacts, showed a reduction of intolerant taxa where major impacts were present.
The IBI was applied only in the dry season, when rivers have low water volume, lower sediment concentrations and higher transparency.Although this limits evaluation of the complete annual cycle, dry season sampling (apart from practical logistic limitations) should be the most appropriate time to apply the IBI because at low flow the parameters evaluated are more sensitive to changes in habitat variables (Pinto et al., 2006).The drastic variation of the rivers during the rainy season reduces the access to the sampling stations and the use of electrofishing standardized in time and space.
The IBI evaluates several attributes of fish communities but its main limitation is the sensitivity of the metrics measured.Some parameters have been shown to be poor or ambiguous predictors of habitat quality (Reynoldson et al., 1997).
In the absence of detailed information on life histories for many species, the baseline values used for with the IBI metrics must be obtained from local streams with low or no human  disturbances.Minimally disturbed rivers should be independently identified, preferably using habitat criteria; however a lack of basic information limits the verification or association with biological indicators.The consideration of local species abundance from protected ecosystems (e.g.national parks) is a useful reference condition for the construction of an IBI.Besides complementing the IBI, multivariate techniques permit incorporation of multiple gradients into the analysis of integrity (Gauch, 1982;Meffe & Sheldon, 1988;Reynoldson et al., 1997;Anderson & Willis, 2003).These provide a better understanding of biotic responses in the fluvial habitats.Many techniques are successfully used to measure or validate the IBI models, such as cluster analysis or discriminant function analysis.The principal component analysis is useful to recognize the environmental gradients but the canonical correspondence analysis is a robust technique to know the relationships between the habitat variables and fish communities.However, the IBI in this situation has the advantage to obtain similar predictive values with less information.For the purpose of natural resource management, IBI has greater utility, lower costs and is more easily appreciated, because the values of the IBI can be used with greater facility in biomonitoring programs.
Among the landscape factors that have been shown to influence the integrity of aquatic ecosystems, land use and riparian forest are good predictors of substrate and water quality (Bojsen & Barriga, 2004;Allan, 2004;Novotny et al., 2005).Among anthropogenic impacts, sedimentation is probably the most detrimental to fluvial habitats.Sutherland et al. (2002) and Casatti (2004) reported decreased fish species richness and changes in local abundance associated with increased sediment loads in rivers.At the local level (tributary or reach) other factors are critical to maintain the fish communities, such as depth and substrate composition (Schlosser, 1991) or the variability and fluctuations in the physical and chemical characteristics (Tejerina-Garro et al., 2005).
Protected areas, such as Yurubí National Park, help maintain flows and substrate heterogeneity but do not guarantee conservation of biota beyond their boundaries.Local disturbances (e.g.dredging, mining, deforestation) have strong and   lasting effects on aquatic ecosystems.Our study sites show this is the case in the Crucito River, which is protected in its headwaters, but highly deteriorated in the lowlands.Bond & Lake (2003) and Lake (2003) demonstrated how stream fish abundance was limited by insufficient habitat availability at a small spatial scale, a situation that applies to rivers that are dredged.The flexibility and robustness of this IBI make it an effective tool for management of hydrobiological resources.IBIs are already used in the regulatory arena in some countries (e.g.U.S. EPA, 2000).The European Community is considering biomonitoring as a fundamental part of management of hydrographic basins to develop a continental conservation strategy (Parlamento Europeo, 2000).In South America, a few preliminary studies using biotic integrity to estimate river conservation status have been published (see for example Rivera & Marrero, 1995;Rodríguez-Olarte & Taphorn, 1995;Hued & Bistoni, 2005;Silveira et al., 2005), but knowledge of taxonomy and ecology of the ichthyofauna is still incomplete.In the same context, several disturbances act on the fluvial ecosystems and their study and understanding are still developing, and the environmental effects on functional metrics of IBI have yet to be completely understood.In time, we can expect biomonitoring to play a greater role in natural resource conservation in the Neotropics.Appendix.Fish attributes considered for the development of IBI.Trophic guilds: carnivores (C), invertivores (I), piscivores (P), omnivores (O) and herbivores (H).Mesohabitats (MH) of occurrence: rapid (R), riffle (Ri) and pool (P).Apparent tolerance: tolerant (T), medium tolerant (MT) and intolerant (I).The abundance (G1/G2) is referred to the variation in means between groups (o: P = 0.01 differences; : P < 0.001; ∅: did not occur in G2 affluent).The species marked with "*" not were included in the IBI model.T. arleoi was not reported in affluents of group G1.

Fig. 1 .
Fig. 1.A:The Aroa River basin.The Yurubí National Park (YNP) is indicated by the shaded polygon.B: each circle indicates one or more sample localities.The white circles denoted tributaries in group 1 (G1) and black circles are tributaries in group 2 (G2).
and 5C), indicating an association between theCr Te Gu Za Ca Cu Tu

Fig. 4 .
Fig. 4. Pattern in fish assemblages distribution among environmental variables based on canonical correspondence analysis ordination.Labeled arrows indicate the direction along which each variable changes most.The arrows were doubled in length for visual purposes.The gradient is dominated by the variations of depth (DE), fine gravel (FG), total dissolved solids (TDS), sand (SA) and width (WI).Each symbol represents a sampling location and degree of disturbance: open circle (low disturbance), gray circle (medium), black circle (high) and triangle (very high).

Fig. 5 .
Fig. 5. Relationships between the PCA 1 and CCA 1 axes with the species richness and rivers stations.The richness (A) was moderately explained, but the multivariate (B) and canonical (C) analyses showed functional relationships with the geographical location of the rivers.

Fig. 6 .
Fig. 6.The correlations of the PCA 1 and CCA 1 axes with the index of biotic integrity explain with moderate correlation and high significance the recognizable variations in the IBI.

Table 1 .
Categories and metrics considered in the model of index of biotic integrity (IBI).Explanations in text.

Table 2 .
Eigenvalues and percent explained variance of the first three principal component axes associated with habitat parameters.

Table 3 .
Canonical correspondence values and percent explained variance associated with the first three axes for habitat and substrate variables measured.PCA analysis for the habitat variables in selected rivers.The graph includes all stations for each river.From environments with greater habitat heterogeneity (Carabobo and Zamuro rivers) the gradient is expressed A) towards rivers with greater sand cover and minor depth (e.g.Crucito: Pt stations) and B) towards rivers with greater total solid dissolved concentration and fine gravel (e.g.Tupe).Stations of a same river showed a varied composition of the substrate with respect to the disturbance gradient (e.g.Guarataro: Gu and St stations).