Hydroclimatological Patterns and Limnological Characteristics of Unique Wetland Systems on the Argentine High Andean Plateau

: High-elevation wetlands in South America are not well described despite their high sensitivity to human impact and unique biodiversity. We describe the hydroclimatological and limnological characteristics of 21 wetlands on the High Andean Plateau of Argentina, synthesizing information gathered over ten years (2010–2020). We collected physical-chemical, phytoplankton, and zooplankton data and counted ﬂamingos in each wetland. We also conducted an extensive analysis of climatic patterns and hydrological responses since 1985. These wetlands are shallow, with a wide range of salinity (from fresh to brine), mostly alkaline, and are dominated by carbonate and gypsum deposits and sodium-chloride waters. They tend to have high nutrient concentrations. Plankton shows a low species richness and moderate to high dominance of taxa. Flamingos are highly dependent on the presence of Bacillariophyta, which appears to be positively linked to silica and soluble reactive phosphorus availability. Climatic conditions show a strong region-wide increase in average air temperature since the mid-1980s and a decrease in precipitation between 1985–1999 and 2000–2020. These high-elevation wetlands are fundamentally sensitive systems; therefore, having baseline information becomes imperative to understanding the impact of climatic changes and other human perturbations. This work attempts to advance the body of scientiﬁc knowledge of these unique wetland systems. that


Introduction
The dry Andean plateau stretches 1800 km along the backbone of the mountain range from southern Peru to northern Argentina between the eastern and western mountains, varies from 350 to 400 km in width and averages 4000 m above sea level (m.a.s.l) [1,2]. Despite the aridity, the plateau is dotted with wetlands that provide essential resources for human activity and habitats for biodiversity highly adapted to extreme environmental conditions [3,4]. Nevertheless, these wetlands have received less attention than other types of water systems on the South American continent. This is likely linked to the difficulty of obtaining samples, the variety of morphological and hydrological characteristics, the complexity of the hydrogeology, and the abundance and variety of basins and water bodies in this environment [3,5]. and Roux [44], Rejas et al. [45], and Frau et al. [46]. These authors also demonstrated the importance of nutrient limitation, especially silica and ionic composition for high-elevation diatom assemblages.
The high Andean wetlands support a diversity of endemic, resident, and migratory birds, and some of them completely dependent on water [47,48]. These wetlands are one of the main drivers of bird spatial distribution at high elevation and have a crucial role in preserving birdlife [49,50]. They also serve as steppingstones during the migration of long-distance migrants, particularly shorebirds [48,51,52]. Among birds, flamingos (Phoenicopteridae) are iconic species in the landscape, moving from one wetland to another, making alternative and complimentary use of these wetlands for feeding and breeding, effectively connecting these wetlands as a network of habitats [53][54][55]. Flamingos are an important component of the increasing nature tourism industry that has developed in the past decade and serve as flagship species for the conservation of high-elevation watersheds. Native mammals that depend on these wetlands to survive include puma (Puma concolor), Andean cat (Leopardus jacobitus), Andean fox (Lycalopex culpaeus), hairy armadillo (Chaetophractus nationi), and the commercially valuable vicuña (Vicugna vicugna) [56].
Despite all the singularities reported above regarding high-elevation wetlands from the Andean Plateau, human activities in combination with global climate change are driving variations in biological and hydrological attributes to a rate faster than we can measure. Baseline information regarding biological conditions at these wetlands, highly sensitive to anthropogenic disturbances [9] and, thus, ideal sentinels of global change [57,58] are becoming necessary to improve our capacity to properly manage the use of these environments [30,59,60]. Catamarca province, located in northwest Argentina, is an area where human activity, mainly mining has increased recently [61]. In Catamarca, we have gathered information over the last ten years to characterize hydroclimatological and limnological aspects, but we lack a holistic perspective of the characteristics and functioning of these systems. We argue that assessing geochemical and limnological characteristics at specific wetland systems in Catamarca can significantly improve the understanding of potential impacts from current and future human activities in these and similar systems across the High Andean Plateau. High-elevation wetlands are considered one of the most comparable ecosystems across continents [57] due to the many ecological features that they share. This makes the information synthesized in this manuscript potentially applicable to many other environments in South America and other arid regions that have not been adequately studied but share similar issues of human modifications and climatic change. In this study, we summarize the findings of a multidisciplinary approach synthesizing information gathered from publicly available data sources and field expeditions, describing the hydrogeological components and relevant biological patterns of wetlands on the High Andean Plateau in northwest Argentina.

Study Area
The study area encompassed 21 wetlands sampled in northwestern Argentina between 25 • -28 • S and 67 • -69 • W and situated between 3000 and 4500 m above sea level (m.a.s.l) ( Figure 1). The area is characterized by predominantly volcanic geomorphology with numerous stratovolcanoes and vast basaltic and pyroclastic deposits. Eight of the wetlands sampled are within the Lagunas Altoandinas y Puneñas de Catamarca Ramsar site, which was designated in 2009. The climate is cold and dry, with temperatures below 0 • C most of the year and average annual rainfall between 50 and 200 mm year −1 [62]. Most of these wetlands are shallow, commonly less than 1 m in depth, and become partially or entirely frozen during winter [10]. They all sit within endorheic basins where water inflow comes mainly from groundwater recharge in the surrounding mountains and local alluvial infiltration and runoff, all of which leave the system through evapotranspiration. comes mainly from groundwater recharge in the surrounding mountains and local alluvial infiltration and runoff, all of which leave the system through evapotranspiration. The surface water occurrence layer depicts the average occurrence of water over the Landsat record on a scale from pale purple meaning very little occurrence (<10%) to dark blue meaning permanent water (100%). This dataset is publicly available through the European Commission's Joint Research Centre [63].

Physical-Chemical Sampling and Laboratory Analyses
Sampling was conducted between 2010 and 2020 during the summer season (January-February) because some of these wetlands are dry or frozen during the rest of the year, over seven consecutive days in a total of 21 wetlands, including peatlands, lagoons, and salars (see Table S1 in the supplementary material for specific sampling dates at each site). At each sampling site, between one and four replicates were taken for the several hydrogeochemical-limnological parameters explained below. It is important to note that The surface water occurrence layer depicts the average occurrence of water over the Landsat record on a scale from pale purple meaning very little occurrence (<10%) to dark blue meaning permanent water (100%). This dataset is publicly available through the European Commission's Joint Research Centre [63].

Physical-Chemical Sampling and Laboratory Analyses
Sampling was conducted between 2010 and 2020 during the summer season (January-February) because some of these wetlands are dry or frozen during the rest of the year, over seven consecutive days in a total of 21 wetlands, including peatlands, lagoons, and salars (see Table S1 in the Supplementary Materials for specific sampling dates at each site). At each sampling site, between one and four replicates were taken for the several hydrogeochemical-limnological parameters explained below. It is important to note that our assessment included 21 wetlands in total; however, no biological or limnological data was collected at Laguna Diamante, meaning there are only 20 wetlands included in those results. For limnological analyses, measurements made in situ were temperature ( • C), pH, and conductivity (mS cm −1 ) using HANNA multiparametric meters; transparency with a Secchi disc (m); and elevation (m.a.s.l) using a GPS altimeter and verified with a digital elevation model in ArcGIS. Water-type categories were derived from the fresh-brackish-saline-brine categorization scheme described by Warren [64].
Water samples were taken for determination of the nutrients nitrate+nitrite (N-NO 3 − +N-NO 2 − ), ammonium (N-NH 4 + ), soluble reactive phosphorus (SRP), and dissolved silica (Si). Samples for determination of N forms and SRP were immediately filtered through membrane filters (0.45 µm pore size). Subsequently, samples for nitrate + nitrite and ammonium determination were acidified to a pH < 2 with concentrated sulfuric acid, whereas samples for SRP determination were preserved by adding 5 mg L −1 of HgCl 2 . All samples were stored in cold and darkness before transporting to the laboratory for analysis. Acidified water samples were neutralized in the laboratory. Samples used for determination of N-NO 3 − +N-NO 2 − were estimated by reducing N-NO 3 − with hydrazine sulfate and subsequent colorimetric determination of N-NO 2 − [65]. N-NH 4 + was determined using the indophenol blue method, SRP using the ascorbic acid method, and Si using the molibdosilicate method. Dissolved inorganic nitrogen (DIN) was calculated as the sum of N-NO 3 − +N-NO 2 − and N-NH 4 + . The methodology proposed in APHA [66] was followed throughout, and nutrient concentrations were expressed as mg L −1 for silica and µg L −1 for DIN and SRP.
Physical and chemical data were explored by using principal component analysis (PCA). All data were log-transformed, except pH, and values obtained were standardized. Neighbor-joining clustering analysis was performed on the physical-chemical data obtained to find common patterns among wetlands sampled.
Water samples analyzed for major ion composition were collected at 10 of the 21 wetlands in February of 2019 and 2020 using a consistent, standardized procedure and repeated sampling from the same location when possible (Supplementary Materials, Table S2). Water samples were filtered through 0.45 µm filters using a plastic 60 mL syringe and were stored in clean HDPE bottles. The concentration of major elements was analyzed using inductively coupled plasma optical emission spectroscopy for major elements (Agilent 5110 ICP-OES) and high-pressure ion chromatography (Dionex Integrion HPIC) for Cl, SO 4 , PO 4 , and NO 4 anions at the University of Massachusetts Amherst. Waters with relatively high TDS were diluted volumetrically before analysis. For ICP OES analysis, samples were acidified to 1% HNO 3 v/v before analysis. Quantification was performed using seven external calibration standards ranging from 0.1 to 100 ppb. Calibration verification standards and blanks were run at every 10th analysis for anions and elements. Elemental analysis was verified with external NIST standard SRM 1643d, and anion analysis was verified with a secondary anion standard (Anion II Std Dionex). Samples that exceeded the calibration by 120% were diluted and reanalyzed. The concentration of HCO 3 in water was determined by titration using an autotitrator. In all cases, the laboratory measurement error was 0.01 mg L −1 . All other major ion data used in the study were obtained from Sureda [67]. A charge balance assessment of all these data was performed, and only samples with less than 10% error were included.

Hydroclimate Characterization
We utilized multiple remotely sensed datasets to assess hydrological and climatological conditions and changes over time across the basins of interest. These include Landsat satellite imagery (spatial resolution of 30 m × 30 m), the TerraClimate climatically aided precipitation interpolation, and the famine early warning systems network land data assimilation system (FLDAS). These data products are publicly available on Google Earth Engine (GEE) platform and were extracted and compiled using tools within the platform. For each wetland system, basin domains were extracted from the HydroSHEDS dataset (Hydrological data and maps based on shuttle elevation derivatives at multiple scales) [68] or a contributing watershed area determined using the Watershed tool in ArcGIS Pro.
Values of total precipitation and mean air temperature were extracted for each satellite pass (~16-day intervals) as spatially aggregated values for each basin or watershed containing the wetland of interest; outlines of these areas are shown in Figure 1. Total living or "green" vegetation extent was also assessed at each wetland complex using the normalized difference vegetation index (NDVI). These data provide a robust assessment of the change in climatic and hydrological conditions from the mid-1980s to the present.
To assess changes in precipitation, we utilized the TerraClimate precipitation product which uses climatically aided interpolation, combining high-spatial-resolution climatological normals from the WorldClim dataset, with coarser spatial resolution but time-varying data from CRU Ts4.0 and the Japanese 55-year reanalysis (JRA55) [69]. The data are provided as monthly precipitation accumulation totals at a resolution of 1/24 • (~4-km) from 1958 to 2021. A time series of annual total precipitation was then created for each basin/watershed. From this time series, we derived a general value of precipitation changes between 1985 and 2020 represented by the percent change from the 1985-1999 mean value to the mean value of the 2000-2020 period.
Mean annual air temperature was determined using the famine early warning systems network (FEWS NET) land data assimilation system (FLDAS) global land surface simulation that provides monthly values of near-surface air temperature at a 0.1 • × 0.1 • resolution [70]. The model is forced by combining the modern-era retrospective analysis for research and applications version 2 (MERRA-2) data and climate hazards group infrared precipitation with station (CHIRPS) 6 hourly rainfall data that have been downscaled using the NASA Land Data Toolkit. Although the precipitation product used in FLDAS is not as reliable for this region, the air temperature data are quite consistent with station data and other datasets and are, therefore, determined to be robust estimates. The value of temperate change from 1985 to 2020 was derived from the linear trend of mean annual air temperature over that period.
To assess changes in vegetation cover, we utilized the NDVI, which is calculated from spectral imagery using the formula: NDVI = (NIR − RED)/(NIR + RED) where NIR is the reflection in the near-infrared spectrum and RED is the reflection in the red range of the spectrum [71]. NDVI is designed to assess the density of vegetation at a given location at a given time. It provides a single band with a range from −1.0 to 1.0 where negative values are clouds, water, and snow; values close to 0 up to 0.1 are from rocks and bare soil. Different types of vegetation are classified in values greater than~0.1. To be conservative and be sure we are capturing only vegetation, we extracted the pixels whose values fell between 0.2 and 0.9. The domains over which each wetland was analyzed were chosen to encompass all vegetation directly associated with the wetland system, at or near each basin floor. Using GEE, we extracted a full series of Landsat 5 and Landsat 7 Satellite images from 1985 to 2020 and determined the number of pixels covered by green vegetation (0.2-0.9 range) from which a total geographic area was then calculated. This provides a time series of the total area covered by living vegetation within each wetland region. From these time series, we derived a general value of change in the vegetated area between 1985 and 2020 represented by the percent change from the 1985-1999 mean value to the mean value of the 2000-2020 period.

Biological Sampling and Analyses
Samples of phytoplankton were collected at the same time as the physical-chemical samples from the subsuperficial area using 120 mL bottles and then were immediately fixed with 1% acidified Lugol solution for preservation. Quantitative sample analyses were carried out following the Utermöhl method [72], and the density obtained was expressed ind mL −1 by accepting a counting error <20% [73]. For estimating zooplankton abundance, between 10 and 30 L of water were filtered through a conventional conical plankton net (50 µm) with a collector and a 2 L bucket. The collected material was fixed in situ with 10% formalin and dyed with erythrosine. The microzooplankton counts (Rotifera and nauplii) were carried out with a conventional optical microscope using Sedgwick-Rafter type chambers with a capacity of 1 mL. Macrozooplankton (Cladocera and Copepoda) counts were performed in a 5 mL Borogov-type counting chamber. Density was expressed as ind L −1 . Taxonomic identification was carried out to the lowest taxonomic possible level. Phytoplankton taxonomic classifications were based on Lee [74]. Taxonomic determinations were made following keys and specific bibliography of each algal group, such as Krammer and Lange-Bertalot [75] for Bacillariophyta, Tell and Conforti [76] for Euglenophyceae, Komárek and Fott [77] for Chlorophyceae, and Komárek and Anagnostidis [78,79] for Cyanobacteria, among other authors and recent revisions. Zooplankton taxonomic identifications were based on Ahlstrom [80], Koste [81], Kořínek [82,83], Korovchinski [84], and Alekseev [85], among others. Flamingos were identified to species level (Phoenicoparrus andinus, P. jamesi, and Phoenicopterus chilensis) using spotting scopes and binoculars and counted from established census points, consistent across years. For large flocks (>500) or groups dispersed throughout the wetland, numbers were estimated using methods described in Marconi [86]. For this study, we used total numbers per flamingo species for each wetland.
Flamingo and plankton abundances were compared among wetlands by using Kruskal-Wallis analyses and the Mann-Whitney post hoc test. Temporal tendencies for flamingo abundances were explored by using Mann-Kendall tests using flamingo counts carried out from 2010 to 2020. Diversity indices, such as Shannon-Wieber, dominance, and richness, were estimated for plankton. Finally, a redundancy analysis (RDA) was performed by using all the environmental variables as predictors and the biological communities (flamingos, phytoplankton, and zooplankton) as response variables to explore the main environmental variables which may influence natural communities' abundance. The RDA analysis was performed after considering the longitude gradient in a detrended correspondence analysis (DCA), and highly correlated variables (VIF > 10) were omitted.

Physical-Chemical and Hydroclimate Characterization
All the wetlands sampled (n = 21) were located above the 3000 m.a.s.l. and range in area from 24 to 1117 ha (determined by surface water extent from the HydroLAKES dataset described by Messager et al. [87] (Table 1). Diamante is not included in Table 1 since we did not collect biological or limnological data there; however, it is one of the largest wetlands in the region and therefore was included in the rest of our analysis. All of them were shallow (<1 m) with high transparency, allowing light to reach the bottom. Alkaline conditions were typical, except for Aguas Dulces showing to be circumneutral (pH = 7). Nutrient concentrations were variable across wetlands. DIN concentrations were above 100 µg L −1 , with the highest concentrations registered in Las Peñas, Carachi Pampa, Grande, and Hedionda and the lowest in Cortaderas, San Francisco, and Archibarca. SRP concentrations followed a similar pattern of high variation among lakes, with the highest concentrations registered in Archibarca and Hedionda. Values registered were above 50 µg L −1 , except in Antofagasta, Azul, and Negra (32.7, 44.1, and 46.1 µg L −1 , respectively). The DIN:SRP ratio showed, in general, low values (<10), indicating DIN limitation in almost all wetlands where information was available. Si concentrations were about one thousand times higher than concentrations of the other inorganic nutrients sampled (>30 mg L −1 for all lakes sampled). Conductivity showed a wide range, with some wetlands having very low values (e.g., Alumbrera, Antofagasta, and Cortaderas had <2.5 mS cm −1 ). About 30% of the wetlands sampled were categorized as freshwater or brackish, all measuring below 10 mS cm −1 . The remaining wetlands (70%) were considered saline, with mean conductivity values between 25 and 110 (mS cm −1 ) being the most saline; Diamante was a large lake composed of brine. (Table 1). Table 1. Median environmental values for each wetland sampled between 2010 and 2020. Elevation (Alt), area (Area), temperature (Temp), pH (pH), water level (WL), Secchi disc (SD), conductivity (Cond), dissolved inorganic nitrogen (DIN), silica (Silica), soluble reactive phosphorus (SRP). ND = no data.

Aguas Dulces
Saline 3325  The PCA ordination analysis explained 49.66% of the total variance in the first two axes. The first axis indicated a gradient of DIN concentration, with maximum values recorded in Carachi Pampa, Grande, and Las Peñas and the lowest value in Cortaderas. The second axis suggested a gradient of SRP concentration and elevation from San Francisco and Archibarca to Alumbrera and Antofagasta. Carachi Pampa, Negra, and Grande stand out for their high conductivity values, Si concentration, and area ( Figure 2). Complementarily, the neighbor-joining clustering analysis suggests four clusters of wetlands. Cluster one (Archibarca, Cabi, Caro, Hedionda, El Peinado, and Las Peñas) had the highest nutrient values, being mostly saline lagoons. Cluster two (Aguas Dulces, Antofagasta, Alumbrera, Carachi Pampa, and Purulla) were those located at the lowest elevation (above 3000 m.a.s.l), had mean temperature values higher than the other wetlands (above 22 • C), and were salinefreshwater lakes. Cluster three (Grande, Verde, Negra, Tres Quebradas, and San Francisco) had the highest elevation, area, and Si concentrations and were saline-brackish. Finally, cluster four (Los Aparejos, Azul, Cortaderas, and Las Tunas) was characterized by having the lowest mean DIN and Si concentration, as was mostly saline-brackish ( Figure 3).
Based on the dissolved ion compositions of 14 wetland surface waters (10 analyzed for this work and 4 from data obtained from Sureda [65], we can describe three distinct groups of wetlands ( Figure 4). Group one: Grande, Diamante, Carachi Pampa, Caro, and Cabi are dominated by Na + , K + , and Cl − with very low concentrations of Ca ++ , Mg ++ , and SO 4 − . Group two: Purulla, El Peinado, Negra, Verde, and San Francisco are also marked by high relative concentrations of Cl − but have higher Ca ++ and Mg ++ concentrations than the first group. The third group of wetlands is marked by substantially higher relative concentrations of HCO 3 and SO 4 − , and include Alumbrera, Antofagasta, Azul, and Los Aparejos. Three of these are the freshest wetlands sampled (<5 mS cm −1 ), and the other is much saltier (~75 mS cm −1 ). Alumbrera and Antofagasta, which comprise the extensive wetland complex at the terminus of the large Punilla River, show very similar characteristics, being dominated by relatively high HCO 3 and SO 4 − but lower relative Cl − , Ca ++ , and Mg ++ concentrations. Azul is like these lagoons but contains higher SO 4 − and much higher Mg ++ . Los Aparejos appears to be more like group two than the fresher lagoons mentioned above but substantially lower in SO 4 − .
wetland complex at the terminus of the large Punilla River, show very similar characteristics, being dominated by relatively high HCO3 and SO4 −− but lower relative Cl − , Ca ++ , and Mg ++ concentrations. Azul is like these lagoons but contains higher SO4 − and much higher Mg ++ . Los Aparejos appears to be more like group two than the fresher lagoons mentioned above but substantially lower in SO4 − .

Figure 2.
Principal components analysis (PCA) considering only those wetlands where full environmental data were available. Abbreviations of variable names as in Table 1.  Table 1.   Our analysis of climatic and hydrological conditions in this region revealed 14 distinct basins or watersheds encompassing the 21 wetlands in this study ( Figure 5). Decadalscale changes in climatic and hydrological conditions within these large internally drained basins provide an effective estimate of the natural drivers of recent hydrological changes in these wetland systems. In the case of Azul, Cortaderas, and Aguas Dulces, the domains used are subwatersheds of larger basins, and in the case of Cabi, a small wetland on the edge of the Grande Basin, its watershed area is too small to assess independently; therefore, we use the Grande basin as a representative of the conditions there. Temperatures since 1985 have increased substantially across the entire region. Every basin has recorded a temperature increase of at least 0.45 °C. Basins in the eastern and northern parts of the study area have seen an increase of between 0.55 and 0.65 °C with the largest increase observed in Archibarca and Caro (Table 2). Precipitation has broadly decreased across all wetland areas, with the smallest decreases observed in the northeastern region of between 9% and 12%. The most significant decreases are in the southern and westernmost regions of between 14% and 22%. Since a major region-wide drought has been identified since the mid-2000s, these changes are likely in large part a reflection of the relative depth of this drought across the study area. In terms of hydrological changes, in this case, represented by the change in vegetation cover, there is a wide distribution of relative changes. The western and southern wetlands (El Peinado, Tres Quebradas, Negra, Azul, and Las Tunas) Our analysis of climatic and hydrological conditions in this region revealed 14 distinct basins or watersheds encompassing the 21 wetlands in this study ( Figure 5). Decadal-scale changes in climatic and hydrological conditions within these large internally drained basins provide an effective estimate of the natural drivers of recent hydrological changes in these wetland systems. In the case of Azul, Cortaderas, and Aguas Dulces, the domains used are subwatersheds of larger basins, and in the case of Cabi, a small wetland on the edge of the Grande Basin, its watershed area is too small to assess independently; therefore, we use the Grande basin as a representative of the conditions there. Temperatures since 1985 have increased substantially across the entire region. Every basin has recorded a temperature increase of at least 0.45 • C. Basins in the eastern and northern parts of the study area have seen an increase of between 0.55 and 0.65 • C with the largest increase observed in Archibarca and Caro (Table 2). Precipitation has broadly decreased across all wetland areas, with the smallest decreases observed in the northeastern region of between 9% and 12%. The most significant decreases are in the southern and westernmost regions of between 14% and 22%. Since a major region-wide drought has been identified since the mid-2000s, these changes are likely in large part a reflection of the relative depth of this drought across the study area. In terms of hydrological changes, in this case, represented by the change in vegetation cover, there is a wide distribution of relative changes. The western and southern wetlands (El Peinado, Tres Quebradas, Negra, Azul, and Las Tunas) as well as the highest elevation wetland, Diamante show large increases in total vegetation cover since the 1980s and 1990s with a range of between 12% and 53%, with the highest values in the Tres Quebrada basin. In contrast to this, many of the wetlands show a substantial decrease in vegetation cover from −10% to −31%. Particularly notable decreases were observed in the northern and eastern wetlands, where Archibarca, Antofagasta, Hedionda, Grande, and Cabi as well as the central basins San Francisco and Aguas Dulces all registered decreases in the vegetated area of >15%.

Biological Communities' Characterization
Three flamingo species, Phoenicoparrus andinus, P. jamesi, and Phoenicopterus chilensis, were recorded in the study area during the sampling periods; however, their abundances were variable among wetlands and across years (Table 3). Flamingo abundance for all

Biological Communities' Characterization
Three flamingo species, Phoenicoparrus andinus, P. jamesi, and Phoenicopterus chilensis, were recorded in the study area during the sampling periods; however, their abundances were variable among wetlands and across years (Table 3). Flamingo abundance for all species among years was variable in some wetlands (e.g., San Francisco, Las Peñas, Aguas Dulces, and Archibarca) while appearing more stable in others (Carachi Pampa, Grande). Phoenicoparrus jamesi was the most abundant in Grande, where numbers were >12,000 in all years. This species was the most abundant in Los Aparejos, Hedionda, and Las Peñas. For P. andinus, numbers fluctuate less across years in San Francisco, Purulla, Negra, and Archibarca. P. chilensis numbers were very low in the study area considering a global population of >500,000 individuals. Las Peñas wetland consistently had P. chilensis present.
The Kruskal-Wallis test showed statistically significant differences among wetlands for the three species (H = 34.66 p < 0.001). These differences were for P. jamesi vs. P. andinus and P. chilensis and between P. andinus vs. P. chilensis (Mann-Whitney test, p < 0.001). When the temporal tendencies of the abundances between 2010 and 2020 were explored with the Mann-Kendall test, a statistically decreasing trend was detected for P. chilensis (Z = 3.42 p < 0.0001), notable since the year 2015. No tendency was seen for the other two flamingo species (p > 0.05) (Supplementary Material, Figure S1. Spearman correlation among all measured environmental variables and P. jamesi (n = 13) in wetlands with high variability among years (Alumbrera, Azul, Los Aparejos, and San Francisco) revealed a positive statistically significant correlation with Bacillariophyta (Rho = 0.71 p = 0.006), which showed a negative correlation with water level (Rho = −0.72 p = 0.02). A statistically significant positive correlation was found between P. andinus abundance and conductivity (n = 11, Rho = 0.87 p = 0.01) which showed a marginally significant correlation with the area of the wetlands (Rho = 0.53 p = 0.06). For phytoplankton species, richness ranged between 3 and 20 species recorded for each wetland. Bacillariophyta species, especially from the genera Navicula, and Nitzschia, were the most frequently recorded. For Chlorophyceae, Monoraphidium and Chlamydomonas were the best-represented genera. For the other groups of algae, Cryptomonas (Cryptophyceae), Lepocinclis, Euglena, Trachelomonas (Euglenophyceae), and Synechocystis, Oscillatoria, Phormidium (Cyanobacteria) were the most representative. The assemblage tends to be dominated by Bacillariophyta, followed by Chlorophyceae and occasionally by Cyanobacteria. Density registered for phytoplankton was highly variable across wetlands with the highest values recorded in Grande, Cabi, and Las Peñas (maximum median values recorded = 18,635 ind mL −1 ) and the lowest in Azul, El Peinado, and Tres Quebradas, with 107 ind mL −1 on average (Figure 6a). The high abundances found were attributed to Bacillariophyta which was the dominant group solely in Grande. The Shannon diversity index showed low values for these wetlands, always <3, and variable dominance values among wetlands (above 0.6 and 1). Alumbrera, Antofagasta, and Aguas Dulces reached the highest values and Hedionda and Las Tunas the lowest (Figure 6b). The Kruskal-Wallis analysis showed differences in total phytoplankton density among lakes (H = 30.49 p = 0.04), showing statistically significant differences between Grande and Antofagasta, Aparejos, and Azul (p < 0.05 in all cases).
index showed low values for these wetlands, always <3, and variable dominance values among wetlands (above 0.6 and 1). Alumbrera, Antofagasta, and Aguas Dulces reached the highest values and Hedionda and Las Tunas the lowest (Figure 6b). The Kruskal-Wallis analysis showed differences in total phytoplankton density among lakes (H = 30.49 p = 0.04), showing statistically significant differences between Grande and Antofagasta, Aparejos, and Azul (p < 0.05 in all cases).  Zooplankton showed a similar pattern to phytoplankton, with between one and five species recorded for each lake and assemblages alternatively dominated by Copepoda, or Rotifera. Alona (Cladocera), Lecane, Keratella, Brachionus, and Boeckella (Copepoda) appeared as the most frequent genera. The density recorded for the different wetlands was, in general, low (median values of 10 ind L −1 ), with Archibarca and Las Peñas being the exceptions (321 and 176 ind L −1 , respectively) (Figure 7a). The Shannon diversity index showed low values (<2) and highly variable values of dominance for the lakes sampled (between 0.2 and 1). Antofagasta, Alumbrera, and Azul had the highest richness values while Verde had the lowest (Figure 7b). The Kruskal-Wallis analysis showed no differences in density among wetlands (H = 20.37 p = 0.11).
Rotifera. Alona (Cladocera), Lecane, Keratella, Brachionus, and Boeckella (Copepoda) appeared as the most frequent genera. The density recorded for the different wetlands was, in general, low (median values of 10 ind L −1 ), with Archibarca and Las Peñas being the exceptions (321 and 176 ind L −1 , respectively) (Figure 7a). The Shannon diversity index showed low values (<2) and highly variable values of dominance for the lakes sampled (between 0.2 and 1). Antofagasta, Alumbrera, and Azul had the highest richness values while Verde had the lowest (Figure 7b). The Kruskal-Wallis analysis showed no differences in density among wetlands (H = 20.37 p = 0.11).  The RDA analysis explained 29.85% of the total variance (F = 1.8 p = 0.01), with the first two axes retaining 71.5% of the total variation. Both P. andinus and P. jamesi were associated with the highest elevation and area wetlands with higher abundances of Bacillariophyta, particularly in Grande, Los Aparejos, Las Peñas, and Negra. Bacillariophyta were positively correlated with higher concentrations of SRP and Si availability. Conversely, P. chilensis appeared associated with Copepoda and Cladocera (zooplankton groups), which were also more associated with other algae groups, such as Cyanobacteria and Chlorophyceae. Euglenophyceae and Cryptophyceae were correlated with higher water temperatures, higher conductivity values, and Secchi transparency while Chlorophyceae tend to be restricted to low-conductivity wetlands and low nutrient concentrations. Rotifera was a group poorly represented in the ordination (Figure 8).
The RDA analysis explained 29.85% of the total variance (F = 1.8 p = 0.01), with the first two axes retaining 71.5% of the total variation. Both P. andinus and P. jamesi were associated with the highest elevation and area wetlands with higher abundances of Bacillariophyta, particularly in Grande, Los Aparejos, Las Peñas, and Negra. Bacillariophyta were positively correlated with higher concentrations of SRP and Si availability. Conversely, P. chilensis appeared associated with Copepoda and Cladocera (zooplankton groups), which were also more associated with other algae groups, such as Cyanobacteria and Chlorophyceae. Euglenophyceae and Cryptophyceae were correlated with higher water temperatures, higher conductivity values, and Secchi transparency while Chlorophyceae tend to be restricted to low-conductivity wetlands and low nutrient concentrations. Rotifera was a group poorly represented in the ordination (Figure 8).  Table 1.

Discussion
Across the Altiplano, most aquatic ecosystems remain poorly described in detail, limiting our ability to understand their ecological importance and sensitivity to human impact. Here, we describe several trophic levels (phytoplankton, zooplankton, and flamingos) along with several limnological and hydroclimatological features for a set of wetlands of a poorly studied region in the High Andean Plateau of Argentina undergoing rapid change from anthropogenic drivers such as climate change and industrial mining development.  Table 1.

Discussion
Across the Altiplano, most aquatic ecosystems remain poorly described in detail, limiting our ability to understand their ecological importance and sensitivity to human impact. Here, we describe several trophic levels (phytoplankton, zooplankton, and flamingos) along with several limnological and hydroclimatological features for a set of wetlands of a poorly studied region in the High Andean Plateau of Argentina undergoing rapid change from anthropogenic drivers such as climate change and industrial mining development.

Geological, Hydrological, and Climatological Features
The ecosystems sampled were highly variable in their surface area, the largest being often co-located with salars (evaporite deposits). Most of them (71%) were located at very high elevations, above 3900 m.a.s.l. Owing to the hyperarid climate conditions on many of the floors of these basins, most of the wetlands (67%) were composed of saline or even brine surface water bodies. These conditions are attributable to the rapid uplift of the plateau and progressive orographic-induced aridification, and its location along the Tropic of Capricorn subtropical high belt [88].
The registered nutrient concentrations also appear highly variable among wetlands sampled, with some lagoons such as Archibarca, Cabi, Caro, Hedionda, Peinado, and Las Peñas registering the highest mean values of DIN and SRP. In a similar system of wetlands in the Atacama Desert (Chile), higher nitrate concentrations were attributed to local windblown salts containing nitrate formed by sun UV ray oxidation of nitrogen in the atmosphere or to dust from other nitrate deposits [89]. In some wetlands, such as Archibarca, Cortaderas, and San Francisco, DIN was a limiting nutrient for phytoplankton with values below 100 ug L −1 , while SRP always registered values above 10 ug L −1 , suggesting it was not limiting [90]. In addition, the DIN/SRP ratio (always < 10) suggested that DIN was the limiting nutrient for those wetlands where data were available during samplings. These results are consistent with those reported for other saline systems around the world e.g., [91][92][93][94]. The high variability observed could also be related to fluctuation in hydrological conditions that induce changes in nutrient inputs; therefore, the variations in hydrological conditions are a critical control in semiarid to hyperarid areas, where surface water levels fluctuate seasonally and interannually [95,96]. The wetlands surveyed in our study tend to be very shallow (<1 m in depth), but they experience high variations in water surface area across seasons and years. Decreases in water volume through evaporation result in increases in salinity and nutrient concentrations due to the accumulation of ions [4].
The high values of Si found in most wetlands, and particularly in Grande, Cortaderas, Las Peñas, San Francisco, and Hedionda, are consistent with the tectonically active region and vast exposures of young volcanic rocks [97,98].
The variability in salinity across wetlands is quite dramatic, even between nearby wetlands. Alumbrera and Antofagasta, the freshest of the lagoons sampled, can be explained because both are fed by the Punilla River, a large river through which recent rainwater and snowmelt are efficiently transported toward the basin floor. However, the other wetlands tend to be mainly fed by groundwater recharge and are rich in bicarbonate, sodium, and potassium, which tend to concentrate and precipitate when they reach the surface through springs (here known as vegas and bofedales). This process is fundamental to the creation and maintenance of salars such as Carachi Pampa or Tres Quebradas, where the formation of thick evaporite sequences primarily of halite (NaCl) and dense brine groundwaters require the persistent concentration of salty water close to the surface e.g., [19,99]. For the other wetlands, such as Negra, Grande, or Purulla, variations in salinity of surface waters can be attributed to several processes, including dissolution of old salts, thermal water contributions, evaporative concentration, and mixing between old saline groundwater, and less saline recently recharged or direct meteoric water [89]. For almost all the wetlands sampled, we found high water transparency, with light reaching the bottom of the water body during sampling, especially in Cortaderas and San Francisco. This pattern emphasizes the shallow water column of these ecosystems (<1 m), which contrasts with lakes from areas of the world at similar high elevations, often of glacial origin [100,101].
Our analysis of major ion compositions in these waters provides valuable insight into the distinctions between these wetland systems based on the interactions of inflow waters (meteoric, streamflow, and groundwater), evapoconcentration, and mineral precipitation. The largest group is dominated by Na + , K + , and Cl − reflecting highly geochemically evolved systems. The surface water in these wetlands is sourced from inflow waters that have undergone substantial evapoconcentration along their flow paths and in the water body itself. This group contains two of the saltiest surface waters measured in this work (100 and 230 mS cm −1 ). The second group is also quite evolved but contains higher proportions of Ca ++ and Mg ++ , indicating a distinct geochemical pathway of inflow waters. The other wetlands display some unique geochemical characteristics but can be described primarily as less evolved versions of the other two groups. Alumbrera and Antofagasta, fresh wetlands upgradient of the massive Carchi Pampa basin floor appear to be diluted versions of the waters that feed Carachi Pampa. The primary difference between these systems is that these fresh, upgradient wetlands contain substantial HCO 3 − and SO4 − . In these hyperarid to arid evaporite-forming systems, as waters progressively evaporate under intense insolation, first carbonate, then gypsum precipitate leaving very little of these ions in the remaining water [64]. At least from a geochemical perspective, these systems represent two end-member wetland systems but appear to be of highly similar source waters and flow paths. A similar case can be made for Azul, Los Aparejos, and Negra, which reside in adjacent basins, where Los Aparejos represents an intermediate composition but appears to have similar hydrogeological characteristics. This is consistent with basins on the Chilean Altiplano and the Salar de Atacama where dense brines can be explained through the intense and progressive evapoconcentration of much fresher inflow waters near the basin floors either as groundwater exiting the ground from long flow paths or ephemeral flow from precipitation events [21,22,102].
Several important observations can be drawn from our hydroclimate analysis. First, there is a substantial temperature increase at every basin over the last several decades, which is consistent with the effects of global climate change predictions for this region [30]. Increasing temperatures have likely increased overall potential evaporation and transpiration in these wetland ecosystems, perhaps contributing to some of the changes seen in vegetation cover over the same period. Basins with the highest temperature increases also show the strongest declines in vegetation cover, and wetlands with the largest increases in vegetation cover generally show lower temperature increases. Previous studies examining carbon exchange in peatlands show that soil temperature and moisture conditions are well coupled to carbon losses (plant respiration and soil decomposition) and carbon uptake (plant productivity) at an ecosystem level [103,104]. It has also been recently suggested that these ecosystems play an important role in controlling microclimates [11,12]. Though quantifying this impact requires further study, this may be a critical tangible and rapid effect that climate change is having on these wetlands. Changing precipitation patterns have likely come along with increasing temperature. Recent research has shown that this may be currently manifesting in South America by shifting the seasonal monsoon southward, leading to an increase in frequency and magnitude of significant precipitation events on the plateau while also increasing interannual variability [105][106][107], resulting in increases in frequency and magnitude of large precipitation events but also more pronounced drought periods. Indeed, a major drought has been documented in the region since the mid-2000s and is likely ongoing [108,109]. It is important to incorporate these two climate-driven phenomena into any hydrological or ecological assessment of the region. The changes in basin-scale hydroclimate we describe here can provide a valuable insight into the characteristics of response to natural conditions in each unique system and by extension, current impacts from or potential vulnerability to non-natural disturbances.
The strong decrease in precipitation across all these wetlands is likely primarily due to the ongoing effects of the drought. However, a pronounced regional distinction exists in the magnitude of this decrease, with basins in the southwest showing about double the decrease as those in the northeast. This may be due to the relative proximity of these basins to the dominant moisture source (from the northeast) and/or differences in the magnitude of large events within these basins over the last 5-10 years. The distribution of these large events may also be an important factor in determining the changes in vegetation area in these wetlands. There is increasing evidence of the importance of wetlands as regulators of the local water balance [110,111]. However, we still have little evidence of the effects that these changes in vegetation, through their effects on infiltration and water retention, may have on flamingo distribution patterns and plankton structure (through water level changes, ion concentration, or even modifications in nutrients availability).
Recent research in the region has shown that large precipitation events, which generally occur over only a few days can have a major and lasting influence on wetland surface waters and vegetation extent [112]. In the Tres Quebradas basin, the largest increase in vegetation may also be partially explained by the contribution of increased meltwater from the only glaciers in the region at Monte Pissis and Ojos del Salado [113,114]. These changes in local precipitation and ephemeral flow paths (days to months' time frame) are likely a dominant factor in whether a wetland has grown or shrunk over the last few decades, especially in higher elevation systems with smaller contributing areas.
Another pattern that can be garnered from the distribution of these hydroclimatic changes is the potential importance of the buffering effect that large-scale, old groundwater flow can have on these wetlands. An assessment of the isotope tritium in several of these lagoon waters made by Moran et al. [115] shows that most water flowing through these basins and sustaining these wetlands is old water (at least 70 years old). In this arid environment, fluxes of recent precipitation at the surface or shallow subsurface are spatially rare, but they become focused on specific locations where water tables remain near the surface such as near large rivers such as the Punilla or small perennial streams, bofadeles at high elevations, or vegas near the margins of wetlands. Although large vegetation changes are observed across the region, most of them can be attributed to the differences in precipitation change from southwest to northeast and the relative amount of the aforementioned old water making its way to these wetlands. Lagoon waters at the floor of the largest and deepest basin (Carachi Pampa) show zero tritium activity, indicating it is sustained nearly entirely by old groundwater. This wetland shows less vegetation cover loss than the wetlands upgradient of it. Since this wetland and others like it (such as Purulla) are likely fed by some of the longest and oldest flow paths in the region, its consistent inflow may buffer this system against the yearly to decadal term climate fluctuation experienced at other wetlands. This concept has important implications for assessing the relative response of individual wetlands in this region to natural change and, therefore, their resilience to potential human impacts such as mining and agricultural activity in the area.

Biological Communities' Characterization
Flamingo abundance and distribution across wetlands were more highly variable than almost all the variables considered in this study. The three species of flamingos make itinerant use of these wetlands, changing their distribution according to food availability. Phoenicopterus chilensis occurs in a wide variety of wetlands, including freshwater and salt lakes, estuaries, and marine coasts from Peru to southern Argentina. It is also distributed from Chile to southern Brazil and Uruguay [55,116,117]. Phoenicoparrus andinus and P. jamesi are primarily restricted to wetlands in the high Andean plateau of Argentina, Bolivia, Chile, and Peru, although a portion of P. andinus population disperses to lowland wetlands in central Argentina during winter, especially when the Andean lakes freeze [53,118]. Phoenicoparrus jamesi have been recorded in lowland wetlands in central Argentina [119,120].
Grande is an important wetland for P. jamesi, where between 12,000 and 18,000 individuals representing >10% of the global population were recorded. Higher densities of Bacillariophyta (diatoms) were also found in Grande across samples. Grande, Hedionda, and Cortaderas had the highest concentration of Si, a primary element to produce the frustule (cell wall) of diatoms [121,122]. Diatoms also appeared to be controlled by SRP availability, a nutrient that can be rapidly mineralized by flamingos through their feces, constituting a perfect cycle [123]. Indeed, several studies indicate that flamingo distribution is influenced mainly by food abundance and quality e.g., [124,125], which would explain the high numbers of P. jamesi in Grande, where the highest density of diatoms has been recorded across years compared with the other wetlands sampled. P. jamesi numbers were highly variable in Los Aparejos, Purulla, Las Tunas, and Las Peñas across years, and positively correlated with diatom abundance, which in turn was negatively correlated with the water level of these wetlands. This suggests that increases in water levels in these wetlands negatively affect diatom abundance and therefore food availability.
Our analysis shows that both species (P. andinus and P. jamesi) abundance is linked to diatoms, and P. chilensis numbers are linked to large zooplankton (Cladocera and Copepoda) and other algae groups. The evidence does not suggest that flamingos are a controlling factor of plankton abundance in these wetlands, as was concluded in Frau et al. [46]. However, an experiment performed by Hurlbert and Chang [126] suggested that P. andinus may influence benthic structure (invertebrates and diatoms), at least locally; therefore, new in situ experiments are necessary to elucidate the grazing impact that flamingos may have on the plankton structure in these kinds of wetlands. While the global population of the Chilean Flamingo appears to be stable or increasing [55], P. chilensis numbers in the wetlands censused during this study decreased. Therefore, further studies to determine if P. chilensis abundance is responding to habitat degradation based on increased human impacts or a consequence of natural cycles in the ecological conditions of these wetlands.
The plankton composition and distribution among the wetlands studies was highly variable. A pattern frequently attributed to salinity and to chemical and temperature gradients that occur within and between these kinds of wetlands [43]. Despite these highly variable patterns registered, the phytoplankton assemblage tended to be dominated mainly by Bacillariophyta (diatoms) as several other studies have reported in similar areas [129][130][131]. The dominance of genera, such as Navicula and Nitzschia, both tychoplanktonic algae (adapted to live in stirred waters), is also consistent with several floristics reports [20,36,132,133] in similar environments from the northern part of Argentina. Diatoms were primarily related to silica and SRP availability and middle to high conductivity values. However, some other studies have reported that diatom abundance may be also related to ionic composition [44,134]. Taxonomic affiliations to the species level are difficult. In the Andes region, numerous endemisms comprising a characteristic flora with frequently new species described in several of these reports e.g., [135,136]. For this reason, making this kind of association between ions composition and diatom species could be difficult to reach.
Chlorophyceae was the second most dominant group of algae, frequently represented by just one or two volvaceans species, linked to low nutrient availability and a broad conductivity spectrum, ranging from very low conductance lakes (Antofagasta) to highly saline ones (Tres Quebradas). Less represented groups, such as Euglenophyceae and Cryptophyceae, were generally restricted to high conductivity wetlands and were not in high densities. Their low representation is consistent with what is reported previously in similar environments [137]. It has been speculated that the relative absence of euplanktonic species might be due to the strong UV pressure on these very high lakes [43], with, however, few reports that corroborate this hypothesis e.g., [138]. Remarkably, the presence of euglenoids in saline lagoons has been reported by other authors [40,139,140], indicating that many Euglenophyceae species could be euryhaline or have a broad salinity tolerance.
Cyanobacteria appeared only occasionally and with middle-low densities in Alumbrera, Las Peñas, and San Francisco. In 2013, we recorded a cyanobacteria bloom in Alumbrera dominated by Chroococcales species, which was consistent with potential wastewater infiltration from a nearby plant. No other cyanobacteria blooms were detected since then. Alumbrera has also been increasingly intensively used for livestock (mainly domestic South American camelids), which congregate near the shore of the lake. In recent decades, an increase in urban and agricultural uses of high-elevation wetlands has become more common in the inter-Andean valleys, although the effects of nutrient enrichment on biological communities have not been deeply studied in arid or semiarid zones [95]. Special attention should be given to Alumbrera, a wetland providing services to the town of Antofagasta de la Sierra, to develop sustainable restoration measures and not repeat the eutrophication problems reported by other authors in similar high-elevation ecosystems [141][142][143].
Our results show that zooplankton was not well represented in species richness or abundance in the Andean plateau wetlands, a pattern that is consistent with the findings of other authors in similar extreme ecosystems, who also conclude that conductivity gradients mainly control zooplankton in these systems [144][145][146][147]. Occasionally, one copepod species, Boechella popoensis appeared dominant in Archibarca and Las Peñas-a pattern already described by other authors in high elevation saline wetlands from South America [4,143,148]. However, the general tendency was for low-density values of all groups, especially for Rotifers, which were well represented in species but not in abundance. The genera found to be dominant, such as Boeckella, Cephalodella, or Brachionus, have been reported as saline tolerant in previous studies [145,149,150], which is consistent with the information reported here.

Conclusions
This study presents a dataset gathered over ten years and a much broader climatological analysis from 1985 that provides a baseline of the physical-chemical, hydroclimatological, and ecological components and their interactions among a group of high-elevation Andean wetlands. We found that these wetlands are highly variable in salinity (from freshwater to brine), reach high nutrient concentrations, and are very shallow. Our analysis of the geochemical composition of surface waters in these wetlands showed strong distinctions between several groups of wetlands and strong connections across large distances, indicating important regional groundwater connections exist. We also found that these long regional flow paths containing very old water may buffer some wetlands against decadal-term climate variations. Regarding flamingos, P. jamesi appears as the dominant species in the study area, primarily because of the high numbers in Grande, and phytoplankton tends to be dominated, in general, by diatoms, while zooplankton is poorly represented. Both plankton groups tend to have low species richness and highly variable abundances among wetlands. The climate is arid, and the general tendency observed in the last decades is increasing air temperatures and decreasing rainfall, with two different trends being observed regarding hydrological responses in the wetlands. One set of wetlands has seen an increase and the other a decrease in the aerial extent of living vegetation. Our data suggest that much of the changes in vegetation patterns are in direct response to the observed climate changes; however, we have no evidence yet of how these changes may affect flamingo and plankton communities. In further studies, we should investigate the effects of these tendencies on the main trophic structure of these wetlands (plankton and flamingos), considering the crucial role that vegetation has in the microclimates of these endorheic basins and, therefore, in the water cycles of these wetlands. In sum, our results suggest that all these wetlands sampled are highly variable and can change rapidly; therefore, further efforts to understand how they function are critical to prevent their degradation and secure the valuable environmental services they provide.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/hydrology8040164/s1, Figure S1: Total number of flamingos from the three species recorded reported on each wetland censured from years 2010 to 2020, Table S1: Wetlands sampled and limnological sampling dates, Table S2: Summary of dissolved major ions in surface waters.