Spatial and Temporal Ecological Uniqueness of Andean Diatom Communities Are Correlated With Climate, Geodiversity and Long-Term Limnological Change

Introduction

Human activities have caused in the past and are currently causing diverse and long-lasting changes in freshwater ecosystems (Vitousek et al., 1997). In mountainous areas, high-elevation lakes are excellent sentinels of current global change and some of the most comparable ecosystems across the world (Catalan and Rondón, 2016). It is widely recognized that predicting how high-elevation lakes will respond to ongoing and future global changes requires a long-term perspective to evaluate recent (last ca. 150 years) human-driven impacts and to characterize background natural variability (Mills et al., 2017; Dubois et al., 2018). Biological assemblages accumulate with lake sediments as natural archives, which can be used to understand temporal dynamics of biodiversity and provide insights into the organization of ecological communities and their responses to natural and human-induced drivers (e.g., habitat loss, human impacts, eutrophication) (Willis et al., 2010; Heino et al., 2016). Because lakes are not isolated in the landscape (rather they form a continuum embedded in a terrestrial matrix), researchers have also examined the role of spatial variables in determining biodiversity patterns using paleolimnological approaches (Castillo-Escrivà et al., 2017; Benito et al., 2019). However, contemporary and paleolimnological studies still remain largely disconnected in biodiversity and environmental change research (Gregory-Eaves and Beisner, 2011). One reason for this may be that the application of the metacommunity concept (i.e., set of local communities potentially connected by dispersal of multiple interacting species, Leibold et al., 2004), has so far seldomly been applied over long time scales. Understanding broad-scale biodiversity patterns is necessary, because many environmental pressures operate at large spatial and long temporal scales but interactions with high-elevation lakes at smaller scales (Catalan et al., 2013) are still mostly unexplored.

Spatial variability in species composition between sites (i.e., spatial beta diversity) is particularly well suited for exploring metacommunity structuring, because its components: species replacement and richness difference (species gain and loss), are often correlated with climatic, geological, and limnological gradients (Winegardner et al., 2017; Castro et al., 2019; Szabó et al., 2019). With an index of local contribution to beta diversity (LCBD) it is possible to examine each site’s disproportionate contribution to maintaining regional diversity. High LCBD values flag ecologically unique sites where regionally rare taxa, exceptional species richness or unusual combinations of taxa are present (Legendre and De Cáceres, 2013). Both niche-based (i.e., species sorting) and dispersal-based community assembly processes may influence beta diversity patterns, either independently or in combination (Heino, 2013). Thus, it is important to include different metrics of beta diversity and their environmental and spatial correlates to assess biodiversity changes across ecological gradients in a more nuanced way.

In aquatic ecosystems, comparisons of temporal beta diversity between localities may be a fingerprint of environmental variation, for instance, as a measure of acidity changes in boreal lakes (Angeler, 2013) or forest cover in United States lakes (Winegardner et al., 2017). Recent meta-analyses have linked beta diversity loss, increasing trends between assemblages’ similarities to reduced environmental heterogeneity due to human actions (Cardinale et al., 2012; McGill et al., 2015). Other studies further indicated that ecosystem function and services, such as carbon storage and food production, could be severely jeopardized by biotic homogenization (Van der Plas et al., 2016). However, patterns prior to intense human occupation remain largely unknown, and a time interval spanning the last 1000 to 2000 years is most appropriate for evaluating the magnitude of recent changes in biodiversity and environmental conditions at regional and global scales (Pages 2k Consortium, 2013).

The tropical Andes is one of the world’s richest biodiversity hotspot (Myers et al., 2000). Lakes are ecologically important regional features and also serve as a crucial source of freshwater for millions of people living in the Andes and the adjacent Amazon lowlands (Buytaert et al., 2006). Historically, humans have been an integral part of Andean lake-catchment systems, shaping cultural landscapes by means of agriculture, pastoralism, and deforestation (Sarmiento, 2002). Andean lakes in tropical South America are valuable model systems for examining spatial and temporal beta diversity patterns for several reasons. First, they are well-defined ecosystems, many of which have persisted over long geological time-scales such as lakes Titicaca and Umayo, allowing evolutionary processes (speciation, extinction) to affect biota (Fritz et al., 2012). Second, despite contrasts in the climatic and evolutionary history of various regions, broadly similar ecoregions and communities are present due to the relative tectonic stability of the Andes throughout much of the Neogene (Baker and Fritz, 2015). Third, they are in a topographically diverse landscape characterized by notable spatial variability in geological, geomorphological, soil, and hydrological features, also known as geodiversity (Killeen et al., 2007; Gray, 2008). In the tropical Andes many geodiversity variables are available for studying their influence on biological communities (Antonelli et al., 2018). However, very little is known about how geodiversity and other macroscale environmental gradients, such as climate, relate to beta diversity and its replacement and richness components in mountain lakes.

Biodiversity studies of Neotropical macroorganisms have historically recognized the roles of environmental, climatic, and geological factors in structuring communities at different spatial and temporal scales (Banda et al., 2016), but their influences on microorganisms are not well characterized (Benito and Fritz, 2020). Diatoms, unicellular siliceous algae, are a very species-rich biological group that disperses widely, responds to local (e.g., limnological) and regional (e.g., climatic) variables, has different traits for resource use and resistance to disturbance, and their remains preserve in lake sediments. In addition to their role as whole-ecosystem indicators, diatoms are also particularly useful for evaluating relationships between spatial and temporal beta diversity (Winegardner et al., 2017). Yet, additional exploration of their use in identifying ecologically unique sites and the mechanisms behind spatial and temporal beta diversity is still needed for tropical Andean lakes. This information is crucial as global change intensifies in the Neotropics (Vuille et al., 2003), including local (e.g., agriculture, fish stocking) and regional (e.g., deforestation, damming, land-use) environmental impacts (Van Colen et al., 2017). Hence, studies comparing modern and paleolimnological records are needed to assess which high-elevation lake ecosystems are most sensitive to environmental changes and to evaluate their resilience.

To start filling these knowledge gaps, we used a database of lake diatoms that spans large gradients of latitude, limnology, climate, and topography in the tropical Andes. We estimated LCBD and applied beta diversity partitioning methods, including species replacement and richness difference and the contributions of sites to these components (LCBDrich and LCBDrepl). LCBDrich and LCBDrepl describe how each individual sample contributes to richness gradients and to replacement gradients (Ruhí et al., 2017). To identify spatial beta diversity trends through time, we analyzed diatom paleolimnological assemblages that span the last ∼2000 years using the same partitioning methods with Generalized Additive Models (GAM). Finally, using diatom beta replacement time series we tested the hypothesis that broad-scale environmental changes have an homogenizing effect across large, disconnected lakes; to provide a regional long-term limnological change perspective in the tropical Andes.

Materials and Methods

Study Area

Our study lakes are distributed across the tropical Andes (8°N–30°S and 58–79°W) (Figure 1). This region encompasses a wide range of physiographic and climatic settings that produce diverse limnological conditions. The investigated lakes occupy a range of high elevations (2500–4500 m a.s.l) and are mostly formed by glacial and/or volcanic processes. The study area has a north-south orientation and is characterized by varied degrees of topographic heterogeneity. Both local and regional climates are influenced by the topographic profile, which creates distinct conditions at both eastern and western flanks (Valencia et al., 2016). Northern Andean lakes in Ecuador and Colombia lie in montane forests, inter-Andean valleys, and páramo ecosystems. Above the tree line (páramo), climatic conditions are characterized by the lack of seasonal changes and cold mean annual temperature. More dry and wet climatic conditions are characteristic of the interandean valleys and montane forests, respectively. In the central Andean Cordillera of Perú and Bolivia, most of the study lakes are in montane grassland or shrubland. In the Altiplano plateau (central Andes), the northern region is characterized by cold and relatively humid conditions. Lakes are mainly freshwater and lie in extensive interconnected hinterland basins (Cohen et al., 2014). The southern Altiplano is drier, and most lakes are isolated and saline due to the basin geology and high evaporation rates (Blanco et al., 2013). Lakes in the south-central Chilean Andes are located at lower elevation (<2500 m a.s.l) and are surrounded by steppe vegetation; climatic conditions are semi-arid to temperate (Carrevedo et al., 2015).

Figure 1. (A) Geographical location of the modern investigated lakes (n = 144) colored by regions within the tropical Andes of South America (see Supplementary Table 1 for individual lakes information in each region); the location of the six lake cores are also shown with different symbols; (B) Proportion of diatom ecological groups for each lake region. Samples represent sediment surface habitats.

Diatom Database

We used diatom abundance data from 144 lakes comprising sediment surface samples (n = 215) from a newly created tropical South American diatom database available in the Dryad repository (Benito et al., 2018b) (Figure 1A and Supplementary Table 1) and GitHub¹. Briefly, the database comprises published and unpublished studies from lentic and lotic environments collected by different authors under different objectives (e.g., paleoclimatic reconstructions, taxonomy, biodiversity). When possible, diatoms were identified to the species level. The samples were collected in the period 1999–2017. Detailed information about sample processing, taxonomic harmonization, and identification of diatom taxa can be found in Benito et al. (2018b).

In addition, we used sedimentary cores from six lakes located in Ecuador (Piñan, Yahuarcocha, Fondococha, Llaviucu), Peru (Umayo), and Bolivia (Titicaca) for temporal beta diversity analyses (Figure 1A). Sediment cores from the Ecuador lakes were collected in summer 2014 (Llaviucu and Fondococha) and July 2017 (Piñan and Yahuarcocha). Cores (mean core length = 61 cm) were retrieved using a UWITEC gravity corer near the center of each lake when possible; lake Piñan was cored at the south-west shallow platform. Sediment cores were sliced in the field at 1 cm intervals. In the laboratory, samples were processed for diatom analyses following standard methods (Battarbee et al., 2002). At least 300 valves were counted per sample and identified using diatom regional floras (Metzeltin and Lange-Bertalot, 1998, 2007; Rumrich et al., 2000; Metzeltin et al., 2005). Core chronologies were established using ²¹⁰Pb and ¹⁴C dating techniques. For Lake Fondococha, details about the ²¹⁰Pb-chronology can be found in Bandowe et al. (2018) and information about the age-depth model is described in Arcusa et al. (2020) and Schneider et al. (2018). Instantly deposited event layers (e.g., tephra layers and flood layers) were masked for the age calculations and reinserted in the combined age-depth model. The list of ¹⁴C and ²¹⁰Pb dates and associated age-depth models for the lakes Piñan, Yahuarcocha, Fondococha and Llaviucu can be found in Supplementary Material. Diatom core assemblages of Lakes Umayo and Titicaca (Ekdahl et al., 2008; Weide et al., 2017) and from Piñan and Fondococha in Ecuador (Luethje, 2020) are published. Altogether, diatom records span the last 1102 (Piñan), 1815 (Yahuarcocha), 2598 (Fondococha), ∼2250 (Llaviucu), ∼6500 (Umayo), and ∼6000 (Titicaca) cal years BP (cal year BP = years before 1950), with a median temporal resolution of 20.8 years (range 8–44 years).

Predictor Variables

We extracted different datasets from several sources to characterize local (limnological) and regional (climatic and geological) environmental characteristics, as well as historical effects of the investigated lakes. Lake water-chemistry variables were collected simultaneously with the modern diatom samples and included water temperature (°C), pH, conductivity (μS/cm), cations (Ca²⁺, Mg²⁺, K⁺, Na⁺; mg/L) and anions (Cl^–, SO₄^2–; mg/L). Nutrient data (N, P) are not included here, because the database has a very high number of missing values. However, our prior analyses indicate that nutrient conditions are strongly correlated with landscape factors—total phosphorus decreases with elevation—and also that geo-climatic factors are more highly correlated with diatom diversity than limnological conditions in most of the Andean lake regions (Benito et al., 2018b). Elevation (m) and lake area (km²) were extracted using ArcGIS from the STRM 90 m Digital Elevation Model (Jarvis et al., 2008) and using ESRI World Imagery layer as a basemap, respectively. In equal grids of 50 km², the Global Lakes and Wetlands Database (GLWD at ∼ 1 km resolution; Lehner and Döll, 2004) was used to extract the surface area occupied by fresh waters. The percentage of aquatic systems in the surface area acts as a proxy for hydrological connectivity. Climatic variables included mean annual air temperature (MAT; °C), mean annual precipitation (MAP; mm), temperature seasonality (SD; °C), and precipitation seasonality (coefficient of variation; mm). These variables were extracted from the WorldClim 1.4 database (Hijmans et al., 2005). WorldClim contains averaged monthly climate data for the period 1950 to 2000 at 1 km of spatial resolution and falls well within the temporal window of the analyzed diatom samples. To account for the effect of geodiversity, for each lake we extracted geological [soil variability (number of soil types per grid cell), long-term erosion (km/Ma)] and topographic (terrain ruggedness index) variables from Antonelli et al. (2018) within 1° × 1° rid cell (∼80 km² at the equator). Soil variability is the number of soil types in each grid cell derived from Hengl et al. (2014), while long-term erosion is derived from termochronometric data using the Herman’s method which accounts for topography and isotopic data to generate maps of averaged erosion rates over time (Fox et al., 2014). Finally, to estimate the human historical footprint in the lakes, we obtained human density and cropland area from the HYDE 3.2 database (Goldewijk et al., 2011). We extracted human density (inhabitants/km² grid cell) and cropland area (km²/ grid cell) values for the last 300 years at 10 years timesteps and averaged over three estimate scenarios (baseline, lower, and upper) for each investigated lake within a spatial resolution of ∼80km².

Statistical Analyses

All statistical analyses were performed using the R software version 3.6.2 (R Development Team, 2016).

Prior to running statistical analyses, predictor variables were checked and transformed accordingly [log10(x + 0.25) or square root] to meet assumptions of linearity and homogeneity of variances. In both modern and fossil diatom matrices, those species having >3% relative abundance in at least one sample were selected.

Hellinger-transformed diatom relative abundances were used for estimating beta diversity indices. First, we calculated beta diversity for each lake and partitioned it into components, namely replacement and richness difference, with the Podani decomposition family of indices using Bray-Curtis dissimilarity (Podani et al., 2013) with the beta.div.comp function of the adespatial package (Dray et al., 2016). The replacement component refers to simultaneous species loss and gain along ecological gradients (in space or time), whereas the species richness component means that one sample contains more unique taxa than another (Podani et al., 2013). Second, to further investigate mechanisms behind beta diversity patterns across space (144 lakes) and time (paleolimnological time series from six lake cores), we estimated the local contribution of each sample (sediment surface and 1-cm sample slice for space and time, respectively) to the total beta diversity (LCBD) using the beta.div function. The significance of each LCBD value was assessed by 999 permutations, and the p-values were corrected for multiple testing using Holm’s procedure (Dray et al., 2016). We also calculated site-specific diatom richness and related with LCBD using Spearman correlation. Third, from the initial richness and replacement matrices, we decomposed LCBD to richness (LCBDrich) and replacement (LCBDrepl) components to assess how each individual sample (in space and time) contributes to richness and to replacement gradients, respectively, in the diatom communities using the LCBD.comp function.

To aid interpretation of local contributions to beta diversity, we calculated mean percentages of diatom ecological groups in each of the regions of study (Figure 1A). We placed diatom species into the following groups: freshwater planktic (e.g., Discostella stelligera, Cyclostephanos andinus), oligosaline planktic (Cyclotella meneghiniana), tychoplanktonic (e.g., Fragilaria capucina, Aulacoseira alpigena), freshwater benthic (e.g., Achnanthidium minutissimum), epiphytic (e.g., Cocconeis placentula), and saline (e.g., Craticula halophila, Navicula salinicola), following regional diatom floristic studies (Steinitz-Kannan et al., 1993; Metzeltin and Lange-Bertalot, 1998, 2007; Weide et al., 2017).

We ran GAM to model the relationship between LCBD and its replacement (LCBDrepl) and richness (LCBDrich) components and local, regional, and historical predictors. GAMs are a non-parametric extension of the Generalized Linear Models and allow fitting linear and non-linear relationships between the response and explanatory variables when there is no a priori reason for choosing a particular function (i.e., linear, quadratic) (Wood, 2017). Only variables that had Variation Inflation Factor (VIF) values <10 were included in the GAMs. We estimated the linear effect of each predictor, accounted for spatial autocorrelation by including smooth splines of geographical coordinates, and included lake region (as in Figure 1) as a random factor (bs = ‘re’). Both statistically significant predictors and the level of complexity of the response shapes to each variable were selected with Restricted Maximum Likelihood (REML) using the mgcv package (Wood, 2011). Finally, we checked residuals for any deviation from normality and linearity using diagnostic plots (Supplementary Figure 2).

To determine if statistically significant change in temporal beta diversity trends at millennial time scales could be identified, we modeled the beta replacement time series using a Hierarchical GAM (HGAM) (Pedersen et al., 2019). Here, we were interested in comparing spatial patterns in temporal beta diversity. We used the vector of beta replacement values that resulted from the beta.div.comp function as a response variable. We tested the hypothesis that variations in temporal beta diversity are homogenous across the landscape (biotic homogenization) or whether temporal changes differed from lake to lake. For this, we built two HGAMs separately: (a) a global smoother and lake-level smoothers having different wiggliness (i.e., individual curves), hence allow for inter-lake variability (HGAM GI model type in Pedersen et al., 2019) and (b) a global smoother and lake-level smoothers that have the same wiggliness (i.e., shared curves), hence do not allow for inter-lake variability (HGAM GS model type in Pedersen et al., 2019). In all models, we accounted for the different amount of time each sediment core sample represents (difference between ages at the top and bottom of each sediment slice) by including these values as weights in the model (Simpson, 2018). Since the beta diversity replacement range from 0 to 1 in the form of relative values, a gaussian link function was applied. We applied two methods for model selection: (a) AIC values using a cut-off level of two units or less from the lowest AIC model (Burnham and Anderson, 2004), and (b) out-of-sample deviance performance, where each model was compared to a null model (intercept-only model with only lake-level random effects intercepts included). We fitted all the models using the gam function in the mgcv package.

Results

Ecological Uniqueness and Its Components

Freshwater benthic diatoms dominate across the study regions, except for Lake Titicaca with dominance of freshwater planktic species (Figure 1B). Freshwater planktic diatoms are the second most abundant ecological group, followed by saline-tolerant and epiphytic benthic species. Among the freshwater benthic taxa, Achnanthidium minutissimum s.l., Amphora veneta, Cymbella cistula, and Staurosira construens var. venter made up the majority of the sediment surface diatom flora (Supplementary Figure 1). These taxa are indicative of a wide range of limnological conditions, ranging from acidic, low nutrient to high conductivity waters. Freshwater planktic diatoms characteristic of low nutrient conditions were dominated by Cyclostephanos andinus, whereas Aulacoseira ambigua, Aulacoseira granulata, Discostella stelligera, and Tabellaria flocculosa strain IV may respond to increased nutrients. Saline-tolerant and epiphytic benthic diatoms included the endemic Amphora carvajaliana, Epithemia adnata, Navicula salinicola, and Cocconeis placentula var. placentula, respectively.

Local contribution to beta diversity did not show statistically significant variation across the study lakes according to corrected p-values for multiple testing. Without correcting for multiple testing, 24 lakes (13% of total) had significant LCBDs, mostly concentrated in the south-central Andes of Chile (Supplementary Figure 2).

There was a significant negative relationship between LCBD and species richness (Spearman rho = −0.58; p < 0.01) (Supplementary Figure 3). GAM results showed that modern LCBD decreased with increasing pH, MAP, and historical footprint (Figure 2A). The full model explained 38.1% of the deviance. When analyzing the replacement and richness components of LCBD, the effects of environmental predictors were generally inverse (Figures 2B,C), as supported by the negative relationship between LCBDrepl and LCBDrich components (Supplementary Figure 4). A consistent set of variables representative of local (limnological), and regional (climatic and geological) conditions explained variation in LCBD indices. LCBDrepl increased with increasing conductivity and seasonality in temperature, and decreased with increasing terrain ruggedness, soil variability, and Na⁺ (full model % deviance explained = 37.7). LCBDrich increased with decreasing conductivity and seasonality in temperature, and increased with higher terrain ruggedness (full model % deviance explained = 31.4).

Figure 2. Effect of predictors on local contributions to beta diversity (LCBD) and to the replacement (LCBDrepl) and richness (LCBDrich) components, estimated as linear coefficients from Generalized Additive Models. Errors bars are ±95% confidence intervals. Colors indicate the significance level α = 0.05.

Temporal Trends of Beta Diversity Components

As expected, the temporal trends of LCBDrepl and LCBDrich components fluctuated over millennial-time scales (Figure 3). In general, LCBDrich fluctuated more and was comparatively higher than LCBDrepl across lakes, especially in the two-deep freshwater Altiplano lakes, Umayo and Titicaca. These two lakes also showed increased trend in LCBDrepl since ca 1000 cal years BP. Replacement and richness fluctuated more similarly in the two remote Ecuadorean páramo lakes (Piñan [Spearman rho: −0.17, p = 0.22] and Fondococha [Spearman rho: 0.51; p < 0.01]) than the two lakes located closer to human settlements (Yahuarcocha and Llaviucu [Spearman rho range: 0.02–0.09, p > 0.05]). Interestingly, LCBDrich and LCBDrepl time series of Piñan, Yahuarcocha and Fondococha showed a peak at ca. 1000 cal years BP. The main diatom stratigraphic changes and dominant taxa for each lake are summarized in the supplementary material (Supplementary Figures 10–15).

Figure 3. Contribution of replacement (LCBDrepl) and richness (LCBDrich) components to beta diversity for the six investigated lake cores (arranged by increased latitude). Isotopic δ¹⁸O measurements from Pumacocha lake (Bird et al., 2011) are interpreted as dry periods with enriched values. Note the two y-axes for replacement and richness components of beta diversity time series.

All the HGAMs models fitted to the beta replacement time series predicted better than the null models (Supplementary Table 2). The best HGAMs included a global smoother plus lake-specific smoothers having different wiggliness (GI model) according to the AIC models. These results suggest that allowing for lake-specific variation explained more variation in beta replacement trends. The shape of the fitted HGAMs differed across lakes. Beta replacement increased in the last ca 500 cal years BP in the deep Altiplano lakes (Umayo and Titicaca) (Figure 4). More coherent trends of beta replacement were observed in the Ecuadorean lakes, which were characterized by slight increases at ca. 1000 and 500 cal years BP and a decrease over the last ca. 200 years (Figure 4).

Figure 4. Temporal trends of Hierarchical GAM-fitted beta replacement for the six Andean lakes investigated in this study (arranged by increased latitude). Model GS (green line) indicates a global curve and inter-lake variation having similar smoothness. Model GI (orange line) indicates a global curve allowing for varied smoothness among lakes. Points are original beta replacement values. Ribbons indicate ±2 standard errors around the mean (green and orange lines for GS and GI models respectively).

Discussion

Previous studies have demonstrated that species composition and taxonomic richness of lake diatom communities in the tropical South America are jointly structured by the local (water chemistry) and regional environmental factors (aquatic connectivity and climate) (Benito et al., 2018b, a). Moreover, biogeographic patterns emerge after determining latitudinal gradients of species richness and estimating the role of dispersal dynamics on diatom community structure (Benito and Fritz, 2020). In the analyses presented here, the calculation of LCBD introduced an additional biodiversity metric for Neotropical diatom metacommunities and biogeography studies. Our results identified a set of local and regional ecological gradients that explained patterns in LCBD, including its replacement and richness components. We found that LCBD was related to pH, MAP, and historical human density. The effect of pH and MAP are not surprising given the relatively high variance displayed across the study lakes (Steinitz-Kannan et al., 1983; Michelutti et al., 2019) and the known direct role of pH in affecting physiological process in diatoms (Van Dam et al., 1994). Precipitation indirectly affects catchment-lake linkages through, for instance, biochemical processes and resource supply (Passy, 2010). Here, we found that high-elevation lakes lying in drier areas are ecologically more unique than lakes receiving more precipitation. Vilmi et al. (2020) showed a distinction between dry/cold and wet/warm conditions related to the assembly processes of high-elevation stream invertebrates and bacteria, reinforcing the role of climate in mountainous aquatic biodiversity patterns.

We found more ecologically unique lakes (i.e., lakes with high LCBD values) in areas with a history of low human impact. Indeed, human impacts can have a homogenizing effect on aquatic communities (Olden et al., 2004). In our case, low LCBD values in lakes with higher human impact may be a result of, for instance, alteration of communities by hydrological modifications and external impacts (e.g., cattle grazing) in the Andean lakes since prehistoric times (Sarmiento, 2002; Van Colen et al., 2018). Although we did not observe statistically significant geographic patterns in the distribution of LCBD when the p-values were corrected for multiple testing, lakes located in the central-south Chilean Andes displayed higher LCBDs (without correcting for multiple testing). A correlation between LCBD and latitude also gave the significant negative relationships (Spearman rho = -0.47, p < 0.05), suggesting a decreasing latitudinal gradient of LCBDs. This may be partially associated with the onset and legacy of historical occupation in the continent of the southern portions of our latitudinal gradient (Gayo et al., 2015; Goldberg et al., 2016).

Partitioning LCBD into replacement and richness components provided further insights into the mechanisms underlying changes in spatial beta diversity of diatoms in South America. Most significant environmental correlates for each of the LCBD indices, i.e., LCBDrepl and LCBDrich, were fundamentally different from the ones observed for total LCBD. For instance, LCBD components are responding inversely to the same water chemistry correlates (e.g., conductivity), supporting previous research on lake communities analyzing turnover and nestedness components (Angeler, 2013). Other studies also highlight that finding consistent predictor variables among beta diversity components is challenging in lentic systems in general, and on tropical aquatic communities in particular, including diatom and invertebrate communities (Jyrkänkallio-Mikkola et al., 2018; Castro et al., 2019). Nonetheless, when the replacement component of LCBD dominates, a regional approach focusing on multiple sites might be needed to conserve ecologically unique diatom metacommunities (Wright and Reeves, 1992). In our study, regional environmental variables for conservation purposes are terrain ruggedness and soil variability. In contrast, a dominating richness component of LCBD suggests the need to focus on a few species-rich lakes and local limnological correlates (Ramos-Jiliberto et al., 2009). In this context, conductivity could be a variable to ensure conservation of diatom metacommunities in certain tropical Andean lakes that are naturally salty such as Southern Altiplano regions (Figure 1 and Supplementary Figure 1).

We found that ecological uniqueness in terms of replacement and richness gradients responded to climate (i.e., seasonality in temperature) and geology (soil variability, terrain ruggedness). The influences of diverse climatic conditions and geodiversity on freshwater biodiversity have recently gained attention (Kärnä et al., 2018; Toivanen et al., 2019). In a study on boreal stream and lake diatoms, Vilmi et al. (2017) found a strong association between LCBD and bed rock, soil, and ecoregion characteristics. Studies relating geodiversity-biodiversity in freshwaters suggest an incipient tight coupling between regional catchment characteristics and local biological dynamics, and are in line with studies from high-elevation lakes (Zaharescu et al., 2016). In the Andes of Ecuador, northern páramo lakes differ in SO₄^2– content compared with their southern páramo lake counterparts, which have much higher Ca²⁺ concentrations (Luethje, 2020). Interestingly, Andean lakes located in high-elevation rugged basins harbor ecologically unique diatom communities in terms of richness. We found a negative relationship between LCBD and species richness (Spearman rho = -0.58, p < 0.01), indicating that lakes with exceptional ecological uniqueness are usually the ones with lower numbers of species. Similar negative relationships have been reported in other contexts as well (Legendre and De Cáceres, 2013; da Silva and Hernández, 2014; Mimouni et al., 2015; Heino et al., 2017). High terrain rugosity promotes lake isolation from the surrounding landscape (Valencia et al., 2016), which may result in more dispersal-limited conditions, even for organisms with high dispersal capabilities, like diatoms (Kristiansen, 1996; Benito et al., 2018b). From a biogeographical perspective, identifying topographically diverse mountain regions that harbor ecologically unique lakes may complement research on evolutionary processes, such as diatom endemism (Spanbauer et al., 2018) and climatic microrefugia (de Novaes Nascimento et al., 2019).

Diatom community structure differs among lake habitats, thereby highlighting the relevance of species sorting due to substrate type (e.g., mud, plants, rocks), and dispersal between different lake habitats (pelagic versus benthic communities) as a result of the fluid aquatic environment (Wetzel et al., 2012; Cantonati and Lowe, 2014). We can assume that each study region operates as a metacommunity (Benito et al., 2018b) and hence hypothesize on mechanisms driving LCBD patterns by considering the variability of lake diatom habitats (Figure 1B). We suggest that it is the diversity from the peripheral (benthic) communities that eventually determines between-lake diversity of diatoms and ultimately the ecological uniqueness of the lakes compared to other sites in the region. For instance, in the Peruvian Andes (Cusco and Wet Puna regions), the presence of heterogenous benthic diatom groups (saline, epiphytic, benthic) may account for the high LCBD values (Supplementary Figure 1). In contrast, the homogenous pattern in terms of lake habitat diatoms in the Sud Lipez and Desaguadero regions of Bolivia could arise because of limited opportunities for dispersal. These Bolivian regions are cold and arid, with a low density and small number of suitable aquatic environments, most of which are shallow hypersaline lakes and wetlands (Servant-Vildary and Roux, 1990). In the Chilean Andean lakes, the higher local contributions to beta diversity could be a function of their relatively poor pelagic diatom community compared with richer-than-average planktic dominated diatom regions, such as Lake Titicaca or the Ecuadorian Andes. Diatom diversity often increases in deeper lake zones, because benthic diatoms are transported from other lake habitats and mixed with pelagic taxa associated with seasonal changes (Pla-Rabés and Catalan, 2018). However, we cannot discard the possibility that the observed relationships between LCBD values and ecological groups in each diatom metacommunity were influenced by the temporal variability in our modern database, as the data do not correspond to the same time point among regions. More research is needed to unveil the effect of benthic area on diatom communities’ structure in deep tropical lakes.

The term metacommunity can also be used to define the diatom community of the whole lake for each sedimentary sequence (Leibold et al., 2004). Sediment samples integrate the local species richness and the beta diversity (replacement) of the lake habitats and the variability in composition among them (Pla-Rabés and Catalan, 2018). Thus, local contributions to beta diversity may differ over time in response to lake habitat changes driven by limnological change. We observed a consistent pattern of more variability in LCBDrich than LCBDrepl time series. This is expected given the role of sediment samples acting as a sink by accumulating different entities (species) from other lake habitats (i.e., source) (Logue et al., 2011). Three lakes in the Andes of Ecuador (Piñan, Yahuarcocha and Fondococha) experienced a coincident peak in LCBDrich, but less in LCBrepl, at ca 1000 cal years BP, likely responding to dry/warm conditions centered around the Medieval Climate Anomaly (MCA). The MCA triggered lower lake levels based on many tropical Andean paleolimnological records (Figure 3; Lüning et al., 2019 and references therein). Despite different conditions in water chemistry and lake depths, our findings suggest that the relatively high synchronous compositional uniqueness at that time may be a fingerprint of regional-scale limnological variation in these three lakes. This is partly supported by the fact that high modern LCBDrich values are explained by low temperature seasonality and conductivity (Figure 2). In this application, however, methodological issues, such as time-averaging processes and the partial representation of the entire population abundance, may introduce bias in beta diversity estimates from paleolimnological assemblages (Birks et al., 2016). To the best of our knowledge, our study is the first to investigate the ecological uniqueness in terms of richness and replacement using sediment diatom assemblages, which hampers comparisons with similar works. Further research in other tropical Andean lakes with available contemporary time series data is necessary to assess the generality of this finding.

Our GAM time series models (HGAMs) further delineate temporal beta diversity patterns in tropical Andean lakes of varied size, limnology, and climatic conditions that can provide a regional, long-term perspective of biodiversity changes (Dornelas et al., 2014; McGill et al., 2015). A decreasing trend in beta diversity over time (biotic homogenization) has been assumed to be a result of increased human impacts (Olden and Rooney, 2006). In our case, we found signs of long-term biotic homogenization of diatom assemblages, as measured by decreasing trends of beta replacement over the past ca. 200 years across the four lakes investigated from Ecuador (Figure 4). Nonetheless, beta replacement values fluctuate around a long-term mean for the whole time series, and no periods with substantial increases or decreases in beta diversity arose. In contrast, biotic differentiation (i.e., increase of beta replacement over time) was found in the two deep freshwater lakes of the Altiplano. Admittedly, the low sample density for the most recent period of these two lakes is not captured well by the GAM models (as in Piñan before ca 900 cal. years BP), resulting in wide confidence intervals, like with any other smooth regression approach (Simpson, 2018).

A high context dependency exists among studies that investigate biotic homogenization and its explanatory factors at varied spatial and temporal scales in aquatic ecosystems. In the case of lake sediment diatoms, Winegardner et al. (2017) did not find patterns of biotic homogenization across the conterminous United States between ca. 150 years ago and modern times. In a study on tropical reservoirs affected by eutrophication, Wengrat et al. (2018) found a decreasing trend of spatial beta diversity over the past 100 years. Eutrophication-driven homogenization was also reported by Salgado et al. (2018) using macrophyte paleoecological assemblages. These observations highlight the usefulness of the HGAM models used here for detecting temporal beta diversity trends across space: as this approach does not assume any specific dynamics in the time series, it is possible to determine if broad-scale environmental change (e.g. climate change) led to uniform diversity patterns across the landscape, or if lake-specific dynamics decouple from the regional signal. The latter could be the case here, indicating no widespread biotic homogenization across the tropical Andes. Additionally, HGAM models allow to circumvent the issue of harmonization among lake samples (e.g., binning), which is problematic given the differences in temporal resolution and length. We suggest that this approach can also be applied to other aquatic systems that may or may not have well-defined boundaries and monitoring time series data but that are subject to strong environmental disturbances (e.g., temporary rivers, wetlands) (Ruhí et al., 2017).

Conclusion

Our approach for investigating ecological uniqueness (i.e., LCBD) has the potential to generate new opportunities to integrate modern ecology and paleolimnology for biodiversity and metacommunity studies. In this context, we emphasize several aspects of our results. First, ecological uniqueness of tropical Andean lakes was linked to local and regional environmental variables and showed an inverse pattern for the replacement and richness components of LCBD. Specifically, mean precipitation and the historical human impact mediated how ecologically unique the lakes were, whereas geodiversity (soil variability, terrain ruggedness), temperature seasonality and conductivity mediated, in an opposite manner, its replacement and richness components. We suggest this finding has different management and conservation measures. For instance, individual lakes characterized by high terrain ruggedness and low conductivity are clear management targets if local diatom richness is a conservation goal; if the goal is to conserve lake-catchment systems within a given spatial context (beta diversity), lakes with low soil variability and low Na⁺ content may be suitable management targets. Second, by identifying diatom richness contributions to beta diversity over time, we were able to observe a shared limnological response to warm/dry climatic changes centered around the MCA driven by decreased temperature seasonality and conductivity. However, diatom sedimentary assemblages should not be interpreted as a unique response to environmental change because, for instance, temperature seasonality might have a stronger effect on planktic communities than on the benthos via thermal structure impacts as a result of warming or drying. If environmental variation mediates the relationship between replacement and richness components of beta diversity, then increasing limnological changes due to climate change and human impacts will likely destabilize long-term metacommunity stability; this could be the case in páramo lakes of Ecuador as seen by higher correlations between both LCBD components over time. Finally, we have provided a broader perspective of aquatic biodiversity change over the Common Era (last 2000 years) with beta replacement trends to test the hypothesis of recent biotic homogenization. Considering diatoms are one of the most sensitive groups of organisms, situated at the base of aquatic food webs, the decreasing trends in temporal beta diversity across lakes in the Ecuadorean Andes, albeit being lake-specific and not unprecedented, may cause unexpected effects in the structure and functioning of these ecological and climatic sensors in vulnerable tropical high-elevation ecosystems.

Data Availability Statement

Diatom and environmental datasets, and R code to perform the analyses are available from GitHub (https://github.com/xbenitogranell/LCBD-diatoms).

Author Contributions

XB and AV conceived the idea. XB constructed the database, performed the statistical analysis, and led the writing. XB, SF, and MLu conducted the fieldwork in Ecuador. MLu performed the diatom analysis for Piñan and Fondococha, and XB for Yahuarcocha. MC performed the diatom analyses from modern Chilean lakes. All authors contributed to writing of the manuscript.

Funding

XB was supported by the National Socio-Environmental Synthesis Center (SESYNC), under funding received from the US NSF DBI-1639145. MLu and SF acknowledge funding from National Geographic #8672-09 and NSF EAR-1338694. MC was supported by a FONDECYT postdoctoral fellowship. AV was supported by the project FRESHABIT LIFE IP (LIFE14/IPE/FI/023).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank D. R. Engstrom, Science Museum of Minnesota, for analysis and interpretation of the 210Pb profiles from the study lakes. We are enormously indebted to all the contributors of the South American diatom database (https://github.com/xbenitogranell/diatoms-biogeography-southamerica), including M. I. Velez, P. Tapia, M. Steinitz-Kannan, B. Valencia, M. Bush, F. Mayle, and J. P. Bradbury. We are greatly thankful to Paul A. Baker, Kewrin Choez, and Maria Arteaga for fieldwork assistance in Imbabura Province, Pablo Mosquera and Henni Hampel for help with Cajas National Park lakes limnology, and ETAPA EP rangers and Simon Quiroz for assistance during fieldwork. We also thank the Ministerio del Ambiente in Ecuador for permitting to conduct this research (151-2017-DPAA/MA and 008-IC-DPACH-MAE-2017).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2020.00260/full#supplementary-material

Footnotes

1. ^ https://github.com/xbenitogranell/diatoms-biogeography-southamerica

References