A modified isolation technique, originally developed for planktonic cyanobacteria, was used, applying a glass capillary under sterile conditions (Zapomělová et al., 2008; Mareš et al., 2015). Colonized areas of the gypsum were scraped off and rehydrated on a microscopy slide in TE buffer (pH = 8), and then mechanically crushed with another microscopy slide. The prepared specimen was observed under the Olympus IX71 inverted microscope and searched for single colonies of selected species. The single colony was separated by repeated washing through a series of droplets of TE buffer, and documented by microphotography. Subsequently, the isolated colony was either directly put into a sterile 0.2 mL Eppendorf tube for further processing (Nostoc sp., Chroococcidiopsis sp.), or gently crushed using a clean cover slip, with five single cells of the target species collected per tube (Gloeocapsa spp., Gloeocapsopsis pleurocapsoides). After being kept frozen overnight at −20°C, total genomic DNA from the isolated cells/colonies was amplified using the Multiple Displacement Amplification (MDA) method utilizing Repli-g Mini Kit (Qiagen). The MDA was performed following the manufacturer’s instructions, which involved cell lysis, DNA amplification by Phi 29 polymerase (16 h at 30°C), and a final denaturation step (5 min at 65°C). The quality of MDA products was tested by gel electrophoresis, performed for 75 min at 80 V using 1.5% agarose gel. Successful MDA products were used as a template for PCR amplification of the 16S rRNA gene and the adjacent internal transcribed spacer (ITS) region using the following primer combinations (forwards/reverse): 16S378F/16S1494R, and 16S27F/23S30R (Taton et al., 2003). The quality of the PCR products was tested by gel electrophoresis; performed for 30 min, at 90 V, using 1.5% agarose gel. PCR products were then purified using an Invisorb^® Fragment CleanUp Kit (STRATEC Molecular GmbH, Germany) or a Monarch PCR and DNA Cleanup Kit (BioLabs^® Inc., New England), and then sent for commercial Sanger sequencing at SeqMe (Dobříš, Czech Republic) using the same primers as for PCR, plus an additional internal sequencing primer (cyano6r) (Mareš et al., 2013).

Coccoid green algae were mostly present in orange subsurface zones of some gypsum samples. Endolithic algae were isolated from crystalline gypsum rock from Monte Gibliscemi (Site 4). The crystalline gypsum from this site was selected due to the highly abundant and easily recognized orange-pigmented algal cells, without the presence of any other autotrophic species.

DNA was extracted from the orange colonization using the commercial kit DNeasy^® PowerLyzer^® PowerSoil^® Kit (Qiagen) as described in the provided protocol. After being briefly vortexed again, the samples were centrifuged at 12,000 g for 2 min, the supernatant was used for testing DNA quality (isolate S14 Sicily) on NanoDrop^®ND–1000 Spectrophotometer (NanoDrop Technologies, Inc.). The 18S rDNA region was then amplified by PCR using the forward primer FC (GGGAGGTAGTGACAAIAAATA), and the reverse primer RF (CCCGTGTTGAGTCAAATTAAG) (Matsuzaki et al., 2015). Each 21.44 μL of PCR reaction for 18S rDNA amplification contained 2 μL of DNA isolates (diluted to a concentration of 5 ng.μL^–1), 4.32 μL 5x MyTaq Red Reaction Buffer (Bioline, Meridian Bioscience, United States), 0.43 μL of each 10 μM primer, 14.04 μL sterile Milli-Q water, and 0.22 μL of 5 U.μL⁻¹ MyTaq HS Red DNA polymerase (Bioline, Meridian Bioscience, United States). Amplification reactions were performed using the following cycle parameters: for 3 min initial denaturation at 95°C, followed by 35 cycles (denaturation for 15 s at 95°C, annealing for 30 s at 56°C, extension for 40 s at 72°C), and final extension for 7 min at 72°C. The quality of the PCR products was tested by gel electrophoresis (45 min, 120 V, 1.5% agarose gel). PCR products were purified and sequenced using an Applied Biosystems automated sequencer (ABI 3730xl) at Macrogen Europe (Amsterdam, Netherlands). Chromatograms of the obtained sequences were examined using the FinchTV 1.4.0 (Geospiza, USA) program. Only sequences showing distinct single peaks in the electropherogram were used.

18S rDNA alignments were constructed for the phylogenetic analyses. The sequences were selected according to Gálvez et al. (2021) to encompass the closest relatives, to isolate S14 Sicily and other lineages in Chlamydomonadales. The 18S rDNA alignment contained 56 sequences (1,682 bp); the representatives of clade Polytominia and Dunaliellinia were selected as the outgroup. The best-fit nucleotide substitution model was estimated by jModeltest 2.0.1(Posada, 2008). Based on the Akaike information criterion, the GTR + I + G model was selected for 18S rDNA. The phylogenetic trees were inferred by Bayesian inference (BI) and maximum likelihood (ML) according to Nedbalová et al. (2017), with a minor modification that the Markov chain Monte Carlo runs were carried out for 3 million generations in BI. Bootstrap analyses and Bayesian posterior probabilities were performed as described by Nedbalová et al. (2017).

All obtained sequences were compared to those in GenBank by BLAST and were submitted to the NIH genetic sequence database GenBank at NCBI during paper revision under accession numbers OQ509475 − OQ509403 and OQ520113.

At this site the major organism identified in the black layer was Nostoc sp., with the dominant pigment scytonemin. Detected Raman bands of scytonemin were 1,712 (weak intensity), 1,633 (medium), 1,600 (strong), 1,558 (m), 1,387 (w), 1,324 (w), 1,174 (ms), 1,098 (mw), 1,023 (w), 755 (w), and 681 (w) cm⁻¹. The most significant Raman band of scytonemin was detected at 1,600 cm⁻¹, corresponding to the ν(CCH) aromatic ring quadrant stretching mode. The band at 1,558 can be assigned to the ν(CCH) p-disubstituted aromatic ring, ν(CCH) p-disubstituted aromatic ring, and the characteristic band at 1,174 cm⁻¹ of the ν(C=C–C=C) system (Edwards et al., 2000).

Gloeocapsa compacta was another major organism in the black layer with the dominant pigment gloeocapsin. Its Raman bands were detected at 1,671 (m), 1,289 (m), and 480 (mw) cm⁻¹. The analysis of Gloeocapsopsis pleurocapsoides cells revealed a strong carotenoid signal at 1,516 (vs), 1,449 (w), 1,388 (w), 1,283 (mw), 1,189 (mw), 1,155 (ms), 1,004 (m), and 957 (w) cm⁻¹. The band at 1,516 cm⁻¹ represents the C=C stretching mode of the conjugated chain in carotenoids. The band at 1,155 cm⁻¹ band represents the C–C stretching mode (coupled with C–H in-plane bending); and the band at 1,004 cm⁻¹ band is attributed to the C–CH₃ rocking mode.

Gloeocapsa compacta was the dominant species in the black layer, and the pigment gloeocapsin was detected as the major pigment of this cyanobacteria. The detected Raman bands were at 1,668 (s), 1,571 (s), 1,424 (ms), 1,340 (m), 1,280 (ms), 1,190 (m), and 465 (m) cm⁻¹. The bands around 1,668 cm⁻¹ can be assigned to ν(C=O) vibration, the bands around 1,571 cm⁻¹ to ν(C=C) aromatic ring vibration, and the bands around 1,276 cm⁻¹ to δ(CH₂) vibration. The signal of carotenoid pigments was also detected in G. compacta, as corroborated by its Raman bands at: 1,513 (vs), 1,283 (w), 1,193 (m), 1,155 (s), 1,003 (m), and 958 (w). The pigment scytonemin was detected as a major pigment in Nostoc sp. cells: 1,711 (w), 1,633 (m), 1,598 (s), 1,555 (m), 1,381 (mw), 1,171 (ms), 573 (w), and 436 (w) cm⁻¹ (see Figure 5B).

At this site the orange layer was colonized by orange-pigmented coccoid cells of Chlorophyta as corroborated by the strong signal of carotenoid pigment that was detected with bands at 1,517 (s), 1,447 (w), 1,389 (w), 1,274 (w), 1,191 (ms), 1,156 (s), 1,005 (m), and 964 (w) cm⁻¹ (see Figure 5C).

We obtained a total of 5,424 unique bacterial ASVs of the 16S rRNA gene, with 100–805 (354 on average) ASVs per sample. Taxonomic assignment of the ASVs to bacterial phyla consistently showed that the dominant group was Cyanobacteria (65.9% on average), followed by Proteobacteria (12.6%), Chloroflexi (4.5%), Bacteroidota (4.2%), Actinobacteriota (3.6%), Acidobacteriota (2.5%), Planctomycetota (2.3%), and Verrucomicrobiota (0.9%)—the latter showing a higher abundance in crystalline gypsum samples (Figure 7). Two families of coccoid cyanobacteria were the most abundant across the majority of the samples: Chroococcidiopsidaceae and Thermosynechococcaceae (Figure 7). Members of Leptolyngbyaceae, Nostocaceae, and Gloeobacteraceae were also generally present, and each of those groups was dominant in several individual samples.

Of the two localities for which sequencing data were obtained, Sicualiana Marina contained both types of the lithic substrate; whereas Santa Elisabetta contained only the selenite substrate. However, both localities included sets of green and black colored endolithic colonizations. Therefore, the data set was only partially representative with regard to the two types of gypsum; and the effect of the substrate type could only partly be separated from any unknown effects of the sampling locality such as, for example, the differences in the nutrient availability. Nevertheless, when we looked at the taxonomic composition between the two localities, the only striking difference was the dominance of Gloeobacteraceae in two (black and green) analyzed selenite samples collected from the Santa Elisabetta site. The two localities did not exhibit a significant difference in average biodiversity (Shannon index; p > 0.1). The samples from Santa Elisabetta were more uniform (Supplementary Figure S1), which was not surprising given their lower number and occurrence of only selenite.

In terms of the type of substrate, selenites on average contained more Nostocaceae, Leptolyngbyaceae, and Gloeobacteraceae (p < 0.05); but less Thermosynechococcaceae (p < 0.01) when compared to crystalline gypsum (Figure 7). The composition of the crystalline gypsum samples was more uniform, and all of them were only dominated by Chroococcidiopsidaceae and Thermosynechococcaceae. However, in the crystalline gypsum samples, Gloeobacteraceae and Coleofasciculaceae sequences were always present, while missing in some samples of the selenites. In summary, these patterns could be expressed as a significantly higher average biodiversity (Shannon index; p < 0.01), but a lower overall variability (Supplementary Figure S2) in crystalline gypsum versus selenite samples.

The dark-pigmented endolithic zones on average contained more Nostocaceae (filamentous cyanobacteria with colored sheaths) compared to the green-pigmented zones (p < 0.05) (Figure 7). On the other hand, Thermosynechococcaceae (coccoid cyanobacteria with colorless envelopes) were more abundant in the green areas (p < 0.01). Chroococcidiopsidaceae (coccoid cyanobacteria with pigmented or non-pigmented envelopes) were present in high abundance within both of these differently colored zones. The biodiversity was not statistically different between the black-pigmented and green zones (Shannon index; p >0.1) however, the samples from green endolithic zones were more homogenous (Supplementary Figure S3).

Significant differences in the overall ASV composition among the endolithic assemblages from selenite vs. white crystalline-gypsum, samples from green vs. dark layers, and samples from the two localities with sequencing data were further supported by NMDS analysis (Supplementary Figures S1–S3); although the variance between the localities was less pronounced. In permutational ANOVA tests, the type of substrate, the color of the endolithic layer, and the locality significantly explained 12% (p < 0.001), 10% (p < 0.001), and 9% (p < 0.01) of the observed variability, respectively.

Phylogenetic analysis allowed us to link the majority of the dominant ASVs discovered by metagenome profiling to cyanobacterial morphotypes observed by light microscopy.

The largest subgroup of dominant ASVs were classified to Thermosynechococcaceae (see Figure 8 for the list of ASVs), clustering in a deep-branching monophyletic lineage together with three of our single cell isolates of the “Chroococcidiopsis” morphotype. Another dominant group of ASVs, classified as Chroococcidiopsidaceae, exhibited considerable taxonomic diversity. It contained ASVs that could be assigned to several species within the lineage of Gloeocapsa spp. and Gloeocapsopsis pleurocapsoides (ASV9, 12, 13, 42, 55, 103, and 1,235), as established by our single cell isolates. Furthermore, it contained two abundant ASVs (8 and 23) clustering in the Chroococcidiopsis sensu stricto clade defined by the sequence of its type species, C. thermalis, two ASVs (36 and 81) in the Sinocapsa lineage, one very abundant “Chroococcidiopsis sp.” of the CC1 lineage (ASV1), and a single ASV (58) related to Myxosarcina. These results demonstrate extensive cryptic diversity of the coccoid endolithic cyanobacteria in the samples, especially regarding the Chroococcidiopsis-like morphotype.

Another very abundant ASV (2), classified as Gloeobacteraceae, was resolved as a sister lineage to Gloeobacter violaceus, in agreement with the observation of G. violaceus by light microscopy. The three most abundant ASVs of Leptolyngbyaceae were resolved in clades corresponding to Scytolyngbya (ASV53), Komarkovaea (ASV29), and “Leptolyngbya” sp. (ASV11); outside of the core Leptolyngbya cluster, showing the relatively broad diversity of simple filamentous cyanobacteria in the samples. These taxa were not captured by microscopic analysis (Table 2). A majority of ASVs classified as Nostocaceae were resolved in the Nostoc sensu stricto clade along with our sequences of Nostoc sp. single filaments identified by light microscopy. The remaining three abundant ASVs of heterocytous cyanobacteria clustered into genera Goleter (ASV16) and Komarekiella (ASV14), which are morphologically cryptic to Nostoc, and to Macrochaete (ASV43), a heteropolar cyanobacterium which was not found by microscopy.

Several taxa identified in the samples using light microscopy, i.e., Petalonema sp., Microcoleus sp., and Anathece sp., were not found among the dominant ASVs; however, sequences classified as Petalonema and Microcoleus were found in the samples as low-abundance ASVs.

Using culture-independent sequencing, the acquired Sanger sequence (385 bp long) matched with 98.96% identity to several closely related strains or environmental sequences that were all assigned to the order Chlamydomonadales (Chlorophyta), clade Stephanosphaerinia: “Chloromonas” sp. KSF063 (MH400034.1), “Chloromonas” sp. KNF032 (MH400032.1), “Chloromonas” sp. KNF030 (MH400031.1), “Chloromonas” sp. 435258601_3_18SR_D09 (MK330880.1), “Chloromonas” sp. 435258601_3_18SF_F09 (MK330879.1), uncultured alga gene Otu011 (LC371433.1), “Chloromonas” sp. KNF0032 (KU886306.1), Chlorominima collina CCAP 6/1 (MW553075), and Chlorominima collina CCCryo 273–06 (HQ404890.1).

Phylogenetic analyses of the 18S rDNA marker placed isolate S14 Sicily in the clade Stephanosphaerinia (Chlorophyta) in a basal position in relationship to the clade containing Balticola, Haematococcus and Stephanosphaera species, and to the clade of Chloromonas-like and Chlorominima species (Figure 9).

Our metagenomic analysis of gypsum from the southern part of Sicily showed the dominance of Cyanobacteria over other phyla such as Proteobacteria, Chloroflexi, Bacteroidota, Actinobacteriota, Planctomycetota, and Verrucomicrobiota. Previous studies on metagenomics showed that cyanobacteria significantly dominate as well, e.g., in gypsum from the Atacama Desert (Crits-Christoph et al., 2016). The other detected phyla have also been described from previous studies in calcite and ignimbrite from the hyperarid Atacama Desert (Wierzchos et al., 2015; Crits-Christoph et al., 2016), sandstone, limestone. and granite from the semiarid regions of the Rocky Mountain region (US) (Walker and Pace, 2007), and interestingly, some of them were also detected in an endolithic community associated with the coral Porites astreoides (Wegley et al., 2007). Cyanobacteria, Proteobacteria, Actinobacteria, and Acidobacteria were also found in gypsum from the Canadian Arctic (Ziolkowski et al., 2013). The similarities of the diversity on the level of phyla from various lithologies and geographical regions of the planet supports the hypothesis of a global metacommunity. This hypothesis says that endolithic ecosystems are comprised of organisms from a relatively limited reservoir of phylogenetic diversity and that these organisms are specifically adapted to their endolithic habitat (Walker and Pace, 2007). Nevertheless, their relative abundance may vary, for example, Chloroflexi and Proteobacteria were dominant over cyanobacteria in gypsum collected from the Lake St. Martin impact crater (Rhind et al., 2014). In halite from the Atacama Desert, Euryarchaeota were also more abundant compared to cyanobacteria (Wierzchos et al., 2018). As already mentioned, the variances of relative abundance can be related to the physical properties and inner architectures of the inhabited substrates. For instance, a higher taxonomic diversity was observed, in calcite than in ignimbrite (Crits-Christoph et al., 2016). This could be explained by the fact that calcite can form larger fissures and cracks compared to ignimbrite, which would lead to more competition for space and nutrient resources (Crits-Christoph et al., 2016). This also results in different metabolic pathways in these two substrates, genes for ABC transporters and osmoregulation are more diverse in the calcite habitat and genes for secondary metabolites are more abundant in the ignimbrite community (Crits-Christoph et al., 2016). The presence of some genes is also related to the nutrient availability in the environment. For example, no genes for diazotrophy have been detected in calcite and ignimbrite from the Atacama Desert where there are high nitrate accumulations from the atmosphere (Crits-Christoph et al., 2016) compare to. Besides that, the greater abundance of cyanobacteria could be related to the lower water availability, which means decreased photosynthesis and thus less CO₂ fixation, which is important for the development of heterotrophs (Wierzchos et al., 2018).

In general, in majority of endolithic ecosystems, cyanobacteria are the main primary producers, and the Chroococcidiopsis is frequently the major genus (Wierzchos et al., 2018). For example, in the dolomite in south China, nine cyanobacterial major phylotypes were found: Chroococcidiopsis sp., Phormidium autumnale, Anabaena oscillarioides, Scytonema sp., Leptolyngbya sp., Calothrix sp., two Nostoc sp., and Brasilonema octagenarum (Tang et al., 2012). In gypsum crust from Chott el Jerid, most cyanobacterial phylotypes were related to the genus Leptolyngbya within the Oscillatoriales and Nostocales (Stivaletta et al., 2010). The molecular studies show that the most abundant genomes found in the calcite and ignimbrite were also related to Chroococcidiopsis and Gloeocapsa (Crits-Christoph et al., 2016). In general, similar compositions have been observed all over the globe in many various substrates (e.g., gypsum crust, Stivaletta et al., 2010), limestone (Ferris and Lowson, 1997; Tang et al., 2012), which could also be explained by the hypothesis of global metacommunity. Also, it has been shown that the overall relative abundance tends to be greater in less extreme environments (Ziolkowski et al., 2013). Nevertheless, the authors acknowledge that it is possible that the higher diversity could be related to the limitations of the molecular methods used (Ziolkowski et al., 2013).

In the current study, we utilized a single cell/filament sequencing technique (Mareš et al., 2015; Kurmayer et al., 2018) to provide a direct link between cyanobacterial morphotypes observed using light microscopy and the results of metagenomic profiling. Such an approach is elaborate, and has only rarely been used in previous studies (Bowers et al., 2022), however, it provides valuable data leading to a better understanding of the microbial communities composition. Another similar, but even more challenging, approach is to link the metagenomic data to barcoded laboratory strains isolated from the same samples (Águila et al., 2021).

The most interesting result was a partial dissection of the cryptic variability in coccoid cyanobacteria that dominated the samples. Surprisingly, the main Chroococcidiopsis-like morphotype (Figure 3D) analysed by single-cell sequencing was resolved as Thermosynechococcaceae (Figure 8), together with a number of major amplicon sequence variants (ASVs). Some of the remaining frequent ASVs could be linked to the Chroococcidiopsis sensu stricto lineage, and even to other lineages of cyanobacteria forming sarcinoid packets of cells (“Chroococcidopsis” lineage CC1, Sinocapsa, “Myxosarcina”). The polyphyly of Chroococcidiopsis-like morphotypes have long been known (Cumbers and Rothschild, 2014); however, these taxa are still, for the most part, awaiting formal taxonomic treatment.

The pattern of Chroococcidiopsis-like phylotypes found in our study matches those previously reported worldwide from both hot and cold deserts (Bahl et al., 2011), except for the dominant clade in Thermosynechococcaceae, for which the morphological information had been missing up to now. Similar organisms were previously sequenced from endolithic limestone communities in the Rocky Mountains and the Mojave Desert, United States (Walker and Pace, 2007; Smith et al., 2014). Another environmental sequence in this clade was erroneously assigned to a true-branching heterocytous cyanobacterium Loriellopsis cavernicola (Lamprinou et al., 2011), probably due to laboratory contamination from the sample of calcite rock taken from a cave in Spain.

A complementary line of evidence was provided by single-cell sequencing of the abundant coccoid cyanobacteria containing dark UV-screening sheath pigments, Gloeocapsa spp. and Gloeocapsopsis pleurocapsoides. All of them clustered in the Chroococcidiopsidaceae clade together with a number of frequent amplicon sequence variants and previously published uncultured “Chroococcidiopsis sp.” sequences (Figure 8). Reliable sequences of these taxa have never been provided before; thus our results will help to distinguish these eco-physiologically different organisms in future metagenomic studies.

Another interesting finding was the relatively high proportion of Gloeobacter in the gypsum samples documented by both light microscopy and amplicon sequencing (Figures 3, 9). While it is known that this basal thylakoid-less cyanobacterium is more common than once thought (Mareš et al., 2013), its global occurrence is still understudied.

The remaining dominant cyanobacteria were classified to Nostocales (Nostocaceae) and Leptolyngbyaceae. Identity of the dominant Nostoc sp. morphotype (Figure 3C) was confirmed by single filament sequencing. However, also in this case, taxa morphologically cryptic to Nostoc (Goleter; Miscoe et al., 2016 and Komarekiella, Johansen et al., 2017) were discovered by the metagenomic analysis along with several other heterocytous cyanobacteria (Figure 8). Surprisingly, the thin filamentous cyanobacteria from Leptolyngbyaceae were completely overlooked by microscopy analysis. This is due to the fact that only limited proportion of the samples can be microscopically observed. On the other hand, some of the taxa seen in the microscope (Microcoleus, Petalonema, Anathece) were very scarce or missing in the amplicon sequencing data (e.g., due to better DNA extraction of some species). These findings illustrate the limitations of both approaches.

Considering the community composition, all samples analysed by metagenomic profiling were generally quite similar, and differed only in the proportion of individual groups (Figure 7). As expected, an increased occurrence of cyanobacteria with dark pigmentation (Nostocales) was confirmed in the macroscopically dark endolithic layers compared to the green layers dominated by Thermosynechococcaceae, in which the sheath pigments are absent.

The orange pigmented coccoid algae inhabiting fine-grained gypsum at Site 4 (isolate S14 Sicily) were found to be genetically related to green algae from the Chlorophyceae family originated from extreme habitats such as: Antarctic or Arctic freshwater (Jung et al., 2016; Cho et al., 2019), red snow on glacier and green snow fields at King George island (Gálvez et al., 2021), snowfields in Ryder Bay, Antarctic Peninsula (Davey et al., 2019), and Yukidori Valley, Eastern Antarctica (Segawa et al., 2018). Metagenomic sequencing of Arctic gypsum endoliths indicated the presence of few green algal sequences closely matching with Trebouxia sp. and Trichosarcina sp., known lichenizing symbionts (Ziolkowski et al., 2013). In other study, coccoid chlorophototrophs were responsible for gypsum layered orange to green cryptoendolithic colonization in Atacama, Northern Chile, but its sequences failed to be amplified by the environmental sequencing; likely due to the primer bias (Wierzchos et al., 2015). Another study on chasmoendoliths in Antarctica in granite and gypsum from the Atacama Desert showed the presence of Trebouxiophyceae (Yung et al., 2014; Ertekin et al., 2020). Further, known alga found in gypsum (spring associated) was the desmid Oocardium stratum (Zygnematophyceae, Streptophyta) (Rott et al., 2010). In our case study, direct DNA extraction isolated from gypsum generated a 18S rDNA sequence corresponding to those from the order of Chlamydomonadales, clade Stephanosphaerinia; mostly to genera of Chloromonas-like morphology and Chlorominima (Figure 9). However, the acquired sequence fragment was too short to allow any closer taxonomic assignment. Some of the related species to our Sicilian algal isolate were shown to form morphologically similar globular-shaped thick-walled cells with cytoplasmic lipid droplets in its life cycle that occupy most of the cell volume. The production of abundant lipid droplets is also correlated with the accumulation of carotenoids (Grünewald et al., 2001), and as mentioned above, carotenoids were also detected in our study by Raman spectroscopy.