Ecosystemic Assessment of Surface Water Quality in the Virilla River: Towards Sanitation Processes in Costa Rica

Water quality information is essential supporting decision making in water management processes. The lack of information restricts, at some point, the implementation of adequate sanitation, which is still scarce in developing countries. In this study, an ecosystemic water quality assessment was conducted in the Virilla river in Costa Rica, in a section of particular interest for future sanitation development. It included the monitoring of physical, chemical, microbiological and benthic macroinvertebrate parameters from 2014 to 2016. Mutivariate statistics and water quality indexes were used for data interpretation. Results indicated that water quality decreased downstream towards more urbanised areas. Particularly, extreme values of phosphorous, nitrogen and E. coli were found. Sample sites were grouped in two clusters, which were consistent with land use. Benthic macroinverterbrates diversity was predominantly represented by Baetidae, Chironomidae, Leptohyphidae, Hydropsychidae, Simuliidae and Physidae. They were mostly influenced by water temperature, nitrite, ammonium, soluble reactive phosphorous, total solids, alkalinity, nitrate and total suspended solids. Three water quality indexes consistently showed the poor condition of the water body. The overall results indicate that the main sources of pollution in the river are likely to be wastewater discharges. Thus, special efforts should be undertaken regarding its regulation in the country.


Introduction
Rivers are an important source of aquatic biota and water for human development. They provide essential resources for recreation, tourism, human water consumption, agriculture, electricity generation and industry [1,2]. However, riverine systems have been constantly modified by human activities [3], causing the alteration of the hydrological cycle and the degradation of their water quality at both a local or regional scale [4]. Whereas natural processes such as erosion, soil mineralisation and meteorological conditions, eventually impact surface water quality; the major impact on its

Study Area
This study was developed in the Virilla river catchment in Costa Rica, between the longitude 84 • 10 48 and 84 • 02 38 , and the latitude 9 • 59 30 and 9 • 57 28 (see Figure 1). The Virilla river flows from the north-east Central Volcanic range of the Central Valley in Costa Rica, to the south-west Pacific region; where it confluences with the Grande river and forms the Grande de Tárcoles river. The covered area includes a section of approximately 20 km of the river from the city of San Miguel to San Antonio, both of the Heredia providence. This section is remarkably important because it will be directly impacted by the implementation of a sewage treatment plant in the upcoming years. In this area, the climate is characterised by a dry season (December to March) and rainy season (May to October) regarding April and November as transition months. During the study period, the average monthly precipitation was 148 mm and air temperature ranged from 17.7 • C to 24.8 • C with an average of 21.4 • C. The elevation gradient was 277 m ranging from 1161 m to 884 m. This area is classified as premontane wet forest zone with irregular relief and soils predominantly vertisols. Land use categorisation (56 km 2 ) included forest (4.9%), pasture (14.9%), arable (21.0%), industry (11.5%) and urban (47.7%), and it was generated by the photo-interpretation technique in ArcGIS 10.4.1 (ESRI) at scale 1:5000, using the satellite images Quick Bird II at 0.6 m resolution distributed by Digital Globe in Google Earth Pro R . In this section of the river, there are wastewater discharges from real estate activities (e.g., sewage disposal from condos), from manufacturing industries of concrete, cement and food products (e.g., coffee, beverages and food preservatives).

Sampling and Methods
Water samples were collected at eight sample sites, every two months, from October 2014 to March 2016 (n = 71). They were collected using high-density polyethylene and glass bottles previously washed with hydrochloric acid 3% m/v and de-ionised water. Samples were then stored at 4 • C and delivered to the laboratory within 6 h of collection. Microbiological samples were collected in 100 mL sterile vessels and stored separately. All samples were analysed using the procedures of the Standard Methods for the Examination of Water and Wastewater [31].
Temperature, dissolved oxygen (DO), pH and conductivity were determined in situ using a field meter Oakton 300 (Vernon Hills, IL, USA) and a Thermo Orion Star A222 (Chelmsford, MA, USA). An aliquot of each water sample was filtered through 0.45 µm pore filter (Advantec R GC-50) for analysis of total dissolved solids (TDS), total suspended solids (TSS), ammonium, nitrite and soluble reactive phosphorous (SRP). Alkalinity was measured by titration using a sulphuric acid standard solution. Total solids (TS), TDS and TSS were determined by gravimetry at 105 • C, 180 • C and 105 • C, respectively. Nephelometry was used to analyse turbidity (Oakton Inst. T100, Vernon Hills, IL, USA). Biochemical oxygen demand (BOD) was determined using the 5-days test and the modified Winkler method. Ammonium was measured using the indophenol blue reaction at 640 nm and nitrite was colorimetric determined at 543 nm. Total phosphorous (TP) and SRP were analysed spectrophotometrically by the stannous chloride method after the application of the persulfate digestion procedure for TP and after filtration for SRP. All spectrophotometric analyses were carried out in a Thermo Aquamate 2000E (Cambridge, UK). Fluoride, chloride, nitrate and sulphate were analysed by ion chromatography (Dionex ICS-5000, Sunnyvale, CA, USA). Finally, total coliform (TC) and Escherichia coli were determined by the multiple tube fermentation technique using the Fluorocult R medium and 24-48 h of incubation at 37 • C. The microbial population was estimated using the table reported by Woomer [32].
The sampling of benthic macroinvertebrate community was undertaken in the water body in safe and accessible sites, on the same day as water sample collection, from February 2015 to March 2016 (n = 49). Sample collection and preservation were achieved according to the methods described by Springer et al. [33]. In brief, composite samples were collected using a D network (500 µm mesh) placed opposite to the water direction flow; while the substrate was gently moved for approximately five minutes in order to collect the maximum amount of macroinvertebrates possible. This procedure was repeated three times in different points around each sample site (e.g., upstream and downstream) for a total of 15 min of effort per site. An initial selection was made in the field where macroinvertebrates were picked live and preserved in ethanol 70% m/v. In addition, samples were stored in plastic bags with ethanol 95% m/v and delivered to the laboratory for a more extensive selection. Macroinvertebrates were processed, classified and quantified as families using a stereoscope according to the taxonomic guidelines established in the Costa Rican surface water legislation [34].

Data Analysis
Data analysis, including non-parametric survival methods and ordination exploratory analysis, were carried out in R statistical package using NADA and vegan libraries [35][36][37]. Analysis of data with values below the quantification limit (<QL) included the calculation of descriptive statistics using the Maximum Likelihood Estimation (MLE) method, and their transformation to tied ranks before performing multivariate non-parametric tests [38]. In order to estimate the degree of association between water quality parameters the non-parametric, Kendall's Tau correlation coefficient was used [39]. An analysis of similarities (ANOSIM) was performed computing 999 permutations to evaluate the difference among sample sites and among sample campaigns. The spatial grouping of the sample sites was defined using hierarchical cluster analysis with Ward's method of association [40] and squared Euclidian distance as a measure of similarity. These last two tests were performed using only the physical, chemical and microbiological data.
Richness, abundance and the Shannon index [41] were calculated as measures of macroinvertebrate's diversity. Redundancy analysis (RDA) was applied to elucidate the relationships between environmental and community data. Based on the percentage of occurrence, the families that were not present more than 24 times were removed for the subsequent analysis. Previous RDA, physical and chemical variables were standardised [42] and biological data were transformed using a Hellinger transformation [43]. Multiple linear regression, variation inflation factors (VIF) and BIOENV function (Best Subset of Environmental Variables with Maximum Correlation with Community Dissimilarities) were employed to select the best subset of physical and chemical variables that influenced the macroinvertebrate data [44][45][46]. A Monte Carlo permutation test (999 permutations) then allowed the statistical validation of the RDA.
Finally, three water quality indexes (WQI) were calculated to evaluate water quality overall. These were the widely used index of the National Sanitation Foundation of the United States (NSF) [47], and the two indexes established in the Costa Rican national legislation; which includes the Dutch index and the Biological Monitoring Working Party modified to Costa Rica (BMWP-CR) [34]. E. coli concentration was used instead of faecal coliform for the calculation of the NSF-WQI. Table 1 presents the summary of the descriptive analysis of the physical, chemical and microbiological parameters in the study area. The spatial trend for some of these parameters is shown in Figure A1. Average water temperature was 19.8 • C, ranging from 15.2 • C to 28.2 • C. Higher temperature values were observed downstream, where riparian vegetation is less abundance. Non-extreme pH values were observed, ranging from 6.69 to 8.75. However, alkalinity concentration increased downstream with values up to 264 mg/L CaCO 3 . These levels were previously associated to vertisol soils presented in the catchment [48], which can also increase conductivity. Average DO concentration was 81.8%, but some critical values were obtained. For instance, the minimum value was 21.6% in sample site V6, but site V5 presented the lowest average DO concentration (69.7%). Turbidity, TS, TDS and TSS concentration also presented the same trend, an increasing towards a more industrialised and urbanised area. Its variation and increase can be related to wastewater discharges into the river as observed at the bottom of the study area. BOD average concentration was 10.15 mg/L O 2 . This is consistent with the results reported by Herrera et al. [48]. However, BOD ranged from <2 mg/L O 2 in sample site V1 to a maximum of 206.8 mg/L O 2 in site V4. Particularly, in the riparian area of this last sample site, landfill was observed. This comes from an informal urban settlement next to the river. In addition, illegal sewage discharges may be present. The concentration of different N-compounds changed with elevation (ammonium, nitrite and nitrate). The first four sample sites presented constant concentration of these compounds, prevailing nitrate. This may be related to the arable land use in this area of the catchment. Since site V3, ammonium concentration increased drastically, being as high as nitrate concentration in some cases.

Physical, Chemical and Microbiological Data
Ammonium, nitrite and nitrate values were found in concentrations that may be toxic for aquatic species [49,50]. Sources of such higher concentration of nitrogen and phosphorus can be either point or diffuse. Some of these are fertilisers, run-off from agricultural fields, industrial effluents and untreated wastewater discharges [51]. Increasing of solids, BDO, ammonium and SRP suggests that the last one may be the main cause of pollution. In contrast, fluoride, chloride and sulphate concentrations were within the standards even for drinking water purposes [52]. This indicates that there are no significant sources of pollution for these compounds.
Total coliform concentration ranged from 2 MPN/100 mL to 1.40 × 10 8 MPN/100 mL, with an average of 7.8 0× 10 6 MPN/100 mL. The lowest and highest average concentrations were found in site T1 (1.83 × 10 5 MPN/100 mL) and site V5 (2.35 × 10 7 MPN/100 mL), respectively. The average E. coli concentration was 4.42 × 10 6 MPN/100 mL and the same pattern was observed for E. coli in accordance to the sites with maximum and minimum concentrations. Coliforms concentration also increased downstream as reported by Leandro et al. [23]; additionally, the results were above the national legislation for wastewater discharges. The high values indicate that there is strong faecal pollution in the river which represents a threat to human health. This includes bacterial gastrointestinal and parasitic infections or harmful virus diseases [53].
Kendall's Tau correlation coefficients of the parameters evaluated in the study area are shown in Water composition among sample sites was significantly different (R = 0.3596, p = 0.001). In addition, cluster analysis generated two clusters at D link /D max < 40 ( Figure A2). The first group is formed by sites T1, T2, T1 and V2; which are located at the top of the study site. The second group included sites V3, V4, V5 and V6. These sites are distributed at the middle and bottom of the study area. Grouping of sample sites was consistent with land use. There is a change from arable-urban land use to a more pasture-arable-industrial use downstream. In addition, riparian land cover changes from forest to pasture-urban as elevation decreases. This change in land use is likely to be influencing water quality [54]. In the other hand, significant differences among sample campaigns were also observed (R = 0.3924, p = 0.001). This suggested that water quality may change seasonally. However, Mena-Rivera et al. [22] reported non-significant differences in water quality in another sub-catchment of the Virilla river. Long-term high-resolution monitoring networks are necessary to track these changes in seasonal patterns.
Despite the slightly improvement showed in sample site V6, the general trend observed with the parameters mentioned above is that water quality in the Virilla river decreased downstream. There is an increase of pollution at some point between sample sites V2 and V3, where the natural condition of the river has been extremely modified. Few meters before site V3, wastewater discharges are more evident, mainly from a coffee processing plant, an oxidation pond and houses. In addition to this condition, riparian landfill as cited previously was observed in site V4. Moreover, there is a reservoir in site V5. There, accumulation of material forms a surface layer, which also affect the minimum flow (ecological flow) observed in site V6. All these conditions are likely to be influencing water quality toward more critical levels.

Benthic Macroinvertebrates
A total of 13,063 benthic macroinvertebrates were collected and classified into 47 families which correspond to 18 orders (Table A1). Insecta was the group with the greatest dominance (n = 11,576). Richness, abundance and Shannon index per sample site are showed in Figure 3. Sample site T1 presented the highest richness with 32 families identified; whereas the lowest richness was observed in site V4 with only 14 families. Average abundance ranged from 467 macroinvertebrates in sample site V1 to 104 in sample site V4. Total Shannon index was H = 2.10 and average diversity ranged from 1.68 (site V1) to 0.90 (site V4). Families with major occurrence were Physidae (65.3%), Simuliidae (69.4%), Leptohyphidae (83.7%), Hydropsychidae (81.6%), Chironomidae (85.7%) and Baetidae (91.5%). In general, diversity and average abundance of macroinvertebrates per family decreased downstream, regardless the slight improvement in richness in sample V6 in comparison to the previous site. This trend is consistent with the physicochemical data previously mentioned.
In tropical rivers, benthic macroinvertebrates community distribution can be affected by the complex dynamic of physical and chemical parameters [55,56]. In this study, families with major occurrence were significantly correlated with some these parameters ( Figure 2). Negative moderate correlations were found between Leptohyphidae and TP (τ = −0.523), Leptohyphidae and nitrite (τ = −0.521), Baetidae and water temperature (τ = −0.507), and Baetidae and TP (τ = −0.506). Low negative significant correlations were presented between Simuliidae and water temperature, turbidity, TS, TSS, BOD, nitrite, TP, SRP and fluoride; Baetidae and turbidity, TS, TSS, BOD, ammonium, nitrite and SRP; Leptohyphidae and water temperature, conductivity, turbidity, TS, TSS, BOD, ammonium, SRP, fluoride, nitrate, sulphate, TC and EC. Nevertheless, Physidae presented low positive correlations with water temperature, conductivity and TDS. In addition, macroinvertebrate indexes also showed negative significant correlations. These were between diversity and ammonium (τ = −0.  According to the RDA (Table 2) approximately 49.1% of the variance of benthic macroinvertebrate data is explained by the physical and chemical variables. The parameters that most influenced macroinvertebrates community were water temperature, nitrite, ammonium, SRP, TS, alkalinity, nitrate and TSS. The RDA biplot (Figure 4) shows that most of the families with the highest occurrence (Simuliidae, Leptohyphidae, Baetidae and Hydropsychidae) presented negative correlations with the degree of pollution. Particularly, there was a high influence of N-compounds in aquatic insects diversity. Their capacity to tolerate high nitrogen concentration can be limited [57,58] and when it increased the insects abundance decreased, even for families that can tolerate some degree of pollution. In summary, diversity of benthic macroinvertebrates community decreased as pollution increased.  Families with less pollution tolerance were found in the upstream samples. For instance, Hydrobiosidae which is usually found under good water quality conditions. However, the families with major occurrence were those who are able to tolerate different levels of pollution. Chironomidae, the most second abundant family in all the samples sites, has the capacity to survive under polluted and anoxic environments. This is likely to be consistent with our results due to, according to the RDA, it was positive correlated with most of the chemical parameters; but negatively correlated with DO. This family can remove up to 70 g of organic matter per m 2 day −1 [59], because of its biomass increases with nutrients concentration [60]. This characteristic and the fact that Costa Rica is a "hotspot" for this family [61], suggest that further investigation have to be done as this group of insects could be relevant for rivers self-depuration in the country. Hydropsychidae and Simuliidae were also present with high occurrence. Both are able to filter fine organic matter and they are usually found in moderate and high discharge, being abundant locally [62]. Despite these characteristics, these families along with Baetidae and Leptohyphidae, were less found where pollution increased. They presented negative correlation with water temperature and TSS. The increase in water temperature between site T1 and V6 (∼5 • C) could be a affecting its distribution. In addition, it has been reported that abrasion and scouring by the increasing of suspended solids can affect macroinvertebrate community drift [63].  Physidae was also found in an important occurrence. This family showed a positive correlation with the temperature and BOD. Temperature seems to strongly influence its life cycle, when warmer environments foster a quicker growth [64]. Physidae also can feed on detritus [65], which could link them to areas with high BOD. This change in environmental conditions like the increase in the temperature and salinity may favour the prevalence of certain species of snails [64]. These pulmonate molluscs are also able to refresh the air from their cavity in the water surface [66]. This might be an advantage to stand in low oxygen levels and under certain pollution conditions.

Axes and Variables Explained Variation (%) Pseudo-F p-Value
No clear pattern was observed between the physicochemical variables and the Hydropsychidae family with the RDA. This could be associated with the taxonomic level employed, since its adequacy to track changes in different conditions depends on the species [67][68][69]. Therefore, it is necessary to study the different groups of macroinvertebrates at the genus or species level in order to better understand their role into water quality assessments. Finally, contrary to the trend of degradation downstream, sample site V6 showed some recovery signs. These were a small increase in diversity (H: 1.16), and a slight improvement in the quality index's score. This may be associated with partial changes in the river structure, which increase the heterogeneity of micro environments [70] along with the DO increasing [57,71].

Water Quality Indexes
Results of the water quality indexes are presented in Table 3. According to the NSF-WQI, the water quality of the Virilla river is classified into two categories, "medium" and "good". Their distribution was consistent with the grouping of the cluster analysis. These results are likely to be considered normal in rivers near the study area. For instance, Leandro et al. [23] and Mena-Rivera et al. [22] reported similar water quality classification using the NSF-WQI in different rivers and years in Costa Rica. In the other hand, the Dutch index showed "incipient" pollution in samples sites located at the top of the study area. Subsequently, from sample site V3 an increase in the index value was observed, resulting in higher levels of pollution. "Incipient", "moderate" and "severe" are typical categories for rivers in areas with low to medium population density in the GAM [25]. BMWP-CR water quality scores also decreased downstream in the Virilla river. The highest average value was 71 in sample site T1, representing a "regular" condition. This condition is presented when the score is between 67 and 95. The lowest score was obtained in sample V4 with an average of 17. This score allowed a water quality classification as "bad" or "highly polluted". BMWP-CR results were similar than NSF-WQI and the Dutch index, showing an almost identical trend in the water quality status. The NSF-WQI and the Dutch index pointed out that sample site V5 presented the worst condition in water quality. Nonetheless, the Dutch index better represented this condition, allowing a change into its classification to "severe". On this site, the collection of macroinvertebrate samples was not possible.
Despite this interpretation, the applicability of water quality indexes could be limited. This is because of the number and characteristics of the parameters incorporated (robustness), and the index's capacity to track changes under different spatial and temporal conditions (sensitivity). For instance, calculation of the NSF index is based on nine parameters including physical, chemical and microbiological indicators; while the Dutch index takes into account only three parameters (DO, BOD and ammonium). Nevertheless, the three water quality indexes employed in this study fairly represent the condition of the river, according to the trends previously mentioned. Nonetheless, the categorical classification used by each index is very different (both phrase and colour) and it may cause misinterpretation of the water quality by decision-makers. For example, Dutch index categorised site T2 as "incipient" pollution while the BMWP-CR as "bad" water quality. The disparate terms may influence the implementation of different mitigation efforts. Thus, the development of an integrative river health index for Costa Rica could lead into a more efficient and accurate monitoring tool for supporting management processes in the country.

Conclusions
Here, we reported useful information to support evidence-based decision making regarding water resources management in the Virilla river catchment in Costa Rica. High values of the physical, chemical and microbiological data were obtained. In addition, the benthic macroinvertebrate community was mainly affected by parameters related to point sources of pollution, such as wastewater discharges. Water quality decreased in the bottom of the study area, where industrialisation and urbanisation increases. Thus, better regulation of these sources of pollution, as well as the implementation of efficient sewerage and drainage systems should be addressed by local authorities.
The results of the water quality indexes established in the national legislation (Dutch index and BMWP-CR) were consistent; considering that there can not be a difference of two classes between the two indexes. However, the development of a more holistic water quality index have to be considered. Finally, the information presented here, as any water quality data in the country, can be useful for first implementations of numerical models, which would improve the efficiency estimation of current (e.g., septic tanks) and future wastewater treatment systems.