Sampling was performed between 24th and 29th August 2023 in the framework of an international expedition with teams from Bolivia, Mexico, Portugal, Spain and Slovenia coordinated by the Universidad Mayor de San Simón and Fundación PROINPA, from Cochabamba, Bolivia.

In order to obtain a robust representation of the microbial communities across the lakes, the sampling strategy was designed to mitigate potential biases by collecting several randomly selected samples per lake and combining them into a pooled sample. A total of eight lagoons were sampled, including four freshwater (salinity < 0.5%) and four saline lagoons (salinity > 0.5%) (Table 1; Figure 1). This sampling strategy ensures that each replicate represents an entire lake rather than a single localized point, minimizing the effects of spatial heterogeneity. These sites were characterized by extreme environmental conditions, severe temperature fluctuations, low precipitation, and altitudes ranging from 3,656 to 4,430 m (Cabrol et al., 2009, 2014). They were specifically selected to investigate the impact of salinity on microbial diversity and functional adaptations. This approach provides critical insights into ecosystem resilience and biogeochemical processes across freshwater and saline ecosystem enhancing our understanding of microbial dynamics in extreme environments.

Approximately, 2 g of the studied Cyanobacterial mats were collected and preserved in Eppendorf tubes containing RNAlater to maintain RNA integrity. Samples were transported to the laboratory and processed immediately upon arrival. For the morphological analysis, the cyanobacterial mats were examined and microphotographed using the LEICA LAS version 4.12.0 image analysis software (Leica Microsystems Limited and CMS GmbH, Switzerland).

Total DNA was extracted from the samples by SGS-Global Biosciences Center (Portugal), using DNeasy PowerSoil Pro Kit, following internally optimized procedures, including quality control, PCR amplification, and library preparation steps. Following internally optimized procedures, the 16S rRNA gene was prepared for each sample along with negative and positive laboratory controls, according to internally optimized procedures. All samples meeting the established quality control thresholds were selected for further processing.

The primer set 338F (5′-ACTCCTACGGGAGGAGGCAGCAG-3′) and 806R (5′ GGACTACHVGGGGTWTCTAAT-3′), tagged with Illumina adapters and linker sequences were used to target the V3-V4 regions of the bacterial 16S rRNA gene. Libraries were sequenced on Illumina NextSeq 2000 platform in 2×300 mode. The sequencing run performance was within the platform specifications. The 16S amplicon sequencing data for this study have been deposited in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB88730.

16S rRNA gene raw sequences from Illumina NextSeq 2000 were pre-processed by removing primers and adapters using QIIME2 platform. Quality filtering, trimming, error models, and dereplication were performed with default settings. Trimmed sequences were denoised to ASV using the dada2 R package. Subsequently, taxonomy at different taxonomic rank levels was assigned using the Naïve Bayes classifier based on SILVA database version 138-NR99.¹

Differences in community composition from freshwater and saline lakes were analyzed using linear discriminant analysis effect size (LEfSe) using a standard threshold of 2.0 logarithmic LDA score and an alpha value of 0.05 for Kruskal-Wallis and paired Wilcoxon test. The core microbiome was assessed using an abundance and occupancy-based method. Briefly, abundance data was transformed into relative abundances via compositional normalization. The definition of core microbiome followed standard thresholds of detection and prevalence, based on the approach outlined by Salonen et al. (2012). Specifically, ASVs were considered part of the core if they were detected at a minimum relative abundance of 0.1% in at least 95% of the samples. The functional analysis was performed using the FAPROTAX software (Louca et al., 2016).

To determine the phylogenetic position of Cyanophyceae sequences from the analyzed samples and obtain a more precise identification at the genus level, the 16S rRNA gene sequences were aligned to those of reference strains of 19 Cyanophyceae orders, according to Strunecký et al. (2023). Alignment was performed using ClustalW, in MEGA11: Molecular Evolutionary Genetics Analysis version 11 (Tamura et al., 2021). The final dataset, comprising Cyanophyceae sequences from this study along with the reference strains, consisted of a total of 520 sequences. The phylogenetic tree was built using maximum likelihood analysis. GTR + G + I evolutionary model was selected by MEGA 11. The robustness of the ML tree was estimated by bootstrap percentages, using 1,000 replications using IQ-Tree online version v1.6.12 (Trifinopoulos et al., 2016). The tree was edited using ITOL.² For the identification of taxa, beyond the phylogeny, 16S rRNA gene identity (p-distance) was used, and sequences were considered to belong to the same genus when presenting >95% identity (Komárek et al., 2014).

Before conducting diversity analyses, samples were rarefied to the same sequencing depth using the R package Microeco (version 1.9.0) (Liu et al., 2021). Taxonomic diversity of the microbiome across both lagoon groups was estimated using the Chao1 and Shannon indices, with analyses performed using the R package vegan (Oksanen et al., 2025). Normality and homogeneity of variances were assessed using the Shapiro–Wilk and Levene’s tests, respectively. When these assumptions were met, analysis of variance (ANOVA) was applied; otherwise, non-parametric comparisons were conducted using the Mann–Whitney U test. All statistical analyses were performed using R version 4.3.2 (R Core Team, 2020, Vienna, Austria). Differences were considered statistically significant at p < 0.05.

The morphological examination of the cyanobacterial mat samples revealed the presence of seven different morphotypes, predominantly dominated by Nostoc Bornet & Flahault, followed by heterocytous and homocytous filamentous forms.

Samples from saline environment of Salar de Uyuni (GBC221, 222), were composed by two different Nostoc morphotypes. In GBC221, trichomes were short, densely arranged, and enveloped by individual mucilaginous sheaths in addition to the colony envelope (Figure 2a). Conversely, in GBC222, trichomes were longer, loosely entangled, and lacked individual sheaths (Figure 2b).

Regarding the freshwater stream samples, in the Catal lake (GCB223) sample, a Nostoc morphotype identical to GBC222 was observed. In the freshwater stream adjacent to Laguna Colorada (GBC224), a homocytous filamentous morphotype was identified. This morphotype presented cylindrical trichomes of 5–6 μm, commonly attenuated toward the ends, within thin colorless sheaths, sometimes forming fascicles and cells shorter than wide. These characters are similar to Microcoleus Gomont and Phormidium-like types (Figure 2c).

In the stream close to Laguna Chalviri (GBC225) a third Nostoc morphotype, resembling that of GBC221, but with less densely arranged trichomes within the colony, was observed (Figure 2d).

The single sample collected from the shoreline of the freshwater Laguna Chojllas Pequeña (GBC226), was composed of a fourth Nostoc morphotype, distinguished by a brownish envelope and denser arrangement of trichomes, compared to the other samples (Figure 2e).

Nostoc morphotypes identified in samples from the shores of saline lakes Laguna Coruto (GBC227), and Laguna Pastos Grandes (GBC229) were identical to those observed in samples GBC226 (freshwater) and GBC221 (saline), respectively. Thus, this suggests a remarkable ecological plasticity within the genus, allowing it to thrive in both freshwater and saline. Sample GBC228 exclusively featured a homocytous morphotype, with parallel arrangement of the trichomes forming fascicles. Trichomes appeared constricted, cylindrical, 4–6 μm wide, with cells longer than wide, and conical apical cells. These characteristics resemble those of the Coleofasciculales (Figure 2f). Samples GBC230 and GBC231 from Laguna Pastos Grandes consisted of a heterocytous filamentous morphotype with parallel arrangement of trichomes, characterized by heteropolar (basal and intercalar heterocytes), attenuated toward the ends, and thick yellow-brownish sheaths. Vegetative cells were shorter than wide, features typical of the Rivulariaceae family (Figures 2g,h).

Amplicon sequencing produced 4,113 16S rRNA gene amplicon sequence variants (ASVs) for the bacterial community. Of these, 310 ASVs were identified as Cyanobacteria, with 134 classified as Cyanophyceae (Strunecký et al., 2023) using a Naïve Bayes Classifier trained on 16S data from the SILVA database. This metabarcoding analysis identified 32 different taxa (Supplementary Table S1). Based on the number of reads, the relative abundance of the class Cyanophyceae was 20%, with bacteria from other classes constituting the remaining 80% of bacterial abundance.

Regarding Cyanophyceae sequences, additional qualitative analyses were conducted to achieve a more precise identification of the cyanobacterial genera present in the samples. For this purpose, a 16S rRNA gene ML phylogenetic analysis (Figure 3) and 16S identity (p-distance) comparisons were performed. These analyses are considered more accurate than SILVA-based metabarcoding approach, as the reference strains included in the phylogenetic tree were manually curated by the authors, ensuring a more precise and reliable taxonomic classification. Supplementary Table S2 provides detailed identification of the Cyanophyceae sequences at the genus level, including the 16S rRNA gene identity (p-distance) with the closest phylogenetically related cyanobacterial genus.

These analyses revealed that among the 134 Cyanophyceae ASVs, 104 could be identified within forty-two different genera, distributed across nine orders. This represents a significantly greater diversity than the morphological analysis, which identified only seven morphotypes predominantly related to Nostoc. Furthermore, the phylogenetic analysis detected an even higher level of diversity compared to the metabarcoding analysis, which identified 32 taxa (Figure 3; Table 1; Supplementary Tables S1, S2).

Comparing genera, according to the phylogenetic analysis, a total of 18 different ASVs were identified as Nostoc, making it the genus with the highest diversity of sequences. Additionally, nine different ASVs were within the Nodosilinea clade, while eight different sequences were identified as Coleofasciculus Siegesmund et al. The genera Cyanobium Rippka & Cohen-Bazire and Mojavia Reháková & Johansen, exhibited five ASVs each (Figure 3, Supplementary Table S2).

The order with the highest diversity of genera was Nostocales, encompassing a total of 11 genera, followed by Chroococcales, with nine, and Oscillatoriales with five (Supplementary Table S2).

In addition to these identified sequences, 30 additional sequences were distributed among 16 distinct monophyletic clades. Interestingly, none of these sequences were closely related to any previously described cyanobacterial genera, exhibiting less than 95% 16S rRNA gene identity with their closest phylogenetic relatives (Figure 3; Supplementary Table S2). These findings suggest the potential existence of at least 16 novel cyanobacterial clades that may warrant formal classification in the future.

Among these clades, the orders with the highest number of unidentified clades were Leptolyngbyales and Chroococcales, with four clades each, followed by Nostocales with three potential new genera. Primarily, these unidentified sequences have originated from the saline lake at Laguna Pastos Grandes (GBC231) and a freshwater stream at Laguna Chalviri (GBC225) (Figure 3).

Regarding the abundance of taxa (number of reads), the limitations of the Naïve Bayes Classifier in identifying taxa at the genus level prevented precise determination of the number of reads for each genus identified by the phylogenetic analysis. Alternatively, the number of reads for “genera complexes” were identified as shown in Supplementary Table S1. Through this approach, the most abundant complex was found to be Mojavia/Nostoc/“unknown genus 6,” with 90,781 reads. This complex was detected in Salar de Uyuni (GBC221, GBC222), Catal Lake (GBC223), Laguna Chalviri (GBC225, GBC226), Laguna Coruto (GBC227), and Laguna Pastos Grandes (GBC229), indicating a remarkable ability of these genera to adapt to various conditions, ranging from freshwater to saline environments. Notably, this complex represented 57% of the total Cyanophyceae reads across all samples. Nostoc-like macroscopic colonies were commonly found thriving in these extreme environments.

The second most abundant genus complex was Coleofasciculus/Pycnacronema Martins & Branco, with 21,366 reads, representing 13% of the total Cyanophyceae across all samples. However, this complex was found exclusively at Laguna Mama Khumu (GBC228), a saline lake (Table 1; Supplementary Table S1). The Coleofasciculus/Pycnacronema mats were observed growing beneath a thin layer of pink bacterial biofilm on soil.

The third and fourth most abundant genera complexes were Rivularia Bornet & Flahault and “unknown genus 1”/Limnoraphis/Microcoleus/Salileptolyngbya, with 19,766 and 17,452 reads, respectively. Rivularia accounted for 12.5% of the total Cyanophyceae reads across all samples. It was identified in Laguna Pastos Grandes (GBC229, GBC230, GBC231) and Laguna Chalviri (GBC225), showcasing its adaptability to both saline and freshwater environments. The complex “unknown genus 1/Limnoraphis/Microcoleus/Salileptolyngbya Zhou & Ling comprised 11% of the total Cyanophyceae reads across all samples and was exclusively sampled from a freshwater stream at Laguna Colorada (GBC224). Additionally, notable findings include Prochlorococcus Chisholm et al., Timaviella Sciuto & Moro, Thermoleptolyngbya Sciuto & Moro and Altericista Averina et al., which are further discussed in the subsequent sections.

Among the samples analyzed, the Laguna Chalviri stream (GBC225) stood out as the most diverse, with 17 taxa identified by metabarcoding, corresponding to 22 genera identified by phylogenetic analysis (Table 1; Supplementary Table S1). Regarding Cyanophyceae abundance across samples, GBC226, collected from the freshwater of Laguna Chojllas, exhibited the highest number of Cyanophyceae reads (26,284), constituting 16.6% of the total reads. Despite its high number of reads, this sample showed low diversity, mainly composed of the Mojavia/Nostoc/“unknown genus 6” complex. Similarly, sample GBC228, from a saline lake in Laguna Mama Khumu, contained 21,719 reads, representing 13.7% of the total reads, with Jaaginema Anagnostidis & Komárek and the Coleofasciculus/Pycnacronema complex as the dominant taxa. The third sample with the highest number of reads was Laguna Chalviri (GBC225) with 18,230 reads, representing 11.5% of the total reads. Notably, as mentioned above, this sample had the highest diversity, comprising 17 different genera complexes (Table 1; Supplementary Table S1).

Metabarcoding analysis revealed a dominance of organisms classified into six major phyla: Actinobacteriota, Bacteroidota, Cyanobacteria, Firmicutes, Proteobacteria, and Verrucomicrobiota (Supplementary Figure S1; Supplementary Table S3). Interestingly, Campilobacterota, Deinococcota, and Halanaerobiaeota appeared to be exclusively represented in the saline samples. Regarding Catal Lake, despite being a freshwater lagoon, this sample exhibited unusually high bacterial viability compared to other freshwater systems (Supplementary Figures S1–S4). Due to this atypical profile, we excluded it from group-level comparisons of bacterial diversity and community composition between freshwater and saline lakes. Instead, the Catal Lake sample was considered separately and used as a reference to provide contextual insight, rather than being included in the comparative statistical analyses.

At the class level, Alphaproteobacteria, Bacteroidia, Gammaproteobacteria, Clostridia, Cyanophyceae, and Campylobacteria were dominant across both sample groups. The bacterial composition of the different samples is illustrated in Figures 4, 5.

Most bacterial groups were present across all samples, although their relative abundances varied. The core microbiome consisted of bacteria from the classes: Acidomicrobia, Actinobacteria, Bacteroidia, Cyanobacteria, Deinococci, Clostridia, Gammatimonadetes, Planctomycetes, Alphaproteobacteria, Gammaproteobacteria, and Verrucomicrobiae.

Samples collected from freshwater reservoirs exhibited the highest bacterial diversity, with Alphaproteobacteria, Bacteroidia, Gammaproteobacteria, Cyanobacteria, and Bacilli being the dominant classes. In contrast, samples from the saline lakes were predominantly composed of Alphaproteobacteria, Bacteroidia, Gammaproteobacteria, Cyanobacteria, and Clostridia (Figure 4).

Figure 5 illustrates the most abundant bacterial genera identified through metabarcoding. Nostoc was prevalent and abundant across all environments. Pseudomonas showed higher prevalence in Laguna Chojllas, particularly associated with Mojavia, Nostoc, “unknown genus 6” and Chamaesiphon Braun, but was scarce in other samples. Weissella soli Magnusson was found in high abundance in Catal Lake, but was absent from the remaining sites. Clostridium Cohn was notably abundant in Laguna Coruto.

When analyzing potential biomarkers for freshwater and saline environments, a differential community composition was detected at all taxonomic levels, including at the class level. Campylobacteria was identified as a biomarker for saline environments, while Blastocatellia appeared to be exclusive to freshwater streams (Supplementary Figure S5). At the family level, eight markers were identified for saline environments and 17 for freshwater streams. As shown in Figure 6, Arcobacteriaceae and Flavobacteria were the most representative families in saline samples, whereas Leuconostocaceae and Beijerinckiaceae were dominant in freshwater streams.

Regarding the relationship between cyanobacterial mats and their respective microbiome, the microbial communities associated with Nostoc mats exhibited significant diversity across different environments (Table 1; Supplementary Figure S4). For example, in Salar de Uyuni (samples GBC221 and GBC222), the Nostoc-associated microbiome comprised similar proportions of Alphaproteobacteria, Bacteroidia, Gammaproteobacteria, Clostridia, Bacilli, and Campylobacteria. In contrast, in Catal Lake (GBC223), the Nostoc microbiome was predominantly composed of Weissella (Bacilli). In Laguna Chojllas Pequeña (GBC226), Pseudomonas (Gammaproteobacteria) dominated the microbial community associated with Nostoc mats, whereas in Laguna Coruto (GBC227) Clostridium (Clostridia) was the dominant genus (Table 1; Supplementary Figure S4). These results suggest that the composition of the microbiome associated with Nostoc mats is likely more influenced by environmental factors than by the Nostoc mat itself.

When comparing samples from freshwater and saline environments, the microbiome associated with cyanobacterial mats showed strong variation driven by abiotic factors. Taken together, these findings underscore the likelihood that the structure of cyanobacteria-associated microbiomes is more regulated by the environmental conditions than by cyanobacterial identity alone.

The inferred metabolic functions of the associated microbiome are presented in Figure 7, while Figure 8 provides a more detailed description of the metabolic characteristics of the most relevant microorganisms. Overall, the major functions of the microbiome encompassed phototrophic and chemoheterotrophic metabolisms, sulfur respiration, nitrogen fixation, nitrate and nitrite metabolism, and fermentation. Notably, nitrogen-fixing cyanobacteria, such as Nostoc, Rivularia, and Coleofasciculus were prevalent, along with bacteria known for their bioremediation potential, including members of the genera Shewanella and Pseudomonas.