Water samples and macroscopic colonies of Nostoc spp. (ca. 10 cm in diameter) were collected from the littoral zone in Lake Chungará (Supplementary Figure S1A). This large lake (22.5 km²) is located in the Andean plateau (18°14′9.67 S, 69°10′53.81 W) at 4520 m above sea level and belongs to the Lauca National Park, a Unesco World Biosphere Reserve (Mühlhauser et al., 1995; Dorador et al., 2003). The littoral zone of this lake has an extensive area of macrophytes (e.g., Miriophyllum elatinoides) providing a habitat for a wide range of organisms including the endemic fish Orestias chungarensis Vila & Pinto and birds (e.g., Fulica gigantea Eydoux & Souleyet) that depends on this area for feeding and breeding (McFarlane, 1975; Vila and Pinto, 1986; Andrew, 1987). The colonies of Nostoc sp. are usually found at the surface of the dense submerged macrophyte belt in the littoral zone.

Samples for molecular analyses were collected from the inner part (homogenous GM) of the colony (Supplementary Figure S1B) using a sterile syringe and scalpel. A sample was obtained from one colony collected during the dry season (DS) in 2013 (sample Chungará_DS2013) and three samples were collected from three different colonies during the DS in 2016 (samples Chungará_DS2016_1, DS2016_2, and DS2016_3). A composite water sample (i.e., same volume pooled from 0.1, 0.6, and 1 m) from the littoral zone was collected with either a 2 L glass bottle (0.1 m depth) or with a 2 L Van Dorn sampler (for the other two depths) during the DS in 2013 (sample DS2013) and the wet season (WS) in 2014 (sample WS2014). Finally, one water sample was collected in 2016 during the WS (sample WS2016) and in triplicate during the DS (samples DS2016_1, DS2016_2 and DS2016_3). In addition, six colonies were collected during the DS in 2016 from the bed of a tributary stream of the Lauca River (Supplementary Figure S1C), located in “Quebrada Culco” (18°34′52.76 S, 69° 3′40.57 W, hereafter referred as to Culco stream) to compare the bacterial community of the inner GM with that of the OL. Due to the difficulty in sampling the two matrices without contamination, the sampling was done in separate colonies. Thus, three colonies were used to obtain the GM with the same methodology described above (samples Culco_DS2016_1, DS2016_2, and DS2016_3) and three other colonies were used to obtain the OL (samples Culco_DS2016_4, DS2016_5, and DS2016_6).

Water samples were kept in cold boxes and afterward (within ca. 2 h) were filtered onto 0.22 μm pore size filters (47 mm, Millipore GPWP), until clogging was observed. Filters and samples from Nostoc were placed in Eppendorf tubes with RNAlater (Qiagen, Germantown, MD). All samples were maintained at -20°C until further analysis.

Genomic DNA was extracted using PowerBiofilm DNA Isolation kit (Mo Bio Laboratories Inc.) following the manufacturer’s protocol. The concentration and quality of DNA were measured with a Nanodrop spectrophotometer (Nanodrop 8000, Thermo Scientific). Illumina Miseq sequencing was used with two different set of primers. First, total DNA from samples in 2013 and 2014 was used as a template for the V4 region amplification of the 16S SSU rRNA with the primers 515F (GTGCCAGCMGCCGCGGTAA) and 806R (GGACTACHVGGGTWTCTAAT) (Caporaso et al., 2011), done at the Research and Testing Laboratory Genomics (Lubbock, Texas, United States). Second, samples from 2016 were used to amplify the V4-V5 region of the 16S SSU rRNA with the primers 515F-Y (GTGYCAGCMGCCGCGGTAA) and 926R (CCGYCAATTYMTTTRAGTTT) (Parada et al., 2016), done at LGC Genomics Gmbh (Berlin, Germany). The 515F-Y/926R primer improves the underestimation of SAR11 clade and the overestimation of Gammaproteobacteria produced by the 515F/806R primer (Parada et al., 2016). Raw amplicons reads were deposited in the sequence read archive (SRA) of NCBI under accession number SRP136789 and SRP136788.

Raw reads from 16 samples were analyzed using Mothur (v. 1.35.1) following the standard operating procedure (Schloss et al., 2011). Briefly, paired-end reads obtained with the primers 515F-Y/926R and 515F/806R were assembled using the USEARCH v.7 (Edgar and Flyvbjerg, 2014) and make.contig command, respectively. After pooling all samples and trimming the reads to the same region (V4), reads were aligned to the SILVA v.132 database using the align.seqs command. Chimeras were detected and removed using UCHIME. The SILVA v132 database was used to classify reads with a confidence threshold of 80%. The remove.lineage command was used to identify and remove mitochondrial, chloroplasts, Archaea, Eukarya, and unknown contaminants. Reads were assigned to operational taxonomic units (OTU) at the 1% level of divergence using the cluster.classic command. All OTUs with less than six reads across all samples were discarded. Samples were normalized by randomly subsampling to the same size according to the sample with the smallest number of reads. After quality control, all samples were normalized to 27157 reads. The final OTU table and sequences are available in https://doi.org/10.6084/m9.figshare.7314875.v2.

To assess alpha diversity of the bacterial communities, the Simpson and Shannon indices that combine measures of richness and abundance were calculated on equal-sized samples using the INEXT package in R (Chao et al., 2014; Hsieh et al., 2016). A Kruskal-Wallis test was made to check for significant differences (P < 0.05) in alpha diversity among samples grouped by sampling site. The VEGAN package (Oksanen et al., 2013) was used to do the ordinations (metaMDS) based on Bray Curtis distance using the Wisconsin square transformation of the OTU relative abundances and to test for significant differences among samples grouped by sampling site in the ordinations using ANOSIM. A maximum likelihood phylogenetic tree, using the general time reversible model (with gamma distribution and bootstrap), was constructed in RAxML v0.6.0 (Kozlov et al., 2018) with the Nostoc reference species (with available 16S rRNA gene) present in the database Taxonomy from the National Center for Biotechnology Information (NCBI) including one sequence from Nostoc sp. (Llayta) reported for the Andean plateau. Then, the short reads (OTUs) belonging to Nostoc spp. were mapped onto the phylogenetic tree. We also repeated the phylogenetic analysis using the database CyanoPhy (cyanobact.000webhostapp.com), however, we obtained the same results (data not shown) and thus, we present only those from the first analysis.

Calculations of the OTUs relative abundance were made excluding the OTUs classified as Nostoc. Samples were grouped according to the littoral water, OL from Culco and inner GM from Lake Chungará and Culco stream. A Venn diagram was created to compare genera among samples grouped by site and sample origin, based on a presence/absence matrix and visualized by the package VennDiagram in R (Chen, 2018).

The functional prediction of the bacterial communities was assessed using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt; Langille et al., 2013) using the galaxy server¹. Briefly, the raw reads were re-analyzed in Mothur (v. 1.35.1) using the Greengenes v.13.5 database. All OTUs were used in downstream analyses, without including those classified as Nostoc. A biom file was produced with the final data and used as input for PICRUSt, where normalization by copy number, metagenome prediction and categorization by function was done. Further, the nearest sequenced taxon index (NSTI) was calculated to quantify the availability of nearby genome representatives for each microbiome sample. A low NSTI value (<0.1) indicates that the samples are highly supported by the reference microbial genome dataset (Langille et al., 2013). A principal component analysis of the predicted functions was made using STAMP v. 2.1.3 (Parks et al., 2014).

In situ measurements of water temperature and pH were done with a portable pH meter and coupled thermometer (HI9126, Hanna Instruments), whereas electrical conductivity was measured with a portable conductivity meter (Orion Star A322, Thermo Scientific). Samples were also collected in parallel for the analysis of major cations (potassium, sodium, calcium, and magnesium) by ion chromatography, and the anion, chloride was performed by the argentometric method, and sulfate by the gravimetric method (Gros, 2003). During 2013 and 2014 in L. Chungará, samples were collected in precombusted (4 h at 450°C) glass bottles for the analysis of dissolved organic carbon (DOC) and dissolved nitrogen (DN). These samples were filtered in situ through two pre-combusted GF/F filters (Whatman). The filtrate was acidified with HCl (pH 2) and analyzed later at the laboratory in Innsbruck, Austria with a Shimadzu TOC-Vc series equipped with a total nitrogen module. The instrument for DOC analysis was calibrated with potassium hydrogen phthalate, while calibration for the DN was done with potassium nitrate. Three to five subsamples were analyzed for each sample and for a consensus reference material (CRM) for DOC (batch 5 FS-2005: 0.57 mg; provided by RSMAS/MAC, University of Miami) that was run in parallel on each occasion. Results differed from the CRM given value by 5%, and the coefficient of variation among subsamples was <2%.

Both alpha diversity metrics (Supplementary Figure S2) were higher for the littoral water (Shannon = 39.3 ± 13.5; Simpson = 13.9 ± 6) and for the OL (Shannon = 25 ± 4; Simpson = 12.7 ± 2.2) than for the inner GM (Lake Chungará: Shannon = 2.8 ± 1.1; Simpson = 1.5 ± 0.05 and Culco stream: Shannon = 3.6 ± 3; Simpson = 1.6 ± 0.8). Further, Shannon and Simpson indexes were significantly different between the GM and the littoral water (Kruskal-Wallis test, P = 0.004), as well as between the GM and the OL (P = 0.011). In contrast, no significant differences in diversity were found when samples of the GM were compared between Lake Chungará and the Culco stream (Kruskal-Wallis test, P > 0.05). An ordination analysis (based on Bray-Curtis dissimilarity) revealed that bacterial communities from the colonies collected in the Culco stream (GM and OL) were grouped close together and that differences among samples from the littoral water and the inner GM of the colonies collected from both sampling sites existed (Figure 1A). Significant differences in OTUs relative abundance (ANOSIM; R² = 0.985, P < 0.001) were found among all samples. This was also true (ANOSIM; R² = 0.875, P < 0.001) even after the OTUs classified as Nostoc were removed from the ordination analysis (Figure 1B). Overall, these results imply that the environment inside the colony selects for the bacteria probably originated from the littoral zone and included in the colony during its morphogenesis. Further, the community in the inner GM of Nostoc sp. from Lake Chungará appeared to be stable over time, at least at the level of taxonomic resolution analyzed (Figure 1), implying that environmental conditions inside the colony are relatively constant.

One environmental factor that is obviously different inside and outside the colony is light intensity and probably also its spectral quality. In fact, in N. sphaeroides, adaptation to low light levels inside the colony is clear when pigments concentration and photosynthetic performance of filaments from the inner and OL are compared (Deng et al., 2008). For example, inner filaments have a lower light saturation point, lower photosynthetic rates and efficiency, but higher chlorophyll a and phycobiliproteins concentrations than those from the OL (Deng et al., 2008). The large size of Nostoc colonies, such as those from Lake Chungará imposes constraints in the uptake of external resources and concentrations of inorganic carbon are probably also limiting inside the colony (Sand-Jensen, 2014). It remains to be tested how these differences affect the associated bacterial community.

A total of 379 bacterial genera (24 phyla) were identified, but only 36 were shared among the littoral water, the OL and GM (Figure 2). Further, the littoral water showed the highest number of unique genera (n = 120) and a high number of shared genera (n = 98) with the GM matrix from Lake Chungará. However, few genera (n = 5) were shared between the GM from L. Chungará and Culco stream supporting the finding that colonies of Nostoc spp. from these ecosystems hold also different bacterial communities (Figure 1). This could be related with the probable existence of two different Nostoc species (see below discussion on identity) or with differences in bacterial community composition between the lake and the stream water included in the colony during morphogenesis, although, this needs to be tested.

Unclassified sequences members within the Rhodobacteraceae (Alphaproteobacteria) and Burkholderiaceae (Betaproteobacteria) families were abundant (3.2–20% of relative abundance) among samples obtained from the GM and OL (Figure 3). Further, all samples from the GM of Lake Chungará showed a high relative abundance of Porphyrobacter sp. (Alphaproteobacteria; 2.6–22.3%), Flavobacterium sp. (Bacteroidetes; 7–25.1%) and Emticicia sp. (Bacteroidetes; 2.6–7.4%).

The GM of Culco samples showed a high relative abundance of Sphingorhabdus sp. (Alphaproteobacteria; 10.6–35%), Rhizorhapis sp. (Alphaproteobacteria; 5.6–26.7%) and unclassified Rhizobiales (Alphaproteobacteria; 2–13%). The Burkholderiaceae Family (abundant in all Nostoc spp. samples), Rhizorhapis sp. and members of Rhizobiales include an extremely diverse group of Betaproteobacteria capable of nitrogen fixation (Coenye, 2014). In addition, some of the bacterial taxa detected in the inner GM have been described for aquatic environments along the Andean plateau. For example, the order Rhodobacteriales, in our study represented by the Family Rhodobacteraceae, has been detected in association with other microorganisms in microbial mats and water samples (Dorador et al., 2013).

The main difference between the GM and OL from Culco colonies, was given by the high relative abundance of the Alphaproteobacteria FukuN57 (4.5–17.8%) in samples from the OL and the high relative abundance of the Alphaproteobacteria UKL13-1 (3.1–4.2%) in the inner GM. Interestingly, the bacterial taxa associated with Nostoc spp. were also common to those found in other cyanobacterial associations. For example, in Microcystis, a bloom-forming genus that produce mucilaginous colonies, Porphyrobacter, Rhodobacterales, Sphingomonadales, and Burkholderiales are also typically associated (Shi et al., 2009, 2012). Further, cyanobacterial associations with Flavobacterium and members of Sphingomonadaceae, Burkholderiales and Rhizobiales have been described from metagenomes of culture collections belonging to different cyanobacterial genera (Cornet et al., 2018). Similarly, Porphyrobacter and some members of the Family Rhodobacteraceae are known from associations with the cyanobacteria Microcoleus sp. (Sánchez et al., 2005), Cylindrospermopsis sp. (Shi et al., 2009) and Oscillatoria brevis Kützing (Hube et al., 2009). One of the groups found in all colonies of Nostoc spp. was an unclassified Burkholderiaceae. This family includes Burkholderia, which is not only found in association with other cyanobacteria, but also with Mimosa species (Angiosperm) in a nitrogen fixation symbiosis (Bontemps et al., 2010).

Despite that in the littoral water from Lake Chungará, there were no clear differences in the main physicochemical variables among seasons (Supplementary Table S1), community composition changed among samplings, though it was always clearly different from that in the inner GM (Figure 1, 3). The Hgcl clade (Actinobacteria; 3.6–13.8%) and Flavobacterium sp. (Bacteroidetes; 1.8–14.5%) were the predominant (high relative abundance) groups in all water samples. In addition, a high relative abundance of Algoriphaghus sp. (Bacteroidetes; 16.7–34%) was observed in water samples from 2016 and unclassified Burkholderiaceae (Betaproteobacteria; 10.3–10.4%) in water samples during 2013 and 2014. Comparing with the bacterial community composition from the pelagic zone of Lake Chungará, the same predominant phyla were found (Aguilar et al., 2018), though differences were observed at the genus level. For example, while Flavobacterium sp. (Bacteroidetes) was the dominant genus in the pelagic zone only during the DS (Aguilar et al., 2018), it was always abundant (1.8–14.3% of relative abundance) in all seasons at the littoral zone. Further, Algoriphagus sp. was the most relative abundant taxa in the littoral zone of Lake Chungará during the DS in 2016, but it has not been detected in the pelagic zone (Aguilar et al., 2018). Only Actinobacteria (hgcl clade) occurred in both pelagic and littoral zones at a high relative abundance.

The PICRUSt analysis showed a low mean NSTI value for all samples (NSTI = 0.07 ± 0.02) indicating that the predicted metabolic functions in our study were highly supported by the reference microbial genome dataset (Langille et al., 2013). The mean predicted metabolic functions of the bacterial community from the littoral zone in Lake Chungará separated well from those of the community in the inner GM (Figure 4). Thus, these results support the idea that bacteria inside and outside the colony of Nostoc spp. differ not only in their composition, but probably also in their physiology. However, the predictions made by PICRUSt has clear limitations (Langille et al., 2013), namely, that they are made based on the comparison between short-reads and reference genomes and thus, their interpretation should be done with caution.

The Otu00001 and Otu00542 found in Lake Chungará were classified as Nostoc commune, whereas Otu00002, Otu00822 and Otu00838 found in the Culco stream were classified as Nostoc sp. indicating that they correspond to different species (Table 1). The Otu00001 and Otu00002 were most abundant (up to 98.5%) in Lake Chungará and Culco stream, respectively (Supplementary Table S2). The existence of two different species is also supported by the phylogenetic analysis, namely, that colonies from Lake Chungará are tentatively assigned to Nostoc commune, whereas those from Culco to N. flagelliforme (Supplementary Figure S3). Surprisingly, the sequence of Nostoc spp. from our study were not related to Nostoc sp. (Llayta), a sequence retrieved from the Andean plateau (no clear isolation source) and reported as the typical species in this area (Galetovic et al., 2017). Although the 16S rDNA gene is valid for phylogenetic studies of Cyanobacteria (Oksanen et al., 2004) and the Nostoc species separate well with other close cyanobacterial taxa (Svenning et al., 2005), we cannot confirm the species only based on a partial short sequence (253 bp).