Biocrusts were sampled in six sites distributed in four preserved areas in the SE region of Brazil (Supplementary Figure S1): Furnas do Bom Jesus São Paulo State Park of (one site), Vassununga São Paulo State Park (one site), Serra do Cipó National Park (two sites: Cipó and Capão), and Serra da Canastra National Park (two sites: Canastra and Zagaia) (Table 1). Ten sequential equidistant samples were collected throughout a 200 m long transect with two parallel transects set at each site, encompassing a total of 120 samples for the study (10 samples/transect × two transects/site × six sites). Soil crusts were collected with a Petri dish (55 mm × 1 cm high) and transported to the laboratory where they were kept dry at −20°C until processing.

Soil grains, rocks, and organic matter (plant roots and leaves) were manually removed, and 1 mg of the sampled biocrust was used for DNA extraction with MoBio Powersoil kit (Laboratories Inc., Carlsbad, CA, USA) according to the manufacturer’s instructions. The region V3-V4 of 16S rRNA gene (420 bp) was amplified with CYA 359 and 781a/b primers as described in Nübel et al. (1997), with an overhang Illumina adapter included in the primers. The PCR reaction contained 10 ng of eDNA, 0.2 μM of each primer, and 1x KAPA HiFi HotStart Ready Mix (KAPA Biosystems) to 25 μl of final volume. The PCR product from the samples 1–3, 4–6, and 7–10 of each transect were pooled and purified using AMPure XP purification kit (Beckman Coulter Inc., Brea, CA, USA). Afterwards, Illumina sequencing adapters and dual-index barcodes were added to the amplicon target using the Nextera XT Index Kit (Illumina, USA), according to the manufacturer’s instructions. The product was purified using AMPure XP purification kit, quantified with Qubit Fluorometric Quantitation (Thermofisher/Life Technology, USA). The samples were normalized and, then, pooled in an equimolar fashion. The preparation of the samples followed the Illumina guidelines for sequencing in a MiSeq platform (Illumina) available at Center of Nuclear Energy and Agriculture (ESALQ/USP) and using MiSeq Reagent kit v3. 2 × 300 cycle.

The 16S rRNA gene forward and reverse sequences were paired using PANDAseq (Masella et al., 2012) and a fastq file was generated. The “QIIME 1.9” [Quantitative Insights into Microbial Ecology (Caporaso et al., 2010)] was used for further analyses. The script “multiple_split_libraries_fasq.py” was run without quality filtering, as quality filtering was done before pairing using Trimmomatic (Bolger et al., 2014). Next, the script “pick_open_reference_otus.py”, clustered reads into 97% self-similar operational taxonomic units (OTUs) using SUMAClust (Schloss, 2016) and SortMeRNA (Kopylova et al., 2012) in combination with the Greengenes 13_8 database (DeSantis et al., 2006). A filter was applied to the OTU table and only the OTUs that appeared in at least two samples were considered in the analyses. The OTU table also was used to do a rarefaction curve using “Chao1” and “goods_coverage” methods, through the script “alpha_diversity.py”, also in “QIIME 1.9”. All sequence datasets are publicly available through NCBI under the project “Diversity and ecology of cyanobacteria of biological soil crusts in Brazilian Savannah” (NCBI identification number: SRP137259).

OTUs were identified by comparison with taxonomic information provided by public databanks. OTUs that presented more than 199 reads were first identified using QIIME (based on SortMeRNA and the Green Genes 13_8 database – DeSantis et al., 2006) and then compared one by one with GenBank data (NCBI nr/nt) using the tool “Basic Local Alignment Search Tool” (Baumann et al., 1990) to refine the initial identification. Cyanobacterial OTUs’ taxonomic assignment at the genus and species level was further informed through phylogenetic placement in the cyanobacterial reference database, Cydrasil (v. rc1, https://github.com/FGPLab/cydrasil). The Cydrasil rc1 database contains 1,161 curated cyanobacterial 16S rRNA gene sequences that are at least 1,100 bp long and includes a phylogenetic tree generated using RAxML 8 (Berger and Stamatakis, 2011). Query cyanobacterial sequences were aligned to the reference alignment with PaPaRa (Berger and Stamatakis, 2011), placed into the reference tree using the RaxML8 Evolutionary Placement Algorithm (Berger and Stamatakis, 2011) without changing the tree topology, and visualized on the iTOL 3 server (Letunic and Bork, 2016). This procedure allowed to relate an OTU sequence to a well-curated specific phylogenetic group or clade in a database that is enriched in cyanobacterial sequences from biocrusts, and thus confirm (or correct) the initial taxonomic identification. In general, OTUs that were 95% similar to an identified sequence (following Yarza et al., 2014) and placed in a highly supported and well-defined clade composed of coherently identified sequences were considered pertaining to the same genus. OTUs that presented less than 94% of similarity to the closest sequence were not identified. OTUs related to a single genus name but distributed in different clades (as the cases of the polyphyletic genera Leptolyngbya and Microcoleus) were considered to be effective genus-level taxa and given provisional identifiers (e.g., Leptolyngbya - Clade I, Leptolyngbya - Clade II…). To facilitate the visualization of the OTUs distribution and the composition of the sites, plots were constructed using the packages “ggplot2,” “scales” (Wickham, 2018a,b) and “reshape2” (Wickham, 2017) written in R language (R Core Team, 2018).

A meta-analysis of a set of 15 libraries of bacterial 16S rRNA gene sequences of biocrusts from North American deserts (Table 2; Velasco Ayuso et al., 2017: NCBI identification numbers PRJNA343817; Fernandes et al., 2018: NCBI identification number PRJNA394792) was used for comparisons. These libraries had been constructed using primers (515F and 806R primers - Caporaso et al., 2012), whereas those in this work used CYA 359 and 781a/b (Nübel et al., 1997). Because of this, the quality-controlled and paired sequence files were merged into a single FASTA file and imported together into the QIIME 22018.2 for analyses. Sequences were clustered at 97% similarity using closed reference OTU picking with VSEARCH (Rognes et al., 2016) with the Green Genes 13_8 database (DeSantis et al., 2006) providing the reference sequences. The resulting OTU table was filtered to only include cyanobacterial sequences. OTUs were then aligned using Mafft (Katoh and Standley, 2013) and a phylogenetic tree was generated using FastTree (Price et al., 2010). Community differences were assessed via permutational multivariate analysis of variance (PERMANOVA) performed on Bray-Curtis distance matrices of relative abundance derived from sequencing and used 9,999 permutations. PERMANOVAs were performed using the function “adonis2” in the <vegan> package (Dixon, 2003) run in “R” (R Core Team, 2018). The <vegan> function “betadispar” was used to test the variances (PERMDISP). A p of 0.05 was set as the significant threshold for all multivariate statistical analyses. Community composition was visualized with NMDS, using 25 restarts and 9,999 iterations.

Climatic data were obtained from the public database available at Center for Weather Forecasting and Climate Research (INPE/MCT; http://bancodedados.cptec.inpe.br/). Average annual precipitation (PRE), average annual high temperature (HT), average annual air humidity (AH), and altitude (ALT) were retrieved. Soil temperature (ST) and pH, which were measured in the field, completed the environmental dataset (Table 1). The packages “vegan” (Oksanen et al., 2016), “ggplot2” (Wickham, 2018b), and “gridBase” (Murrell, 2014), written in R language (R Core Team, 2016) were used to relate environmental data to the distribution and abundance of OTUs. After normalizing the OTU table, a redundancy analysis (RDA) was run as a constrained ordination method to search for possible spatial patterns in the cyanobacterial database.

Altogether, we detected 14,465 cyanobacterial OTUs, of which the 600 most frequently represented comprised 70% of total reads. Rarefaction curves showed that most samples (except CA1, ZA3, and ZB1) reached a plateau, and therefore most of diversity was accessed (Supplementary Figure S2). The samples that did not reach the plateau were still included in the subsequent analyses because they were similar in community composition with other samples from the same localities that did.

With our taxonomic identification constrained at the genus level, a large portion of the cyanobacterial diversity remained unassigned (Canastra 34.5%, Capão 19.8%, Cipó 46.5%, Furnas 9.4%, Vassununga 26.9%, and Zagaia 35.9%). But a majority of the unassigned OTUs also failed to align to sequences found within the NCBI database. Even among publicly available environmental sequences, such OTUs did not align with greater than 95% sequence identity, speaking for the presence of a significant level of biodiversity novelty in our biocrusts. At this level of resolution, community composition was relatively homogeneous among samples from the same site (Figure 1). OTUs assignable to various clades of Leptolyngbya had the most reads. They were dominant in Furnas, Cipó, Capão, and Canastra. In some Capão samples, sequences allied to Pycnacronema were also abundant. The most compositionally divergent communities within our set were those from Zagaia and Vassununga, where biocrusts were dominated by sequences assignable to Porphyrosiphon notarisii Kützing ex Gomont.

On the basis of OTUs (i.e., independently of taxonomic assignments and including unassigned OTUs), communities also differed significantly in composition (Supplementary Figure S3) among localities but not within them (PERMANOVA, Pseudo-F = 5.6904, df = 5/30, p < 0.0001), each locality being statistically different from each other in pairwise comparisons (p < 0.05). Again here, Zagaia and Vassununga seemed to be the most divergent cerrado localities along a compositional continuum. These differences, however, were much less marked than those found between North American desert communities and those from the cerrado as a whole (Figure 2; PERMANOVA, Pseudo-F = 16.442, df = 1/50, p < 0.0001).

The detailed phylogenetic placement using Cydrasil revealed the occurrence of several important clades of OTUs, which we describe in more detail here. The placement of major cerrado biocrusts inhabitants within the cyanobacterial radiation is given in Figure 3. Details for the different taxa are in Supplementary Figures S4–S6.

Leptolyngbya Clade I fell within the poorly defined, polyphyletic complex of sequences belonging to thin filamentous cyanobacteria that are usually attributed to this morphogenus and was formed by OTUs related to the public sequence of Leptolyngbya frigida (Fritsch) Anagnostidis and Komárek. It was clearly polyphyletic to the clade containing the sequence for the generic type species (L. boryana Anagnostidis and Komárek), so it likely represents a new generic unit. A second clade in this complex, Leptolyngbya clade II, was composed of OTUs similar to the sequence of Leptolyngbya sp. 7FUR (MF109116), but they were also, as Leptolyngbya Clade I, polyphyletic to true Leptolyngbya. Finally, OTUs placed within Leptolyngbya Clade III gathered around the sequence of the type species and are likely Leptolyngbya sensu stricto. A tree for the “Leptolyngbya complex” with assignments can be found in the Supplementary Figure S4.

A set of abundant OTUs with high representation in some of our localities was affiliated with Porphyrosiphon notarisii as judged by close similarity with several new sequences derived from bona fide cultures obtained from Furnas biocrusts (now permanently added to the Cydrasil database). This clade of sequences is well separated phylogenetically from sequences of other large-celled members of the Phormidiaceae, indicating that they do indeed represent a differentiated generic entity, as the traditional taxonomy would predict (Supplementary Figure S5).

A second complex of biocrust cyanobacteria corresponds to the epithet “Microcoleus steenstrupii” which has been recognized as a supra-generic entity in need of re-evaluation, not only because its members are not related to M. vaginatus (Vaucher) Gomont ex Gomont (the type species for the genus Microcoleus, which falls within the Complex Filamentous clade of Figure 3), but also because it encompasses a variety of separated clades with apparently diverging ecological traits (Fernandes et al., 2018). Sensu lato, the complex also includes sequences aligned with cultured strains assigned to M. paludosus Gomont, and the recently described Pycnacronema savannensis Martins, Machado-de-Lima, and Branco isolated and described from Brazilian savanna biocrusts (Martins et al., 2019). All these cyanobacteria are morphologically indistinguishable from Phormidium, except for the fact that they often form trichome bundles. Many OTUs in our crusts clustered within this complex, and given the difficulty in systematics, we have maintained the epithet “Microcoleus steenstrupii complex” to refer to them, except for OTUs clearly affiliated with P. savannensis, which conformed a major component in some of our crusts (Supplementary Figure S5). Several new full sequences of P. savannensis from bona fide cultures were included in the database before analyses to ensure that this choice was correct. OTUs in this complex not affiliated with Pycnacronema were not very numerous in terms of total reads.

OTUs most similar to the sequences of Microcoleus vaginatus, the type species for the genus Microcoleus in the Oscillatoriaceae, and likely the most common terrestrial cyanobacterium globally (Garcia-Pichel, 2009), while not very common here, constituted another clear assignment.

Within the heterocytous cyanobacteria, many OTUs were phylogenetically cognate with members of Brasilonema sp. (Supplementary Figure S6), previously unreported from soil crusts. Among the three major heterocytous types found in North American arid land crusts (Yeager et al., 2007), only Scytonema found significant representation in the cerrado crusts.

Our RDA analyses (Figure 4) showed that the combined effect of all environmental variables considered (Supplementary Table S1) could explained 58% of the total variability in cyanobacterial community composition among Brazilian biocrusts. The two statistical axes represented in Figure 4 explained 40% of this variability. The opposing vectors of high annual temperature (HT) and precipitation (PRE)/air humidity (AH) explained the majority of the variation in community composition among locations, particularly separating Zagaia and Vassununga (related to hotter, drier climate) from the rest. The arid end of this continuum was determined significantly by the strong contributions of Porphyrosiphon, the opposite end by those of Leptolyngbya Clade I. Variability among Canastra, Cipó, Capão, and Furnas communities seemed to be driven by soil pH (Figure 4).

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