Central Andean Highlands are a vast territory shared by Argentina, Bolivia and Chile, with an altitude above the 3500 m asl. They include two ecoregions, namely Puna (an upraised flatland) and High Andes (with steepest slopes). In Argentina, these ecological regions span over the western of Jujuy, Salta, Tucumán, Catamarca, La Rioja, San Juan and Mendoza provinces. The area is characterized by extreme climatic conditions, like low temperatures, very high insolation, strong winds, little rainfall (mostly restricted to the summer season), and wide daily thermal amplitude around 30°C. Water is the main limiting ecological factor and thus wetlands are key functional units supporting a high proportion of the biodiversity.
The Andean wetlands include shallow and relatively small lakes and ponds, saline flats (known as salares), and peatbogs (Maldonado et al., 2011). These aquatic systems can be permanent or temporary with different degrees of salinity. Particularly, peatbogs are placed in topographic depressions characterized by high soil moisture, low rates of decomposition and low temperatures (Izquierdo et al., 2015a). The combined action of these factors promotes their peculiar vegetation and soils that contrast with the surrounding landscape, and make them the main sources of water and biomass to sustain the biota as well as human populations in the region (Izquierdo et al., 2015a).
The high biodiversity includes emblematic native camelids, flamingos, and a rich variety of endemic plants (Morello et al., 2012). Although less known than other faunal elements, the aquatic macroinvertebrates also include unique species, even genera not found elsewhere. This aquatic fauna plays a key role in the ecosystem either recycling nutrients and processing organic matter or providing food for many animals such as fishes, amphibians (Giorgi and Tiraboschi, 1999; Domínguez et al., 2006) and birds like flamingos (Tobar et al., 2014). Among the few regional-scale efforts to assess patterns of macroinvertebrates diversity in the Andean highlands, there are the works of Roldán-Pérez (1988) and Zúñiga et al. (2004) in Colombia, Jacobsen (2003, 2008a) and Jacobsen et al. (1997, 2003) in Ecuador, Maldonado et al. (2007) and Molina et al. (2008) in Bolivia. In Argentina, those efforts are of narrow taxonomic scope and refer to certain groups present in mid to low elevation Andean streams such as studies about the spatio-temporal distribution of larvae of Chironomidae performed by Scheibler et al. (2008) in Mendoza and Tejerina and Malizia (2012) in Tucumán. Recently, Nieto et al. (2016) described the community of seven peatbogs in the Andean region of Salta.
Beta-diversity broadly reflects the differences in community composition between sites. The term was introduced by R.H. Whittaker and defined as “the extent of change in community composition, or degree of community differentiation, in relation to a complex-gradient of environment, or a pattern of environments” (Whittaker, 1960). The beta diversity does not only account for the relationship between local and regional diversity, but also informs about the degree of differentiation among biological communities (Baselga, 2010). It may reflect two different phenomena: nestedness and spatial turnover derived from processes of loss and replacement of species (Harrison et al., 1992; Baselga et al., 2007). The knowledge of patterns of beta diversity is a major step towards understanding how and why biodiversity is distributed in the way it is (Baselga, 2010).
This poorly surveyed extreme environment is thought to harbor a homogeneous pool of organisms. In this context, the study of beta diversity results appropriate to answer questions around this topic. The beta diversity accounts for the compositional dissimilarity between communities. The main purpose of the present work is to analyze and to predict spatial patterns of beta diversity associated with assemblages of aquatic organisms occurring at high altitude. To do so, we used generalized dissimilarity models (GDM). First proposed by Ferrier et al. (2007), GDM relies on a statistical technique to analyze and predict spatial patterns of beta diversity that relates biological distance to ecological distance. This method builds a single linear-predictor axis, and delivers a plot showing the curvilinear relationship between ecological distance and observed compositional dissimilarity (Ferrier et al., 2007). This approach is an extension of matrix regression and is specifically designed to accommodate the nonlinearity commonly found in largescale ecological datasets. The GDM method has been used for terrestrial ecological community modelling and biodiversity assessment (Ferrier et al., 2007), conservation planning (Ferrier, 2002; Marsh et al., 2010), regional scale survey design (Ashcroft et al., 2010; Leaper et al., 2011; Laidlaw et al., 2015), and river (Leathwick et al., 2011) and marine environmental classifications (Lasram et al., 2015).
More specifically, this study aims at i) measuring beta diversity through pairwise comparison of abundance profiles but contextualized to the overall table; ii) identifying patterns of beta diversity associated with macroinvertebrate communities from peatbogs in Northwestern Argentina; iii) searching the environmental factors that best predict the composition dissimilarity; iv) classifying the sampling sites into groups based on the ecological distance and characterize them biologically.
We sampled 13 peatbogs (Fig. 1) located in the mountains of the Puna region of the Salta province, Argentina (Punean and High Andean ecoregions) between 3780 and 4750 m asl. The small streams draining this area converge in endorreic saline depressions forming lagoons or peatbogs. The water of this region comes from rainfall, aquifer and snowmelts (Ruthsatz, 2012; Alzérreca et al., 2001). Tocomar (P1), Olacapato Chico (P9) and Catua Abajo (P11) peatbogs belong to the sub-basin Cauchari-Olaroz. Chorrillos (P2), Incachule (P3) and Acay Norte (P13) belong to the sub-basin Salinas Grandes. Aguas Calientes (P4), Santa Rosa Pastos Grandes (P5), Salinas Pastos Grandes (P6) and Condor Huasi (P7) belong to the subbasin Pastos Grandes. Quirón (P8) belongs to the subbasin Pocitos/Quirón; Catua Arriba (P10) belongs to the sub-basin Rincón. Finally, Acay Sur (P12) belongs to the sub-basin Calchaquí, the only one draining into a non-endorreic watershed (i.e., Paraná river).
The region is cold and arid being rainfalls mainly concentrated around the months of austral summer (November to April). According to the set of global climate layers WorldClim, annual average rain is 100 mm, with a minimum of 91 mm (P7) and a maximum of 111 mm (P3). The annual average temperature is 5.8°C, with a minimum of 4.1°C (P3) and a maximum of 7.5°C (P6). Three field trips were carried out during the rainy season (Tab. 1). Peatbogs 1 and 2 were sampled in December 2013, peatbogs 3-7 in February 2014 and peatbogs 8-13 in April 2015.
Data collection and sampling
In each site, different samples of macroinvertebrates and water were taken, depending on the size and the heterogeneity of each site. All samples were taken from 10:00 AM until 05:00 PM. This time frame did not affect significantly either the physical and chemical characteristics of the water nor the macroinvertebrates fauna, as reported Wasson et al. (1989) from Bolivia and Jacobsen (2008a) from Ecuador.
To obtain a quantitative and comparable measurement of the abundance of macroinvertebrates, the sampling effort was standardized in time. The time window spanned 30 min and biological samples were obtained through a D-frame net (mesh size of 300 μm). Depending on the size and environmental heterogeneity of each site, the number of sampling units ranged from two to six (in total 40). Within a single site, sampling units were taken as far as possible aiming at maximize local differences in hydraulic habitat and vegetational structure from the immediate surroundings. For that reason, sampling units portray valuable information about the within-site variation and we decide to treat them separately in the subsequent analyses. The material was conserved in vials of 15 ml with 96% ethyl alcohol. In the laboratory, each sample was identified and quantified with the aid of a stereomicroscope. Regional taxonomic keys for benthic macroinvertebrates were used (Domínguez and Fernández, 2009) and identification was performed to the lowest possible level. The material examined is housed at the IBN (Instituto de Biodiversidad Neotropical, CONICET-UNT), Tucumán, Argentina.
In each sampling point, altitude was recorded and different physical and chemical characteristics of water measured with Horiba Multiparameter Water Quality Meters U-52. The variables related to the water column were: water temperature (°C), oxidation reduction potential (ORP mV), conductivity (μS/cm), turbidity (NTU), dissolved oxygen (%). Geographical distance among sites was also calculated for further analysis.
Measurement of beta diversity
The difference in the composition of assemblages must be calculated before running GDM. A coefficient based on the ordering properties of frequency values from the perspective of both sites and taxa was used to accomplish this task. The logic behind its calculation can be stated as follows: whenever two communities are closer along a sorted array induced by taxon occurrences, the respective compositional dissimilarity should decrease. Unlike other metrics, the novel proposal is erected around the notion of a double reading of the samples-by-taxa or community matrix, comprising thus a horizontal or rowwise and vertical or column-wise reading. It is aimed to capture the resemblance between samples but i) considering the regional context from where they have been drawn, and ii) comparing the structural relevance of consisting taxa. Although Euclidean distance and Bray-Curtis (Bray and Curtis, 1957) index are two dissimilarity indexes commonly used in this kind of analysis (Nollet and De Gelder, 2014), the coefficient employed here is preferred over them because double absences are taken into account (a matter of concern when communities under comparison have been sufficiently sampled and involved organisms can be considered a subset of a common regional pool). Specifically, calculation of Bray-Curtis is based on absolute differences and in the sums of the abundances of each i-species in the two samples. It is strongly influenced by the dominant species and very little influenced by rare species (Valentin, 2012). Since this index disregards the influence of multiple matching in the occurrence of rare species and it strongly depends on the disparities in size of count data, a loss of information occurs in comparing the structure of communities.
Following, we describe the measure for computing compositional dissimilarity. This is sensitive to the location of sampling units in an ordered vector after sorting them by incidences of taxa. In doing so, all elements in the lists of taxa incidences, including joint absences, contribute to the compositional dissimilarity between two assemblages. The measure depends upon the context of a given dataset and is based on the ordering properties of values of abundance across rows and columns of the samples- by- taxa matrix. A formal definition is given below.
Let A = [aij] being the table accounting for the abundance aij of species j in sampling site i.
Let R = [rij] being the table accounting for the relative incidence rij of species j in site i. In other words, relative incidences are obtained within rows (sites). Compositional similarity between two assemblages (CS), measured in the context of a given data set, is obtained by computing the following Iverson brackets applied on both raw abundance and relative incidence profiles:
where s is the number of sampling sites and n the total number of sampled taxa. The Iverson bracket is a notation that denotes a number that is 1 if the condition in square brackets is satisfied and 0 otherwise. More exactly, [P] = 1 if the statement P is true and [P] = 0 if P is false. In the equation above, logical expressions ask whether focal value aij (or equivalently rij) are outside of the range defined by ayj and azj (or equivalently ryj and rzj). A normalizing factor has been introduced to the left of summation so that the output is scaled between 0 and 1. Given a dataset, the theoretical minimum between two community samples is achieved whenever they both are in the opposite extremes of abundance and relative incidence across sampled taxa. On the contrary, the theoretical maximum occurs whenever the profiles of two samples allow defining abundance and relative incidence ranges differentiated from the rest. The R script in the Supplementary Material calculates the beta diversity of assemblages, implicit in eq. 2, using as input the respective sample-by-taxa matrix.
Generalized dissimilarity modeling
GDM (Ferrier et al., 2007), a nonlinear statistical regression technique, is used here to reveal patterns of beta diversity and to predict the compositional dissimilarity (quantified with the previously described measure) between sampled peat bogs (macroinvertebrates are the descriptor variables). This approach aims to estimate the magnitude and rate of taxa turnover along environmental/ geographical gradients, representing the dissimilarity between pairs of sampling units as a function of environmental differences and geographical distance (Fitzpatrick et al., 2013). Inputs are site-by-taxa and site-by-predictors matrices. The advantages of the GDM over classical linear matrix regression include: i) the curvilinear relationship between compositional dissimilarity and increasing environmental/ geographical distance between sites, and ii) the flexibility in the rate of compositional turnover (non-stationarity) at different positions along a given gradient (Ferrier, 2002). To deal with non-stationarity, GDM uses I-splines to transform each of the predictor variables and provide the best supported relationship between intersite environmental/geographical separation and compositional dissimilarity. I-splines are appropriate for fitting positively monotonic functions of flexible shape (Ramsay, 1988). GDM can be implemented via the toolkit package gdm (Manion et al., 2016) for R software. This package currently provides basic functions to fit, plot, predict and summarize generalized dissimilarity models. The default option of three I-spline basis functions per predictor was used. Spatial autocorrelation was accounted for by including the geographic distance between pairs of sites as a predictor variable.
The I-spline associated to each variable describes the relationship between beta diversity and that gradient. The total amount of compositional turnover associated with each variable, holding all other variables constant, can be inferred from the maximum height of the I-spline associated to it and is quantified through summation of its coefficients (each spline has three coefficients) (Fitzpatrick et al., 2013). So, the I-splines are partial regression fits that are indicative of the contribution of each variable. Additionally, changes in the rate of species turnover can be tracked by observing the first derivative of I-spline at any point along the respective gradient. Importantly, the difference in height between any two sites along the I-spline corresponds to the modeled contribution of that predictor variable to the total ecological distance between those sites (Ferrier et al., 2007; Fitzpatrick et al., 2013). Overall model significance was tested by using Monte Carlo permutation on the input matrix of predictors followed by calculation of per cent deviance explained. To evaluate the unique contributions of variables in explaining beta diversity, the respective GDMs were run after removing predictors one by one. The explained deviance of the full GDM was compared with each of the partial models to evaluate the imprint of each variable in the difference between benthic communities.
Hierarchical clustering of sites
The previously fitted model generates a matrix of predicted ecological distances between the sampling units. A hierarchical cluster analysis was performed over such matrix using the complete linkage method (Sneath and Sokal, 1973). For this procedure, the distance between the merged pair and the others correspond to the maximum of pairwise distances. Because of this, all units clustered below a particular level of distance will have inter-sample distance less than that level. At a later instance, the final dendrogram was cut at two different heights in order to extract the most inclusive and nested groups in the hierarchical arrangement of samples. All statistical analyses and graphics were done with R (R Core Team, 2015).
An IndVal analysis (Dufrêne and Legendre, 1997) was carried out to biologically characterize the different clusters of sampling sites classified by their distance along the composite ecological gradient. This analysis identified the taxa with the highest association values to the groups of samples and depends upon two factors: i) fidelity or relative frequency (fraction of sites occupied by the target taxon) and ii) exclusivity or relative average abundance of taxa within a particular group. The calculations were performed with the R package labdsv (Roberts and Roberts, 2015) which also provides the statistical significance of the values through a permutation test. A given taxon was selected as indicator whenever its probability of obtaining an indicator value as high as observed over 1000 iterations was lower than 0.05.
A total of 4117 aquatic macroinvertebrates were collected. 28 taxa were identified to the lowest feasible taxonomic level: hexapods (11 families, 3 subfamilies, 5 genera and 3 species), arachnids (1 subclass), crustaceans (1 genus), mollusks (2 classes), and annelids (1 class and 1 family). With 948 specimens, Hyallela was the dominant taxon across the great majority of peatbogs. The next taxa in abundance were Austrelmis (adults and larvae) with 630 specimens, Simuliidae with 583 specimens, Andesiops peruvianus with 527 specimens, and Protallagma titicacae with 324 specimens. All these taxa occurred in more than half of the studied peatbogs. Importantly, Glossiphonidae was detected in 24 out of 40 sampling points, even if it is usually infrequent in studies conducted in other ecoregions (Mesa et al., 2013). Insects were the dominant item in the composition across the totality of studied communities.
Peatbog 4 (Aguas Calientes) was the richest one with 21 taxa recorded, followed by peatbog 5 (Santa Rosa de Pastos Grandes) with 16 taxa and peatbogs 1 (Tocomar) and 6 (Salinas Pastos Grandes) with 15 taxa recognized. Peatbog 3 (Incachule) was the poorest one with only 6 taxa reported. The daytime temperature of the water fluctuated between 4°C (P12) and 25°C (P4). The conductivity achieved a minimum of 30 μS cm–1 (P1) and a maximum of 1300 μS cm–1 (P10). The total dissolved solids (TDS) varied from 0.08 mg L–1 (P2) to 0.8 mg L–1 (P10). Most of the peatbogs had pH neutral or alkaline with a maximum value of 8 (P10), but in P1 an acidic pH (= 3.5) was recorded. Turbidity ranged from 1 (P1 and P4) to 70 (P11); redox potential (ORP) from 40 (P3) to 190 (P9); and the dissolved oxygen from 60% (P12) to 145% (P1).
The full GDM, including both environmental variables and geographical distance, explained around 33% of the deviance in observed benthos dissimilarities (Tab. 2). This percentage was always significantly higher than those calculated on environmental tables with rows randomly permutated. When the deviance explained from GDM was partitioned, the unique contributions of altitude (17%) and conductivity (7%) were the most relevant. They also were the most important predictors since they gave the highest scores for the magnitude of effects (measured by summing the coefficients of the I-splines from GDM). Thus, altitude and conductivity contributed jointly to 50% of the composite ecological gradient. The relationship between observed compositional dissimilarity of each pair of samples and the linear predictor of the regression equation from GDM (predicted ecological distance between sample pairs) was increasingly monotonic as expected (Fig. 2). The shape of all the fitted functions (GDM-fitted I-splines) for variables associated with benthos beta diversity was consistently nonlinear suggesting a complex behavior for the rate of compositional turnover along the gradient (Fig. 3). Some predictors reached a plateau when a given threshold is surpassed, namely ORP, temperature and geographical distance, meanwhile the remaining predictors showed a non-asymptotic behavior. The peatbogs are sufficiently separated each other, as discrete units throughout the landscape, beyond the distance of 20 km, explaining thus the observed plateau in the geographical distance I-spline. The concavity of the turbidity curve means that the turnover comes slowly in clear water, but dissimilarity increase rapidly as water becomes more turbid. Contrarily, the concavity of conductivity curve means that compositional dissimilarity increases quickly at low values of conductivity (hyperfresh condition), but it increases more slowly as mineralization occurs. The distinctive cubic polynomial shape of DO I-spline clearly reflects the similar imprint of sub-saturation as well as hyper-saturation in the structure of benthos. The altitude curve increases abruptly around 4000 m asl, suggesting that elevation changes make a deep effect on community structure when they operate at higher altitudes.
A dendrogram (Fig. 4) was constructed from the predicted ecological distance between sampling sites. The hierarchical arrangement showed two main clusters, the first group included peatbogs 1-3, 12-13 (hereinafter upland), and the second group the remaining ones (hereinafter lowland). These two groups can be differentiated from each other by altitude alone, being the height 4200 m asl a hard threshold of classification (Fig. 5). These main clusters can be further split into two subgroups and distinguished between them by conductivity. The critical common threshold was around the very low value of 135 μS cm–1 in terms of salinity (Fig. 5), suggesting a distinction between demineralized and slightly more mineralized waters instead. Importantly, the same cutoff has been established as delimiter of a hyperfresh condition in which no endogenic minerals precipitate in lacustrine sediments from the Tibetan Plateau (Wang et al., 2012). Consequently, the resulting groups will be referred as Fresh-lowland, Freshupland, Hyperfresh-lowland and Hyperfresh-upland.
The IndVal recovered the indicator taxa for each of the previous clusters (Tab. 3). From these results, an account about community composition turnover associated with beta diversity can be performed. For lowland peatbogs, the list includes Austrelmis larvae and adults, Protallagma titicacae, Simuliidae and Massartellopsis. The upland peatbogs were characterized by Podonominae, Chironominae, Oligochaeta and Limnephilidae (Antarctoecia). Focusing on the pool of lowland peatbogs, the dipteran families Dolichopodidae and Stratyomidae can distinguish the hyperfresh from fresh habitats. In the pool of upland samples Dolichopodidae was again restricted to hyperfresh sites meanwhile Trichocorixa and Podonominae were distinctively present in the more mineralized subset. Hyalella, the most conspicuous taxon across the study area, was less frequent in the Hyperfresh-upland cluster.
In this work, the main abiotic factors behind the patterns of beta diversity in peatbogs from Central Andean Highlands were identified via statistical modeling. Beta diversity was measured through a coefficient sensitive to the rank order of taxa in the composition of communities sampled from a common regional pool of organisms. Among other aspects, absences can be valuable to assess community dissimilarities however the Bray Curtis index does not account for them. Therefore, in the context of our study, Bray-Curtis can be a misleading choice. The relevance of guidelines (coincidence, complementarity, relative invariance, independence of joint absence, localization and dependence on totals) are certainly contextdependent and should not be accepted by any ecologists at all circumstances (Clarke et al., 2006). With the index here presented the influence of a few dominant items and disparities in count data is down-weighed. To assess quantitatively the compositional similarity between a pair of samples, the novel procedure calculates first their propinquity in an overall array of samples ordered by taxon occurrence frequency (both relative and absolute frequency in a row-wise and column-wise reading of the sample-bytaxa matrix, respectively), and averages then such propinquity across all taxa of the regional pool. We know there are available in the current literature better alternatives for the Bray-Curtis index such as the adjusted indices introduced by Chao et al. (2006), but they still rely on the logic calculations around classic Jaccard and Sorensen indices already proved mathematically to be misleading (Dos Santos and Deutsch, 2010).
Beta diversity is defined as compositional dissimilarity between communities and it is the result of different processes occurring at multiple scales. The analysis and the prediction of spatial patterns are crucial to understand the local communities. Altitude and conductivity turned out the main drivers of beta diversity for the benthos of peatbogs from the Puna region in Salta province, Argentina. Critical values of altitude and conductivity have been set to 4200 m asl and 135 μS cm–1, respectively. The altitude is considered the main critical variable for the biota in the rivers (Vannote et al., 1980; Rahbek, 1995). The altitude is correlated with atmospheric pressure, oxygen concentration, riparian vegetation, availability of hydraulic habitats along the streambed and concentration of organic matter (Vinson and Hawkins, 1998; Jacobsen, 2003; Jacobsen et al., 2003). Altitude shapes the diversity of plants and animals (Rahbek, 1995, 2005; McCoy, 1990; Brehm et al., 2007); as a general rule, diversity decreases with increasing altitude. However, Jacobsen (2008a, 2008b) sustained that this pattern is much less pronounced when aquatic macroinvertebrates are considered. In a previous study (Nieto et al., 2016), we found the same conclusion, although we presumed that it is because of the short range considered therein (700 m). Here, the communities are clearly affected by the altitude; some taxa frequently present in the ‘Lowland’ disappeared in the ‘Upland’, and vice versa.
The altitudinal limit of 4200 m asl could correspond with the limit between Punean and High Andean ecoregions, which is consistent with Burkart et al. (1999), who proposed this limit to be between 4300-4500 m. However, altitude alone is not the best and single classifier to distinguish these two ecoregions. Reboratti (2005) sustained that Puna has a relatively flat relief whereas High Andes exhibits steepest slopes. The characteristics of the relief certainly impacts on water flow affecting the aquatic communities. Despite both ecoregions share many features attributable to their common geological origin (Morello et al., 2012), they should be considered two different units. Curiously, Alzérreca et al. (2001), establish the altitude of 4100 m as the critical limits to distinguish different types of peatbogs in Bolivia.
Electric conductivity resulted to be the second main characteristic behind compositional changes of benthos. The conductivity is a variable of primary interest since it is correlated with salinity which is in turn an indirect measure of geology and concentration of dissolved ions in the water body (González-Achem et al., 2015). It has also been found to be a synthetic variable for understanding the spatio-temporal dynamics of benthos from Yungas in Northwestern Argentina (Reynaga and Dos Santos, 2012). Critical limits have been established but they seem to depend on the study region. Some authors consider the level of 1500 μS/cm as critical for explaining differences in the structure of macroinvertebrate communities (Hart et al., 1991; Marshall and Bailey, 2004). However, the critical cutoff was established at 800 μS/cm by Horrigan et al. (2005) and Gonzalez-Achem et al. (2015). In terms of salinity as a predictor of community structure, the region of Puna sensu lato has not been sufficiently explored yet. Our critical threshold of conductivity fixed at 135 μS/cm is low in comparison to the previous values, and probably it is because of the highly oligotrophic and low order nature of streams (if not headwaters) from which peat bogs analyzed here originate. Harrel and Dorris (1968) found that conductivity increased with stream order because of increased leaching of soluble minerals with greater stream flow and depth of stream channels. On the other hand, the more mineralized peatbogs are close to the Calama-Olacapato-El Toro fault. Many stratovolcanoes are aligned along or near to this system (Norini et al., 2013) which could in turn have an influence on the water conductivity through volcanic ash (Stewart et al., 2006). So instead, the low cutoff of 135 μS/cm marks the distinction between virtually demineralized and more mineralized waters. Strikingly, Wang et al. (2012) studied the deposition of minerals dissolved in lakes from Tibetan Plateau and coined the term hyperfresh condition whenever conductivity is less than 0.13 mS/cm.
The relatively good agreement between compositional similarity and hierarchical arrangement of samples in function of their ecological distance, calls for a hypothetical model of ‘ecological filters’ (sensu Poff, 1997) acting on the study area. Filters are essentially environmental factors belonging to different hierarchical levels of a landscape configuration that work together by a successive selective process. In that way, filters limit the development of some organisms and are associated with species constraints, thus explaining the structure of local communities (Myers and Harms, 2009). Heino et al. (2007) stresses the importance of this approach (multi-scale regulation) to understand the organization of freshwater macroinvertebrates. Here, the first group division is completely explained by altitude alone at an unambiguous threshold of 4200 m asl between the lowland and upland groups. The next splitting of elements (fresh and hyperfresh subgroups) is highly concordant with a classification rule established on the ground of conductivity set at 135 μS cm–1. Altitude can be interpreted as a proxy for a filter acting at a regional scale (correlated with distinction between High Andean and Punean ecoregions), while conductivity represents a filter at a local scale. The perspective of filters needs to be explored more in depth, separating effects on taxonomic and functional aspects (Heino et al., 2007). Future research will be oriented in this direction.
Indicator taxa were projected onto the classification of sites from the ecological distances predicted by GDM model. They helped to reveal the biological side or compositional turnover behind the beta diversity patterns. Taxa pointed out here as characteristics from a given class have already been reported in other surveys from comparable regions. This is the case of Podonominae, which was positively related to the altitude (Tejerina and Molineri, 2007; Scheibler et al., 2008) and was considered as a cold stenothermal group from Andean rivers (Acosta Rivas, 2009; Epele et al., 2012; Medina et al., 2008; Scheibler et al., 2008, 2014; Tejerina and Molineri, 2007; Rodríguez Garay et al., 2015). Dolichopodidae has been consistently recovered as indicator of hyperfresh environments. This finding is remarkable since the family has been reported in hyper-saline inland waters (Velasco et al., 2006). Certainly, much effort to improve natural history knowledge into this family would shed light about specific preferences to aquatic habitats. Hyalella was not recovered as an indicator taxon of a specific cluster of sites. In fact, this amphipod is dominant throughout the area and can be considered a typical component of the region. Given the large amount of relative biomass, it is likely the main energetic subside for waterbirds. Moreover, the presence of Hyalella in distant and isolated peatbogs could be explained by dispersal me- diated by waterbirds (Swanson, 1984). Another remarkable taxon is the leech family of Glossiphonidae. Leeches were collected from a relatively large percentage of peatbogs (60%) averaging around 5 specimens when present. This taxon was collected in the Acay Sur peatbog at an altitude of 4731 m and this record possibly constitutes the highest collecting locality ever for leeches. Weber’s description of four Helobdella species in 1916 (as cited in Siddall, 2001) reported purportedly the maximum altitudinal record for leeches (i.e., 5140 m) corresponding to lakes near Cerro of Pasco from Perú. However, we have checked the altitude of such lakes (e.g., Huaron, Naticocha) and never overpass the height of 4700 m. So, to our best knowledge, specimens collected in the Acay peatbog represent the upper bound in the altitudinal range of leeches worldwide. Studying a similar environment but in the Bolivian Andes, Siddall (2001) described four species of Hirudinea and pointed out that this kind of ecosystem at high altitude is particularly favorable to the diversity of the group.
Lastly, we acknowledge the importance of the taxonomic resolution in our results. The invertebrates were identified to the lowest feasible taxonomic level. In an ideal scenario, analyses should be based on taxonomic units identified to the species level. With benthic fauna, this goal is hard to reach since organisms are often immature and not identifiable to species level. The taxonomic impediment is here enhanced by the poor knowledge of the fauna associated with this remote area. Eventually, analyses could rely on a higher taxonomic resolution, like family level, for all groups. However, this last procedure would entail a loss of information to outline patterns of compositional dissimilarities. Anyway, our results are not expected to vary significantly if the great majority of taxa were monotypic.
Overall, according to our results, benthic communities can be classified based on their ecological requirements, first by altitude (cutoff = 4200 m) and then by conductivity (cutoff = 135 μS cm–1), suggesting a system of ecological filters acting on them. Altitude is taken as the surrogate of a regional factor (distinction between Puna and High Andean ecoregions) whereas conductivity the proxy for a more local factor (distinction between hyperfresh and more mineralized habitats). By being near the edges of life, high elevation peatbog communities are potentially very sensitive to climatic and land use changes. The region is expected to experience substantial increases in temperature in relation to global atmospheric changes (Izquierdo et al., 2015b); the importance of altitude, a geographic variable highly related to temperature, suggests that significant climatic change can result in significant biotic changes. During the past decades, the region has experienced a significant decrease in rainfall (Carilla et al., 2013; Morales et al., 2015). While our results do not show a strong association between community composition and rainfall, trends in rainfall patterns may result in changes in water ion concentration, which was shown to be a major environmental control. Mining and urbanization in the region could also affect water characteristics (Izquierdo et al., 2015b), and this needs to be studied to explore scenarios of environmental change in the region.