Analysis of molecular diversity within single cyanobacterial colonies from environmental samples

Attached or floating macroscopic cyanobacteria can be found in shallow waters and can be easily hand-collected, but their identification is often challenging due to their high morphological variability. In addition, many members of environmental samples lose their morphological adaptations under controlled conditions, making the integration of analyses of field populations and derived isolated cultures necessary in order to evaluate phenotypic plasticity for identification purposes. Therefore, in this study, twenty-nine macroscopic field samples were analyzed by Illumina sequencing and parallel optical microscopy. Some colonies showed the typical morphological characteristics of Rivularia biasolettiana, and others showed those of Rivularia haematites. However, other Rivularia-like colonies showed ambiguous morphologies, and some of them showed the phenotypic features of the new genus Cyanomargarita, which is virtually indistinguishable from Rivularia in the field. In all of the colonies, phylotype composition was highly heterogeneous, with abundances varying depending on the analyzed sample. Some colonies were dominated (97–99%) by a single phylotype, while in others, the percentage of the dominant phylotype decreased to approximately 50–60%. Surprisingly, the same dominant phylotype was found in R. biasolettiana and R. haematites colonies. The relationships between environmental and/or biological factors and morphological variability in these colonies are discussed.

Cyanobacteria are a very special group of gram-negative prokaryotes. Because they developed the ability to use water as an electron donor in oxygenic photosynthesis, approximately 3.6 billion years ago 1 , they played an important role in the evolution of life on Earth, the subsequent generation of molecular O 2, and the oxygenation of the atmosphere 2 . Cyanobacteria are also the ancestors of chloroplasts 3,4 . In addition, many cyanobacteria are capable of fixing atmospheric N 2 and thus play an important role in nitrogen cycling 5 .
It is easy to understand why cyanobacteria were previously called blue-green algae since, in addition to performing plant-like oxygenic photosynthesis, macroscopic green forms similar to algae can be easily found and observed by the naked eye, with a cosmopolitan distribution ranging from hot hyperarid deserts to polar aquatic environments 6 . These alga-like forms, in which individuals are microscopic but exhibit macroscopic growth, have been described as colonies, thalli, mats, tufts, thin coatings, soft gelatinous covers, subspherical globes, tightly attached felts, etc. [7][8][9] .
The importance of characterizing and identifying natural cyanobacterial populations in order to compare them with corresponding cultures, which are further used, e.g., in phylogenetic analysis or potential applications in biotechnology, has been emphasized [10][11][12][13] . However, unfortunately, the majority of recent studies were based only on isolated strains and morphological descriptions corresponding to culture media, and controlled conditions were often different from the characteristics found in the environments where highly diverse cyanobacteria live 7,13 . In a previous study, we compared phenotypic characteristics of natural populations and isolated cultures from environmental samples and found large differences between them, which can lead to misidentification of strains 12 .
Rivulariaceae includes all the heterocystous cyanobacteria with tapered trichomes, with a clear distinction from the base to the apex. A terminal heterocyst (heterocyte) occurs at the broad basal end of the trichomes, and a sheath encloses everything but the heterocyst 7 . These tapering cyanobacteria include many visually conspicuous

Results
Morphological characterization of environmental samples. Twenty-eight collected samples (Table 1) showed a macroscopic Rivularia-like morphology in which hemispherical or slightly irregular-hemispherical lobate colonies approximately 0.5-3 cm in diameter/length could be observed (Figs. 1, 2, 3, 4 and 5). After microscopic evaluation, the Rivularia-like colonies could be separated according to their specific features.
Seven other colonies (GOR11, GOR12, OSI3, ARA4, GOR2, END1 and END8) were also hemispherical, some of them soft and gelatinous and others highly calcified with sections showing zonation ( Fig. 3a-d). The trichomes also showed the typical tapering of those in Rivularia colonies and other similar morphological characteristics ( Fig. 3a-h). However, the main difference between these colonies and the previously analyzed colonies was the size of the cells ( Table 2, Fig. 3g,h), which were clearly wider in these colonies (significant differences, p < 0.05). The mean trichome diameter at the filament base ranged from 7.8 to 10.1 μm. The heterocysts were also clearly larger (Table 2). These features corresponded with those of the new genus Cyanomargarita 15 .
However, four Rivularia-like colonies (MU4, GDL1, BAT2 and BAT4) showed ambiguous morphology, varying in the degree of calcification and zonation ( Fig. 4a-d) as well as displaying high heterogeneity in the trichomes ( Fig. 4e-h). For instance, GDL1 was a highly calcified colony with a certain degree of zonation, and most of the filaments were of the Rivularia type, but there was a considerable proportion of other morphotypes within the colony, such as thinner and nontapered trichomes, without heterocysts, that were found between the Rivularia filaments and sometimes emerging outside the colony (Fig. 4a,e,g). The MU4, BAT2 and BAT4 colonies also presented filaments and trichomes with very different sizes and shapes in the same colony ( Fig. 4b-d,f,h). In BAT4, clear zonation could be observed (Fig. 4d,f), as well as a large number of filaments without heterocysts (Fig. 4f).
Microscopic examination of the hemispherical and calcified BAT13 colony (Fig. 5a) showed features different from those of the previously analyzed colonies (Fig. 5b-i). Most of the filaments did not taper, although this was due to the presence of a wide and colored sheath (Fig. 5b-i). When observed in detail, the trichomes narrowed from the base, which presented a heterocyst, to the apex ( Fig. 5c-g). Although this colony appeared macroscopically as a Rivularia colony (Fig. 5a), the trichome arrangement resembled that in the genus Dichothrix 7 , such as subdichotomic falsely branched filaments, where secondary trichomes remained within the 'mother' sheath ( Fig. 5c,d,f). Some filaments showed a bunch of trichomes, open at the end with a funnel-like ending (Fig. 5c).  www.nature.com/scientificreports/ quadratic or shorter than wide at the end, tapering into a hyaline hair ( Fig. 5e-g). The heterocysts were basal and mainly conical (Fig. 5d,f), but some intercalary and cylindrical heterocysts could be observed (Fig. 5e). Filaments of Calothrix were also found inside the sample (Fig. 5h). www.nature.com/scientificreports/ Finally, sample MUD1 showed the typical characteristics of the genus Dichothrix 7 , such as macroscopic brush-like fasciculated tufts encrusted by calcareous precipitate (Fig. 6a-d) and filaments in a characteristic dichotomic arrangement that were repeatedly falsely branched (Fig. 6a,b). The filaments were tapered and had strong lateral false branching adjacent to heterocysts (Fig. 6c,d). The branches formed individual sheaths and diverged from the basal filament in parallel, and the new trichome shared the sheath with the old trichome in the www.nature.com/scientificreports/ basal part (Fig. 6c,d). The sheaths were firm, yellow to brown, and lamellated in only some filaments (Fig. 6d). The trichomes gradually narrowed towards the apex, with barrel-shaped cells (Fig. 6c,d). No hairs were observed. The basal heterocysts were hemispherical or barrel-shaped with a characteristic intense blue-green color (Fig. 6d). www.nature.com/scientificreports/  hypervariable region of the 16S rRNA gene from all the collected samples revealed in total 25 OTUs present at an abundance of 1% or more in at least one of the colonies. OTUs were numbered in order of decreasing abundance, taking into account all the colonies. The first two OTUs accounted for 81% and the first 10 accounted for 95% of all the reads, considering all the colonies together. Dedicated phylogenetic trees constructed with our sequences from cultures and field samples and those downloaded from the NCBI database ( Fig. 7) allowed us the taxonomic assignment of these OTUs, as previously carried out 27,28 . Fourteen OTUs were mapped to the heterocystous subtree (Fig. 7a). The most abundant OTU (OTU1) fell into a very well-supported cluster with other Rivularia sequences. The percent identity between OTU1 and the other sequences in this cluster ranged from 97.6 to 100%, with the OTUs being 100% similar to an environmental colony collected from the Alharabe River in southwestern Spain 11 and to other environmental samples isolated from the Baltic Sea 29 . Thus, OTU1 was assigned to Rivularia sp. OTUs 2 and 13 also fell in a very well-supported cluster with Cyanomargarita melechinii and Cyanomargarita calcarea. OTU2 was more similar to C. melechinii (99.52%), while OTU13 was more similar to C. calcarea (98.07%). Both OTUs were assigned to Cyanomargarita sp. OTUs 3, 11, 14, and 26 clustered with other Calothrix sequences; thus, they were assigned to this genus. OTUs 20 and 4 were included in a cluster with Macrochaete strains and were assigned to this genus; OTU4 was 98.3% similar to Macrochaete lichenoides, and OTU20 was 96.4-96.63% similar to Macrochaete psychrophila. OTU6 was mapped to a cluster that included other Calothrix sequences and Scytonematopsis contorta strains, with similarities ranging from 95.9 to 98.3%, as well as other environmental sequences obtained from Dichothrix tufts from the Muga River, which were 99.5% similar; therefore, it was assigned to Dichothrix sp. OTU8 was assigned to Scytonema sp. since it mapped to a well-supported clade with other Scytonema sequences. OTU19 clustered with other sequences of Nostoc and was assigned to this genus, while OTU15 was mapped to a sister clade of Anabaena, Nostoc and Tolypothrix strains and was therefore assigned only to the Nostocaceae family. Finally, OTU9 mapped to a cluster with an uncultured bacterium, and the maximum similarities found in the NCBI database were 95.22 with an uncultured bacterium and 93.78 with Calothrix sp. CAL3363 29 ; therefore, this OTU was assigned to the Nostocales order.
Regarding the nonheterocystous cyanobacteria, eleven OTUs mapped to this subtree (Fig. 7b). The most abundant OTUs, OTU5 and OTU7, fell in two well-supported clusters with Oculatella strains and Phormidium strains, respectively. OTUs 16, 17, 18, 24 and 27, assigned to Leptolyngbya sp., mapped to clades with other Leptolyngbya sequences. OTU28 clustered with other Schizothrix strains. OTU10 and OTU21 could be assigned only to Leptolyngbyaceae and Oscillatoriales, respectively, due to their low percent identity with other sequences in the databases. Finally, OTU12 remained unassigned since no match was found in the databases for this sequence, and there was no clustering with any close relatives in the phylogenetic tree (Fig. 7b).  www.nature.com/scientificreports/ In six of the seven colonies of the R. haematites type, the abundance of OTU1 ranged from 97 to 99.8%, and in one colony, the abundance was 85.5% (Fig. 8). OTU3, corresponding to Calothrix sp., and OTU5, corresponding to Oculatella sp., were other phylotypes found in these colonies, with abundances ranging from 1.4-8.6% for Calothrix and 1.5-20.9% for Oculatella (Fig. 8).
In the seven Cyanomargarita colonies, OTU2 dominated, ranging from 97.1 to 99.2% in 4 colonies, decreasing to 86.6% in another colony, and further decreasing to 67.6 and 57.28% in the other two colonies. Regarding the other phylotypes found in these colonies, OTU7, corresponding to Phormidium sp., was the most abundant and present in all the colonies, with abundances of 12% and 13.5% in two of the colonies (Fig. 8). A phylotype corresponding to Calothrix sp., found in previously analyzed Rivularia colonies (OTU3), was also found in a colony of Cyanomargarita (GOR2), with an abundance of 25%. OTU13, which corresponded to another phylotype of Cyanomargarita, was present in only one colony (END1), with a 13.5% abundance.

Discussion
Genotypic heterogeneity in single Rivularia-like colonies. Rivularia-like colonies have a global distribution, occurring in marine or freshwater habitats, where they are usually attached to a rocky substrate; however, many studies have reported that Rivularia spp. are associated with unpolluted environments 14 . In addition, the relationships between some morphological or physiological features and the environment make these species excellent environmental indicators of changes in running water quality, mainly related to eutrophication processes 14,30 . Therefore, they have been included in biomonitoring programs 21,31,32 . On the other hand, because Rivularia colonies sometimes persist for very long periods, avoiding grazing, the toxicity of these colonies is being investigated 33 . It is undoubted that in all of these approaches, where genera and species must be strictly identified from environmental samples, accurate cyanobacterial characterization is essential.
Traditional identification of cyanobacteria involves assigning a colony to a morphospecies, and conventionally, a bacterial colony is defined as a visible mass of clonal microorganisms, all of which originated from a single cell. However, the results from the present study show that the majority of the analyzed colonies consist of different clones growing together. Among the 28 Rivularia-like colonies, the phylotype corresponding to Rivularia sp. was present in 19 colonies, with abundances ranging from 59.4 to 99.8% depending on the studied colony. Nevertheless, it should also be noted that in most of the colonies, this phylotype dominated, whereby in 14 colonies, it presented an abundance of ≥ 90% (and within 7 of these colonies, the abundance was close to 99%). However, in three colonies, the abundance ranged from 72 to 85%, and in two of them, the abundance decreased to approximately 60%. The other highly abundant phylotypes found in these colonies, which reached abundances up to approximately 21%, corresponded to Calothrix sp. and Oculatella sp., the latter a genus morphologically similar to Leptolyngbya but separated from it because of genetic differences 34 . These results indicated great variability in the abundance of the phylotype corresponding to Rivularia depending on the analyzed colony, as well as variation in the other phylotypes and their abundances found in these colonies.
One of the surprising findings was that among the twenty-eight analyzed Rivularia-like colonies, seven corresponded to the new, recently described genus Cyanomargarita, which as the authors described, is virtually indistinguishable from Rivularia in field samples 15 . In these colonies, genotypic heterogeneity was also found, in which the abundance of the phylotype corresponding to Cyanomargarita varied from 57,28% in a colony with clear lamination resembling R. haematites (see Fig. 3b) to 99.2% in a soft colony resembling R. biasolettiana. Interestingly, in these colonies, Phormidium sp. was the dominant nonheterocystous cyanobacterium instead of Oculatella from Rivularia colonies, but the phylotype corresponding to Calothrix was also found.
Furthermore, phylotypes corresponding to Cyanomargarita and Rivularia were never found together in the same colony, although both types of colonies coexisted in the same rivers (e.g., Gordale Beck and Endrinales). Allelopathic effects could explain these results, as previously suggested for other cyanobacteria 35 . In fact, García-Espín et al. 33 showed that extracts obtained from Rivularia colonies affected the photosynthetic activity of several diatoms and a red alga. Further experiments with extracts from both colonies would confirm this possible effect.
Another very surprising finding was that two Rivularia-like colonies did not present any phylotypes corresponding to Rivularia or Cyanomargarita (or contained them at an abundance ≤ 0.7%). In one of these colonies (colony BAT4), five different phylotypes were found at similar abundances (approximately 15-20%), of which three corresponded to different Calothix spp. and the others corresponded to other Nostocaceae and Leptolyngbyaceae. In the other colony (BAT13), the dominant phylotype corresponded to the new genus Macrochaete 16 . This genus has been described only from cultures, so to the best of our knowledge, this is the first report in which a natural population is morphologically and genetically characterized. Nevertheless, it is noteworthy that the morphological characteristics of filaments and trichomes in this environmental sample were different from those reported in the description of this new genus, in which the phenotypic features resembled those of Calothrix. However, these features corresponded only to isolated strains, which are known to exhibit morphological variability and differences from natural populations 7,12,13 .

R. biasolettiana vs R. haematites.
However, what was very interesting and deserves to be highlighted is that when we tried to differentiate the two typical Rivularia colonies found in calcareous streams, R. biasolettiana and R. haematites, we did not find genetic differences, at least at the studied level, the 16S rRNA gene.
16S rRNA is the most widely used marker gene 36,37 , which fits the criteria of ubiquity, regions of strong conservation, and regions of hypervariability 38,39 . This gene is supported by reference databases containing over a million full-length 16S rRNA sequences, therefore spanning a broad phylogenetic spectrum 40 . The 16S rRNA gene has served as the general framework and as the benchmark for the taxonomy of prokaryotes 41 . Advances in high-throughput sequencing technologies have enabled almost comprehensive descriptions of bacterial diversity through 16S rRNA gene amplicons, which have been used in surveys of microbial communities to characterize the composition of microorganisms present in environments worldwide [42][43][44][45] . Although some issues have been raised, such as identification of metabolic or other functional capabilities of microorganisms when studies focus only on this gene, recent studies have shown that the phylogenetic information contained in 16S marker gene sequences is sufficiently well correlated with genomic content to yield accurate predictions when related reference genomes are available [46][47][48][49] . Therefore, the 16S rRNA gene continues to be the mainstay of sequence-based bacterial analysis, vastly expanding our understanding of the microbial world 50  www.nature.com/scientificreports/ In particular, in cyanobacteria, as in other prokaryotes, the 16S rDNA gene is currently the most commonly used marker for molecular and phylogenetic studies 51,52 . The information obtained from 16S rDNA gene phylogenetic reconstructions, together with morphological, ultrastructural, and ecological data, led Komárek et al. 53 to propose the current accepted classification of cyanobacteria. There have also been specific studies by this group concerning the problems associated with single-gene phylogenies, in which robust phylogenomic trees of cyanobacteria derived from multiple conserved proteins have also shown congruence between the multilocus and 16S rRNA gene phylogenies, which once again demonstrates the considerable strength of the 16S rRNA gene for phylogenetic inference and evaluation of prokaryote diversity [54][55][56][57] .
In this study, in contrast to the genetic identity found in R. biasolettiana and R. haematites colonies, showing a dominance of OTU1, the remainder of the studied representatives of Rivulariaceae showed a wide range of variation in the 16S rDNA sequences and with OTU1. Sequence identity between OTU1 and the remaining OTUs belonging to this family was as low as approximately 90%, ranging from 90.73 to 93.41%, and when it was compared with other Rivulariaceae from the databases, in the different clusters of the phylogenetic tree, this value ranged from 87.12 to 93.90%. A large difference between the sequences of this gene was also found in other studies on Rivulariaceae [15][16][17]29,58 . In fact, several new genera are emerging on the basis of these differences [15][16][17] . Comparisons of phylogenies using other markers, such as the phycocyanin operon and the intervening intergenic spacer (cpcBA-IGS) with the 16S rRNA gene in previous studies in Rivulariaceae, have shown largely consistent results, with a high level of divergence between the components of this family 11 .
In addition, the results of the present study showed correlations between morphological characteristics and the analyzed genes in all the cyanobacterial colonies/tufts, except for those of R. biasolettiana and R. haematites. In these two cyanobacteria, only distinct macroscopic phenotypic features were observed due to zonation and different degrees of calcification since no significant differences were found in the size measurements or other microscopic characteristics.
Therefore, although the remainder of the genome has not been studied in these populations, the genetic identity of the studied marker, phenotypic features, together with environmental preferences point out that R. biasolettiana and R. haematites are ecotypes of the same species, as previously suggested 59 .
R. biasolettiana and R. haematites have very similar morphotypes, and traditional taxonomical classification and studies have distinguished them primarily by their degrees of calcification. R. biasolettiana-type colonies are described as more gelatinous and less calcified, and the crystals are disseminated; however, R. haematites colonies are very hard and exhibit extensive calcification in concentric zones, which leads to clear lamination 24,25,60,61 . Because of its heavy mineralization, R. haematites is a model for stromatolite-binding organisms 25,26 .
Microscopic observations from this study showed that some colonies presented typical R. haematites morphology with concentric bands of intense calcification (see Fig. 2a,b), and others were soft and less calcified, such as R. biasolettiana, although all of them presented the same dominant phylotype. Many others with this dominant phylotype have also shown ambiguous morphology with no clear lamination, although some dark/light zones could be observed (see, e.g., Fig. 4b,d, f). Even in Cyanomargarita colonies, whose genotype was clearly separated from that of Rivularia, concentric zones and extensive calcification could be observed (see, e.g., Fig. 3b,d). These results suggested that these phenotypic features are not diagnostic characteristics for further identification.
In a two-year study, Obenlüneschloss and Schneider 61 found that not all analyzed R. haematites colonies showed distinct concentric calcification layers. In the stromatolites of both types of Rivularia, the same lamination was observed, and the differences in calcification appeared later 60 . Pentecost and Franke 26 compared populations of R. biasolettiana and R. haematites and argued that although both could be distinguished by their form of calcification and their trichome diameter, some populations of R. biasolettiana were more intensely calcified than others, suggesting that a continuity of forms may exist, even within the same stream, and therefore, a continuum of colony forms probably occurs between these taxa.
Differences in the calcification pattern have been attributed to seasonality and cyanobacterial activity, in particular to photosynthesis 24,26,62 . The calcification in R. haematites occurred in concentric bands, which varied in thickness and the density of crystals. Since characteristic zonation is formed by filaments of different successive generations, the thickness will vary depending on the growth rate, while crystal density will depend on the rate of calcification. Calcification is the result of photosynthesis (with a maximum of 14%) and evaporation during the warmer seasons, while it is entirely abiogenic during winter as a result of CO 2 evasion 63 . Therefore, dense calcified bands similar to those formed in winter have been described that are caused by a reduction in trichome growth and EPS production, allowing the development of abiotic surface precipitate, and less calcified layers are formed during spring and summer, when calcification is associated with photosynthesis in zones of growth with cell division 24,26 . Thus, differences in climatic conditions and/or biological activity seem to lead to differences in the degrees of zonation and calcification.
The growth of Rivularia colonies is seasonal and strongly correlated with water temperature 24,26 . The colony growth rates were 12-14 µm/day in summer and 2 µm/day in winter 24 . The occurrence of R. biasolettiana was more closely related to high temperatures than that of R. haematites 21 . Moreover, colonies of R. haematites were generally collected under temperatures below 15 °C in mountain running waters 64 , and R. haematites stromatolites have been described as preferentially developed in wet periods, particularly in autumn and winter 60 . Our own field observations during the sampling for this and previous studies were that the gelatinous and weakly calcified R. biasolettiana type was more abundant in warmer locations, and in contrast, R. haematites was dominant in cold locations (data not shown).
One possible explanation for the results found in this study could be related to these differences in the degree of zonation and calcification in relation to climate, which could include microclimatic conditions. In warmer sites or climatic conditions, when growth is rapid, the number of filaments will increase, moving towards the surface in a weakly dense and unaligned arrangement, on which calcite crystals spread, providing a lighter and less calcified structure. www.nature.com/scientificreports/ colder conditions, such as in winter, or microclimatic conditions, when growth slows down for other reasons, such as low light, filaments become more densely packed, allowing the development of extensive precipitates and leading to a dark band. When these conditions change, e.g., in the spring and summer, increases in temperatures and/or light will result in increasing and faster growth, leading to a less calcified new layer, and successive seasonal and/or microenvironmental changes will result in the typical lamination of R. haematites. Therefore, warmer places with high temperatures and/or light will allow the occurrence of the R. biasolettiana type, while in colder sites and/or sites with alternating environmental conditions, the R. haematites type will develop. Shaded colonies and colonies that lie in the supratidal spraywater zone often contain small, irregular and more densely packed crystals 61 .
Cyanobacteria are known to modify EPS production, pigments, and morphology under environmental stimuli 6 . The production of EPS also varies depending on the cyanobacteria, whereby Rivularia has shown a well-developed exopolymer layer 65 , which is of great importance for this epilithic cyanobacteria, as it acts as an adhesive that allows cells to stick to the stones in the running waters, and it holds the filaments together, minimizing cell damage during intermittent drying exposure to the air and evaporation in the warmer seasons 66 . The C/N ratio is an important parameter for the variation in EPS production since high amounts of fixed C compared to N levels drive EPS synthesis to store excess C 67,68 . Therefore, Rivularia colonies that are exposed, in spring and summer, to high light intensities and temperatures will increase their photosynthetic rates and therefore the amount of EPS, as shown by the R. biasolettiana morphotype. In addition, most of the analyzed populations were dark in color, probably in relation to the accumulation of the yellow-brown scytonemin pigment in the sheaths or EPS, as previously observed in shallow and clear oligotrophic ecosystems, where water clarity allows UV radiation to penetrate well, protecting the cells from the damaging effects of this radiation 69,70 .
In conclusion, environmental factors can lead to differences in macroscopic phenotypic features, such as those found in the Rivularia colonies studied here. However, further sampling under different climatic conditions and/ or microenvironmental conditions or of Rivularia cultures grown under distinct temperature and/or illumination conditions, as well as analysis of other genes, are needed to confirm this hypothesis.

Methods
Sampling. The locations of the sampling sites from which samples were collected and their codes are shown in Table 1. All the sampled rivers or streams were characterized by highly calcareous waters. The Muga, Guarga, Osia, and Arás rivers are located in northeastern Spain, the Hoyas and Endrinales streams and Bogarra River are located in southeastern Spain, and the Guadiela River is located in central Spain. The Gordale Beck sampling site was located near Malham, West Yorkshire, northern England, from which Rivularia haematites has been reported 24,26 . From these locations, twenty-one Rivularia-like colonies were sampled in Spain, and seven in the United Kingdom. In addition, a Dichothrix tuft was collected from the Muga River, Spain. The samples were kept cold until reaching the laboratory, where they were divided into two parts, after washing them and checking by microscopy that there were no cyanobacterial contaminants outside the colonies. One part from each sample was used for morphological characterization and the other frozen and stored at − 20 °C until DNA extraction for genetic characterization. The colonies or tuft were named after the river or stream where they were collected, followed by a number (Table 1).

Morphological characterization.
Samples were inspected under a Leica MZ12.5 dissecting stereomicroscope (Leica, Leica Microsystems, Wetzler, Germany) equipped with epifluorescence and video camera systems and by an Olympus BH2-RFCA photomicroscope (Olympus, Tokyo, Japan) equipped with phase-contrast, epifluorescence, and video camera systems (Leica DC Camera, Leica Microsystems). To better observe the samples via microscopy, several colonies were decalcified with 0.5 M EDTA. Representative colonies of Rivularia-like samples were analyzed, except those with ambiguous morphologies, for which the variability of the morphotypes found within the colony precluded analysis. Size measurements were analyzed using SigmaScan Pro 5.0.0 software. Kruskal-Wallis one-way analysis of variance of ranks was used to determine the significance of differences between size measurements. Genetic characterization. Genomic DNA was extracted from individual colonies using a Power Biofilm DNA Extraction Kit (Mo Bio, Carlsbad, CA, United States) following the manufacturer's instructions, with a modification at the beginning of the protocol as previously described 28 to improve the lysis of the cells.
The phylotype composition of the cyanobacterial colonies was assessed by amplicon metagenomics targeting the hypervariable V3-V4 region of the 16S rRNA gene, using an Illumina MiSeq sequencer (Illumina Inc., San Diego, CA, USA) at the Genomic Service of the Universidad Autónoma de Madrid. First, PCR was performed using the cyanobacterial-specific primers CYA359F and 781Ra/781Rb 71 in separate reactions for each reverse primer as suggested by Boutte et al. 72 plus a targeted sequence. PCR was carried out with 1 ng of DNA, Q5 Hot Start High-Fidelity DNA Polymerase (New England Biolabs), and 100 nM primers in a 25 μl volume, and the cycling conditions were 1 × 98 °C for 30 s, 23 × 98 °C for 10 s, 54 °C for 20 s, and 72 °C for 20 s, and 1 × 72 °C for 2 min. The second PCR using barcoded primers for each colony was performed in a 20 μl final volume using the same DNA polymerase, 1 μl of the first PCR product, and 400 nM primers with the following cycling conditions: 1 × 98 °C for 30 s, 7 × 98 °C for 10 s, 60 °C for 20 s, and 72 °C for 20 s, and 1 × 72 °C for 2 min.
PCR products were quantified, pooled in equimolar amounts and purified using AMPure Beads (Beckman Coulter) prior to sequencing with a MiSeq sequencer (Illumina Inc., San Diego, CA, USA) at a read length of 2 × 300 bp. At least 100,000 sequences were obtained for each amplicon.
Sequence data were processed using QIIME v.1.9.0 73 as previously described 28 . The first taxonomical assignment of the OTUs was performed against the Greengenes database 74 by the classifier method of the Ribosomal Scientific Reports | (2020) 10:18453 | https://doi.org/10.1038/s41598-020-75303-2 www.nature.com/scientificreports/ Database Project with a confidence value of 0.8 75 , which allowed the removal of OTUs identified as chloroplasts or noncyanobacterial organisms. Afterwards, OTUs were matched against our sequence dataset of isolated cultures and field samples by employing the 'uparse ref command' in USEARCH 76 . In addition, a previously sampled Dichothrix field tuft from the Muga River was analyzed by cloning almost the entire 16S rRNA gene, as previously described 16 . Then, all this information was compared with the taxonomic assignments made against the Greengenes database 74 as mentioned above and against the SilvaMod database 77 using the lowest common ancestor (LCA) algorithm implemented in CREST 78 . In addition, the representative OTU sequences were compared with those in the NCBI database using a BLAST search (https ://www.ncbi.nlm.nih.gov/blast ), and sequences similar to them were downloaded to construct phylogenetic trees as previously described 27,28 . The first tree was constructed with only full-length 16S rDNA sequences (approximately 1500 bp), which were aligned using ClustalW in BioEdit 7.0.5.3 79 . After this, the OTU sequences (approximately 420 bp) were aligned against the alignment of the first tree, and a new tree was constructed and compared with the previous tree to ensure that the clades were maintained. Phylogenetic trees were computed using the neighbor-joining method 80 in Mega 7 81 with 1000 bootstrap replicates. Evolutionary distances were computed using the Tajima-Nei method 82 , and pairwise deletion was used to account for sequence-length variation and gaps. The percent identity between sequences was determined as (1-p-distance)*100. Because OTU1 appeared at markedly high abundance, amplicon sequence variant (ASV) determination was performed using the DADA2 pipeline in QIIME 2 83 , showing that ASVs assigned to this OTU differed by only one nucleotide and were therefore indistinctly distributed between the different samples.
Alpha diversity indices (Chao1 estimator, Good's coverage and observed OTUs) were calculated using QIIME. Good's coverage estimates reached 100% in all the samples, indicating that the large majority of the cyanobacterial diversity was captured.
The OTU sequences have been deposited in the GenBank database under accession numbers MT335702-MT335726. Raw sequencing data have been deposited in the NCBI Sequencing Read Archive under accession number PRJNA648107.