Spatial Characteristics and Temporal Evolution of Chemical and Biological Freshwater Status as Baseline Assessment on the Tropical Island San Cristóbal (Galapagos, Ecuador)

Abstract

The fragility of geographically isolated islands stresses the necessity of evaluating the current situation, identifying temporal trends and suggesting appropriate conservation measures. To support this, we assessed the freshwater quality of three stream basins on San Cristóbal (Galapagos) for two consecutive years. Abiotic conditions generally complied with existing guidelines, except for the pH in the Cerro Gato basin (<6.5) and orthophosphate concentrations in 2016 (>0.12 mg P L⁻¹). Macroinvertebrate communities were characterized by low family richness (3–11) and were dominated by Atyidae or Chironomidae, thereby causing low diversity (0.33–1.65). Spatial analysis via principal component analysis (PCA) illustrated that abiotic differences between stream basins were mostly related to turbidity, pH, dissolved oxygen (DO), and conductivity. Biotic differences were less clear due to masking by anthropogenic disturbances and dispersal limitations, yet indicated a negative effect of reduced pH and DO on Atyidae presence. In 2017, significantly narrower ranges were found for turbidity, temperature, pH, and diversity (p < 0.01), suggesting a decrease in habitat variability and a need for conservation measures, including mitigating measures related to dam construction for water extraction. As such, further follow-up is highly recommended for the sustainable development and environmental protection of this unique archipelago.

1. Introduction

Freshwater systems located on geographically isolated islands represent sensitive environments prone to irreversible changes when experiencing an unexpected disturbance [1]. Human colonization of these remote areas has occurred throughout history, relying on ground- and/or rainwater to provide for daily freshwater needs. Currently, population growth, tourism, and climate change are major factors affecting both freshwater quantity and quality [2,3,4,5]. For instance, land conversion for human purposes limits rainwater interception and replenishment of water resources, while allowing evaporation to occur at higher rates, hence altering water quantity [3,6]. Water quality, on the other hand, is affected by domestic wastewater, fertilizers, and pesticides, with groundwater levels receiving additional pressure from rising sea levels [2]. As quantity and chemical quality change, ecological quality changes, too, though with a higher degree of irreversibility on isolated island(s) due to the oceanic migration barrier and potentially high degrees of endemism [3,7,8]. Consequently, freshwater diversity on islands is often unbalanced in terms of macroinvertebrate communities and characteristically different from the mainland, especially when considering volcanic islands (e.g., Macaronesia, Galapagos) that have never been in physical contact with the mainland [1,9].

Monitoring freshwater systems returns information on prevailing conditions and supports spatiotemporal analyses, given that appropriate baseline information exists. A myriad of multivariate statistical analyses is available to improve interpretations of spatial analysis based on highly dimensional monitoring data, including correspondence analysis (CA), hierarchical clustering (HC), principal component analysis (PCA), and linear discriminant analysis (LDA). For instance, Li et al. [10] applied PCA and HC as part of a spatial analysis within a lake during and following a period of drought, observing that spatial differences did occur. Similarly, by applying PCA and HC, Rakotondrabe et al. [11] found that sites were affected differently by gold mining, while Ferrier et al. [12] obtained separate clusters of urban rivers on the one hand and highland rivers on the other hand. Moreover, Ferrier et al. [12] complemented their study with a temporal assessment via smoothing splines, obtaining a clear increase in nitrate concentrations. In addition to these indirect gradient analyses, direct gradient analyses are often performed to extract patterns between a set of explanatory variables and one (or more) response variables. For instance, Forio et al. [13] applied a combination of canonical correspondence analysis (CCA) and a multivariate linear regression model (LRM) to infer that velocity, sediment composition, conductivity, and dissolved oxygen affect the prevailing macroinvertebrate community. In contrast, when studying the macroinvertebrate-based ecological water quality in the Guayas basin (Ecuador), Damanik-Ambarita et al. [14] derived the relatively high importance of flow velocity, sediment type, and conductivity from CA following sample clustering based on the observed biological water quality.

Tropical islands are generally believed to host limited macroinvertebrate diversity while facing increasingly growing pressures due to anthropogenic activities, causing biological invasions, wildlife exploitation, and cultivation [5,9,15,16]. Nevertheless, studies on related effects remain limited, despite their potential in steering the development of management plans [5]. A particular example is the Galapagos archipelago, which was placed on the “In Danger” list of the UNESCO World Heritage Sites in 2007 due to overdevelopment, followed by removal from that list in 2010 [17,18]. The archipelago has been subject to a boost in tourism, with about 1000 tourists per year in 1960 increasing to 241,800 in 2017, causing people to migrate from the mainland to the islands to profit from the higher life standards there [4,19]. As population and tourism increase, higher pressures are exerted on local ecosystems, leading to an elevated interest in researchers to study processes, pressures, and impacts within the archipelago, driven by (i) the ecosystems’ high level of unicity [20] and/or (ii) the need for information to support adequate decision-making by local authorities. In contrast to the vast amount of research related to marine [21,22,23] or terrestrial ecosystems [24,25], recent information about inland water systems remains scarce [7,26,27]. For instance, Colinvaux [26] even stated that, when searching for freshwater systems five decades ago, no permanent streams existed on the islands. Moreover, the demand for agricultural products to support the increasing local population requires the application of fertilizers and pesticides, which act as additional potential pressures on the available fresh water and highlight the need for a baseline assessment [2,8]. As such, the collected information on the insular freshwater system allows for a general description of the prevailing conditions and an assessment of spatial patterns. Moreover, following temporal changes in tourism, the climate, and related threats, a temporal evolution in these conditions is expected.

This research provides an assessment of the ecological water quality of the prevailing freshwater systems to provide local authorities on the island of San Cristóbal (the only island with significant permanent fresh water [28]) with baseline information to support adequate management. Within this assessment, we aim to (1) provide an overview of the freshwater systems both chemically and biologically, (2) determine spatial patterns between stream basins, and (3) identify temporal differences over two consecutive years. These findings contribute to a larger study project that aims to understand spatiotemporal variability in water quality, enriching the current limited knowledge on freshwater conditions within tropical islands in general and in the Galapagos in particular. Furthermore, they support local water management by providing a basis for developing stream management plans to avoid potential future degradation.

2. Materials and Methods

2.1. Site Description and Data Collection

The Galapagos archipelago is about 3 to 5 million years old and consists of 13 major and 7 minor islands, consisting of a land area of about 8000 km² [4,28]. Only four islands are permanently inhabited: San Cristóbal, Santa Cruz, Isabela, and Floreana, with the first being the oldest island and containing the only significant permanent source of fresh water within the archipelago [28]. Additionally, small brooks can be found year-round on San Cristóbal, running from the highlands toward the ocean.

Sampling campaigns were organized on San Cristóbal in the dry season (running from June until December, with the main wet season from January until May [29]) of October 2016 and July 2017 to collect chemical and biological data from 10 sampling sites located in three groundwater-fed basins: Toma de la Policía (P), Chino basin (C), and Cerro Gato basin (G). Sampling sites were selected along the stream paths based on accessibility, the surrounding land use, and sampling permit specifications, with the lowest site numbering assigned to the most upstream location (see Figure 1). The P-basin is characterized by a main water extraction site between sites P2 and P3, with an additional barrier having been constructed directly upstream of site P4 between 2016 and 2017. Similarly, within the C-basin, a small extraction site is situated directly downstream of site C2, and a small waterfall (±10 m) is present directly upstream of site C3. Within the G-basin, a vegetated wetland occurs upstream of site G1. All methods were carried out in accordance with the relevant guidelines and regulations of the Galapagos National Park Directorate under research permit PC-02-19. All experimental protocols were reviewed and approved by a Galapagos National Park Directorate committee, which assesses animal care in research activities.

Water temperature (T, °C), acidity (pH, -), conductivity (EC, µS cm⁻¹), and dissolved oxygen (DO, mg O₂ L⁻¹) were measured with handheld multiprobes (2016: YSI (9620v2-2 version 1.2, Elscolab Belgium, YSI Xylem Inc., Yellow Springs, Ohio, USA); 2017: Aquaread (AP5000 version 4.07, Eijckelkamp The Netherlands, Aquaread Ltd, Broadstairs, Kent, UK)). Turbidity (nephelometric turbidity units (NTUs)) was measured with a 2100Q portable turbidimeter (Hach; 2016) (Hach USA, Loveland, Colorado, USA) and via multiprobe measurement (Aquaread AP5000; 2017). Calibration was performed prior to sampling and was supplemented with a consistency check after sampling. Chemical oxygen demand (COD, mg O₂ L⁻¹) and nutrient status (ammonium (NH₄⁺, mg N L⁻¹), nitrite (NO₂⁻, mg N L⁻¹), nitrate (NO₃⁻, mg N L⁻¹), and orthophosphate (o-PO₄³⁻, mg P L⁻¹)) were determined spectrophotometrically (COD: Hach DR3900 (version 1.15, Hach Belgium, Hach Lange GmbH, Berlin, Germany); nutrients: ThermoFisher Genesys 10S UV-VIS (ThermoFisher Scientific, Madison, Wisconsin, USA) by using low-range Hach (LCK, Hach Belgium, Hach Lange GmbH, Düsseldorf, Germany) and Merck test kits (Spectroquant, VWR Belgium, Merck KGaA, Darmstadt, Germany), respectively (Appendix A,

Table A1. List of physicochemical variables measured in 2016 and 2017 and the associated method.

Variable	Technique	Remark(s)
Temperature	Probe	2016: YSI V6920; 2017: Aquaread AP5000
Conductivity	Probe	2016: YSI V6920; 2017: Aquaread AP5000
pH	Probe	2016: YSI V6920; 2017: Aquaread AP5000
Dissolved oxygen	Probe	2016: YSI V6920; 2017: Aquaread AP5000
Turbidity	2016: Turbidimeter 2017: Probe	2016: Hach 2017: Aquaread AP5000
Chemical oxygen demand	Test kit code 2415815	Hach, limits: 0.7–40.0 mgO2 L−1; follows EPA 5220 D
Ammonia-nitrogen	Test kit code 114752	Merck, limits: 0.010–3.00 mg N L−1; follows EPA 350.1, APHA 4500-NH3 F, ISO 7150-1, and DIN 38406-5
Nitrite-nitrogen	Test kit code 114776	Merck, limits: 0.002–1.00 mg N L−1; follows EPA 354.1, APHA 4500-NO2-B, and DIN EN 26 777
Nitrate-nitrogen	Test kit code 109713	Merck, limits: 0.1–25.0 mg N L−1; follows DIN 38405-9
Orthophosphate-phosphorus	Test kit code 114848	Merck, limits: 0.0025–5.00 mg P L−1; follows EPA 365.2 + 3, APHA 4500-P E, and DIN EN ISO 6878

). Biological samples were collected by executing the kick-net sampling technique according to De Pauw and Vanhooren [30]. This technique consists of walking upstream in a backward way while kicking sediments loose and capturing the dislodged macroinvertebrates in a hand net (0.15 m × 0.15 m metal frame, 500-µm mesh size) within a 5-min timeframe and proportionally covering all available habitats. Collected macroinvertebrates were live-sorted (after sample-sieving using a 500-µm mesh size) and put in >70% ethanol to prevent degradation prior to identification up to the family level (referred to as “taxon”). Subsequently, the Biological Monitoring Working Party score adapted for Colombia (BMWP_Col) [14,31,32], overall richness, and Shannon diversity (Equation (1)) were calculated [33,34]. Data collection methods were similar between sampling sites and years, allowing for comparisons of the obtained results.

$$H^{\prime }= -{\displaystyle \sum }_i{p}_i\times \ln p_i$$

(1)

Here, H’ is the Shannon entropy representing diversity, i is the taxon, and p_i is the fraction of all individuals represented by taxon i. Low values represent low evenness and dominance by a subset of species present.

2.2. Data Exploration and Analysis

Data exploration was performed in MS Excel^® and RStudio [35,36]. Prior to analyses, data were checked for the observed ranges, and values below the detection limit were replaced by the lowest detection limit (i.e., a conservative approach). Principal component analysis (PCA) was performed on the normalized abiotic data [37] to identify potential differences between the three basins and the two sampling years. Similarly, PCA was applied to the macroinvertebrate community composition after being transformed according to Equation (2) (Hellinger transformation):

$$y_{i,j}^{\prime }=\sqrt{\frac{y_{i,j}}{{{\displaystyle \sum }}_j{y}_{i,j}}}$$

(2)

where y is the original value, y’ is the transformed value, i is the sample site, and j is the species.

Redundancy analysis (RDA) was performed on the combined data (original abiotic data and transformed biotic data) to link the observed communities with prevailing environmental conditions. Permutation Analysis of Variance (ANOVA) was run to test the significance of the RDA axes. PCA and RDA were preferred over correspondence analysis (CA) and canonical CA (CCA), as (1) a preliminary detrended correspondence analysis (DCA) resulted in an axis length shorter than 3 units, suggesting the use of a linear ordination [38]; (2) the Hellinger distance provides a better trade-off between resolution and linearity, compared to chi-squared-based techniques [37]; and (3) CCA is not recommended to be used if data richness covers a wide range [39]. The transformation of the biological data reduces the skewed abundance distribution between sampling sites among sampling sites [37] while allowing for a Euclidean-based distance analysis to be performed, with similar results as with CCA [39].

The abiotic-based PCA biplot was investigated for clustered sampling sites with specific attention on clustering based on (1) basin (P, C, and G) or (2) sampling year (2016 and 2017). When shifts in year-based clusters were identified along one of the axes, the correlated variables (>0.5) were selected for further analysis. For each of the selected variables, the difference in observed values between both sampling years (∆₁₇₋₁₆) was determined for each location and subsequently linked with the value observed in 2016 via bivariate linear regression (see Equation (3)):

$$\Delta {\mathit{X}}_{17-16}=\alpha +\beta ·{\mathit{X}}_{16}+\mathbf{\epsilon }$$

(3)

where X is the selected variable (list of 10 observations), α the intercept, β the slope, and ε the residuals of the linear regression.

Significant slope values (β, p < 0.05) indicated an uneven change in the variable’s value, reflecting (1) a convergence when slope values were negative and (2) a divergence when slope values were positive. The intensity of this behavior was reflected by the absolute slope value. Similarly, the biotic-based PCA and overall RDA were analyzed for the presence of clusters. When clusters or related shifts could not be clearly identified, no further analysis was performed and differences were considered to have occurred by chance. Additionally, richness, BMWP_Col, and Shannon diversity were analyzed via linear regression to determine potential converging or diverging behavior. Whenever the slope was identified to be nonsignificant (p > 0.05) and generally negative or positive differences were found, a paired t-test was performed. Residuals were plotted to assess the potential presence of patterns excluded from the linear model. The R packages used for data analysis were vegan [40] and BiodiversityR [41]. It should be noted that this approach did not represent a complete temporal analysis, as it only considered two consecutive years and excluded additional temporal variation (e.g., diel, monthly, and seasonal patterns). Nevertheless, it could provide a first signal for any temporal evolution being present, allowing for the appropriate response (e.g., early warnings requesting increased monitoring and/or improved protection).

3. Results

3.1. Baseline Data and Spatial Analysis

3.1.1. Abiotic Conditions

Relatively low conductivity values were observed, ranging from 69 up to 156 µS cm⁻¹ in 2016 and between 34 and 130 µS cm⁻¹ in 2017, while pH fluctuated around 7.5 on the National Institute of Standards and Technology (NIST) scale (overall range between 6.06 and 8.51). The lowest dissolved oxygen (DO) concentration (4.62 mg O₂ L⁻¹) was observed in site G1 in 2017, but it was mostly higher than 6 mg O₂ L⁻¹. Turbidity ranged from 2.6 up to 46.4 NTU in 2016, but presented a narrower range in 2017 (4.1 up to 23.7 NTU). In contrast, COD was characterized by a wider range in 2017 (5.8–21.2 mg O₂ L⁻¹) compared to 2016 (2.7–9.2 mg O₂ L⁻¹). Ammonium levels were often below the detection limit (0.010 mg N L⁻¹; five sites in 2016, four sites in 2017), as were nitrate levels in 2016 (0.1 mg N L⁻¹; nine sites) and orthophosphate levels in 2017 (0.0025 mg P L⁻¹; nine sites). Consequently, nitrate and orthophosphate were excluded from the data, which were subsequently normalized prior to principal component analysis (PCA). The resulting primary and secondary axes explained 32% and 20% of the overall variance, respectively, with pH, DO, and turbidity negatively correlated with the first axis (−0.57, −0.51, and −0.46, respectively), while conductivity and temperature were negatively correlated with the second axis (−0.62 and −0.55, respectively). A summary of the observed range and mean values is provided in

Table 1. Physicochemical variables measured in streams on San Cristóbal island (Galapagos). The observed range (Min–Max), mean, standard deviation (sd), and number of sites below the detection limit (n_BDL) are reported for the campaigns in 2016 and 2017, analyzing the same sites (N = 10). Probe measurements were never below the operational limits (reported as “-” for n_BDL). NA: Not available due to a high number of sites below the detection limit; NTU: Nephelometric turbidity unit.

Variable (Unit)	Basin	2016					2017
		Mean	sd	Min	Max	nBDL	Mean	sd	Min	Max	nBDL
Temperature (°C)	P	18.6	0.1	18.6	18.7	-	20.5	0.1	20.3	20.6	-
	C	20.5	1.1	19.3	21.5	-	21.0	0.6	20.5	21.6	-
	G	18.7	0.4	18.4	19.0	-	19.7	0.5	19.3	20.0	-
	Mean	19.2	1.0	18.4	21.5	-	20.5	0.6	19.3	21.6	-
Conductivity (µS cm−1)	P	84	21	72	122	-	63	7	58	72	-
	C	139	15	129	156	-	126	3	124	130	-
	G	77	11	69	84	-	38	5	34	42	-
	Mean	99	32	69	156	-	77	36	34	130	-
pH (-)	P	8.32	0.17	8.15	8.51	-	7.59	0.09	7.48	7.72	-
	C	7.59	0.31	7.30	7.91	-	7.41	0.21	7.23	7.64	-
	G	6.14	0.11	6.06	6.21	-	6.43	0.09	6.37	6.49	-
	Mean	7.67	0.90	6.06	8.51	-	7.30	0.48	6.37	7.64	-
Dissolved oxygen (mg O2 L−1)	P	8.66	0.23	8.39	8.92	-	8.69	0.09	8.54	8.77	-
	C	7.15	0.97	6.36	8.23	-	7.77	0.84	6.82	8.39	-
	G	6.88	0.35	6.63	7.13	-	5.30	0.96	4.62	5.98	-
	Mean	7.85	0.99	6.36	8.92	-	7.74	1.45	4.62	8.77	-
Turbidity (NTU)	P	25.3	15.6	10.2	46.4	-	18.4	4.1	14.8	23.7	-
	C	6.4	3.6	2.7	9.9	-	8.5	5.8	4.1	15.1	-
	G	5.3	3.7	2.6	7.9	-	5.8	1.0	5.1	6.5	-
	Mean	15.6	14.7	2.6	46.4	-	12.9	7.1	4.1	23.7	-
COD (mg O2 L−1)	P	6.8	2.5	2.7	9.2	0	7.7	1.1	5.8	8.9	0
	C	4.8	2.2	2.8	7.1	0	16.7	4.0	13.6	21.2	0
	G	3.8	0.1	3.7	3.9	0	14.0	2.5	12.2	15.7	0
	Mean	5.6	2.3	2.7	9.2	0	11.6	4.8	5.8	21.2	0
Ammonium (mg N L−1)	P	0.076	NA	0.076	0.076	4	0.066	0.046	0.014	0.100	2
	C	0.048	0.012	0.039	0.056	1	0.074	0.072	0.023	0.125	1
	G	0.038	0.016	0.027	0.049	0	0.064	NA	0.064	0.064	1
	Mean	0.049	0.020	0.027	0.076	5	0.069	0.044	0.014	0.125	4
Nitrite (mg N L−1)	P	0.026	0.006	0.020	0.037	0	0.006	0.003	0.002	0.010	0
	C	0.017	0.002	0.016	0.019	0	0.020	0.018	0.006	0.040	0
	G	0.012	0.001	0.011	0.013	0	0.008	0.001	0.007	0.009	0
	Mean	0.021	0.010	0.011	0.037	0	0.011	0.011	0.002	0.040	0
Nitrate (mg N L−1)	P	NA	NA	NA	NA	5	0.52	0.39	0.24	0.96	2
	C	0.67	NA	0.67	0.67	2	0.88	0.25	0.67	1.16	0
	G	NA	NA	NA	NA	2	0.59	0.02	0.57	0.60	0
	Mean	0.67	NA	0.67	0.67	9	0.67	0.30	0.24	1.16	2
Orthophosphate (mg P L−1)	P	0.68	0.67	0.16	1.43	0	0.021	NA	0.021	0.021	4
	C	0.29	0.21	0.15	0.53	0	NA	NA	NA	NA	3
	G	0.12	NA	0.12	0.12	1	NA	NA	NA	NA	2
	Mean	0.49	0.54	0.12	1.43	1	0.021	NA	0.021	0.021	9

, while the abiotic-based PCA can be found in Figure 2A.

Six different clusters were identified, representing the two sampling years and three different basins (Figure 2A). Abiotic differences were greater between than within basins and showed an additional temporal influence. For instance, conductivity was generally higher within the C-basin compared to the P- and G-basins in 2016 (139 ± 15 µS cm⁻¹ versus 84 ± 21 µS cm⁻¹ and 77 ± 11 µs cm⁻¹, respectively) and 2017 (126 ± 3 µS cm⁻¹ versus 63 ± 7 µS cm⁻¹ and 38 ± 5 µS cm⁻¹, respectively), with variability between sites being lower in 2017 (Table 1, Figure 2A). Similarly, turbidity scores were higher in the P-basin compared to the C- and G-basins in 2016 (25.3 ± 15.6 µS cm⁻¹ versus 6.4 ± 3.6 µS cm⁻¹ and 5.3 ± 3.7 µs cm⁻¹, respectively) and 2017 (18.4 ± 4.1 µS cm⁻¹ versus 8.5 ± 5.8 µS cm⁻¹ and 5.8 ± 1.0 µs cm⁻¹, respectively), additionally showing reduced variability in 2017 (Table 1, Figure 2A). Moreover, both pH and DO were lower within the G-basin compared to the P- and C-basins, with the highest values in the P-basin: In 2016, the average pH was 6.14 ± 0.11, lower than that observed in the P-basin (8.32 ± 0.17) and C-basin (7.59 ± 0.31), with an analogous pattern in 2017 (6.43 ± 0.09 versus 7.59 ± 0.09 and 7.41 ± 0.21, respectively). Similarly, dissolved oxygen was lower in 2016 (6.88 ± 0.35 mg O₂ L⁻¹ versus 8.66 ± 0.23 mg O₂ L⁻¹ and 7.15 ± 0.97 mg O₂ L⁻¹, respectively) and 2017 (5.30 ± 0.96 mg O₂ L⁻¹ versus 8.69 ± 0.09 mg O₂ L⁻¹ and 7.77 ± 0.84 mg O₂ L⁻¹, respectively). At the level of sampling year, some minor shifts reflected that, for instance, P-sites shifted in the positive direction of PC1, linked with lower pH, DO, and/or turbidity values. In contrast, C-sites shifted in the negative direction of PC2, linked with increased conductivity and/or temperature, while G-sites showed a combined shift in the positive direction of PC1 and negative direction of PC2.

3.1.2. Biotic Conditions

In total, 4381 organisms belonging to 19 taxa were collected in 2016, the majority of which were Atyidae (2197 individuals, 50%) and Chironomidae (1730 individuals, 39%). Yet, only 2950 organisms were collected in 2017, belonging to 18 taxa, with Atyidae (1279 individuals, 43%) and Chironomidae (1019 individuals, 35%) being again the most abundant. The remaining 11% (2016) and 22% (2017) of the overall macroinvertebrate count were occupied by other taxa, including Dytiscidae (2016: 1.3%; 2017: 8.7%), Hyalellidae (2016: 3.3%; 2017: 2.3%), and Veliidae (2016: 1.4%; 2017: 1.3%), with abundances being relatively low compared to Atyidae and Chironomidae, yet no EPT (Ephemeroptera, Plecoptera, and Trichoptera) taxa were observed. Chironomidae were encountered in each sampling site, while Atyidae were only found in the P- and C-basins and were absent in the G-basin (Appendix B,

Table A2. List of macroinvertebrate taxa, their tolerance score, and abundance at each of the 10 locations sampled in 2016 and 2017. Tolerance scores (TS) were based on Roldán [31] and Álvarez Arango [32].

Taxon	TS	2016										2017
		P1	P2	P3	P4	P5	C1	C2	C3	G1	G2	P1	P2	P3	P4	P5	C1	C2	C3	G1	G2
Atyidae	8	253	349	366	191	313	88	1	636			221	341	595	8	2	13	3	96
Lymnaeidae	8									2										1
Palaemonidae	8																		1
Planorbidae	8						1					1				8	2	3
Coenagrionidae	7			2		2				15	6					3				23
Hyalellidae	7	98	46									7	60
Simuliidae	7	1					4						2	1	4	50	2
Veliidae	7						11		44		6						10		25	3
Aeshnidae	6	5		3	2	10	1				2			1	1	15	9	12		2	1
Dugesiidae	6	1
Pleidae	6									2
Staphylinidae	6					1
Ceratopogonidae	5															1
Gyrinidae	5	12	14	31	9	1	3						3	3		2
Mesoveliidae	5						3
Notonectidae	5							15										25	1
Hydrophilidae	3				1	1				5	3	1			6	6	5		2	2
Chironomidae	2	346	1	159	56	198	145	13	9	683	120	84	27	57	7	327	175	115	22	168	37
Tubificidae	1	1					1	1	2									6		7
Dytiscidae	-	1	1			1				29	26		7		3	3	1	5		230	7
Lumbricidae	-											5	6	10			10
Lumbriculidae	-	9		8		4	7					11	9	9		12			7

). In contrast, higher counts of Dytiscidae were observed within the G-basin (up to 30 times higher) compared to the P- and C-basins, both containing sites without Dytiscidae. Several taxa only occurred in a single sampling site, including Dugesiidae (2016, P1), Mesoveliidae (2016, C1), Pleidae (2016, G1), Staphylinidae (2016, P5), Ceratopogonidae (2017, P5), and Palaemonidae (2017, C3), while others only occurred within a specific basin (e.g., Hyalellidae in the P-basin, Notonectidae in the C-basin, and Lymnaeidae in the G-basin) or were absent within a specific basin (e.g., no Veliidae in the P-basin, no Coenagrionidae in the C-basin, and no Simuliidae in the G-basin). A detailed overview of the site-specific community composition can be found in Appendix B, Table A2.

The principal component analysis of the biotic data showed clusters and shifts that were less clear compared to the abiotic analysis, with the primary axis explaining 61% and the secondary axis 11% of the overall variance. The first axis was characterized by a negative correlation with Atyidae and a positive correlation with Chironomidae (−0.79 and 0.51, respectively), while the second axis represented a positive correlation with Notonectidae (0.70) and a negative correlation with Dytiscidae (−0.59) (Figure 2B). Absolute correlation values of the remaining taxa were lower than 0.5 and were excluded from the graphical representation to ease interpretation. The G-sites were distinctly situated to the right, reflecting the absence of Atyidae, while additionally being linked with elevated Dytiscidae presence (Figure 2B). In contrast, the P-sites were mostly situated to the left of the biplot (Figure 2B), suggesting Atyidae-dominated sites, with only limited influence by other taxa, except for sites P4 and P5 in 2017. The C-sites were scattered within the biplot (Figure 2B), with C2 clearly in the top-right corner due to relatively high Notonectidae counts in a Chironomidae-dominated community. In contrast, C3 was clearly Atyidae-dominated (left side of the biplot, Figure 2B), while C1 shifted from an almost balanced situation in 2016 to being Chironomidae-dominated in 2017.

The lowest taxa richness was observed at site G2 in 2017 (3 taxa), resulting in a BMWP_Col score of 8 (very bad quality), while the highest richness was observed at site P5 in 2017 (11 taxa), with an accompanying BMWP_Col score of 51 (poor quality). Similarly, these sites scored relatively low (0.53) and high (0.94) regarding Shannon diversity, but did not represent the minimum and maximum of the diversity range, respectively. The lowest diversity scores were obtained for sites C3 and G1 in 2016 (respectively 0.33 and 0.34), while multiple sites scored higher than P5 in 2017, including sites P4 in 2017 (1.65) and P1 in 2016 (1.18). The average basin-specific metric scores were lowest for the G-basin, although clear differences could not be derived. Moreover, within the G-basin, BMWP_Col and richness declined from 2016 to 2017 (

Table 2. Summarizing metrics of macroinvertebrate community within streams on San Cristóbal island (Galapagos). Metrics include richness (number of taxa), the BMWP_Col score (taxa-based tolerance score), and diversity (combining richness and abundance). The observed range (Min–Max) and arithmetic mean (mean) for each metric–year combination as well as the mean metric value for each basin are reported.

Basin	Code	Richness		BMWP_Col		Diversity
		2016	2017	2016	2017	2016	2017
Toma de la Policía	P1	10	7	42	28	1.18	0.91
	P2	5	8	22	29	0.53	0.91
	P3	6	7	28	28	0.91	0.48
	P4	5	6	24	26	0.73	1.65
	P5	9	11	37	51	0.86	0.94
	Mean	7	7.8	30.6	32.4	0.842	0.978
Chino	C1	10	9	49	41	1.15	0.96
	C2	4	7	16	30	0.94	1.10
	C3	4	7	18	33	0.33	1.13
	Mean	6	7.7	27.7	34.7	0.807	1.063
Cerro Gato	G1	6	8	26	34	0.34	1.02
	G2	6	3	25	8	0.89	0.53
	Mean	6	5.5	25.5	21	0.615	0.775
Total	Min	4	3	16	8	0.33	0.48
	Max	10	11	49	51	1.18	1.65
	Mean	6.5	7.3	28.7	30.8	0.786	0.964

), with only diversity increasing due to a decrease in dominance by Chironomidae (i.e., 683 to 186 individuals at site G1). In contrast, all metric scores for the P- and C-basins increased from 2016 to 2017, supporting an overall average increase between the two consecutive years (Table 2), but none of the metrics showed significant differences (all p-values > 0.05).

The combination of both datasets (abiotic and biotic) showed (after applying a redundancy analysis) a pattern similar to Figure 2B, with overlaid environmental gradients allowing for the inference of potential relationships between prevailing macroinvertebrate dominance and abiotic conditions. For instance, an increase in DO and/or pH had a positive effect on the Atyidae count and was, consequently, negatively correlated with the Chironomidae count (Figure 3). Similarly, an increase in conductivity and/or temperature had a positive effect on Notonectidae presence, while in contrast, reduced values of DO (and/or pH) and conductivity (and/or temperature) provided better conditions for Dytiscidae (Figure 3). The presence of Atyidae was linked to higher DO and pH values, while Chironomidae was linked to increased values of ammonium as well as reduced values of DO and pH (Figure 3). The primary axis of this constrained analysis explained 47% of the overall variance, while the secondary axis explained an additional 8%, and yet only the first axis was found to be significantly important (ANOVA permutation, p = 0.004), suggesting a relatively high importance of DO in defining the prevailing macroinvertebrate community within the sampled sites.

3.2. Temporal Evolution

The overall rightward shift of the sampling sites along the primary axis (Figure 2A) had its origin in generally lower pH values (8 out of 10 sites) and partially lower turbidity values (4 out of 10 sites). The former was characterized by a narrower range in 2017, with initially high pH values decreasing and vice versa regarding low pH values (Table 1, Figure 4A). Similarly, high turbidity values in 2016 were not observed in 2017, reflecting a reduction in the more extreme conditions (Figure 4B). For both, a significant decreasing behavior could be observed (slope of linear regression significantly different from 0,

Table 3. Results of linear analyses relating the observed differences between 2016 and 2017 (∆₁₇₋₁₆), with the observations made in 2016 as represented in Figure 4, Figure 5 and Figure 6. Slope values are given with their accompanying p-value and 95% confidence interval (CI). The R² value describes the model fit to the observations. The significance of the slope value is indicated for p < 0.01 (**) and for p < 0.001 (***).

Variable	Slope	p-Value	R2	CI 2.5%	CI 97.5%
pH	−0.49 ***	0.00002	0.90	−0.61	−0.36
Turbidity	−0.66 ***	0.0005	0.79	−0.93	−0.38
DO	0.085	0.8	0.01	−0.71	0.88
Temperature	−0.59 **	0.003	0.68	−0.92	−0.26
Conductivity	−0.032	0.9	0.00	−0.48	0.42
Richness	−0.62	0.06	0.37	−1.30	0.03
BMWP_Col	−0.57	0.1	0.27	−1.30	0.19
Diversity	−1.30 **	0.008	0.61	−2.1	−0.44

). In contrast, no pattern was observed for dissolved oxygen (DO), which showed similar oxygen concentrations in 2017 compared to 2016 (Figure 4C, Table 3).

C- and G-sites showed a slight downward shift within the abiotic PCA (Figure 2A), suggesting an increase in temperature and/or conductivity. Indeed, increased water temperatures were observed in 2017, with a similar maximum value but a higher minimum value (Table 1 and Figure 5A). The significant slope coefficient (−0.59, p = 0.003; Table 3) supported the observation of a narrowing range, with low values increasing and high values remaining almost constant. However, no pattern in conductivity values was observed, with generally lower values in 2017 (Figure 5B). These values were confirmed to be significantly lower based on the performed paired t-test (t = 3.95, df = 9, p = 0.002).

No general clusters or patterns of stream basins were observed when investigating species composition of the different sites (Figure 2B), with some sites shifting upward along the secondary axis and others moving down. Sites shifting away from the primary axis reflected a more complex community composition than sites situated on or near the primary axis. Richness, BMWP_Col, and diversity each showed a higher value in 2017 for six sites, while the remaining four sites showed a decrease. A regression analysis of richness (Figure 6A), BMWP_Col (Figure 6B), and diversity (Figure 6C) showed that biological convergence (i.e., the simultaneous increase of initially lower and decrease of initially higher biotic metric values) occurred on the island. This convergence was only significant for the Shannon diversity index (slope = −1.3, p = 0.008; Table 3) and was marginally nonsignificant for overall richness (slope = −0.62, p = 0.06; Table 3). In contrast, the BMWP_Col score was not observed to have changed significantly (slope = −0.57, p = 0.1; Table 3).

4. Discussion

4.1. Baseline Data and Spatial Analysis

The observed temperature and dissolved oxygen complied with the water quality criteria set by the United States Environmental Protection Agency (USEPA) [42]. COD values remained relatively low due to the limited distance covered by the sampled streams, causing typical mid- and downstream sections (and, hence, the plausibility of reduced quality) to be absent [13]. Nitrate and nitrite concentrations remained well below the limits for drinking water (i.e., 10 mg N L⁻¹ and 1 mg N L⁻¹, respectively) [43]. In contrast, orthophosphate concentrations in 2016 exceeded the suggested threshold to avoid eutrophication (0.1 mg P L⁻¹), while pH within the G-basin was also below the quality guideline (i.e., 6.5–9.0) [42]. The latter was related to the presence of a small, vegetated wetland upstream of the sampling sites, causing decomposition of organic matter in low-flow conditions and related reduced pH and oxygen concentrations. The elevated phosphate concentrations in 2016 increased the risk of eutrophication, but the simultaneous low levels of nitrate suggested a nitrogen deficiency in 2016. However, the indication of a potential phosphorus deficiency in 2017 hinted at an underlying temporal nutrient dynamic that requires further study.

Sampling sites were observed to be mainly situated along an Atyidae/Chironomidae axis, explaining up to 61% of the overall variability in community composition among sites. Atyidae are common inhabitants of littoral drainages within the South American continent and on tropical islands [7,9,13], and hence their survival within the sampled waters is expected and adds to the observations of Peck [7], who reported the overall presence of 52 crustacean species within the archipelago, 2 of which belonged to the Atyidae family (i.e., Typlatya galapagensis and Potimirim glabra cf., the former being endemic). Their absence from the G-basin was linked to the prevailing reduced water quality (low pH and dissolved oxygen), as Atyidae are relatively sensitive to pollution (i.e., tolerance score of 8, with 0 reflecting high tolerance and 10 reflecting high sensitivity, see Appendix B, Table A2). In contrast, Chironomidae are highly tolerant (i.e., tolerance score of 2) and are widely observed within insular and mainland freshwater systems [13,14,44], hence their occurrence in each sampling site. Dytiscidae presence within the Galapagos archipelago was reported three decades ago by Peck and Kukalová-Peck [45], along with 53 other beetle families (including Gyrinidae and Hydrophilidae). Additionally, the presence of both Coenagrionidae (damselfly) and Aeshnidae (dragonfly) is also supported by Peck [46], who reported the presence of seven native and one endemic species in the Galapagos. Yet, diversity remained low, which was clearly reflected in the low richness and diversity scores at C3 (4 and 0.33, respectively) and G1 (6 and 0.34, respectively). Both locations were characterized by an extreme dominance of Atyidae (92% of all organisms) and Chironomidae (93% of all organisms), respectively, with the latter representing the overall highest number of individuals found within a single site (683, Table A2). G-sites were clearly Chironomidae-dominated as a consequence of reduced water quality (low pH and DO), while C-sites shifted longitudinally from being Chironomidae-dominated to Atyidae-dominated between C2 and C3. Site C2 had a slightly lower pH and DO, which in combination with a small dam and waterfall acted as a migration barrier for upstream migration of Atyidae, while Chironomidae could benefit from their nonaquatic life stage to disperse. Similarly, a temporal shift was observed for sites P4 and P5, moving from Atyidae-dominated (characteristic of the P-basin, Figure 2B) in 2016 to a balanced (P4) and Chironomidae-dominated (P5) situation in 2017, illustrating the effects of human disturbance. Indeed, between the campaigns, minor construction works were performed directly upstream of P4 (which consisted of substrate removal and the creation of a small dam), thereby decimating the Atyidae population in downstream sites.

Despite the abiotic differences between basins, no clear biotic-based clusters could be defined due to the general dominance of Atyidae or Chironomidae. These taxa occurred in each sampling site except for the G-basin (Atyidae) due to (i) reduced water quality and (ii) dispersal limitations. The former was supported by observations in the C-basin showing decreased Atyidae presence in site C2 compared to C1 and C3, as this site was also characterized by slightly decreased water quality (see sites G1 and G2). In addition, Atyidae have a completely aquatic life cycle, which impedes easy dispersal to a nonconnected basin (i.e., from P or C to G). The relatively high explanatory power of the first principal axis of the biotic PCA (Figure 2B) suggested that differences in insular macroinvertebrate communities were mainly influenced by biotic processes including (1) dispersal, (2) competition, and (3) stochastic extinctions [1,8]. Additionally, (historical) disturbances could have exacerbated these effects by causing sudden alterations in the community composition [9].

With recent (a)biotic data of freshwater systems on San Cristóbal island being limited, observations were compared with the Portoviejo and Guayas rivers, two main river basins located in the coastal zone of Ecuador. Both river basins were assessed recently, following the same assessment methodology and were assumed to be exposed to similar climatic conditions as the basins in San Cristóbal. In general, the observed abiotic ranges reported in Table 1 were narrower than the ranges reported for the Portoviejo river [44], except for temperature and conductivity, which were consistently lower on San Cristóbal (i.e., 25.6–31.3 °C and 164–2447 µS cm⁻¹ within the Portoviejo river). Similarly, broader ranges were obtained in the Guayas river [14], with reported ranges largely covering the observed ranges in Table 1, including temperature and conductivity (i.e., minima of 19.0 °C and 37 µS cm⁻¹). In general, the observed range of the BMWP_Col index was low compared to the overall range found in the Guayas river (ranging from 0 up to 168) and the Portoviejo river (ranging from 16 up to 140) [14,44]. These discrepancies were inherently linked to the spatial extent of both the Portoviejo and Guayas rivers, as the influences of surrounding land use accumulated toward the mouth, extending the range of most variables [13,14,44]. Hence, a comparison to prevailing conditions within mainland basins could be performed, though care should be taken when inferring conclusions. Assuming that the least disturbed conditions within the Portoviejo and Guayas basins were found in upstream sites, e.g., higher than 250 m a.s.l. [14], reflecting conditions similar to the streams on San Cristóbal (limited anthropogenic influence, relatively high shading, and higher flow velocities), no specific differences between mainland and insular freshwater systems in terms of abiotic conditions were observed. Still, there was a clear difference that occurred at a biological level, as most BMWP_Col scores were higher than 100 [14]. Consequently, based on the observed index scores, it could be inferred that, in comparison to mainland river basins, biological quality in freshwater streams was very bad and highly degraded. Yet, as the development of currently existing assessment techniques and biotic indices focused on continental water bodies, both techniques and indices are considered to be only limitedly type-specific, neglecting the uniqueness of insular systems [1]. For instance, the observed low diversity and nearly complete absence of vulnerable species did not necessarily reflect a highly degraded system, but represented a system that was underdeveloped compared to mainland freshwater systems due to the absence of a physical connection with the mainland, the presence of a vast migration barrier, and the relatively small size and elevation of the considered insular system [9,20,46]. Nevertheless, the mismatch between the uniqueness of island systems and the continental-focused indices stresses the need for (i) an appropriate baseline assessment and (ii) the development of more specific indices to summarize prevailing conditions and determine discrepancies between potential and current conditions. With a freshwater biodiversity crisis raging around the world and with insular biodiversity being threatened by biological invasions, wildlife exploitation and cultivation [15,47], the development of relevant monitoring, and conservation programs have never been more urgent, requiring (local) adaptations of assessment schemes and indices developed for mainland systems. Hence, a translation of the obtained BMWP_Col score into a general condition as part of a baseline assessment (e.g., “very bad and highly degraded”) is to be avoided until specific reference conditions have been defined.

4.2. Temporal Evolution and Implications

Among the steering abiotic variables, no significant converging behavior was observed for DO and conductivity. Conductivity was shown to be significantly lower in 2017 as a potential artifact of seasonality (sampling performed in October 2016 versus July 2017): Sampling in 2017 was performed earlier in the dry season (July), and hence lower conductivity values were expected due to reported increases in conductivity toward the end of the dry season [48]. In contrast, the average values of the remaining variables did not clearly change, and yet average richness, BMWP_Col, and diversity increased from 2016 to 2017. Still, site-specific differences suggested an overall convergence of both chemical and biological freshwater status. The discrepancy between the reported ranges in Table 1, which were wider in 2017 compared to 2016, and the observed temporal evolution in Figure 4, Figure 5 and Figure 6 suggesting a narrowing pattern was attributable to outliers. For instance, diversity in 2016 ranged from 0.33 up to 1.18 (Δ = 0.85) while ranging from 0.48 up to 1.65 (Δ = 1.17) in 2017. However, cutting out the outermost values decreased this difference to 0.82 (2016) and 0.65 (2017), confirming that “extreme” values masked the overall observed range narrowing. The evolution toward similar conditions among freshwater sites, situated in different stream basins, decreased overall variability within the island, hence limiting the presence of specific conditions and related specialized organisms. Nevertheless, apart from the converging behavior, an improvement in average biological water quality status occurred, reflected by an increase in average richness, BMWP_Col, and diversity.

With tourism, agriculture, and related anthropogenic activities continuously increasing and exacerbating the effects of climate change [1,4], it is urgent to develop freshwater management plans to preserve (and improve) current conditions. For instance, limiting surface water extraction by relying more on the desalination of seawater (or brackish groundwater) for human consumption, recycling domestic wastewater (after proper treatment), or training water resource management personnel are only a few options to be considered [2]. In contrast, the construction of dams to support small-scale water extraction activities leads to local disturbances that have the potential to wipe out the established community and make future dispersal harder [8,49]. Moreover, the subsequent reduction in flow velocity alters the biotic and abiotic conditions both at upstream and downstream sites, including the prevailing macroinvertebrate assemblages [13,14]. This was illustrated by the observed decrease in turbidity within the P-basin, which was being used for small-scale water extraction via dams. Following these disturbances, species diversity (alpha, beta, and gamma diversity) is expected to decrease as opportunistic and disturbance-tolerant taxa (e.g., high numbers of Atyidae and Chironomidae) are highly likely to disperse and colonize, subsequently outcompeting less-tolerant species and potentially driving them to (unnoticed) extinction [8].

5. Conclusions

Abiotic conditions indicated differences between stream basins but were not strong enough to uniformly affect the prevailing macroinvertebrate community within each basin. Anthropogenic disturbances and dispersal constraints helped in explaining the prevalence of the most common macroinvertebrates, yet for most macroinvertebrates, presence was a mere consequence of stochastic processes and/or (historical) disturbance impeding high biotic diversity. Overall, biotic water quality improved within the considered basins, yet masked the converging behavior of both abiotic conditions and biotic metrics. Further follow-up of abiotic conditions, nutrient status, and macroinvertebrate communities by means of a continuous monitoring scheme is required and recommended to obtain a synthesis of the underlying temporal pattern and to conserve limited freshwater resources in the Galapagos and on islands facing similar pressures (i.e., increased tourism, demographics, and agriculture), not only at the quantitative but also at the qualitative level. Moreover, as specific climatic, meteorological, and hydrographic information is lacking, installment of a meteorological station as well as mapping the different basins is highly recommended to (1) update the data described in this paper and (2) support the implementation, evaluation, and improvement of management actions.