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

Abstract

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.

1. 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 degradation is due to anthropogenic activities such as the excessive use of pesticides and fertilizers, alteration of land use and untreated wastewater discharges from industries and houses [5,6].

Untreated wastewater discharges are the main source of surface water pollution in urban areas [7,8,9,10]. They reduce the quality of the water body and stimulate the proliferation of pathogenic organisms that can cause severe diseases concerning public health [11,12,13]. Its collection, transport, treatment, disposal and reuse are taken into account in the framework of sanitation according to the United Nations and the World Health Organisation. In addition, the provision of adequate treatment for the disposal of urine and human faeces is considered a human right [14,15]. Despite this, there remain countries, in particular developing ones, where access to sanitation is limited [16,17,18].

Central American countries face substantial challenges regarding water and sanitation. In Costa Rica, for example, 76% of houses use septic tanks and only 22% are connected to sewerage; but only 8.2% of the wastewater is treated before being discharged into the rivers [19]. The effect of this lack of wastewater treatment in the riverine ecosystem could be even worse in areas with accelerated socio-economic development [2]. One of these areas is the Virilla river catchment, which drains the Greater Metropolitan Area (GAM) and receives approximately 67% of the wastewater discharges of the country [15]. The GAM corresponds to only 4% of the Costa Rican territory but its population density is nearly 1200 hab km${}^{-2}$. It also contains the largest sector of industrial and commercial activities in the country [15]. To reduce the potential impact of these activities in the environment, implementation of new legislation and environmental programs have been undertaken during the last years [19,20]. However, the success of these policies is difficult to estimate, due to the functions overlapping across different institutions and the lack of accurate information about the ecological status of the water bodies [17,19,20]. Notably, information about the surface water bodies is essential regarding control of pollution and the application of mitigation strategies. These networks are useful in providing reliable information to better understand temporal and spatial changes in water quality and supporting integrated water management processes [17,21]. But only a few studies have been reported in the country and water monitoring networks are scarce [22,23,24,25,26].

Water quality in rivers can be estimated using physical and chemical characteristics or macro and micro biological indicators [22,23,27]. The category and number of such indicators are not generally uniform among similar studies and its selection depends on the objectives and financial resources of the monitoring program. This makes difficult the comparison among different studies and different disciplines are rarely integrated together. One approach oriented to simplify this is the water quality index; which takes into account the value of certain parameters in an overall water quality score, through the estimation of a relative weight of each parameter [28,29]. However, there is still uncertainty related to the selection and weight of the parameters. Thus, there is an interest in more holistic approaches, which will provide better and integrative perspectives of the riverine ecosystem status [30].

This study presents results of an holistic integrative approach to surface water quality assessment in the Virilla river catchment in Costa Rica, where sanitation mitigation measures will be developed. Nevertheless, there is not enough background information about the water quality of this area in general. The ecosystemic approach presented here includes physical, chemical, microbiological and benthic macroinvertebrate data; which are interpreted using land use information, multivariate statistics and water quality indexes. This information is useful to better understand the water quality status in the Virilla river catchment; and therefore, to generate evidence-based mitigation strategies in the near future. Finally, this information will provide a surface water quality background in Costa Rica for further implementation of sanitation processes in the country.

2. Materials and Methods

2.1. Study Area

This study was developed in the Virilla river catchment in Costa Rica, between the longitude 84${}^{\circ }$10${}^{\prime }$48${}^{″}$ and 84${}^{\circ }$02${}^{\prime }$38${}^{″}$, and the latitude 9${}^{\circ }$59${}^{\prime }$30${}^{″}$ and 9${}^{\circ }$57${}^{\prime }$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 ${}^{\circ }$C to 24.8 ${}^{\circ }$C with an average of 21.4 ${}^{\circ }$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 ^®. 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).

2.2. 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 ${}^{\circ }$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 $\mathsf{\mu }$m pore filter (Advantec^® 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 ${}^{\circ }$C, 180 ${}^{\circ }$C and 105 ${}^{\circ }$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^® medium and 24–48 h of incubation at 37 ${}^{\circ }$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 $\mathsf{\mu }$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].

2.3. 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.

3. Results and Discussion

3.1. Physical, Chemical and Microbiological Data

Table 1. Summary of the physical, chemical and microbiological parameters in the Virilla river.

Parameter	Units	Average	Minimum	Maximum	SD
Temperature	${}^{\circ }$C	19.8	15.2	28.2	2.85
pH	-	7.69	6.69	8.75	0.35
Conductivity	$\mathsf{\mu }$S/cm	183	37	408	94.9
DO	%	81.8	21.6	113.1	1.35
Turbidity	NTU	10.59	1.41	66.53	13.8
TS	mg/L	170	48	543	89.4
TDS	mg/L	137	11	312	64.5
TSS *	mg/L	30	<1	382	223
Alkalinity	mg/L CaCO${}_3$	81.5	12.2	264	57.2
BOD *	mg/L O${}_2$	10.15	<2	206.8	8.06
Ammonium *	mg/L N-NH${}_4$${}^{+}$	1.050	<0.07	9.300	2.88
Nitrite *	mg/L N-NO${}_2$${}^{-}$	0.1243	<0.012	0.9260	0.59
Nitrate	mg/L N-NO${}_3$${}^{-}$	1.954	0.047	8.546	1.78
TP *	mg/L	0.435	<0.02	3.520	0.71
SRP *	mg/L	0.350	<0.02	3.706	0.69
Fluoride *	mg/L	0.106	<0.043	0.270	0.07
Chloride	mg/L	7.75	0.95	25.45	5.79
Sulphate	mg/L	9.55	1.02	22.43	5.39
Total coliform	MPN/100 mL	7.80 × 10${}^6$	2	1.40 × 10${}^8$	2.21 × 10${}^7$
Escherichia coli	MPN/100 mL	4.42 × 10${}^6$	2	1.07 × 10${}^8$	1.41 × 10${}^7$

DO: dissolved oxygen; TS: total solids; TDS: total dissolved solids; TSS: total suspended solids; BOD: biochemical oxygen demand; TP: total phosphate; SRP: soluble reactive phosphorus; NTU: nephelometric turbidity units; MPN: most probable number. *: calculated using the Maximum Likelihood Estimation method.

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 ${}^{\circ }$C, ranging from 15.2 ${}^{\circ }$C to 28.2 ${}^{\circ }$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. 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 Figure 2. Of the 406 possible correlations, 273 were significant (p < 0.05). TDS and conductivity ($\tau$ = 0.776), and sulphate and chloride ($\tau$ = 0.721) presented high positive significant correlations. Moderate correlations (0.5–0.7) were found in 37 cases, where 33 were classified as positive correlations and 4 as negative. Low correlation distribution (0.3–0.5) included 84 positive and 38 negative, for a total of 122. Finally, negligible correlations (<0.3) were presented between 122 variables, which 55 were positive and 57 negative. In general, most of the physical and chemical parameters were positively correlated each other. However, PSO presented a negative correlation with most of the parameters as expected, whereas the positive correlations previously mentioned indicate that the level of pollution increased.

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.

3.2. Benthic Macroinvertebrates

A total of 13,063 benthic macroinvertebrates were collected and classified into 47 families which correspond to 18 orders (

Table A1. Cumulative of families per sample site and occurrence of the macroinvertebrate data.

Order	Family	Sample Site							Occurrence (%)
		T1	T2	V1	V2	V3	V4	V6
Annelida	Oligochaeta	3	7	0	13	4	53	10	30.6
Arhynchobdellida	Erpobdellidae	8	8	8	3	36	27	37	49.0
Arhynchobdellida	Salifidae	3	5	1	6	3	0	6	22.4
Basommatophora	Lymnaeidae	0	0	0	0	0	0	1	2.0
Basommatophora	Planorbidae	0	0	0	1	0	0	1	4.1
Basommatophora	Physidae	9	234	6	262	328	83	94	65.3
Blattaria	Blaberidae	1	0	0	0	0	0	0	2.0
Coleoptera	Curculionidae	0	0	0	1	0	0	1	4.1
Coleoptera	Dytiscidae	1	0	0	0	0	0	1	4.1
Coleoptera	Dryopidae	2	0	4	1	0	0	0	6.1
Coleoptera	Elmidae	1	2	0	0	0	0	0	4.1
Coleoptera	Hydraenidae	1	0	0	0	0	0	0	2.0
Coleoptera	Hydrophilidae	0	0	0	0	0	1	1	4.1
Coleoptera	Lampyridae	1	0	0	0	0	0	0	2.0
Coleoptera	Staphylinidae	6	4	9	5	1	0	4	32.7
Coleoptera	Psephenidae	0	0	1	11	0	0	0	4.1
Diptera	Ceratopogonidae	4	0	0	0	0	0	0	6.1
Diptera	Chironomidae	343	588	324	688	493	505	430	85.7
Diptera	Empididae	1	1	0	2	0	0	0	6.1
Diptera	Muscidae	5	11	1	0	1	0	16	16.3
Diptera	Psychodidae	11	1	3	9	1	2	29	26.5
Diptera	Simuliidae	991	595	383	104	28	8	26	69.4
Diptera	Stratiomyidae	1	0	0	0	0	0	0	30.6
Diptera	Tipulidae	5	0	8	1	1	0	0	8.2
Entomobryomorpha	Entomobryidae	0	0	1	0	0	0	1	4.1
Ephemeroptera	Baetidae	370	378	409	319	133	23	21	91.8
Ephemeroptera	Leptohyphidae	599	926	700	334	23	8	12	83.7
Ephemeroptera	Leptophlebiidae	176	8	57	11	2	1	0	57.1
Hemiptera	Belostomatidae	0	5	0	0	0	0	0	4.1
Hemiptera	Cicadellidae	0	3	0	0	0	0	0	4.1
Hemiptera	Guerridae	1	0	0	0	0	0	0	2.0
Hemiptera	Veliidae	1	14	0	2	0	0	3	10.2
Isopoda	Asellota	0	0	1	0	0	0	0	2.0
Lepidoptera	Pyralidae	4	8	2	1	0	0	0	18.4
Littorinimorpha	Hydrobiidae	0	0	0	0	2	0	1	4.1
Odonata	Calopterygidae	12	10	9	3	2	1	0	32.7
Odonata	Coenagrionidae	3	3	4	0	0	0	3	22.4
Odonata	Libellulidae	4	0	3	2	0	0	0	16.3
Rhynchobdellida	Glossiphoniidae	0	66	5	29	45	10	52	26.5
Trichoptera	Glossosomatidae	37	53	60	16	0	0	0	36.7
Trichoptera	Helicopsychidae	0	0	1	1	0	0	0	4.1
Trichoptera	Hydrobiosidae	22	18	3	6	0	0	0	28.6
Trichoptera	Hydropsychidae	156	293	166	109	24	6	275	81.6
Trichoptera	Hydroptilidae	2	9	2	4	0	0	1	4.1
Tricladida	Planariidae	2	0	2	0	0	0	0	4.1
Trombidiformes	Hydrachnidia	0	1	4	1	0	0	0	4.1
Veneroida	Sphaeriidae	0	19	0	1	1	1	1	12.2
Total of invertebrates identified		2786	3270	2177	1946	1128	729	1027
Total of families identified		33	28	29	30	18	14	25

). 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 ($\tau$ = −0.523), Leptohyphidae and nitrite ($\tau$ = −0.521), Baetidae and water temperature ($\tau$ = −0.507), and Baetidae and TP ($\tau$ = −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 ($\tau$ = −0.450), diversity and nitrite ($\tau$ = −0.363), diversity and TP ($\tau$ = −0.368), richness and turbidity ($\tau$ = −0.362), richness and TSS ($\tau$ = −0.302), richness and TP ($\tau$ = −0.335), richness and SRP ($\tau$ = −0.302), abundance and turbidity ($\tau$ = −0.508), abundance and TSS ($\tau$ = −0.385), and abundance and TP ($\tau$ = −0.330).

According to the RDA (

Table 2. Redundancy analysis (RDA) results showing the variance of the macroinvertebrate data explained by canonical axes and explanatory variables.

Axes and Variables	Explained Variation (%)	Pseudo-F	p-Value
Canonical axes
First axis	32.3	26.6	0.001 *
Second axis	10.7	8.8	0.001 *
All axes	49.1	4.2	0.001 *
Explanatory variables
Temperature	−70.9	18.5	0.005 **
Nitrite	−64.9	17.2	0.005 **
Ammonium	−64.3	13.0	0.005 **
SRP	−63.3	9.9	0.005 **
TS	−62.6	9.8	0.005 **
Alkalinity	−59.8	7.8	0.005 **
Nitrate	−61.1	5.2	0.005 **
TSS	−57.8	3.4	0.005 **

* p < 0.01; ** p < 0.05.

) 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 ${}^{\circ }$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.

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].

3.3. Water Quality Indexes

Results of the water quality indexes are presented in

Table 3. Average of water quality indexes per sample site in the Virilla river (standard error in parenthesis).

Site	NSF-WQI		Dutch Index		BMWP-CR
	Average	Classification *	Average	Classification **	Average	Classification *
T1	73 (0.8)	Good	5 (0.3)	Incipient	71 (7.2)	Regular
T2	73 (1.1)	Good	5 (0.2)	Incipient	54 (6.5)	Bad
V1	74 (1.3)	Good	5 (0.3)	Incipient	55 (7.2)	Bad
V2	71 (0.8)	Good	6 (0.4)	Incipient	50 (10)	Bad
V3	69 (1.2)	Medium	7 (0.5)	Moderate	28 (3.4)	Bad, highly polluted
V4	64 (2.7)	Medium	9 (1.0)	Moderate	17 (3.5)	Bad, highly polluted
V5	63 (2.2)	Medium	10 (0.8)	Severe	-	-
V6	65 (3.4)	Medium	9 (0.5)	Moderate	28 (4.1)	Bad, highly polluted

NFS-WQI: National Sanitation Foundation water quality index; BMWP-CR: Biological Monitoring Working Party modified to Costa Rica. * in terms of water quality; ** in terms of pollution.

. 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.

4. 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.