Taxonomic and functional insights into rock varnish microbiome using shotgun metagenomics

ABSTRACT

Rock varnish is a microbial habitat, characterised by thin (5–500 μm) and shiny coatings of iron (Fe) and manganese (Mn) oxides associated with clay minerals. This structure is well studied by geologists, and recently there have been reports about the taxonomical composition of its microbiome. In this study, we investigated the rock varnish microbiome using shotgun metagenomics together with analyses of elemental composition, lipid and small molecule biomarkers, and rock surface analyses to explore the biogeography of microbial communities and their functional features. We report taxa and encoded functions represented in metagenomes retrieved from varnish or non-varnish samples, additionally, eight nearly complete genomes have been reconstructed spanning four phyla (Acidobacteria, Actinobacteria, Chloroflexi and TM7). The functional and taxonomic analyses presented in this study provide new insights into the ecosystem dynamics and survival strategies of microbial communities inhabiting varnish and non-varnish rock surfaces.

INTRODUCTION

The functional capabilities of microorganisms make them an integral component of geochemical processes in specific environments including rock varnish formation. Rock varnish is a microbial habitat, characterised by a thin (5–500 μm), shiny coatings of iron (Fe) and manganese (Mn) oxides associated with clay minerals (Dorn 2013). Varnish accretes on the rock surface in a time frame ranging between 25 and 300 000 years (Engel and Sharp 1958; Knauss and Ku 1980; Liu and Broecker 2000; Dorn 2007) and the rate of varnish formation varies with location, time and micro-climate conditions (Hayden 2006). The wide range of rock varnish implications in astrobiology, paleoclimatology, archaeology and quaternary studies in arid and semi-arid environments (Krinsley, Dorn and DiGregorio 2009; Eggleston and Parkinson 2012) have attracted scientists to investigate the genesis of such rock coating in a great detail (López-García et al. 2005; Abreu et al. 2011; Williams et al. 2013). Several theories have been established to explain both abiotic and biotic genesis of rock varnish (Hooke, Yang and Weiblen 1969; Allen et al. 2004; Mergelov et al. 2012). The theory of abiotic genesis assumes that the enrichment of Mn and Fe oxides is largely governed by factors such as pH fluctuation and abiotic oxidation (Hooke, Yang and Weiblen 1969). Alternative biotic theories suggest that the genesis of rock varnish is associated with the metabolism of microorganisms and their physiologies (Perry and Kolb 2004; Krinsley, Dorn and DiGregorio 2009; Northup et al. 2010a; Wang et al. 2011; Esposito et al. 2015). In the thin section of rock varnish, biofilms composed of microbial communities and extracellular polymeric substances (EPS) are arranged in finely repetitive layers that are wavy and discontinuous (Perry and Kolb 2004). Each cellular structure, metabolism-independent process and EPS associated with the microbial cell interact with the surrounding lithic substrate to form the rock varnish (Friedmann 1982; Amy et al. 1992; Wierzchos, Ascaso and McKay 2006; Kuhlman et al. 2006a; Tang, Beckage and Smith 2012). According to the biogenic origin of rock varnish, its formation starts when microorganisms mediate the accumulation of Mn and Fe oxides followed by cementation of microbial remains in the varnish layers (Dorn and Oberlander 1981; Krumbein and Jens 1981; Palmer et al. 1986; Drake, Heydeman and White 1993; Krinsley 1998).

While the microbial taxonomic diversity of rock varnish is well-studied (Kuhlman et al. 2006b; Northup et al. 2010b; Parchert et al. 2012; Zhang et al. 2012; Esposito et al. 2015; Lang-Yona et al. 2018), there is general lack of knowledge about the functional diversity (Lang-Yona et al. 2018). Rock varnish has been found in different environments, including recently deglaciated alpine areas where barren rock can present areas with varnish in close proximity with areas without varnish (Esposito et al. 2015). Previous studies suggested that endolithic microbes have an important role in pedogenesis (Mergelov et al. 2018). In this regard, we employed an integrated approach involving high-throughput metagenomic sequencing, analyses of elemental composition, lipid and small biomolecule biomarkers, and rock surface analyses to explore the in situ patterns of microbial diversity and potential functionality. By sampling varnish and non-varnish rock surfaces from the same bulk rock, with the same mineralogy and the assumption that inputs to the system and environmental conditions are roughly the same, we were able to investigate how varnish and non-varnish rock surfaces in close proximity may accommodate different microbial communities and biogeochemical processes.

MATERIALS AND METHODS

Study area and sampling

The study area is located in Matsch Valley (South Tyrol, Italy) (Fig. 1), a southeast-oriented valley with a continental climate. This valley has lower total annual precipitation compared to other mountain areas at the same elevation with a mean annual precipitation of 800–1000 mm at 2000 m above sea level (a.s.l), typically as snowfall from November to late April and rainfall in summer. The mean annual air temperature is 6.6°C at 1500 m a.s.l. (Penna et al. 2014). The geological substrate of the Matsch Valley consists of red schist and gneiss (Habler, Thöni and Grasemann 2009). The most common soil types are leptosols, regosols and umbrisols with a mean pH of 4.3. The whole study area is located above the tree line, at an altitude between 2400 and 2500 m a.s.l. where a rapidly retreating glacier has left >3 km² of foreland (Knoll and Kerschner 2010). According to photographical evidences and phytosociological analysis, the area where the sampling sites are located should have been left uncovered by the glacier between 1850 and 1926.

Rocks were sampled on June, 2013 from two sites (Fig. 1), where varnish surfaces which characterised by a brown to black shiny colour, were in close proximity to non-varnish surfaces that characterised by a red colour and absence of visible biota (mosses, lichens and visible biofilms). The distance between the two sampling sites is approximately 300 m. The geographic coordinates of the two sites are as following: site #1 (46° 46′ 38.85 N 10° 41′ 50.31 E) and site #2 (46° 46′ 35.12 N 10° 41′ 53.67 E). At each site, two rock chips of varnish and two rock chips of non-varnish were collected. The four rock chips per each site lay in an area of about 1.5 m². Each rock chip had a surface area of 0.2–0.3 cm². To ensure adequate sterile conditions in rock chip sampling necessary for molecular analysis, these chips were taken using flame-sterilized tweezers and immediately placed in 1.5 ml sterile tubes containing the first buffer solution and the crushing beads of the Ultra-Clean DNA extraction kit (MO BIO, Arcore, Italy). Tubes were directly put at −20°C until laboratory analysis. Additional rock chips of the same size as above from varnish and non-varnish surfaces were collected and stored at 4°C for geochemical analyses.

DNA extraction, metagenomic library construction and sequencing

DNA extraction was performed on the eight samples using the Ultra-Clean kit (MO BIO Laboratories, a Qiagen Company, Carlsbad, Ca, USA) according to the protocol of the manufacturer with minor modifications in order to maximize the yield and the quality of the extraction. The tubes containing the rock samples, the crushing beads and the buffer solution were placed into the grinding jars of a Retsch mill after a first round of crushing (30 s at max speed), we added once more the Ultra-Clean beads buffer (without the beads). This procedure was implemented to improve the quantity of the extract (Feinstein, Sul and Blackwood 2009). The following DNA extraction steps were executed in an UV-sterilized extractor hood. Initial test of the integrity of DNA extraction was done by loading 2 μL of eluted DNA in a 0.8% agarose gel with lambda marker (Thermo Fisher Scientific, Milan, Italy). DNA quantification was done using a Qubit fluorometer (Invitrogen, Carlsbad, CA). For metagenomic library construction from each sample, 2–10 ng of DNA was prepared with the Rubicon ThruPlex kit (Rubicon Genomics, Ann Arbor, Michigan, USA) according to the instructions of the manufacturer. Cleaning steps were performed using MyOne™ carboxylic acid-coated superparamagnetic beads (Invitrogen, Carlsbad, CA, USA). Constructed libraries were sequenced at the National Genomics Infrastructure (NGI; Solna, Sweden; https://portal.scilifelab.se/genomics/) on a MiSeq instrument (Illumina Inc., San Diego, CA, USA) with paired-end reads of 300 bp in length. On average, 27 million paired-end reads were generated per sample. These sequence data have been submitted to MG-RAST (Glass and Meyer 2011) with the sample number mgs 388 395–388 423.

Bioinformatics analyses

For taxonomical and functional assignment, we used MG-RAST with default settings for reads pre-processing, quality checking and homology searching. Taxonomy was assigned against Genbank database, whereas function were assigned using the SEED subsystems database (Overbeek, Disz and Stevens 2004). Tabular outputs, containing information on reads mapping to taxonomically or functionally assigned sequences were exported and parsed using the LEfSe (linear discriminant analysis effect size) algorithm. To investigate taxonomic and functional biomarkers discriminating between varnish and non-varnish rocks, an LDA score of 2.5 was used (Segata et al. 2011). Box-plots of biomarker taxa and Venn diagram of specialist and ubiquitous taxa were created using R software (R Core Team 2012).

Reads were co-assembled using MEGAHIT (Li et al. 2016), and the resulting contigs imported into the Anvi'o genomics visualisation platform for subsequent analysis (Eren et al. 2015). Briefly: (i) the contig database containing information about the contigs such as k-mer frequency, ORF and taxonomy were inferred by centrifuge (Kim et al. 2016), (ii) original raw reads were mapped on the assembled contigs using Bowtie2 (Langmead and Salzberg 2012) and (iii) separate profiles were merged and contigs binned using CONCOCT (Alneberg et al. 2014). After binning, each bin was manually checked to ensure that the completion level (calculated as the number of single copy genes occurring only once) was >90% and redundancy level was below 10% to improve the quality of the metagenome-assembled genomes (MAGs), as described in the MIMAG guidelines (Bowers et al. 2017). A putative taxonomy was assigned to the bins using PhyloPhlAn (Segata et al. 2013), considering the average nucleotide identity (ANI) of the genome bin with the closest relative. Functions encoded in the MAGs were classified according to the COG categories using eggnog-mapper (Huerta-Cepas et al. 2017).

Analysis of elemental composition

Varnish and non-varnish rock chips were ground with a Retsch mill and one gram of the ground rock was digested with HNO₃ and diluted using MilliQ water (final acid concentration 2% v/v) to avoid precipitation of dissolved elements. All the samples were analyzed by inductively coupled plasma atomic emission spectroscopy (ICP-OES) (Varian Vista AX, SPS-5, USA) to determine the concentration of dissolved elements. To test for difference in elemental composition between varnish and non-varnish we run Mann–Whitney test using R.

Rock surface analysis

The varnish and non-varnish rock samples were fixed and dehydrated prior to scanning electron microscopy (SEM) analysis using the modified method of Echlin, (2009). For fixation, the rock samples were incubated in 2.5% glutaraldehyde (diluted in 0.1 M phosphate buffer, pH 7.4) for 20–40 min then rinsed twice in 0.1 M phosphate buffer (pH 7.4). The sample surface was washed by rinsing in distilled water twice for 1 min. Rock samples were then dehydrated in a series of ethanol/water dilutions (70%, 90% for 10 min and 100% for 15 min) and incubated in acetone twice for 5 min. Finally, the rock was immersed in 100% hexamethyldisilazane for 3 min and were dried in a desiccator. All samples were coated with gold, and the surface was analyzed by SEM (JEOL JSM-7000F, Japan). The surface elemental composition of the un-coated rocks was analyzed using energy dispersive X-ray spectroscopy (EDS) (JEOL JSM-7000F, Japan).

Lipid extraction and analysis

For lipid and small molecule biomarker analysis in the varnish and non-varnish surfaces, rock surfaces were scraped using a low-speed diamond saw coupled with a continuous flow of de-ionised water. Resulting rock suspensions were transferred to glass vials and placed in a drying cabinet for 3 h at 40˚C.

Dried ground varnish and non-varnish surfaces were weighted (approx. 1 g) and extracted using 3x Methanol, 3x (1:1) Methanol: Dichloromethane (DCM) and 3x DCM. Dried extracts were subsequently methylated by heating with methanol and HCl (Ichihara and Fukubayashi 2009) then silylated by adding 20 µL pyridine and 20 µL BSTFA and heating at 60˚C for 20 min. Ethyl acetate (25 µL) was added as the injection solvent and samples were analysed using a Shimadzu QP 2010 Ultra gas chromatography mass spectrometer (GC/MS) fitted with a Zebron ZB-5HT column (30 m x 0.25 mm x 0.10 µm) using 1.5 mL per minute helium carrier gas flow. Samples were injected splitless onto the column at 60˚C, ramping to 180˚C at a rate of 15˚C/min, followed by ramping to 325˚C at a rate of 4˚C/min with a final hold time of 20 min, the total run time was 60 min. The MS was set to scan from 50 to 800 m/z with an event time of 0.7 s and a scan speed of 1111 u/s. Contamination was not introduced into the work up blank samples which was comparable to the instrument blank.

RESULTS AND DISCUSSION

The elemental composition among the varnish and non-varnished areas did not display significant differences (Fig. S1, Supporting Information). Although not significant (P < 0.05), trends were visible for elements such as Co, Cr and Ni, suggesting that the sampling site had a major effect on elemental composition, thus masking the differences between varnish and non-varnish. The relatively small sample size and the fact that varnished areas could have similar composition but different elemental oxidation state to non-varnish could further bias the statistical inference. The SEM-EDS analysis showed hyphal structures diffused all over the rock surface, other shapes (such as rods, cocci or spirilli) were found more rarely (Fig. 2). The size of such hyphae (∼3 μm) suggest that it could belong to fungal species. It was not possible, however, to determine clear differences in the morphology of microorganisms from varnish and non-varnish (Supplementary Files 1). The lipid biomarkers analysis did not either show large differences between varnished and non-varnished surfaces, beside the observation that isoprenoids, which are lipid markers for eukaryotic organisms such as algae and plants, were detected only on non-varnish surfaces (Table S1, Supporting Information). The saturated and unsaturated fatty acids found are chemotaxonomically non-specific and sterols and ergosterol (specific fungal biomarker lipids) were found at both varnish and non-varnish sites. In addition to lipid, small molecules were observed in the GC/MS run. The organic acid oxalate was only present in varnish samples.

The sequencing run resulted in 25604579 ± 3101080 reads per sample and the 10.7% ± 1.6% of reads mapped on ribosomal RNA, the 30.5% ± 2.4% mapped on genes with known function, the 58.5% ± 2.5% mapped on genes with unknown function (Table 1). The raw reads assembled into 287661 contigs, with length ranging 1000–166761 bp for a total length of a little more than 754.5 Mbp and an N50 of 53448 bp (Fig. S2, Supporting Information). There were 952315 genes detected, and according to the available sets of single copy genes in previous studies (Campbell et al. 2013; Waterhouse et al. 2018), it was estimated that there were 81 bacterial genomes and nine eukaryotic genomes (Fig. S3, Supporting Information).

There were 4666 Taxa detected spanning all kingdoms of life, most of them were ubiquitous, with a core microbiome of 2578 Taxa (55.2%, Fig. S4, Supporting Information). LEfSe analysis revealed that some bacterial taxa were significantly more abundant in varnish environment compared to non-varnish. Those results were consistent with our previous study based on 16S (i.e. Chloroflexi, represented mostly by the species Ktedonobacter racemifer, Esposito et al. 2015), However, this study revealed different structure of the abundant species that may have been overlooked by the 16S rRNA metabarcoding approach (Fig. 3 and Fig. S5, Supporting Information). For example, in the present study we found among the varnish biomarkers Acidiphilium cryptum, a heterotrophic alpha-proteobacterium previously isolated from acidic rocks and Rubrobacter xylanophilus, which was consistently found in previous studies on rock varnish microbiome (Kuhlman et al. 2006a; Zhang et al. 2012; Lang-Yona et al. 2018). Bacterial species more abundant in non-varnish included Terriglobus saanensis and Acidobacterium capsulatum (along with several other uncharacterized Acidobacterial taxa, Fig. 3).

Interestingly, no photosynthetic eukaryotes were found as varnish biomarkers, whereas non-varnish rocks hosted two photosynthetic species, namely, a Chlorella-like algae and Tetraselmis maculata (Fig. 3). This is consistent with the hypothesis formulated in Esposito et al. 2015 that oxygenic photosynthesis is disadvantaged in such environment.

Archaeal species with higher abundance in varnish samples were Sulfolobus solfataricus, Methanospirillum hungatei and Halomicrobium mukohatei. Conversely, in non-varnish rocks Methanosphaerula palustris, Haloarcula marismortui and Halobacterium salinarum were prevalent (Fig. 3 and Fig. S6, Supporting Information).

Mycosphaerella graminicola (class: Dothideomycetes) (Fig. S7, Supporting Information) was the only fungal biomarker taxa for varnish samples. Lang-Yona et al. (2018) also found species belonging to the same class (Dothideomycetes) in varnish samples.

In non-varnish samples, Lodderomyces elongisporus, Candida albicans and Rhizoplaca chrisoleuca were more abundant. Viral diversity appeared to reflect the host distribution, i.e. Bacillus phage SP-beta had a higher occurrence in varnish, where species belonging to the Bacillus Genus were more abundant and, conversely the Enterobacteria phage N4 was more abundant in non-varnish where Enterobacteria were more numerous (Fig. S8, Supporting Information). Considering the absence of physical dispersal barriers between the two environments, the results described above suggest that environmental selection, rather than dispersal, determine the diversity of microbial communities on the rock surfaces under study.

Reads mapping on functionally annotated genes provided strong evidence for a chemolithoautotrophic lifestyle of varnish microbiome. CO₂-fixation pathways associated genes were detected, revealing the existence of a microbiome with diverse autotrophic potential. CO₂-concentrating carboxysomes (bacterial organelles associated to photosyntheis), glyoxylate bypass and photorespiration were significantly more abundant in varnish surfaces (Fig. 4 and Fig. S9, Supporting Information). The abundance of taxa belonging to the Chloroflexi Phylum suggests that another C-fixation process may operate via the Wood–Ljungdahl pathway, which is encoded in uncultured Chloroflexi genomes (Hug et al. 2013). Fungi like Mycosphaerella graminicola contain genes encoding oxalate oxidase within the glyoxylate and dicarboxylate pathway. They require high amount of Mn and Fe (that is actively recruited from the external environment) for the production of organic acids (Gadd 1999).

Besides, varnish samples displayed a wide variety of genes associated with heterotrophic lifestyle, including genes for central carbohydrate metabolism, biotin biosynthesis, methionine biosynthesis, sucrose, purine and ketogluconate utilisation (Fig. 4).

Genes for denitrification and nitrogen acquisition using allantoin utilisation pathway were prevalently detected in varnish metagenomes (Fig. 4). With respect to micronutrient uptake, the varnish samples featured genes for siderophore biosynthesis and energy-coupling factor (ECF) transporters (Fig. 4). These ECF transporters are a subgroup of ABC (ATP-binding cassette) transporters which mediate the translocation of element-organic compound complexes across the cytoplasmic membrane (Köster 2001).

Non-varnish metagenomes biomarkers included genes for copper homeostasis, cobalt (Co), heme and siroheme biosynthesis, molybdopterin cytosine dinucleotide and molybdopterin oxidoreductase (Fig. 4).

Both varnish and non-varnish samples contained a number of stress response pathways. We found a higher occurrence of oxidative stress genes (like glutathione analogs mycothiol, Fig. 4) in varnish metagenomes. This could be an adaptation to oxidative stress due to the exposure to ionizing radiations (Fredrickson et al. 2008; Webb and DiRuggiero 2013). Another survival strategy is represented by the high abundance of DNA repair genes (bacterial RecBCD and MutL-MutS system, Fig. 4) in the varnish metagenomes, which are key parts of DNA protection against desiccation, radiation and oxidative damage (da Costa, Santos and Galinski 1998; Smith 2012; Holzinger and Karsten 2013).

On the other hand, genes for EPS, along with genes encoding for secretory systems, were more abundant in the non-varnish metagenomes (Fig. 4). A possible reason is that non-varnish surfaces retain more liquid water and EPS might play a significant role in water retention during dry seasons (Hu et al. 2012). In contrast, the varnish surfaces have lower water retention and therefore EPS biosynthesis would not be an effective survival strategy. In this regard, we explored also several genes for microbial motility in the non-varnish metagenomes but not in the varnish metagenomes (Fig. 4).

Other stress response pathways were also detected. For instance, in varnish samples, we found genes encoding for host defense and restriction modification systems, including clustered regularly interspaced short palindromic repeats (CRISPR) genes and type I restriction modification genes that confers resistance to foreign genetic elements (Clokie et al. 2011) (Fig. S9, Supporting Information). Differently, in non-varnish metagenomes, we found several genes involved in detoxification, metal resistance, host defense, heat shock (like, macromolecular synthesis operon) and osmoregulation (like, synthesis of osmoregulated periplasmic glucans, Fig. 4).

Binning with CONCOCT resulted into 209 bins, accounting for >754 Mbp, however, after manual checking only eight genomes had completion level above 90% and redundancy lower than 10%. Five of the 90% completed genomes were assembled mainly from non-varnish surfaces and three were assembled mainly from varnished surfaces (Table 2). None of the genome bins could be classified at the genus or species level, meaning every bin had an ANI below 85% to the known closest relative. Sequences of the genome bins are available as Supplementary Files 2. Most proteins encoded in the MAGs genomes were classified with unknown functions (Fig. 5), despite this, proteins connected to amino-acid and carbohydrate metabolism and transport, proteins involved in replication and repair machinery and outer structures were represented with an abundance >5% in all genomes. No significant differences were found in the abundance of proteins classified in any given COG category in genomes retrieved from varnish or non-varnish. This could be due to the low N (5 versus 3), but also to the fact that the different functions could be performed by the organisms at low abundance (which could not be assembled due to low coverage).

CONCLUSION

The present work elucidates the functional diversity of microbes inhabiting varnish rock surfaces and non-varnished rock using whole-metagenome shotgun approach. We found new insights into the ecology and survival strategies of microbial communities colonising varnish and non-varnish rock surfaces. The main conclusion is that varnish and non-varnish surfaces present differential middle-scale functional patchiness and adaptations to stresses. Functional biomarkers involved in the scavenging and metabolism of transition metals (namely, NiFe hydrogenase, ECF class transporters and siderophores), and the lack of enzymes specific for manganese oxidation, raise a puzzling question regarding the role of microbiome in varnish formation. The absence of specific for manganese oxidating enzyme was consistent with a previous study (Lang-Yona et al. 2018). The presence of oxalate (inferred by the abundance of Micosphaerella graminicola, and detected by the GC/MS analysis) and the siderophore encoding genes suggest that process of Fe and Mn accretion could be biological catalyzed by distinct microbial communities. The role of the microbiome in varnish formation may go through the local alteration of pH that fractionates manganese leading to the increase in Mn/Fe ratio (Thiagarajan and Lee 2004).

Such a wealth of manganese is useful in mechanisms of resistance to oxidative stress (e.g. it acts as a co-factor for enzymes such as superoxide dismutase), whereas iron can be scavenged using siderophores.

Organisms performing oxygenic photosynthesis are less abundant on varnish. However, varnish surfaces host autotrophic organisms with alternative CO₂ fixation pathways. It remains to be elucidated why oxygenic photosynthesis is disadvantaged on varnish.

ACKNOWLEDGEMENTS

We acknowledge the help of Kjell Jansson for his suggestions for SEM analysis and Boyke Bunk for bioinformatics. The authors would like to acknowledge support from Science for Life Laboratory, the National Genomics Infrastructure (NGI) and Uppmax for providing assistance in massive parallel sequencing and computational infrastructure. The present study was supported by the Erich-Ritter and the Herzog-Sellenberg Foundations within the Stifterverband für die Deutsche Wissenschaft within the EMERGE project (CUP I41J11000490007) to Brusetti L., and grants from Liljevalchs fund and Bergs foundation to Ahmed E.

Conflicts of interest. None declared.

REFERENCES

Author notes