Introduction
Glacierised catchments harbor streams from different water sources, such as glacier-fed (kryal), rainfall/snowmelt-fed (rhithral) and groundwater-fed (krenal) streams (Füreder, 1999; Brown, Hannah & Milner, 2003), whose confluence creates a complex mosaic of reaches with varying levels of glacial influence and environmental conditions. This, in turn, influences the assemblage of local aquatic communities and produces complex patterns of species diversity and composition at the catchment scale. Glacier-fed streams, in particular those closer to glaciers, are expected to produce strong niche selection, because harsh environmental conditions filter non-tolerant species from the regional pool (Chase, 2007). This should produce communities with a more predictable composition and higher proportions of rare taxa, than those from less extreme habitats located farther away from the glacier snout (Leibold et al., 2004; Chase, 2007; Jacobsen et al., 2012). Several studies on glacier-fed streams around the world have found that taxonomic and functional diversity decrease with increasing glacial influence (Milner & Petts, 1994; Milner et al., 2001a; Jacobsen et al., 2010; Jacobsen & Dangles, 2012; Crespo-Pérez et al., 2020), probably due to lower water temperature, and higher turbidity, hydrologic instability, and streambed disturbance closer to the glaciers (Milner et al., 2001a; Crespo-Pérez et al., 2020). Other studies have revealed peaks in local richness at intermediate levels of glacier influence, (Jacobsen et al., 2012; Cauvy-Fraunié et al., 2014a; Quenta et al., 2016), and more diversified communities at the confluences of streams with different water sources (e.g., kryal-rhithral, rhithral-krenal) (Lencioni, Rossaro & Maiolini, 2007).
Commonness, rarity and taxa turnover rates have been the central focus of various ecological studies around the world (e.g., Kunin & Gaston, 1993; Cao, Larsen & Thorne, 2001; Alahuhta et al., 2014). Commonness is defined as the antagonistic definition of uncommon or rare, where rarity can be categorized in two different ways: (1) by total abundances and (2) by occurrence or occupancy in the samples (Josefson, 2009). Lately, with ongoing climate change, and glacial retreat, there has been increased interest in the study of levels of commonness and rarity in invertebrate and microbiota communities in glacierised catchments (e.g., Wilhelm et al., 2014; Robinson et al., 2016; Lencioni, 2018; Alther et al., 2019). Although common species are thought to be more habitat generalists and occur in higher abundances, there has been no consensus in the underlying processes that make species common or rare (Alahuhta et al., 2014). Some studies suggest that there are no strong differences in the mechanisms that make species rare or common at a regional scale (Heino & Soininen, 2010; Siqueira et al., 2012). Others propose that commonness and rarity depend on niche differentiation (Rabinowitz, Rapp & Dixon, 1984; Cucherousset et al., 2008), whereas another group of studies propose that rarity is associated with dispersal abilities and colonization dynamics among species (Kunin & Gaston, 1993; Resh, Bêche & McElravy, 2005). In the case of aquatic macroinvertebrates, and particularly the highly adapted dipteran family Chironomidae, Robinson et al. (2016) found highly diverse and spatially-structured assemblages with higher levels of rarity in mountain catchments. There is speculation that the insularity of high mountain catchments and the harsh environmental conditions associated with glacial influence may be the cause of this high degree of rarity often found in headwater streams (Lencioni, 2018).
In this study we examined local (α) and among site (β) diversity of benthic macroinvertebrate communities, as well as environmental instream conditions in 51 stream site types in a glacierised catchment of the tropical Andes. We used the same macroinvertebrate database as Cauvy-Fraunié et al. (2015) to answer the following questions: (1) are macroinvertebrate diversity, rarity, commonness and spatial distribution patterns different among site types with different water sources? (2) which environmental variables influence the density and presence of macroinvertebrate taxa, and in particular of the subfamilies of the ubiquitous chironomids? We predicted that, compared to the other stream site types, glacier-fed sites would shelter a higher proportion of rare and specialized taxa adapted to the harsh local conditions. We also expected that environmental variables associated with glacial meltwater, such as turbidity, water temperature, conductivity, and physical stability would influence the distribution and density of freshwater macroinvertebrates and of chironomids in particular. Understanding such diversity patterns could yield important information, useful for designing effective management and/or conservation strategies for these highly threatened ecosystems.
Materials and methods
Study area: The study was conducted at 51 sites located in the western slopes of Mt. Antisana, located in the eastern cordillera of the Ecuadorian Andes (0°28’S, 78°09’W, 5 760 m.a.s.l.), c. 50 km south of the equator. All sites are tributaries of the Río Antisana, headwater of the Napo River, a main tributary of the upper Amazon River. All study sites were located between 3 886 m.a.s.l. and 4 835 m.a.s.l., at distances from 15 m to 15.2 km from the glacier snouts of Mt. Antisana (Fig. 1). For additional information on the study area see Cauvy-Fraunié et al. (2015).
Following Brown et al. (2003), stream sites were first divided into three main groups, depending on their origin: kryal (glacier melt-dominated streams), krenal (groundwater-fed streams), and rhithral (seasonal rain and snowmelt-fed streams). We then applied a sub-classification - based on the percentage of glacier cover in the catchment (GCC) (Füreder, 2007; Milner et al., 2009) - to all sites with glacier melt water influence, and to those resulting from the confluence of streams with different water sources (see Cauvy-Fraunié et al., 2014a for details on GCC calculation). Kryal sites originated from three glaciers on the western side of Mt. Antisana: glacier 12, also called Crespo, with an extension of 1.82 km², in 2010; glacier 14, with an extension of 1.24 km², in 2011; and glacier 15, with an extension of 0.60 km², in 2010. We divided kryal sites into two groups: kryal 1 (K1) sites were located between 4 520-4 835 m.a.s.l. and had the highest percentages of GCC (mean = 68%, range = 39-93%). Kryal 2 (K2) sites consisted of five sites between 4 109-4 332 m.a.s.l. and had a GCC of 17 to 27% (mean = 23%). We included 18 rhithral and six krenal sites, which were located between 3 917-4 368 and 4 006-4 124 m.a.s.l., respectively, and had no glacial influence. The mixed sites (M), included 16 sites located between 3 930-4 246 m.a.s.l., with 1 to 18% GCC (mean = 8 %) (Fig. 1, Table 1).
Stream Site Types | Number of sites per site type | Stream Order | Mean Altitude (m.a.s.l.) | Mean Glacier Cover in Catchment (%) | Mean Distance from Glacier (m) | Unique Taxa per site type | Total Common Taxa |
K1 | 6 | 1 | 4 648 (4 520-4 835) | 68 (39-93) | 927 (15-1 847) | 2 (2.4%) | 13 (15.9%) |
K2 | 5 | 1, 2 | 4 206 (4 109-4 332) | 23 (17-27) | 5 573 (3 715-7 189) | 2 (2.4%) | |
M | 16 | 1, 2, 3 | 4 072.4 (3 930-4 246) | 8 (1-18) | 9 048 (3 900-15 183) | 8 (9.8%) | |
Rhithral | 18 | 1, 2, 3 | 4 080.4 (3 917-4 368) | - | - | 8 (9.8%) | |
Krenal | 6 | 1 | 4 063.4 (4 006-4 124) | - | - | 4 (4.9%) |
Environmental variables: In all sites, on the same dates of macroinvertebrate sampling (see below), water temperature (°C), conductivity (at 25 °C), pH, turbidity and the channel bottom component of the Pfankuch index (used to quantify the physical stability of the stream) were measured by Cauvy-Fraunié et al. (2015, see details therein). Autochthonous and allochthonous food resources available for macroinvertebrates were measured by quantifying chlorophyll a (including phaeopigments) concentration in epilithic algae on randomly collected pebbles from each site and by calculating the weight of benthic detritus obtained in each Surber sample after the sorting of the animals (see Cauvy-Fraunié et al., 2015). Current velocity, stream slope, depth and width were measured by Crespo-Pérez et al. (2020, see details therein). For this study, on the same dates as measurements performed by Crespo-Pérez et al. (2020), we also characterized the stream bottom at three transects at each site. We estimated the number of substrate types, according to the Wentworth Scale (Giller & Malmqvist, 1998): silt (0.004-0.062 mm), sand (0.063-2 mm), gravel (2-16 mm), pebble (16-64 mm), cobble (64-256 mm), boulder (>256 mm), and algae, macrophytes and/or moss (Supp. Table S1).
Macrobenthos sampling: Macroinvertebrate sampling was carried out at each site between May and October 2009 (i.e., during the dry season). At each site, five quantitative Surber samples (500 cm2; mesh size 200 µm) were collected randomly from pebble-cobble substratum in riffle/run habitats. Invertebrates were identified to morphospecies, genus, subfamily (in the case of Chironomidae) and family according to South American macroinvertebrate keys (Domínguez & Fernández, 2009). For further details in macroinvertebrate sampling see Cauvy-Fraunié et al. (2015).
Data analysis: To describe biodiversity patterns across all five stream site types, we calculated two types of diversity levels: α and β. Alpha diversity (local diversity) was calculated as: (1) taxa richness (number of taxa); (2) the Shannon-Wiener (H) diversity index (based on the number of taxa corrected by their abundance), providing more information about community structure (including evenness) than a simple measure of richness (Magurran, 2013); (3) rarefied richness, using the minimum abundance across all sites (gives the expected species richness in random subsamples to correct for the density effect related to sampling imbalances); and (4) true diversity index (exponential Shannon) to convert the non-linearity characteristic of H to a linearized form, thereby allowing meaningful interpretations of the effective number of taxa (Jost, 2006). Shannon-Wiener, rarefied richness and true diversity indices were calculated in R software using the “vegan” package (version 3.6.3, R Development Core Team)
We used a multivariate approach to assess the spatial variation in community composition (i.e., Beta diversity). First, we tested similarity among two or more within-group dispersion matrices (Borcard et al., 2018) with two indices: Sorensen (presence-absence data) and Bray-Curtis (raw abundance data). To further understand how community composition changed among site types, we partitioned our multiple-site overall beta diversity (calculated with the Sorensen index) into two components: taxa replacement (i.e., taxa turnover) and nestedness (i.e., species loss or gain), with the poorest assemblage being a strict subset of the richest one (Baselga & Orme, 2012). Finally, we evaluated differences in beta diversity components among site types, using analysis of variance, followed by Tukey’s HSD post-hoc tests. These analyses were performed using the “betapart” package of R software (version 3.6.3, R Development Core Team). Differences in macroinvertebrate taxa relative abundance among site types were further assessed with an analysis of similarities (ANOSIM) (Chapman & Underwood, 1999), which tested the null hypothesis that within-group similarity was equal to between-group similarity. To determine which macroinvertebrate taxa contributed the most to differentiating site types, we performed a similarity percentage analysis (SIMPER). ANOSIM and SIMPER analyses were ran using the Bray-Curtis distance measure and 9 999 permutations, with the free software PAST (PAleontological STatistics, version 3.20).
To evaluate the influence of environmental variables on community structure we performed a Redundancy Analysis (RDA) and generalized linear models (GLM). Prior to these analyses, we excluded rare taxa (those occurring in only one sample), as well as those that accounted for less than 0.3% of the total species abundance. However, taxa with abundances less than 0.3% but commonly distributed (i.e., > 7 sites) were retained. Taxa with abundances exceeding 0.3%, but occurring only in one or two sites were also excluded, in order to avoid bias caused by rare taxa (Legendre & Legendre, 1998; Cao et al., 2001), or the ‘dilution’ of taxa that are uncommon and may influence statistical analysis. We also determined “rarity in abundance” by counting taxa with only one or two individuals at each site and “rarity in distribution” by counting taxa that were only present in one or two sites. Commonness was calculated as the number of taxa that were shared between all our site types. Before running RDAs we tested if our count data were suitable for linear ordination methods by testing their homogeneity with Detrended Correspondence Analysis (DCA) (Lepš & Šmilauer, 2003). RDA computes axes that are linear combinations of the exploratory variables, in our case environmental variables, and preserves Euclidean distances among objects (Legendre & Legendre, 1998). All correlated parameters (r ≥ 0.7) were excluded from the RDA, to reduce collinearity (Supp. Table S2). We ran the RDA on log10 (x+1) transformed environmental variables and on Hellinger transformed taxa density data, to reduce the effect of taxa with low counts and many zeroes. RDA analysis was performed with the “vegan” package of R (version 3.6.3, R Development Core Team).
The effect of instream environmental conditions on community diversity patterns (richness and density) was further examined with a generalized linear modeling (GLM) approach. We included 12 environmental predictors as candidate variables, after the removal of highly correlated variables (r ≥ 0.7). Depending on the response variables, we fitted models with a different structure. Poisson models were used for non-overdispersed count data, ‘‘quasi-Poisson’’ models were used for responses with an overdispersion lower than 15, negative binomial models were used in cases of overdispersion higher than 15, and binomial models were used for presence/absence response data (Tagliaferro & Pascual, 2017). The best fitting models were chosen based on the lowest corrected Akaike’s Information Criterion (AICc) value. Quasi-likelihood methods for small sample sizes (QAICc) were used with overdispersed data (Burnham & Anderson, 2002). Note that we also ran GLMs for Chironomidae density, and its subfamilies’ density and presence-absence. Dispersion tests, negative binomial models, multiple models and multi-model inference were carried out using “AER”, “MASS”, “glmulti” and “MuMIn” packages, respectively, in R Software (version 3.6.3, R Development Core Team).
Results
Community structure and diversity distribution among stream site types: We collected a total of 90 078 specimens belonging to 82 taxa and 15 macroinvertebrate orders from the 51 sampled sites (including rare taxa). Of these, 34 taxa belonged to Diptera (10 of which were from the family Chironomidae), 13 to Trichoptera, 11 to Coleoptera, 5 to Lepidoptera, 4 to Rhynchobdellida, 3 to Haplotaxida, 2 to Hemiptera, and 1 to Lumbriculida, Tricladida, Gordioidea, Hydracarina, Basommatophora, Podocopida, Amphipoda, Veneroida, Plecoptera, and Ephemeroptera. Removing rare taxa (see data analysis section) left us with 89 767 specimens and 40 taxa, representing 99.7% of total abundance and 48.8% of total richness, respectively.
Kryal (K1 and K2) and mixed (M) site types had an important proportion of rare taxa (singletons and doubletons) (K1 = 9, K2 = 13 and M = 19, Fig. 2). These accounted for 42.9, 35.1 and 26.4% of all the taxa collected in those three site types, respectively. In contrast, rhithral and krenal stream site types, which were dominated by non-rare taxa, contributed to more than 80 % of the total number of collected individuals (Fig. 2). Rarity in distribution revealed that the site type with the highest proportion of taxa found in a single locality were K2 sites, with 22 taxa of a total of 41 (54%, Fig. 3). The other site types had lower proportions of taxa in single localities, ranging from 31% in rhithral sites to 37% in M sites (Fig. 3). In terms of uniqueness per site type, M and rhithral sites had the highest number of taxa present in only these types of sites (8 taxa each). In contrast, K1 and K2 sites harbored only two taxa unique to these types of sites (Table 1, Supp. Table S3). Commonness (i.e., shared species among all site types) was 15.9%, meaning that 13 out 82 taxa were present at all types of sites (Table 1, Supp. Table S3).
Rhithral sites harbored the highest abundance of macroinvertebrates, whereas M and K1 sites had the lowest numbers of individuals (Table 2). Mean richness was highest in the M sites and lowest in K1 sites (Table 1&Table 2). Rarefied richness (subsample of 44 individuals) showed that K2 sites were the richest (7.05) and K1 sites, the poorest (4.05, Table 2.). We found similar patterns for Shannon-Wiener (H) and true diversity (exponential Shannon) (Table 2, Fig. 4a).
Stream Site Types | Individuals | Taxa Richness | Rarefied Richness | Shannon- Wiener (H) | True Diversity (exp(H)) | Bray-Curtis’s Dissimilarity | Sorensen’s Dissimilarity | Sorensen’s Turnover | Sorensen’s Nestedness |
Kryal 1 (6) | |||||||||
Mean | 159.7 | 8.83 | 4.05 | 1.16 | 3.38 | 0.68 | 0.52 | 0.14 | 0.37 |
SD | 122.2 | 5.64 | 1.30 | 0.38 | 1.21 | 0.18 | 0.15 | 0.15 | 0.24 |
Range | 44-340 | 2-16 | 2-5.48 | 0.66-1.61 | 1.93-5.03 | 0.40-0.91 | 0.23-0.78 | 0.23-0.78 | 0.03-0.78 |
Kryal 2 (5) | |||||||||
Mean | 292.2 | 15.8 | 7.05 | 1.82 | 6.73 | 0.74 | 0.45 | 0.35 | 0.10 |
SD | 300.5 | 3.70 | 1.73 | 0.52 | 2.66 | 0.19 | 0.09 | 0.12 | 0.09 |
Range | 58-809 | 10-20 | 4.24-8.95 | 0.94-2.23 | 2.56-9.28 | 0.27-0.89 | 0.35-0.63 | 0.1-0.53 | 0.00-0.3 |
Mixed (16) | |||||||||
Mean | 1 576.7 | 29.2 | 6.68 | 1.68 | 6.64 | 0.64 | 0.29 | 0.21 | 0.08 |
SD | 925.1 | 4.90 | 1.60 | 0.58 | 2.41 | 0.19 | 0.07 | 0.08 | 0.06 |
Range | 226-2 969 | 20-36 | 2.83-8.94 | 0.58-2.39 | 1.79-10.93 | 0.06-0.94 | 0.15-0.46 | 0.04-0.45 | 0.00-0.26 |
Rhithral (18) | |||||||||
Mean | 2 833.3 | 25.7 | 5.68 | 1.54 | 5.23 | 0.70 | 0.36 | 0.27 | 0.09 |
SD | 2 238.6 | 5.65 | 1.61 | 0.53 | 2.36 | 0.18 | 0.10 | 0.10 | 0.08 |
Range | 232-10 050 | 14-37 | 2.44-8.76 | 0.46-2.36 | 1.59-10.61 | 0.24-0.99 | 0.12-0.67 | 0.05-0.55 | 0.00-0.36 |
Krenal (6) | |||||||||
Mean | 1 905.5 | 24.8 | 5.7 | 1.63 | 5.21 | 0.55 | 0.39 | 0.24 | 0.15 |
SD | 767.2 | 7.28 | 0.67 | 0.19 | 1.01 | 0.14 | 0.09 | 0.14 | 0.10 |
Range | 555-2 797 | 16-34 | 4.89-6.43 | 1.42-1.89 | 4.15-6.59 | 0.32-0.81 | 0.20-0.55 | 0.00-0.46 | 0.03-0.36 |
Beta diversity (using Sorensen and Bray-Curtis similarity indices) varied across site types, with kryal sites (K1 and K2) presenting the highest values for Sorensen (0.52 and 0.45, respectively) and K2 and rhithral, for Bray-Curtis (0.74 and 0.70, respectively, Fig. 4b). On the other hand, the nestedness component of the Sorensen index differed among site types (ANOVA, F4,46 = 6.05, p< 0.001), whereas the turnover component presented non-significant differences (ANOVA, F4,46 = 1.6, p> 0.05). Post hoc HSD Tukey test showed that K1 sites differed from all the other site types in the nestedness component, but differed only from M sites in the overall component (Table 2, Fig. 4c).
Community composition patterns: Our analyses of similarity (ANOSIM) revealed significant differences between communities of the five types of sites (R = 0.283, p < 0.01, average dissimilarity 67.9%). SIMPER analysis showed that Orthocladiinae (Chironomidae, Diptera), Hyalella sp. (Hyallelidae, Amphipoda) and Andesiops sp. (Baetidae, Ephemeroptera) contributed the most to this dissimilarity (˃ 45% of cumulative contribution). K1 sites differed significantly from all the other site types (ANOSIM: R > 0.6, p < 0.05, average dissimilarity of 84.3%), with Podonominae Type 1, Type 3 (Chironomidae, Diptera) and Orthocladiinae as the major contributors to those differences (Table 3). Hyalella sp. and Andesiops sp. were completely absent from K1 sites. Finally, K2 sites differed significantly from krenal sites (ANOSIM: R = 0.55, p < 0.05, average dissimilarity of 71.9%), with Hyalella sp., Orthocladiinae and Neoelmis sp. (Elmidae, Coleoptera) as the major contributors to those differences, and Hyalella sp. presenting low densities (0.002 Ind m-2) in K2 sites (Table 3).
Sites | SIMPER | ANOSIM | ||||
Most discriminating taxa | Contribution (%) | Cumulative percentage | Overall Average Dissimilarity | R | p | |
K1 vs. K2 | Podonominae Type 1 | 19.8 | 19.8 | 82.6 | 0.627 | 0.042 |
Podonominae Type 3 | 16.4 | 36.3 | ||||
Orthocladiinae | 12.5 | 48.7 | ||||
K1 vs. M | Podonominae Type 1 | 20.4 | 20.4 | 83 | 0.769 | 0.001 |
Orthocladiinae | 16.9 | 37.3 | ||||
Podonominae Type 3 | 16.4 | 53.7 | ||||
K1 vs. Rhithral | Podonominae Type 1 | 20.4 | 20.4 | 84.7 | 0.675 | 0.002 |
Podonominae Type 3 | 16.3 | 36.7 | ||||
Orthocladinae | 15.7 | 52.4 | ||||
K1 vs. Krenal | Hyalella sp. | 20.6 | 20.6 | 87 | 0.848 | 0.027 |
Podonominae Type 1 | 19.9 | 40.5 | ||||
Podonominae Type 3 | 15.9 | 56.5 | ||||
K2 vs. M | Orthocladiinae | 23.2 | 23.2 | 62.9 | 0.254 | 0.47 |
Neoelmis_sp. | 15.4 | 38.6 | ||||
Andesiops sp. | 10.5 | 49.1 | ||||
K2 vs. Rhithral | Orthocladiinae | 19.9 | 19.9 | 69.7 | 0.241 | 0.71 |
Neoelmis_sp. | 13.9 | 33.9 | ||||
Hyalella sp. | 13.4 | 47.3 | ||||
K2 vs. Krenal | Hyalella sp. | 24.7 | 24.7 | 71.9 | 0.547 | 0.02 |
Orthocladiinae | 14.6 | 39.3 | ||||
Neoelmis_sp. | 13.4 | 52.7 | ||||
M vs. Rhithral | Orthocladinae | 24.3 | 24.3 | 59.9 | 0.0503 | 1 |
Andesiops sp. | 16 | 40.3 | ||||
Hyalella sp. | 15.6 | 55.9 | ||||
M vs. Krenal | Hyalella sp. | 27 | 27 | 59.8 | 0.227 | 0.286 |
Orthocladiinae | 21.8 | 48.8 | ||||
Andesiops sp. | 13.1 | 61.9 | ||||
Rhithral vs. Krenal | Hyalella sp. | 24.5 | 24.4 | 56.7 | -0.098 | 1 |
Orthocladiinae | 21.3 | 45.8 | ||||
Andesiops sp. | 15 | 60.8 |
Environmental gradients and driving factors: In our RDA, the first two canonical axes explained 60.2% of the variance in the macroinvertebrate-environment interaction (axis 1: 39.82% and axis 2: 20.36%). Axis 1 was significantly defined by water temperature, CPOM, slope and turbidity, while axis 2 by current velocity, chlorophyll a and stream width. The eigenvalues of axis 1 and 2 were 0.09 and 0.05, respectively. Based on environmental vector lengths: water depth, turbidity, water temperature, current velocity and CPOM were the most important environmental factors (Fig. 5a). K1 sites were associated with higher values of GCC and turbidity and higher densities of Podonominae type 1, 3 and 4 (Fig. 5b). Water temperature, amount of CPOM and densities of Andesiops sp., Anomalocosmoecus sp. and Hyalella sp. were negatively associated with Kryal sites. Hyalella sp. and Ostracoda presented an association with higher pH and Pfankuch stability, whereas Orthocladinae was associated with higher current velocity and chlorophyll a concentration.
GLMs showed that the main variables influencing macroinvertebrate richness were CPOM, current velocity (both with a positive relationship), the Pfankuch index and GCC (negatively related) (Table 4). Turbidity and GCC were the best variables to explain macroinvertebrate density with a negative relationship (Table 4). In the case of chironomids density and presence/absence, GCC appeared in most of the best-fitted models. The increase in GCC was related to high densities of Orthocladiinae, Diamesinae, Tanypodinae and Chironominae, and to the presence/absence of Orthocladiinae and Tanypodinae. Water temperature had a negative effect on the densities of Podonominae, Diamesinae, Tanypodinae and Chironominae (Table 4).
Dependent Variable | Model | Variables in model | Df | logLik | (Q)AICc | Δ (Q)AICc | Distribution |
Richness | Best | CPOM(0.01) + V(0.52) + Pfk(-0.01) + GCC(-0.017) | 5 | -153.872 | 250.6 | 0 | Poisson |
Full | All | 13 | -151.706 | 273 | 22.43 | Poisson | |
Null | ― | 1 | -219.018 | 341.1 | 90.54 | Poisson | |
Density (Ind/m2) | Best | Trb(-0.002) + GCC(-0.03) | 4 | -481.72 | 800.8 | 0 | Bin.Neg. |
Full | All | 14 | -475.285 | 822.6 | 21.84 | Bin.Neg. | |
Null | ― | 2 | -503.299 | 831.3 | 30.54 | Bin.Neg. | |
Chironomidae Only (Ind/ m2) | Best | T(-0.23) + Sbs(-0.16) + Pfk(-0.04) + GCC(-0.05) | 6 | -446.591 | 820.8 | 0 | Bin.Neg. |
Full | All | 14 | -459.162 | 840.8 | 19.97 | Bin.Neg. | |
Null | ― | 2 | -442.624 | 833.3 | 12.54 | Bin.Neg. | |
Podonominae Density (Ind/ m2) | Best | Chlo(0.02) + CPOM(-0.07) + T(-0.2) + Slp(-0.01) + Wdt(-0.01) + Dpt(-0.04) + Sbs(-0.14) | 9 | -272.875 | 568.1 | 0 | Bin.Neg. |
Full | All | 14 | -271.931 | 583.5 | 15.39 | Bin.Neg. | |
Null | ― | 2 | -287.73 | 579.7 | 11.57 | Bin.Neg. | |
Podonominae Presence/Absence | Best | T(-0.52)+Wdt(-0.005) | 3 | -16.311 | 39.1 | 0 | Binomial |
Full | All | 13 | -14.709 | 65.3 | 26.12 | Binomial | |
Null | ― | 1 | -20.397 | 42.9 | 3.74 | Binomial | |
Orthocladiinae Density (Ind/ m2) | Best | Wdt(0.003) + Pfk(-0.06) + GCC(-0.06) | 5 | -427.299 | 781.8 | 0 | Bin.Neg. |
Full | All | 14 | -424.871 | 807.2 | 25.44 | Bin.Neg. | |
Null | ― | 2 | -443.605 | 803.7 | 21.9 | Bin.Neg. | |
Orthocladiinae Presence/Absence | Best | Chlo(2.4) + Trb(-0.08) + GCC(0.22) | 4 | 0 | 8.9 | 0 | Binomial |
Full | All | 13 | -4.922 | 35.8 | 26.97 | Binomial | |
Null | ― | 1 | 0 | 11.9 | 3.06 | Binomial | |
Diamesinae Density (Ind/ m2) | Best | Chlo(0.02) + pH(1.1) + T(-0.48) + Pfk(-0.11) + GCC(-0.06) | 7 | -208.214 | 433 | 0 | Bin.Neg. |
Full | All | 14 | -207.266 | 454.2 | 21.16 | Bin.Neg. | |
Null | ― | 2 | -219.352 | 443 | 9.92 | Bin.Neg. | |
Diamesinae Presence/Absence | Best | CPOM(0.08) + Trb(-0.004) + Slp(-0.01) | 4 | -27.828 | 64.5 | 0 | Binomial |
Full | All | 13 | -25.899 | 87.6 | 23.11 | Binomial | |
Null | ― | 1 | -33.112 | 68.3 | 3.78 | Binomial | |
Tanypodinae Density (Ind/ m2) | Best | CPOM(0.1) + pH(2.1) + T(-0.26) + Pfk(-0.14) + GCC(-0.26) | 7 | -161.772 | 340.1 | 0 | Bin.Neg. |
Full | All | 14 | -160.791 | 361.2 | 21.04 | Bin.Neg. | |
Null | ― | 2 | -184.12 | 372.5 | 32.34 | Bin.Neg. | |
Tanypodinae Presence/Absence | Best | Chlo(0.03) + pH(2.1) + Slp(0.02) + Dpt(0.08) + GCC(-0.09) | 6 | -20.671 | 55.3 | 0 | Binomial |
Full | All | 13 | -16.819 | 69.5 | 14.23 | Binomial | |
Null | ― | 1 | -35.105 | 72.3 | 17.04 | Binomial | |
Chironominae Density (Ind/ m2) | Best | CPOM(0.09) + T(-0.44) + Slp(-0.01) + Pfk(-0.06) + GCC(-0.15) | 7 | -198.306 | 413.2 | 0 | Bin.Neg. |
Full | All | 14 | -196.618 | 432.9 | 19.68 | Bin.Neg. | |
Null | ― | 2 | -219.875 | 444 | 30.78 | Bin.Neg. | |
Chironominae Presence/Absence | Best | Trb(-0.01) + Wdt(-0.01) + Dpt (0.26) | 4 | -13.016 | 34.9 | 0 | Binomial |
Full | All | 13 | -9.334 | 54.5 | 19.6 | Binomial | |
Null | ― | 1 | -30.896 | 63.9 | 28.97 | Binomial |
Discussion
Previous studies from the Antisana area have focused on longitudinal, spatial and temporal patterns of taxa composition and functional diversity, on community and metacommunity structure and on experimental ecological responses to environmental drivers, physical forces, ecosystemfunctioning and climate change induced reduction in water flow (Jacobsen et al., 2010; Kuhn et al., 2011; Fugère et al., 2012; Cauvy-Fraunié et al., 2014b; Cauvy-Fraunié et al., 2015; Cauvy-Fraunié et al., 2016; Tiegs et al., 2019; Crespo-Pérez et al., 2020). The novelty of the present study lies on the study of α and β diversity (local diversity and taxa turnover, respectively), on the sub-classification of stream site types (i.e., kryal 1, kryal 2, mixed, krenal and rhithral) and on the analysis of rare and common taxa between these site types. Our a priori grouping of sites was confirmed to be acceptable according to the results obtained in the ANOSIM and RDA analysis, and supports the hypothesis that distinct benthic macroinvertebrate assemblages would inhabit streams of different origins (Brown, Milner & Hannah, 2006). Our study also revealed how diversity patterns (α and β) varied along environmental gradients (i.e., GCC, temperature, turbidity, etc.) and among different site types at the catchment scale.
Commonness, rarity and diversity gradients: Understanding how the distribution of rare and common taxa contributes to the emergence of biodiversity patterns is a key concern in conservation biology (Lennon et al., 2004; Pearman & Weber, 2007). In this study, we found differences in the composition of aquatic communities between sites with different water sources. We found unique taxa in all the site types (i.e., taxa that were present only in a particular site type), with rhithral and M sites having 12 and 9 singletons, respectively, which corresponded to 13 and 10% of the total richness found in the study (Fig. 2). In kryal sites (2, 8 and 9) we found two unique taxa, Dytiscidae and Staphylinidae type 1 (Coleoptera), that were not found in any of the other site types. This could represent only chance encounters, or could reveal recent colonization of these harsh environments by these predatory taxa. This result could also be related to the reduced number of K1, K2 and krenal sites, compared to mixed and rhithral sites. More data, especially from kryal and krenal sites, are needed to confirm this finding. However, some Staphylinidae have indeed been found to be adapted for life in subaquatic conditions near the stream channel (Domínguez & Fernández, 2009; Lancaster & Downes, 2013). Podonominae type 4 (Diptera) were present in all the sites that were influenced by glacial melt water (K1, K2 and M), but not in rhithral or krenal sites, which could mean that this chironomid is highly adapted to conditions of glacier-fed streams (Cauvy-Fraunié et al., 2015). On the other hand, there were 13 taxa that were present in all five stream site types [e.g., Atopsyche sp. (Trichoptera), Cailloma sp. (Trichoptera), Simulium sp. (Diptera), Hemerodromia sp. (Diptera), Molophilus sp. (Diptera), Podonominae type 3 (Diptera), etc.], but more analyses (e.g., genetic barcodes) are necessary to determine if each taxon corresponds only to one or to several species. If the former were true this would mean that these species have plastic habitat requirements and/or high dispersal capacity (all of them have an adult flying stage), and adaptations (e.g., morphological and physiological) for life under all these high altitude conditions. As observed by Alther et al. (2019) in the Swiss Alps, high altitude glacierised catchments harbor a surplus of generalist species, with a few rare species, compared to other ecosystems.
As found in previous studies (Jacobsen et al., 2012; Cauvy-Fraunié et al., 2014a), our results revealed a peak in richness at intermediate levels of glacier influence. As stated by these authors, this could be explained by the intermediate disturbance hypothesis (Connell, 1978). Indeed, intermediate disturbance at mixed sites could favor the presence of disturbance tolerant taxa, while reducing the abundance of competitive, disturbance averse taxa (Crespo-Pérez et al., 2020). Another possible explanation to the higher diversity of mixed sites, could be that these sites receive more individual colonists, both from upstream reaches (e.g., from kryal and rhithral sites) and from lateral movements between branches (Crespo-Pérez et al., 2020, Brown and Swan, 2010). Our results based on Sorensen’s measures of dissimilarity showed that nestedness (sites with lower richness are subsets of sites with higher richness) (Ulrich & Gotelli, 2007) better explains the difference between our communities. Sites with stronger glacier influence (K1) differed from all the other types of sites according to RDA, ANOSIM and Tukey’s HSD post hoc tests (Fig. 4c), which could reflect a non-random process of species loss across these site types that promotes these disaggregated assemblages (Gaston & Blackburn, 2000). As expected, and in accordance to previous studies, Sorensen’s overall dissimilarity increased with increasing glacial influence (which is significantly correlated with altitude, Supp. Table S2) (Finn et al, 2013; Castro et al., 2019). This process could be explained by an increment of environmental heterogeneity between sites, which in turn, increases regional diversity (Cauvy-Fraunié et al., 2015; Crespo-Pérez et al., 2020).
Community composition, environmental characteristics and driving factors: Hyalella sp. (Amphipoda) and Andesiops sp. (Ephemeroptera) were distributed in all our site types except for the K1 sites. Anomalocosmoecus sp. (Trichoptera) was present in all sites but not in all kryal sites. Hyalellidae amphipods have been found to present some degree of trophic plasticity (Wantzen & Wagner, 2006) and successful adaptation to a variety of environmental conditions (Acosta & Prat, 2011), but they were found to have high mortality rates when transplanted to high altitude streams (c.a. 4 500 m.a.s.l.) (Madsen et al., 2015). On the other hand, Andesiops ephemeropterans, are commonly found and are relatively abundant in hig altitude ecosystems of Ecuador, and as proposed by Finn et al. (2016), this genus presents recent population growth, in evolutionary terms, at the upstream catchment scale. These populations of Andesiops sp. and Anomalocosmoecus sp. proved to be very sensitive to temperature in transfer experiments, with high mortality rates above 4 200 m.a.s.l. (Madsen et al., 2015). Therefore, these three taxa seem to be competitive and dominant under relatively benign conditions, but are excluded by the extreme environmental conditions of kryal sites.
Macroinvertebrate communities seemed to be affected primarily by environmental factors such as turbidity, water temperature, CPOM, current velocity and chlorophyll a. Certain macroinvertebrate taxa (e.g., Diamesinae, Baetidae, Leptophlebiidae and Podonominae) have been shown by previous studies to have greater tolerance to glacial influence than other less tolerant taxa (Milner et al., 2001b, Ilg & Castella, 2006; Milner et al., 2009). GCC influences the stream thermal regime (Hood & Berner, 2009), with temperature increasing farther downstream as glacial influence dampens; this temperature gradient is expected to decline with warming climates and glacial recession (Slemmons, Saros & Simon, 2013). Channel stability and water temperature are key physicochemical drivers influencing longitudinal trends in macroinvertebrate assemblages along proglacial rivers (Milner et al., 2001a). In our study, channel stability seemed to be an antagonistic force to current velocity and enhanced the density of certain taxa like Ostracoda and Hyalella sp. Several studies have found current velocity to be an important driver of stream macroinvertebrate assemblages (e.g., Hieber et al., 2005; Cauvy-Fraunié et al., 2014a; b), and in our case, positively affected the density of certain rheophilic macroinvertebrates such as Blepharicera sp. (Diptera) (Fig. 5b).
Macroinvertebrate richness was positively influenced by CPOM and stream current velocity. CPOM has been previously linked to the number of taxa occurring in headwater streams (Fenoglio et al., 2005; Straka, Syrovátka & Helešic, 2012); although, a previous study limited to nine glacier-fed sites in the same region, did not find a relationship between the fauna and CPOM levels (Kuhn et al., 2011). On the other hand, we found that macroinvertebrate density was negatively related to turbidity and GCC, with both environmental factors being strongly correlated with water temperature and conductivity (Table 4).
According to previous studies (e.g., Jacobsen et al., 2010; 2014; Kuhn et al., 2011) chironomids are numerically dominant in our study region. In temperate European streams, as one gets closer to the glacier, Diamesa sp. (Diamesinae, Chironomidae) becomes the only genera present (Milner et al., 2001a; Clitherow, Carrivick & Brown, 2013); a pattern observed with Podonominae in Ecuadorian tropical glacierised catchments (Jacobsen et al., 2010). More specifically, we found that Podonominae type 1 and 3 were abundant at kryal sites, but Podonominae type 4 was the only taxon found exclusively at the upper-most sites. The latter taxon probably belongs to the genus Parochlus, but further taxonomic identification is mandatory in order to confirm its endemic status in these Ecuadorian Andes streams. Diamesa species found in a harsh glacier-fed stream in the Central Austrian Alps were found to have an omnivorous diet rich in detritus and diatoms (Clitherow, Carrivick & Brown, 2013). Likewise, according to our GLM analysis, the relative abundance of the subfamily Podonominae was positively influenced by CPOM quantity and chlorophyll a concentration (but note that CPOM and chlorophyll a were negatively related to GCC), and their presence was negatively influenced by water temperature and stream width (Table 4, Supp. Table S2). All this suggests that the ecological niche necessary for the presence and abundance of certain chironomid species in harsh glacier conditions are similar for temperate and tropical glacier-fed streams.
In conclusion, given the relatively little knowledge about mountaintop invertebrates, and the increased threat posed by climate change, more research and monitoring is urgently needed to predict and adapt to the effects of taxa extirpation on ecosystem integrity and functioning (Muhlfeld et al., 2011). There is an urgent need to expand the spatial range of tropical glacierised catchment studies, to apply state-of-the-art research techniques (e.g., next generation sequencing, stable isotope analyses, species distribution modeling) and to expand research topics (e.g., functional diversity, genetic diversity) in glacier-fed catchments in order to discover new endemic species, unveil patterns of genetic structuring and key functional and ecological roles in these highly fragile and understudied tropical ecosystems. Finally, as proposed by Jacobsen et al. (2012), the potential threats to biodiversity suggest that strategic conservation should, from now on, take a holistic approach that includes both invertebrate and vertebrate aquatic species. Detailed analyses of benthic communities and the development of databases are key for conservation strategies. In general, water management municipalities and/or enterprises tend to focus on mountain streams’ flows, but this study highlights that water managers should also consider water quality and stream site types for more sustainable management of these important ecosystems.
Ethical statement: authors declare that they all agree with this publication and made significant contributions; that there is no conflict of interest of any kind; and that we followed all pertinent ethical and legal procedures and requirements. All financial sources are fully and clearly stated in the acknowledgements section. A signed document has been filed in the journal archives.