Environmental and ecological measurements of the four studied systems
At each zone (1-4), we selected one stream site of the glacial stream derived from Carihuairazo’s glacier. Invertebrate community was estimated based on five Surber samples (0.05 m2; mesh size 200μm, and preserved in 70% ethanol). This sampling was done at six different dates from May 2015 and May 2017 for sites 2-4 while site 1 was only sampled on May 2017 since the glacier was still covering this zone. In addition, at each sampling site and date, fifteen pebbles were collected in five groups of three stones, and conserved in 96% ethanol in the dark to assess chlorophyll A concentration in order to quantify the available resource.
In the laboratory, Surber samples were rinsed through a 200μm sieve and sorted thoroughly by hand in a standardized manner, without use of magnification. Macroinvertebrates were identified under a microscope at 4.5x magnification range to different taxonomic levels, mostly to family level using a dichotomous identification key of benthic invertebrates of South America (Fernandez and Dominguez, 2001). The benthic detritus left on each Surber was collected and dried in an oven for 24 hours at 80°C and then burnt to quantify the course particulate organic matter (CPOM), an additional measure of available resources for benthic macroinvertebrates. The ethanol of the stones samples were transferred to a spectrophotometer and absorption was measured at 665 and 750 nm. The ChloA (μg.cm-2) was then calculated using the method proposed by Københavns Universitet (1989) that uses the absorbance, the volume of ethanol used and the pebble surface in order to calculate the ChloA. The pebbles surface was calculated according to the formula proposed by Graham et al. (1988) that uses the length (L), width (W) and height (H) of the pebbles. In addition, physicochemical measurements were conducted for every taxa sampling. Water temperature (°C) and pH were measured with a pH meter. Conductivity (μS.cm-2) was measured using a conductivity meter. Oxygen saturation (%) and dissolved oxygen (mg.L-1) were measured using an oxygen meter. Turbidity (NTU) was measured using a turbidity meter
Terrestrial system – Flora
At each zone (1-4), terrestrial plants were observed between September and October 2017 in ten quadrants (1 m2) along 2 transects (five quadrants per transect). Quadrants were placed on both sides of the river for each site. Quadrants were subsampled in smaller quadrants of 10×10 cm in order to be more precise on the estimations. For each quadrant, the cover area (%) was determined for every flora taxon or substrate present. Flora was determined to a species level when possible. The total cover of each substrate for each smaller quadrant was then summed and divided by 100 in order to obtain a total cover for each quadrant.
Terrestrial system – Carabidae
At each zone (1-4), Carabidae were sampled using pitfall traps. Six pitfalls traps were placed during 12 or 14 days in three dates between 2016 and 2017 using wine vinegar as attractive solution. Taxa was then collected and preserved in ethanol. In the laboratory, Carabidae were identified under a microscope at 4.5x magnification rate to species level.
Terrestrial system – Soil DNA
Five soil sites were selected for each zone (2-5). Sampling was done using a shovel that was sterilized for every site so DNA would not mix up between sampling sites. Four soil samples were extracted per site and preserved in a sterilized container for further analysis. Total DNA was extracted using the PowerSoilVR DNA isolation kit (MoBio Laboratories) and two amplicon libraries were prepared for bacteria and eukaryotes. Libraries were sequenced on the Ion Torrent Personal Genome Machine. DNA sequences were filtered using the OBITOOLS software (http://metabarcoding.org/obitools), and were assigned to the relevant taxa using ecoTag program that finds highly similar sequence(s) in a suitable database (Ficetola et al. 2010).
Terrestrial system – environmental measures
At each zone (1-4), relative humidity and temperature was measured every 30 minutes using relative humidity/temperature data loggers for Carabidae or only temperature data loggers for soil DNA (Hobo relative humidity and temperature loggers, Onset Computer Corp) during two pitfall sampling duration on each site and one month period for soil DNA. Data were then used to calculate minimal, maximal, mean and standard deviation temperature as well as mean relative humidity for the measurement period.
Beta-diversity and physicochemical/spatial measurements effects on it
β-Diversity can be defined as a variation in species composition among sites which provides a link between α-diversity (local diversity) and γ-Diversity (regional diversity) (Anderson et al. 2011). To analyse the dissimilarity among the study zones, we calculated Bray-Curtis dissimilarity index (β-Bray, based on abundance data) (Bray & Curtis, 1957) and the Sørensen dissimilarity index (β-Sør, based on presence-absence data) (Sørensen, 1948). For each system, β-Bray was used to perform a non-metric multidimensional scaling (NMDS) in order to assess dissimilarities in community composition among all the sites for each system. Physicochemical and spatial measurements were also fitted into the NMDS plots using the envfit() function of the vegan package in R. This function allows fitting environmental variables onto the ordination in the form of vectors which will have the maximum possible correlation with the environmental variable used, indicating thus which environmental variable has an effect on the community variation among sites.
The variation of taxa assemblage can be partitioned into two components: dissimilarity due to turnover and dissimilarity due to nestedness where turnover represents the replacement of some species by other species and where nestedness represents the loss of species (the poorest assemblage being a subset of the richest assemblage) (Baselga, 2010). This concept of partitioning β-Diversity was first applied by Baselga (2010) on the Sørensen dissimilarity index where β-Sørturn represents the turnover and β-Sørnest represents the nestedness. It was then applied by Baselga (2013) to partition β-Bray into the balanced variation in species abundances, β-Brayturn, representing a replacement of a number of individuals for the same number of individuals of another species from one site to the other; and the abundance gradient, β-Braynest representing the loss of some individuals from one site to the other. In our study both β-Sør and β-Bray were calculated, partitioned and used to perform dissimilarity matrices between sites. These analyses were done in R (R Development Core Team 2013, version 3.3.3), using the R package betapart. Distance matrices were then performed for environmental and spatial variables for each system comparing all measurements done per variable with each the other measurements done. β-Sør, β-Bray and their partitions were then plotted against the differences of physicochemical and spatial variables then fitted to linear models for each system. Multiple regression of distance matrices (MRM) was used following Lichten’s (2007) method in order to assess statistical significance.
Interactions between studied systems
Because available data from the studied systems did not have the same length, simulations were done for the systems with less data to test the effect of biotic system between each other. All data was then fit to linear models where the factor was the studied site and the linear response corresponded to every possible combination of α-Diversity and taxa count (density, cover, relative abundance; depending on the system) from each system. An anova was then performed to the models in order to assess possible interactions between the biological systems and also to assess if interaction between the site and the variable had a role in the response. Combinations of variables having a significant p-value in the anova test were then plotted.
All analyses and figures produced were done in R (R Development Core Team 2013, version 3.3.3) using packages vegan, PMCMR, betapart, cowplot and ecodistSite differences on α-Diversity, abundance, Shannon and Pielou’s evenness per system.
Communities composition and physicochemical/spatial variables explaining this composition
The NMDS done had all good enough values of stress (0.17; 0.14; 0.006; 0.1 for aquatic, flora, Carabidae and soil DNA respectively). On each one of the systems NMDS revealed us possible differences on community compositions. As we can see for benthic macroinvertebrates (Fig. 4a), site 1 (that has a single point) community composition appears to be very distant from the rest of the sites community composition. On this system, site 2 seemingly is explained by the altitude and conductivity and is different from sites 3 and 4. Sites 3 and 4 have much more similarities than the rest of the sites and are rather explained by the temperature, the chlorophyll A content and the distance to glacier. For the flora system (Fig. 4b) we can see that site 1 remains very different from the rest of sites but unlike benthic macroinvertebrates sites 2 and 3 are more similar in terms of community composition meanwhile site 4 appears to be the much more different from these two sites. The community composition of site 2 is explained by the altitude while for site 4 community composition is explained by distance to glacier. For Carabidae (Fig. 4c) we can see that all sites appear to have different community composition with a slight similarity between sites 3 and 4. The composition of the community on site 1 appears to be explained by the distance to the stream. Finally, for soil DNA (Fig. 4d) we can see that all sites have different community compositions. This finding is rather interesting given the fact that on the previous section no differences were shown on calculated indices between sites. The composition of site 2 appears to be explained by altitude and max temperature while the composition of site 4 appears to be explained by distance to stream and standard deviation temperature.
P-values of MRM tests with significant values (p-value <0.05, Annex 6) showed us that variation in β-Diversity (Bray-Curtis with its partitions and Sørensen with its partitions) could be explained by the variation of spatial or physicochemical measures. For benthic macroinvertebrates ΔpH affected β-Sør and β-Sørnest. ΔTemperature (ΔT) affected on β-Brayturn and β-Sør. ΔChloA affected both β-Brayturn and β-Sørturn. Finally, the variation of distance to the glacier ΔDistance to Glacier (ΔDG) affected both β-Sør and β-Bray (Fig. 5a).
Succession patterns along the deglaciation gradient
Results showed an increasing gradient of α-Diversity and abundance with increasing deglaciation age among sites for both flora and aquatic invertebrates as found in Milner et al. (2000) for aquatic invertebrates. Increase of α-Diversity was more pronounced on late colonization stages for the flora and aquatic invertebrates. The increase in both flora cover and α-Diversity was accompanied by an increase of existing litter particularly in the latest stage of primary succession that may be due to the establishment of new vascular plants along the deglaciation gradient. The litter increase could be the responsible of the observed increase of Enchytraeidae at the same stage of succession. Surprisingly, α-Diversity of Carabidae didn’t change with along the deglaciation gradient, which is opposed to the findings of Greben-Krenn et al. (2011) in the Alps where sites of later age since deglaciation presented higher values of both α-Diversity and abundance. In the Carihuairazo only Carabidae abundance increased along deglaciation gradient for studied sites. The first colonizers for flora and aquatic invertebrates were well-adapted taxa to cold temperature and freezing such as Chironomidae (Irons, Miller & Oswood, 1993) or well adapted taxa to desiccation such as Bryophyta (Proctor, Ligrone & Duckett, 2007), which have more facility to establish. Such was the case of our stream where Chironomidae were the dominant taxa in all sites that decreased with the establishment of new taxonomic groups. This dominance decreased along the stream as new taxa colonized it. This increase in taxa could be translated as an increase in present functional groups along the stream. Most abundant Chironomidae was Orthocladiinae followed by Podonominae. There was an absence of Diamesinae along the stream, which is opposed to the findings of Fureder et al. (2001) where Diamesinae were the dominant taxon in the Alps. The finding of high content of chlorophyll A (2.937 μg.cm-2) among studied zones in earliest colonization stage gives us also an idea on the earliest colonization in streams by primary producers assembled in biofilms that could be one of the consumed resources by Chironomidae. The most homogeneous evenness values were found in a mid-late stage of primary succession (site 3), which could mean a stabilization of community structure before the colonization of a second wave of arriving taxa. Studied sites presented overall differences of community composition for all systems, of which studied sites 3 and 4 were more similar for aquatic invertebrates and sites 2 and 3 were more similar for flora. Thus, the community composition along deglaciation gradient varies with the establishment of new taxa or individuals with increasing deglaciation age as well as the species or individual turnover. Primary succession after glacial retreat along deglaciation gradient appeared to follow the following patterns for both flora and taxa: An earliest colonization (site 1) of high-tolerant taxa as pioneers of colonization; an early colonization (site 2) characterized by the arrival of new taxa and an increase in both functional and taxonomic diversity that could imply the establishment of more complex interactions within studied communities. A rather stable mid-late colonization (site 3) with arrival of a few more taxa, high values of evenness and the presence of high vascular plant DNA of taxa that possibly could not yet establish in this site because of time-lag between glacial retreat and upward migration in the tropical Andes (Zimmer et al. 2018). Finally a late colonization (site 4) characterized by a sharp increase in present taxa and abundance where possibly a second wave of colonizers arrived due to more favorable conditions and interactions that allowed them to establish. Studied Carabidae didn’t follow these patterns as all species were found in all sites and only abundance increased along deglaciation gradient. This could mean that there are differences in primary succession that could depend on taxonomic systems even though we lack on information of other terrestrial invertebrates such as spiders to see if terrestrial invertebrates have a similar pattern as the one observed for studied Carabidae.
Table of contents :
II. Material and Methods
1) Study site
2) Environmental and ecological measurements of the four studied systems
a) Aquatic system
b) Terrestrial system – Flora
c) Terrestrial system – Carabidae
d) Terrestrial system – Soil DNA
e) Terrestrial system – environmental measures
3) Data analyses
a) Comparison between study sites
b) Beta-diversity and physicochemical/spatial measurements effects on it
c) Interactions between studied systems
1) Site differences on α-Diversity, abundance, Shannon and Pielou’s evenness per system.
2) Communities composition and physicochemical/spatial variables explaining this composition
3) Interaction between biological systems
1) Succession patterns along the deglaciation gradient
2) Environmental and spatial filtering primary succession
3) Biotic interactions