Geographical and physical characteristics of each sampling site at the Dos Novillos River (1-5) and the Mercedes Stream (6), Guácimo, Province of Limón, Costa Rica.
Costa Rica is an ideal reference point for global tropical ecology. It has an abundance of tropical forests, wetlands, rivers, estuaries, and active volcanoes. It supports one of the highest known species density (number of species per unit area) [1, 2] on the planet and possesses about 4 % of the world´s total species diversity . Because of its tropical setting, it also serves as an important location for agricultural production, including cultivars such as coffee, bananas, palm hearts, and pineapples. The country has also attracted more ecotourists and adventure travelers per square kilometer than any other country in the world .
The agrorural frontier on the Caribbean side of Costa Rica started to spread during the 1970s, especially in its northeastern area. Migrations of land-poor people from the Pacific and mountain areas of the country started to colonize the land that the government had made available [5, 6]. These waves of immigrants tended to establish themselves along river systems. In this way, towns, small to medium-scale family farming, ranching, and plantation agriculture began to base themselves along the main river systems. It was during this time that the human settlements originated along the Dos Novillos River .
Residual waters produced by all of the aforementioned human activities are at present discharged into the river systems in the Costa Rican Caribbean area. Households not situated in the neighborhood of rivers will use septic tanks; homesteads situated along riverbanks will discharge their effluents directly into the rivers. Other activities like the production of residual waters from dairy farms, pigsties, banana packing plants, plantations´ excess fertilization, etc. will drain eventually into a river. The Dos Novillos River is no exception. Water sewage systems in this part of Costa Rica are almost non-existent.
Aquatic biomonitoring in Latin America started at the end of the last century, commencing in Colombia [8-10] and then spreading to other Latin American countries. Revisions of the use of aquatic biomonitoring indices in Latin America are given in de la Lanza Espino et al. , Prat et al. , and Springer . In the case of Costa Rica, several studies based on the importance of macroinvertebrates for biomonitoring water quality and community structure, and function in banana and conventional and organic rice systems [14-22], as well as macroinvertebrate field-guides [23, 24] have been published. Some studies of macroinvertebrate assemblage structures also have been reported for rivers on the southern Caribbean coast of Costa Rica [25-34].
Although ecosystem studies have been used in Costa Rica for evaluating possible impacts caused by crop activities on river water quality at specific points, almost nothing is known about how a mix of other human activities can impact macroinvertebrate biodiversity and ecosystem structure along the length of a river. The aim of this study was to use macroinvertebrate biodiversity under the influence of different human activities along the length of a river in order to describe their impact on community structure and function in tropical agrorural environments.
2. Materials and methods
2.1. Study area
This study was conducted at the Dos Novillos River (Figure 1, Table 1) in the province of Limón, Costa Rica. Samples were collected two times a month from January 2005 to March 2005 and monthly from April 2005 to January 2006. The Dos Novillos River drains from the Central Cordillera at an elevation of approximately 2380 masl towards the Caribbean lowlands of the province of Limón. This river is part of the 2950 km2 Parismina River watershed in a premontane wet forest and tropical moist forest region . The underlying geology is represented by quaternary sedimentary and volcanic rocks under the influence of nearby volcanoes, with a flat to undulating topography and poorly drained alluvial soils susceptible to flooding . Banana plantations have been developed on the lower reaches of this watershed. The study area is characterized by a humid tropical climate with a mean temperature of 25.8 °C, an average annual relative humidity of 87 %, and an annual precipitation average of 3460 mm ± 750 mm without a pronounced dry season. The sampling area runs in a straight line through the towns and localities of La Argentina, Pocora, and EARTH University in the Province of Limón. Each fore mentioned landmark is separated by approximately 4 km. The total sampling area runs along a length of 13.6 km, across an area with a mixture of a premontane wet forest, pastureland, small town, riparian tropical moist forest, and banana agricultural areas.
Six sampling sites were located along the Dos Novillos River (Figs. 1-7; Table 1) where macroinvertebrates were sampled. The first sampling site (Figure 2) (site 1, “Don Eladio”) served as a reference site, being part of the rhithral region, located upstream of the first anthropogenic disturbance (pastureland). This site is surrounded by tropical rain forest; therefore, natural good water conditions were expected, as well as high taxa richness and an assemblage composition dominated by pollution-sensitive organisms. Site 2, “La Argentina” (Figure 3), was located approximately 5.5 km downstream from site 1. This site was selected to examine the possible extent that small livestock farming might have on river water quality. Site 3, “Chiquitín” (Figure 4), was located approximately 8 km downstream from site 1. High anthropogenic influence was expected at this site because the houses situated at the riverfront discharge their grey and black waters directly into the river. Site 4, “Puente La Hamaca” (Figure 5), was located within the property of EARTH University, approximately 2 km downstream from site 3. As the intervening river length between sites 3 and 4 runs through forest areas, site 4 was selected in order to examine if water quality was improved by a forest filtering processes. Site 5, “Desembocadura” (Figure 6), is within the EARTH University campus and is located approximately 500 m upstream from the confluence of the Dos Novillos and Parismina Rivers. Site 5 was selected to analyze the impact of banana plantations on river water quality. Site 6, “Quebrada Mercedes” (Figure 7), is one tributary of the Dos Novillos River flowing through a forested area within the EARTH University campus, approximately 2 km downstream from site 3. Site 6 was chosen to examine if the intermittent discharge of a small drain of water used to wash bananas in a packing plant had any effect on the stream. Table 1 indicates the exact location, depth, width, and current conditions for each site.
|Longitude (N)||10° 07´ 09.7´´||10° 09´ 14.3´´||10° 10´ 40.8´´||10° 13´ 00.9´´||10° 14´||10° 12´|
|Latitude (W)||83° 39´ 15.2´´||83° 37´ 24.7´´||83° 36´ 10.6´´||83° 35´ 18.4´´||83° 34´||83° 35´|
|Depth (m)||0.3-1.3||0.25-0.8||0.25-0.8||0.2-0.86||0.2 – 1.2||0.15-0.3|
A plastic strainer with a diameter of 20 cm and 0.5 mm mesh size, and tweezers were used for directly collecting macroinvertebrates. The main criterion for this semi-quantitative collecting method is time; there were no defined sampling areas. All types of microhabitats present at a particular site were examined equally for the macroinvertebrates. Collected organisms were fixed immediately in 70 % ethanol at the time of sampling. Exact details of the sampling methodology can be found in Stein et al. .
A presampling was carried out in order to determine sampling time . Out of the achieved results, an accumulated taxa curve was elaborated, and 120 min were determined to be a representative sampling time per site. Sampling took place during early morning and under normal current situations in order to avoid the negative effects of flooding and high water conditions.
The government of Costa Rica has officially suggested this method under the water quality monitoring regulation . This method is coupled with the use of a modified Biological Monitoring Working Party index for Costa Rica (BMWP´-CR), a biotic index utilized to define different levels of water quality. Each family of macroinvertebrates has a sensitivity value ranging from 1 to 10, reflecting tolerance to pollution based on the knowledge of distribution and abundance. The values for each family are then summed up independently from abundance and generic or species diversity. Sensitivity scores higher than 120 points indicate undisturbed aquatic ecosystems, while low values indicate serious contamination (mostly organic) of the environment [13, 23, 38,].
2.3. Data analysis
The analyzed data comprised the values of the physical-chemical water quality variables (Table 2) and macroinvertebrate abundances (Table 3) during a collecting period of 13 months (January 2005 to January 2006).
The following physical-chemical variables were measured: pH, temperature, O2, O2 saturation, suspended solids, turbidity, conductivity, NO3-, NH4+, PO4+, BOD, and COD. The water samples also were analyzed for the following agrochemicals associated with banana production using gas chromatography-MS and liquid chromatography-PDA: Chlorpyrifos, Diazinon, Dimethoate, Edifenfos, Etoprofos, Fenamifos, Malathion, Parathion-methyl, Parathion-ehtyl, Terbufos, Difenoconazol, Propiconazol, Imazalil, Ametrine, Atrazine, Hexazinone, Terbutylazine, Bromacil, Bitertanol, Chlorothalonil, and Thiabendazole. However, no traces of them could be detected in the river water samples. This does not come as a surprise because in order to monitor pesticides very frequent sampling would be required to detect peak concentrations during pesticide application periods , whereas low concentrations are very difficult to detect.
Following the suggestion of Ramírez and Gutiérrez-Fonseca , this study is also undertaking an ecosystem process analysis of the functional feeding groups (FGG) of the aquatic macroinvertebrates. This sort of analysis is based on two key aspects of macroinvertebrates: morphological characteristics related to the obtainment of food resources (e.g., mouthparts and related structures) and behavioural mechanisms (e.g., feeding behaviour). FGG is a very useful tool that provides valuable information on ecosystem functioning, facilitating stream ecosystem comparisons, and avoiding the traps of gut content analysis, which is more appropriate for assigning trophic guilds .
In this ecosystem study, the use of different parameters of the structure and composition of macroinvertebrate assemblages are presented: total and relative abundances, taxa richness, and functional feeding groups. Also, correlations of different genera and functional feeding groups with environmental variables were analyzed.
2.4. Statistical analysis
The model comparisons between physical-chemical variables, macroinvertebrate abundances, and the BMWP´-CR index at different collecting sites was done by performing an analysis of variance (ANOVA; α=0.05). For all three cases, the proposed hypothesis is to test the existence of significant variable differences between sites. Abundances were square root transformed in order to comply with error normality. Evaluation of the best model (homocedastic or heterocedastic) for each variable was performed using the Akaike information criterion (AIC), which is one of the benchmarks of mixed models based on penalized likelihood [41, 42]. When the model detected significant differences, a DGC (Di Rienzo – González – Casanoves) statistical test was performed for the comparison of means .
On the other hand, taxonomic groupings and FFG relative frequencies analyses were done using a Chi-square test in order to assay for statistically significant differences between sites. The Chi-square analysis is testing for independence between the sites and the studied variables. Any p value below 0.05 shows the existence of an association between the site and the studied variable.
PLS regression is a technique that combines Principal Component Analysis and Linear Regression . It is applied when it is desired to predict a set of dependent variables (y), in this case the abundance of macroinvertebrate genera and FFG abundances and the BMWP index values, from a set of predictor variables (x), in this case physical-chemical variables. To represent the results obtained from the PLS analysis, a Triplot graph was superimposed on a Biplot graph , thus correlating all variables. Then, the observations appear ordered in a Triplot graph (sites), depending on the values of the dependent variables (macroinvertebrate and FGG abundances and BMWP´-CR index) and their correlation with the predictor variables (physical-chemical water-quality variables). For the macroinvertebrate genera PLS analysis, out of the 127 collected taxa (of which, 123 could be identified to the genera level and their different developmental stages: larva, pupa, adult), the multimetric analysis included only 58 taxa, which were chosen using a PCA (Principal Component Analysis). The rest were characterized for repeating the same information. These 58 taxa, composed of 15 688 individuals, showed high projection values on the first two principal components. All statistical analyses were done using the InfoStat program .
In order to correlate the abundance of macroinvertebrates, FFG, and the BMWP´-CR index values of different sites, these variables were correlated with physical-chemical variables using the Spearman rank correlation coefficient. In the present case, the hypothesis tries to establish if one variable can be effectively substituted by another one, due to the existence of a significant correlation. The Spearman correlation coefficient was selected, versus Pearson, because its use is recommended in the case of having a small sample.
Finally, in order to evince the consistency of the sites´ congruences arranged in one plane unto physicochemical variables and FFG and macroinvertebrate genera abundances, a Generalized Procrustes Analysis was performed. This analysis is used for harmonizing multivariate configurations obtained on the same set of observations with different types of variables or time points . Alignment is performed through a series of steps including normalization, rotation, reflection, and scaling of data to obtain a consensus array between groups of variables. This series of steps should maintain the distances between individuals from the individual configurations and minimize the distance between similar points . The result of this multivariate method is to present a graph that displays the configurations arrived at by each variable type and the consensus configuration. A percentage consensus analysis was also undertaken. A high consensus indicates that any group of variables characterizes the different sites in the same way; therefore, using any group of variables is indistinct for site characterization.
3.1. Physical-chemical parameters
Table 2 shows the results of the physical-chemical analysis (ANOVA, DGC-test, p>0, 05). All sites presented neutral pH and high dissolved oxygen levels and saturation. Temperatures varied in a statistically significant way, site 1 being the coolest place and site 5 the warmest. Conductivity was quite low at all sites, but statistically significantly lower in sites 5 and 6; whereas, turbidity was statistically significantly higher in site 6. NO3 was statistically significantly higher in sites 5 and 6, and lowest in sites 1, 2, and 3.
3.2. Diversity and composition of macroinvertebrate assemblages
The study collected a total of 17 163 specimens, distributed in the following 6 classes (number of total amount of specimens in parentheses): Clitellata (2), Turbellaria (10), Gastropoda (403), Arachnida (21), Malacostraca (40), and Insecta (16 706) with the following orders: Ephemeroptera (6299), Coleoptera (3150), Trichoptera (3010), Diptera (2868), Plecoptera (683), Odonata (435), Hemiptera (101), Megaloptera (91), Lepidoptera (64), Blattodea (2) (Table 3). The total abundance analysis (Figure 8) shows that Ephemeroptera, Coleoptera, Trichoptera, and Diptera were the most abundant groups, comprising 89.2 % of the collected macroinvertebrates.
|Site 3 Chiquitín||Site 4 Hamaca||Site 5|
|Trombidiformes||Various families||Various genera||Pr||12||2||2||3||0||2||21|
3.3. Comparison between sampling sites with different agrorural influence
If we consider the mean total number of collected individuals per site, clear abundance differences are evinced (Figure 9). The greatest mean total numbers, with almost 300 individuals collected per sampling date, are found in the two least impacted sites (ANOVA, DGC-test; p>0.05),the reference (site 1) and the livestock-pasture site (site 2), whereas the sites under more intense human influence present statistically significantly reduced mean total values (around 150 individuals collected per month).
Table 4 presents a family, genera, and EPT richness analysis of the different collecting sites, as well as the mean BMWP-CR values and resulting water quality. The total family and genera richness does not show significant frequency differences between sampling sites, although the highest number of genera were found at the first two sites, with over 60 genera. The site with the lowest taxonomic richness, both on family and genus level, was site 5, close to the river mouth. The EPT taxa richness was highest at site 1 (14 EPT taxa), decreased downstream towards only nine EPT taxa at sites 5 and 6, although no statistically significant frequency differences were found using the Chi-square test. The BMWP’-CR index shows statistically significant differences (ANOVA, DGC-test; p>0.05), where the index values diminish according to the following site groupings: (sites 1-2)-(site 6)-(sites 3-4)-(site 5). The water quality, indicated by the mean BMWP’-CR index falls into the categories “good quality” (sites 1-2) and “regular quality” (sites 3-6).
|Site 5 Desembocadura||Site 6|
On the other hand, a very clear change in the structure of macroinvertebrate assemblages can be observed along the different sampling sites of the Dos Novillos River (Figure 10). At the first site, the undisturbed sampling point, beetles are significantly more abundant than at the other sites (Chi-square, p<0.0001), which are all under the influence of human impact. The last site (site 6) shows a significantly greater abundance of Gastropoda, Odonata, and Plecoptera (Chi-square, p<0.0001). Here the banana packing plant discharges effluents, carrying small banana pieces and other suspended organic material, into the water of the stream. Also, a tendency of an increase in mayfly (Ephemeroptera) abundance can be observed towards sites 4 and 5, while caddisfly abundance is relatively steady throughout all sampling sites.
3.4. Functional feeding group analysis
The FFG analysis is presented in Figure 11 (Chi-square; p<0.0001). Collector-gatherers are the dominant group at all sites along the main river, with around 50% of all individuals collected; at Quebrada Mercedes (site 6), they are the second largest group after the filter-feeders. Filter-feeders show also high percentages at each site, with exception of site 1 (reference site), where they have statistically significant lower relative frequencies. One the other hand, piercers have a statistically significant greater relative frequency at this site, compared to the rest of the sampling sites. Finally, site 6 has a statistically significant lower relative frequency of collector-gatherers and greater relative frequencies of predators, scrapers, and shredders.
The FFG abundance PLS (Figure 12) presents high variation explanatory values, with factor 1 explaining 39.1 % of the total variance and factor 2 explaining 25.4 %. The variables´ temperature, as well as the piercer and predator abundances, mainly separates the different sites on the horizontal projection, as does the shredder abundance on the vertical projection. The FFG most closely allied to the reference site, is piercer abundance. On the other hand, the most closely allied variables to the agriculturally impacted areas, sites 5 and 6, are the variables turbidity, conductivity, and NO3.
3.5. Relationship between macroinvertebrate assemblages, physical-chemical parameters and sampling sites
The macroinvertebrate abundance partial least square analysis (Figure 13) presents high variation explanatory values, with factor 1 explaining 42.0 % of the total variance and factor 2 explaining 29.9%. The variable NO3, and to a lesser degree conductivity and turbidity, separate mainly the different sites on the horizontal projection, as does suspended solids on the vertical one. The genera most closely allied to the reference site are: Brysopterix, Heterelmis, and Paltostoma. On the other hand, the most closely allied genera to the agriculturally impacted areas, sites 5 and 6, are larvae of the genera Caenis, Dythemis, and Hexatoma. Interestingly, the Spearman rank analysis has resulted in a plethora of correlations. The strongest physical-chemical/biological correlations are presented in Table 5.
|Variable (1)||Variable (2)||Spearman rank||p-value|
The Procrustes analysis produced a good site ordination consensus based on macroinvertebrate abundances, FFG, and the physical-chemical variables. The first axis explains 45.8 % of the variance and the second axis explains 33.8 % (Figure 14). The proportional consensus percentages are also good, ranging from 76 % to 92.4 %, with a mean value of 82.9% (Table 6). One can conclude that the site ordination has a good congruence with the biological and physical-chemical data sets.
An important consideration for the present study is the ability to distinguish between natural variability and human impacts . In the present case, all sites are located along the same river, spanning only a small distance of 10 km (Figure 1). Moreover, elevation, stream size, and surface geology are relatively similar (Table 1) in order to properly assess human impacts and reduce natural gradients.
4.1. Diversity and composition of macroinvertebrate assemblages
The present study recorded the existence of 16 macroinvertebrate orders. In an analysis undertaken by Castillo et al.  very near to the present collecting localities, the authors reported the existence of 15-16 macroinvertebrate orders in their reference sites and 12-16 orders in the banana plantation sites. The results of the present study fall within the order range for similar studies in nearby regions. The present analysis resulted in 53 collected families (Table 3). This number compares well with the number of families collected by Lorion and Kennedy  and O’Callaghan and Kelly-Quinn  in Costa Rican and Honduran neotropical rivers, who reported 56 and 60 families, respectively. The present study also identified a total of 98 genera (Table 3). Montoya Moreno et al.  reported 69 genera for the Negro River in Colombia; whereas, Sánchez Argüello et al.  reported 96 genera in their study in Panama. Considering that some groups, such as water mites, were not identified to genus level in the present study, the total number of genera present in the Dos Novillos river and its tributaries is likely to be over 100 genera, reflecting a very high taxa richness.
It is interesting to compare the dominance results obtained in this study with analyses made in neighbouring areas under similar ecological conditions. Ramírez et al.  sampled the Carbón and Gandoca Rivers and found the following dominance gradient for total abundances: Ephemeroptera-Diptera-Trichoptera-Odonata. Castillo et al. , sampling on sites very near to the present ones, found the following order dominance in their reference sites: Ephemeroptera-Trichoptera-Coleoptera and the following one in banana plantation sites: Ephemeroptera-Diptera-Coleoptera-Gastropoda-Trichoptera. Lorion and Kennedy  studied several streams in the Sixaola River Valley. Considering the total number of individuals, they found a diminishing abundance gradient as follows: Ephemeroptera-Diptera-Coleoptera-Odonata. Gutiérrez-Fonseca and Ramírez  have reported at La Selva Biological Station, in a 15-year study, the following dominance sequence in unpolluted streams: Diptera-Trichoptera-Odonata. The present study (Figure 8) has found the following abundance dominance order: Ephemeroptera-Coleoptera-Trichoptera-Diptera. It would seem from these results that dominance sequences might vary depending on several factors dependent on the collecting site, such as substrate, current, and water quality, but also can be a result of different sampling device and mesh size . However, Ephemeroptera would appear to be the most constant dominant group in most lowland rivers and streams in Costa Rica.
4.2. Comparison between sampling sites with different agrorural influence
The analysis of abundance differences demonstrates that less impacted sites clearly show statistically significantly higher abundances than sites under stronger human influence (Figure 9). Similarly, Paaby et al.  and Lorion and Kennedy  detected greater abundances in forested areas versus pastures under neotropical conditions; whereas, the study by Ramírez et al.  detects greater abundances in Costa Rican tropical riffle habitats than in other habitats.
Taxonomic composition also varies in a clear way along the sampling sites under the influence of different ecological impacts (Figure 10). The reference site, arguably not under the influence of human impact, presented a significantly greater relative abundance of Coleoptera. This is due to the high amount of individuals collected from the riffle beetle family (Elmidae), which are especially diverse and abundant in well oxygenized rivers and streams in forested areas (Springer, unpubl). A study of rivers in the Guanacaste area in Costa Rica (Kohlmann, unpubl.) has also found very high numbers of Coleoptera in unpolluted rivers. At the other end of the scale, the banana packing plant discharge, carrying banana debris and showing statistically significantly higher values of conductivity, turbidity, and NO3, presents statistically significantly high numbers of Gastropoda, Odonata, and Plecoptera. This is very interesting because stoneflies (Plecoptera) have always been considered as good indicators of oxygenated, clean, and cool running waters [51-54]. Additionally, the triplot (Figure 12) and Spearman rank correlation analyses show a strong negative correlation between stonefly (Anacroneuria) abundance and temperature values (Spearman rank correlation value=-0.89, p=0.0476). Gutiérrez and Springer  reported the widespread species A. holzenthali from coffee plantations in Costa Rica. Tomanova and Tedesco , as well as Thorp and Covich  also indicated that stonefly presence is not necessarily a sure sign of water cleanliness. The results from these analyses seem to support this claim.
The EPT taxa richness (Table 4) showed a steady decline along the collecting sites, following an increasing anthropogenic impact trend, similar to what Lorion and Kennedy  reported following a forest-forest buffer-pasture gradient. The number of families per site (Table 4) did not vary much (33-40) along the collecting transect in this study (Table 4). An analysis by Castillo et al.  reported a family number variation going from 39 to 47 families. Similarly, their family numbers did not vary much between their reference (46-47) and their banana plantation collecting sites (39-46). Some families (Table 3) were restricted to only a single site (site number in parentheses), like (1): Anomalopsychidae, Calamoceratidae; (3): Ancylidae, Lumbriculidae, Physidae; (4): Corduliidae; (5): Gyrinidae; (6): Ampullaridae, Chironomidae, and Tabanidae. Certain genera (Table 3) appeared only once in a specific site: (1): Chloronia, Contulma, Dythemis, Hemerodromia, Ochterus, Palaemnema, and Polycentropus; (2): Hebrus, Leptohyphes, Mayobaetis, and Psephenops; (3) Stenonema; (4): Anchitrichia and Neocylloepus; (5): Haplohyphes; and (6): Chironomus, Chrysops. These taxa would seem to be closely associated with the ecological conditions of each site, e.g. Ampullaridae, Chironomidae, and Tabanidae in a banana packing-plant effluent with a high organic waste discharge versus Anomalopsychidae and Calamoceratidae in a forested undisturbed condition. Interestingly, the number of families and genera do not show the same trend among sampling sites. The reference and slightly disturbed sites (site 1 and 2, respectively) show a higher number of genera (over 60), in comparison to the more influenced sites (with around 40 to 50 genera), even though the number of families were similar at the latter, and even higher in one case (site 3). These results suggest that family richness is not necessarily an adequate indicator for biomonitoring, and generic identification is necessary to achieve results that are more reliable.
The BMWP index in its different variations has been popularly employed in Latin America. These studies usually have found this index to be satisfactory for reflecting water quality [10, 12, 22, 60-63], especially the Costa Rican adaptation . Sánchez Argüello et al.  undertook in their Panamanian study a comparison between the Colombian and Costa Rican adaptations of the BMWP index and found the latter to be more unforgiving in its water quality evaluation. Rizo-Patrón et al.  undertook an analysis using the BMWP index modified for Costa Rica, studying the environmental impact caused by conventional and organic-irrigated rice fields on the macroinvertebrate communities. Their BMWP´-CR results show that the index values were greater in the organically irrigated rice fields. On the other hand, Fenoglio et al.  recommend, from their experience in Nicaragua, the use of the Indice Biotico Esteso  because of its ease of use and low cost. However, these studies have mostly assessed the comparative performance of the various indices; no attempts were made to correlate the BMWP index to specific physical-chemical variables.
The results of the BMWP’-CR index for the present analysis (Table 2) do indeed show a discriminating capacity of the index following a diminishing environmental quality site trend, especially under agricultural impact conditions (site 5), but it also shows a tendency of reporting a higher value when in river waters with high organic pollution (site 6). Interestingly, there is a statistically significant negative correlation of the BMWP’-CR index with temperature, not with pollution variables, as one could expect from the general assumption that the BMWP index reflects organic water pollution quality. These results generate some doubts about the reliability of the BMWP’-CR index as an environmentally representative tool, as the following studies indicate. Sermeño Chicas et al.  tried to implement the BMWP-CR index in El Salvador where rivers showed consistently high organic pollution conditions that were not reflected by the BMWP index . In their analysis of selected macroinvertebrate-based biotic indices in Honduras, O´Callaghan and Kelly-Quinn  found that a BMWP-CR-based version of the ASPT index performed much better than the aforementioned index. Without doubt, more studies will be necessary in order to adjust the biotic indices used for aquatic biomonitoring in Costa Rica and Central America according to the different ecoregions.
4.3. Functional feeding group analysis
FFG relative abundances also change significantly depending on the human impact conditions on the quality of river water. It would seem that under undisturbed conditions filterers´ relative abundances tend to be minimal, their increase at disturbed sites might be a result of higher dissolved organic matter. In this study, under conditions of high organic pollution, shredders, scrapers, and predators tend to have maximal relative abundances while collector-gatherers tend to have minimal values. Finally, filter feeders seem to react positively to high concentrations of dissolved O2 (Table 5), which is positively correlated with fast flowing waters, a condition that also favors the feeding mechanism of filterers.
The taxonomic grouping triplot analysis (Figure 13) suggests a correlation between the reference site and the piercers, and it would appear that the first axis is characterized by a piercers’ abundance gradient, diminishing from the reference community (site 1) to the high organic waste discharge sites (site 6). In the present study, piercers are mainly represented by the caddisfly family Hydroptilidae, which is especially abundant in the splash zone of big rocks in riffle areas, which were characteristic for the reference sampling site. The second axis appears to be characterized by filterers, arranging the different agrorural ecosystems along a diminishing abundance gradient.
The FFG triplot analysis (Figure 12) supports a strong agrorural ecosystem ordination process mediated by the abundance of piercers along the main axis as suggested in the previous triplot analysis; however, as this analysis benefits from having more information (127 taxa versus 58 taxa), it also evinces the importance of predators and shredders as relevant ecosystem ordinating biological variables. The strong negative correlation between the piercer’s abundance and NO3 and PO4 values (Table 5) stresses again the importance of the piercers as an ecosystem characterizing variable, although the presence of suitable microhabitats might be another important factor to consider.
4.4. Relationship between macroinvertebrate assemblages, physical-chemical parameters and sampling sites
The taxonomic grouping triplot analysis (Figure 13) and the Spearman rank correlation analysis (Table 5) show the existence of several genera and FFG that are highly correlated with physical-chemical variables and that possibly could be used as surrogates for these variables. The larvae of Maruina showed one of the strongest correlations with NO3, although other genera like Farrodes and Cora, and the piercers´ functional feeding group were also significantly correlated with this chemical variable. Species of the genus Caenis have been found regularly in organically enriched streams . Of the statistically significant genus list (Table 5), only Heterelmis and Farrodes already had been cited before as good quality bioindicators for toxicity and pollution-sensitivity testing by Castillo et al.  and Rizo-Patrón et al. , respectively, in a similar type of analysis.
Finally, the Procrustes analysis allows an assessment of the goodness of fit of the taxonomic, FFG, and physical-chemical analyses (Figure 14, Table 6). The analysis resulted in a good match of the collecting sites (landmarks) with the values derived from the three blocks of variables, indicating that any one of the three is describing in a similar way the ecology of each study site. The consensus values indicated in Table 6 show a very good concordance in this study between the different variable blocs and the consensus values, where generic abundances seem to generate a better ecological ordination. An environmental impact gradient also becomes very apparent on the first ordination axis, ranging from the undisturbed reference site inside the forest on the right of the graph to the banana packing-plant effluent on the left.
The present study clearly shows that tropical river macroinvertebrate diversity changes and at the same time characterizes and defines different river ecosystem conditions under various agrorural impacts. Changes in taxonomic composition and functional feeding group structure are very indicative of ecosystem function.
Ephemeroptera seem to be, in general, a rather constant and numerous group, present in the great majority of collections in neotropical rivers. However, high relative abundances of Coleoptera, especially from the family Elmidae, seem to be indicative of unpolluted conditions in tropical rivers; whereas, high relative abundances of Gastropoda, Odonata, and Plecoptera show up in sites with relatively high (plant-derived) organic pollution, although with well-oxygenized waters and forested stream margins. Addtionally, high numbers of individuals were found in unpolluted or slightly polluted sites; whereas, lower abundances were found in sites under human impact (town and fruit packing plant discharges, agricultural plantations).
Especially illustrative is the change in structure of functional feeding groups. Piercers showed the highest relative abundance in unpolluted sites and seem to be especially sensitive to human impacts because they quickly disappear under altered conditions, most probably reacting to the decrease of their microhabitat and food items. Filterers have the lowest relative abundance under unpolluted conditions and quickly become relatively more abundant under human impact conditions, most probably reflecting an increase in particles in the water column. At the other end of the spectrum, under high (plant-derived) organic pollution, shredders and scrapers show their highest relative abundance concomitantly with an increase in particulate organic matter, and probably as a response to it. Predators here show their highest relative abundance as well and at the same time, collector-gatherers here show their lowest relative abundance. It would seem then that predation, as well as scraping and shredding, increases significantly in river areas with high (plant-derived) organic pollution.
The biomonitoring analysis presented in this chapter used an adaptation of the BMWP index to Costa Rica. In all cases, the method revealed or indicated the existence of an anthropogenic/agricultural impact gradient going from the unpolluted site to the most perturbed locality. ANOVA tests evidenced the fact that the BMWP index has enough sensitivity and discriminating power to detect changes in macroinvertebrate biodiversity, which can be translated into statistically significant differences between sampling sites. The method determined the existence of changes in macroinvertebrate communities associated with agricultural areas, even when analytical methods could not detect the presence of pesticides in river water. A previous analysis of the BMWP score by Pinder and Farr [68, 69] found that it was significantly negatively correlated only with dissolved organic carbon. The present study detects a strong negative correlation between the BMWP score and temperature (Figure 12, Table 5). Due to its simplicity, speed of use, efficiency, and cheapness, this method shall undoubtedly continue to be a very popular one in the future. The BMWP is considered an extremely successful index according to Spellerberg . However, there are some doubts about its suitability as a tool for detecting organic pollution in some regions of Latin America, especially if one considers that in this case the index was more sensitive to the impact of a banana plantation than the one caused by relatively high (plant-derived) organic pollution. The PLS/Procrustes analysis seemed, on this occasion, to be a more suitable method for describing and evaluating anthropogenic/agricultural environmental impacts.
The use of multivariate ordination for environmental studies is becoming more and more common. In particular, the PLS/Procrustes analyses represent a very powerful combination tool as they not only perform a site ordination but different taxa and environmental variables can be correlated at the same time. This makes this method extremely useful for taxa, physical-chemical, and FFG variables’ correlations as specific environmental variable surrogates. Castillo et al.  and Rizo-Patrón et al.  found similar promising results for the study of agricultural ecosystem analyses. Castillo et al.  indicated, based on their results, that multivariate analyses are more sensitive in distinguishing pesticide effects than toxicity tests. Therefore, multivariate analyses should be incorporated as an approach for future ecosystem/biodiversity analyses.
We would like to thank Mildred Linkimer for her help in the design of some figures; as well as Roxana Araya for her help in the composition of tables and lists. Finally yet importantly, we are also indebted to Jane Yeomans and Kent McLeod for critically revising the text and to Pablo Gutiérrez for useful suggestions in organizing the results. This material is based upon work supported by the Department of Energy [National Nuclear Security Administration] under Award Number DE FG02-04ER 63856.