Introduction

During recent decades, aquatic macroinvertebrates have been widely used as indicators of impacts on freshwater ecosystems worldwide and increasingly in Latin America (^{Armitage et al., 1983}; ^{Hellawell, 1986}; ^{Alba-Tercedor & Sánchez-Ortega, 1988}; ^{Rosenberg & Resh, 1993}; ^{Weigel et al., 2002}; ^{Roldán, 2003}; ^{Juárez & Ibáñez, 2003}; ^{Acosta et al., 2009}; ^{Ferreira et al. 2011}; ^{Sánchez et al., 2012}; ^{Rios-Touma et al., 2014}; ^{Ramírez & Gutiérrez-Fonseca, 2014a}). Although several studies have dealt with the ecology of aquatic macroinvertebrates in the Central American streams, including ^{Wolda and Flowers (1985)}, ^{Flowers and Pringle (1995)}, ^{Astorga et al. (1997)}, ^{Ramírez and Pringle (1998)}, ^{Paaby et al. (1998)}, ^{Umaña (1998)}, ^{Boyero and Bosch (2002)}, ^{Ramírez et al. (2006)}, ^{Stein et al. (2008)}, ^{Vásquez et al. (2009)}, ^{Chaves-Ulloa et al. (2014)}, more research is needed to obtain the required level of knowledge about ecology and taxonomy (^{Springer, 2008}; ^{Ramírez & Gutiérrez-Fonseca, 2014a}) to assess the effects of global change scenarios on freshwater ecosystems in the region.

One of the main drivers of aquatic communities is climate (^{Hynes, 1970}). Many tropical areas are from moderately to highly seasonal with respect to rainfall. Rainfall seasonality occurs both in short term cycles, with strong events within the year and in long-termsones, such as “El Niño Southern Oscillation” (ENSO) phenomena. Macroinvertebrate assemblages have adapted to such variations through natural history. Some studies have shown that aquatic insects exhibit seasonal fluctuations even in non-seasonal tropical areas (^{Wolda & Flowers, 1985}). These patterns are increasingly affected by the alteration of natural flow regimes resulting from dams and hydropower production. Several studies show the effects of hydropower production on Neotropical rivers (^{Pringle et al., 2000}), highlighting the importance of using multiple measures of macroinvertebrate assemblage structure for assessing this type of environmental impact (^{Chaves-Ulloa et al., 2014}).

On the other hand, land use changes in mountain areas, such as clearing of native vegetation, lead to changes in magnitude, frequency, duration, and seasonality of flows, as well as increased sediment input to streams, which in turn change aquatic and riparian habitats (^{Wohl, 2006}; ^{Mancilla et al., 2009}). Cattle areas, especially dairy-intensive farms, produce forceful impacts on mountain tropical streams, with higher dissolved solids, lower dissolved oxygen and, thus, pollution-tolerant aquatic communities (^{Giraldo et al., 2014}). Understanding changes in tropical streams in response to land use impacts is a priority for tropical ecosystems conservation and management (^{Boyero et al., 2009}).

However, there are few data regarding the composition and structure of aquatic communities in tropical volcanic highlands. In this study, a two-year sampling effort led to the characterization of aquatic macroinvertebrate assemblages in the Birrís River Basin (Costa Rica), located at the southern slope of the Irazú Volcano. The hypothesis that volcanic headwater communities are highly vulnerable to human impacts was tested, as volcanic ash generates highly permeable soils and, thus, low baseflow stream discharge in these areas.

The main hypothesis is that volcanic headwater communities are highly vulnerable to human impacts, so remaining taxa are resilient to them. Thus, objectives of this research are: (1) test the ecological responses of macroinvertebrates assemblages under human pressure in highland tropical headwaters; (2) evaluate the resilient responses of macroinvertebrate assemblages to degradation and rehabilitation gradients in these environments; and (3) analyze the recovery potential of these ecosystems.

Material and methods

Study area. The Birrís River Basin is located in the highlands and southern slopes of the Irazú Volcano, 3400 m.a.s.l., an active volcanic cone located in the Central Volcanic Range of Costa Rica (Fig. 1). The basin was selected as it is representative of the high-priority water producing mountain headwaters in the tropics. This study is focused on the two main streams, the Birrís River, with a gradient of 2139 m over its 15.4 km distance, and its tributary Pacayas Creek (or Quebrada Pacayas) with a gradient of 1750 m along a distance of 11.2 km. Such high gradients are the origin of high slope values in fluvial reaches ranging from 7.7% to 11.1%. Annual rainfall averages 2300 mm caused by cold polar fronts (October to January) and convective storms (May to October) with a relatively dry season from January to May. During wet season flash-floods occur, due to the combination of steep slopes and heavy rainfall.

Figure 1 Spatial location of the studied area and sampling stations in the River Birrís Basin (Costa Rica). 

Nine sampling sites were located in this basin to assess the responses of aquatic communities to increasing human pressure balanced with higher self-purification capacity. Two stations were located in the less impaired Birrís River and Quebrada Pacayas headwaters (B1 and P1). The rest (B1’, B2, B3, B4, P2, P3, P4) were located in an increasing pressure design (Fig. 1, Table 1). Sampling campaigns started in February 2002 (1^st dry season) and were repeated bimonthly until October 2003 (2^nd wet season), with a total of 11 campaigns. Upper reaches (B1, P1) have fine silty and sandy substrates, because of the constant input of volcanic ash from previous eruptions (1963-1965). Forest cover and riparian forest corridor decrease downstream in both the Pacayas and Birrís streams (Table 1). Fieldwork evidence indicated that the main sources of contaminants in the Birrís River Basin were the intense agricultural and cattle-raising activity as well as the discharge of human sewage. The substrate was characterized in the field following the methodology proposed by Wolman (1954), consisting of 100 random surveys in the channel.

For each sampling station, the hydrologic drainage basin was determined using ArcSwat software. Relevant parameters, such as the percentage of land cover types (forest, pasture, crops, volcanic ash, and urban), upstream population, and cattle density were calculated using ArcView GIS for each sub-basin.

Macroinvertebrate sampling. Benthic macroinvertebrates were collected by aggregating two samples: river bed SURBER sampling (1-square foot sampling area, three samples per station and campaign) and multi-habitat sampling (15 minutes of total collecting time, using a hand-net). Ethylic alcohol (70%) was used for preservation. Benthic invertebrate samples were sorted in the laboratory and identified, when possible, to the genus level using taxonomic keys for Neotropical (^{Roldán, 1988}) and Costa Rican fauna (^{Springer et al., 2010}). As required by Costa Rican law, samples were deposited at the Zoological Museum of the University of Costa Rica. Taxa were assigned to functional groups: shredders, scrapers, collector-gatherers, collector-filterers, and predators-parasites (^{Merritt & Cummins, 1996}; ^{Ramírez & Pringle, 1998}; ^{Ramírez & Gutiérrez-Fonseca, 2014b}). SURBER data were used for quantitative indexes and density was shown as individuals per square meter. Shannon’s diversity (H’), richness, and percentages of indicator taxa (Ephemeroptera and Trichoptera, summarizing the EPT (Ephemeroptera-Plecoptera-Tricopthera parameter), since there were no Plecoptera) were calculated based on the total of sampled individuals (SURBER + multi-habitat samples). Along with invertebrate sampling, water quality (pH, temperature, dissolved oxygen, biological oxygen demand-BOD, and fecal coliforms), and quantity (flow discharge) were also measured. Analyses were done in the National Water Laboratory (AyA Costa Rica) following the “Standard Methods for Examination of Water & Wastewater” by the ^{APHA, AWWA, and WPCF (1975)}.

Statistical analysis. Non-parametric tests were chosen to deal with a wide range of variables, scales, and units, since these could affect the normality of the selected variables. Two MDS analyses were conducted in order to characterize the responses of macroinvertebrate communities to spatial and temporal gradients. For this purpose, a global matrix was built integrating the above-mentioned climatic, land use, physicochemical, biotic, and geomorphological variables for the sampling stations over the studied period. The spatial matrix included values for each sampling station of the available non-redundant variables, while a temporal matrix with campaign-data was used for the temporal analysis. Statistical analyses were carried out using the commercial software package Statgraphics (Statistical Graphics Corp.).

Results

Physicochemical water properties. Physicochemical water characteristics changed along basin gradients as shown in Table 2. Statios B1’ was affected by dairy-farm waste waters, while P3 showed the aggregated impacts of agricultural diffuse pollution and untreated sewage discharge from the town of Pacayas. Water temperature ranged from 10 to 21 ºC along the spatial gradient. Organically polluted sites showed peaks of fecal coliforms, phosphates, BOD, and COD, especially B1’, located in the headwaters of the Birrís River among dairy farms. Both streams had better values in downstream stations B2 and P4, via self-purification and dilution processes.

Composition and structure of the macroinvertebrate assemblages. A total of 368 samples were taken and sorted, with 66,182 individuals identified up to a family or genus level (Table 3). These individuals were classified into 16 orders, 50 families, and 77 genera, similar to other studies in Central American streams (^{Lorion & Kennedy, 2009}; ^{Kohlmann et al., 2015}). Global composition was dominated by Diptera (72% of all individuals), mainly Chironomidae and Simuliidae, and, secondly, by Ephemeroptera (mainly Baetidae) (17%) (Fig. 2). Glossosomatidae (Trichoptera) and Oligochaeta accounted respectively for 3% and 4%, while the rest of the taxa made up 4% of the population, in decreasing order of weight: Amphipoda, Gastropoda, Coleoptera, Hemiptera, Collembola, Hydrachnidia, Lepidoptera, Odonata, Hirudinea, Megaloptera, Mollusca, and Tricladida. Some individuals belonging to the Diptera and Coleoptera orders could not be identified for lack of detailed regional taxonomy. The unidentified individuals were classified morphologically, accounting for 23 different Coleoptera morphotaxa, 3 Diptera Psychodidae morphotaxa, 8 Diptera Tipulidae morphotaxa, and 41 Diptera morphotaxa belonging to families that could not be identified. Additionally, two Diptera Dixidae taxa were collected, one identified as Dixella and the other one classified as undetermined (cf. Meringodixa), in the upper station of the Birrís River (B1), at 2940 m.a.s.l. very close to the Irazú volcano crater.

Figure 2 Composition of the Birrís River (Costa Rica) macroinvertebrate assemblages at family level. undet. = undetermined. 

Three families accounted for 84% of the population. The highest abundance was recorded for the Chironomidae (34%), Simuliidae (33%), and Baetidae (17%) families (Fig. 2). At a genus level, Simulium represented 33% of the population, whereas Camelobaetidius, Baetodes, and Fallceon (Ephemeroptera, Baetidae) accounted for 6%, 6%, and 4% of the population, respectively. Glossosomatidae (Trichoptera) represented 3%. The rest of the population included an undetermined Oligochaeta (4%) and less than 2% belonging to other taxa. None of the analyzed samples had individuals belonging to the Plecoptera order.

Spatial gradients and responses. The composition of macroinvertebrate assemblages varied along the spatial gradient (Fig. 3). Diptera dominated throughout the entire basin, especially in organically polluted reaches, such as the P3 (sewage) and B1’ (cattle) stations. Trichoptera were more abundant in clean headwaters such as P1 and B1, while Ephemeroptera dominated at stations P1, P2, and P4. The functional structure varied also at the spatial scale. Shredders were more abundant in forested clean headwaters, being almost absent at the remaining sites. Collector gatherers were dominant at impacted sites, whereas collector filterers increased in re-oxygenated reaches after self-purification processes.

Figures 3a-h Average values and standard deviation for relevant variables of composition and diversity of macroinvertebrate assemblages along the spatial gradient in the Birrís River (Costa Rica). a) Frequency of Ephemeroptera + Plecoptera + Trichoptera; b) Frequency of Chironomidae; c) Density; d) Ratio EPT/(Diptera + Oligochaeta); e) Shannon index (Familes); f) Richness (Families); g) EPT Richness (Families); h) Richness of Diptera (Families). 

Organically polluted reaches showed the minimum values for EPT%. However, EPT values tended to increase as self-depuration processes operated downstream. Chironomidae frequency varied oppositely to EPT patterns of variation (Fig. 3). Differences among stations (Kruskal-Wallis median test) showed that they were significant for EPT average frequencies (p < 0.0001), and Chironomidae frequencies (p <0.0001).

Analyzed data show that density tended to increase downstream in both channels. Stations affected by organic pollution (B1’ and P3) had higher densities than less-affected sites nearby (B1, P2). Differences in density were statistically significant among sites (Kruskal-Wallis Test, p =0.04). Integrity of aquatic communities was evaluated using the ratio between EPT and Diptera+Oligochaeta frequencies, showing a minimum at the polluted sites, and recovering downstream from these sites, possibly by self-depuration (p <0.0001) (Fig. 3). Spatial trends in diversity and richness were also analyzed. The analysis showed that organically polluted reaches had lower diversity and richness values than the rest.

Diversity (p <0.0001) and richness (p =0.04) were significantly different between sites in accordance with physicochemical values (Fig. 3). EPT taxa richness decreased in organically polluted sites (B1’ and P3), where Diptera richness was medium. Reaches such as P2 (relatively clean waters) and B4 (self-depurated waters) showed higher values for EPT diversity. Diptera showed higher diversity in clean Birrís River headwaters (B1) at the highest elevation for the basin. Differences amongst sites were significant according to Kruskal-Wallis test for EPT richness (p <0.0001) but not for Diptera richness (p =0.31) (Fig. 3).

The spatial MDS identified four main dimensions. The first one was linked to organic pollution contrasting less impaired sites (structure and diversity in positive coordinates) and organically polluted ones (negative coordinates). The latter represents the self-depuration gradient from impacted sites (positive coordinates) to regenerated reaches (negative coordinates) through re-oxygenation and motion (waterfalls, high gradient reaches, and others). Figure 4 shows the plot of the sampling stations in the two main dimensions identified in the spatial MDS analysis.

Figure 4 Plot of the location of the nine sampling stations in the first two dimensions of the MDS spatial variability analysis of Birrís River (circles) and Pacayas Creek (diamonds) (Costa Rica). 

An exploratory PCA analysis was carried out, taking into account the recommendations by ^{Box and Andersen (1955)} regarding non-linearity. This analysis estimated that the first and second axes accounted for 59% of total variance (71% if a third dimension was included).

Temporal gradients and responses. The highest EPT abundance was detected during the dry season, with minimum frequencies of Chironomidae (Fig. 5). Data show that density increased in the dry season (February-April), with statistical differences recorded between months (Kruskal-Wallis Test, p <0.0001).

Figures 5a-h Average values and standard deviation (February, 2002-October, 2003) for relevant variables of composition and diversity of macroinvertebrate assemblages along the temporal gradient in the in the Birrís River (Costa Rica). a) Frequency of Ephemeroptera+ Plecoptera+ Trichoptera; b) Frequency of Chironomidae; c) Density; d) Ratio EPT/(Diptera+Oligochaeta); e) Shannon index (Familes); f) Richness (Families); g) EPT Richness (Families); h) Richness of Diptera (Families). 

The ratio EPT/ Diptera+Oligochaeta varied following the same pattern, but differences between months were not significant at 95% (Kruskal-Wallis Test, p = 0.07). Diversity indexes at the family level reached their maximum values by the end of the rainy season (December).

The pattern of variation for diversity and richness indexes was significant at 95% (p <0.0001). Maximum EPT and Diptera richness values were achieved in the transition from wet to dry season. EPT variation was not significant between months (p =0.22) whereas Diptera variation was (p <0.0001).

The temporal MDS analysis identified another four main dimensions. The first was related to the gradient between dry and wet seasons. The second dimension was related to the difference between the two sampled years, 2002 and 2003. The former year had higher rainfall values than the latter, which showed an abnormal dry period between July and August. Figure 6 shows the plot of the sampling stations in the two main dimensions identified in the temporal MDS analysis.

Figure 6 Plot of the location of the eleven field campaigns (2002-triangles; 2003-crosses) in the first two dimensions of the MDS temporal variability analysis of the Birrís River Basin (Costa Rica). 

Discussion

Land use and macroinvertebrate assemblages. The Birrís River Basin is highly affected by its headwaters, as are many other volcanic basins in the tropics. At the top of the basin, organic waste from dairy farms affects small headwater streams and this impact is amplified because of its natural low baseflow discharge, which results in high vulnerability. Pressure from human and cattle-raising activity increases downstream, affecting medium and low reaches. Lower reaches recover quicker than headwaters from organic effluents because of their higher flows and gradient, which allow water self-purification at numerous steps, i.e., rapids, riffles, and waterfalls as reported by ^{Ramírez and Pringle (2001)}. More tolerant taxa dominate in impaired reaches (^{Prygiel et al., 1999}; ^{Del Rosario et al., 2002}), whereas more sensitive ones occur in headwaters and self-purified reaches. We note that Plecoptera are absent, even in the less affected sites in headwaters, despite the fact that this order is present in other reference sites in the Reventazón Basin (^{Astorga et al., 1997}).

Assemblages in organically polluted sites were dominated by collectors such as Chironomidae and Oligochaeta (tolerant taxa, as found by ^{Prygiel et al., 1999}; ^{Del Rosario et al., 2002}; ^{Figueroa et al., 2003}; ^{Córdova et al., 2009}; ^{Ruiz-Picos et al. 2016}), with less non-tolerant taxa such as Ephemeroptera and Trichoptera (^{Resh, 1993}; ^{Chang et al., 2014}). Dominance of tolerant taxa led to low diversity and EPT richness in these sites. However, the slope of river channels and the abundance of steps and cascades allow quick self-purification of water downstream of organically polluted sites. Consequently, Ephemeroptera and Trichoptera abundance recovered shortly after impaired sites, as reported by ^{Astorga et al. (1997)} for the Reventazón River Basin, in which the Birrís Basin is included. Macroinvertebrates taxa display adaptations for high resilience and resistance, as do other communities in the tropics (^{Longo et al., 2010}).

Spatial and temporal responses. The spatial variability of the assemblages is determined by natural gradients, mainly elevation, from 2900 to 1200 m.a.s.l., as well as human impacts. The main dimensions of spatial variability are pollution and self-purification. Most intolerant taxa belonging to the Ephemeroptera and Trichoptera orders decreased in density, frequency, and richness at the polluted sites. Tolerant taxa, such as Chironomidae and Oligochaeta, followed opposite trends, prevailing in the polluted sites. Oligochaeta tended to dominate at higher elevations and low-temperature polluted sites with fine substrates (sand and silt derived from volcanic ash), while Diptera taxa dominated in medium and lowland polluted sites. However, a single maximum of Diptera richness was found at the top of the volcano in clean headwaters, at 2940 m.a.s.l., where many undetermined Diptera taxa (i.e. one unidentified genus of Dixidae) were collected. These results coincide with the conclusions of ^{Chaverri-Sánchez and Borkent (2007)}, suggesting that further collecting, especially at high elevations, will reveal further Dixidae species, taking into account that 50% of the named Dixella species in Costa Rica were collected at over 1200 m. Height and volcanic activity at the top of the Irazú cone lead to special ecological conditions in headwaters and stream communities characterized by high sulphate/metal pulses, cool temperatures, slow growth cycles, and acidic pH to which singular macroinvertebrate taxa would be adapted.

Temporal variability is ruled by inter and intra annual climatic variability, which were identified as the main dimensions of temporal variability. Annual diversity maximums were reached in the transition between the rainy and dry seasons (December), as seen in other Neotropical assemblages (^{Sánchez-Argüello et al., 2010}).

Rainfall, high flows and macroinvertebrates assemblages. Macroinvertebrate density was reduced during rainy periods, as reported by ^{Rincon and Cressa (2000)}, increasing rapidly at the beginning of the dry season and again in the short dry period between June and July. ^{Ramírez and Pringle (1998)} observed biomass peaks for Ephemeroptera and Diptera in riffle habitats in the dry period (November-April). High flows affected most Ephemeroptera and Trichoptera taxa, reducing their richness, while Diptera richness increased at the end of the rainy season. ^{Ramírez and Pringle (1998)} found that high flows reduced abundance affecting mainly collector-gatherers, predators, and shredders. Filterers and scrapers, as Simuliidae (Diptera) are more adapted to high flows, and thus, were not that affected in this study. In other seasonal tropical streams with less stream power, Baetidae dominate during the rainy season (^{Ribeiro & Uieda, 2005}), thanks to their morphological adaptations. Additionally, ^{Boyero and Bosch (2004)} found that the recolonization process in a Costa Rican stream on cobble substrates usually takes less than 20 days and the most abundant families on recolonizing stones are Simuliidae and Baetidae, being also the most abundant families in drift samples (^{Boyero & Bosch, 2002}).

Thus, we conclude that aquatic macroinvertebrate assemblages in the Birrís River Basin are adapted to shifting spatial and temporal conditions that include unpredictable severe floods, steep channels, active geomorphology, moderate to high organic pollution, self-purification processes, and hydrological seasonality. Diversity generally decreased when impacts were more severe, while tolerant taxa seemed to be favored by increased impacts.

The composition of aquatic macroinvertebrate assemblages in the Birrís River Basin is dominated by tolerant taxa, especially Chironomidae, Simuliidae, and Baetidae, adapted to moderate organic pollution and high stream power. However, composition varied spatially and temporally, with a resilient response linked to self-purification processes. Further taxonomic and ecological research should be undertaken to study assemblages in higher volcanic reaches, dominated mainly by Diptera and Coleoptera undetermined taxa, in order to elucidate the taxonomy and functioning of these assemblages.

The structure of aquatic communities varied spatially, according to environmental pressures, and temporally, responding to seasonality and inter-annual variations. Maximum diversity was reported at the end of the rainy season (beginning of the dry season), when most taxa showed mature larval development.

Special attention should be paid to the high vulnerability of headwater courses in volcanic basins, with low base flows that amplify the impacts of cattle and agriculture. Special management schemes should be applied to protect and restore these headwaters, which are strategic sources of freshwater for surrounding populations. Assemblages in medium and lower reaches showed greater resilience, thanks to the higher flows and gradients that triggered intense self-purification processes. However, the impact of human settlements is still important in these reaches, because of the lack of adequate sewage treatment facilities.

Further research should be focused on the taxonomy and ecology of aquatic communities in highlands volcanic streams. Maintenance of adequate ecological status in these streams is very important to ensure the conservation of strategic freshwater resources in highly populated tropical regions.