The origin and role of biological rock crusts in rocky desert weathering

In drylands, microbes that colonize rock surfaces have been linked to erosion because water scarcity excludes traditional weathering mechanisms. We studied the origin and role of rock biofilms in geomorphic processes of hard lime and dolomitic rocks that feature comparable weathering morphologies, although these two rock types originate from arid and hyperarid environments, respectively. We hypothesized that weathering patterns are fashioned by salt erosion and mediated by the rock biofilms that originate from the adjacent soil and dust. We used a combination of microbial and geological techniques to characterize rock morphologies and the origin and diversity of their biofilms. Amplicon sequencing of the SSU rRNA gene suggested that bacterial diversity is low and dominated by Proteobacteria and Actinobacteria. These phyla only formed laminar biofilms on rock surfaces that were exposed to the atmosphere and burrowed up to 6 mm beneath the surface, protected by sedimentary deposits. Unexpectedly, the microbial composition of the biofilms differed between the two rock types and was also distinct from the communities identified in the adjacent soil and settled dust, showing a habitat-specific filtering effect. Moreover, the rock bacterial communities were shown to secrete extracellular polymeric substances (EPSs) that form an evaporation barrier, reducing water loss rates by 65 %–75 %. The reduced water transport rates through the rock also limit salt transport and its crystallization in surface pores, which is thought to be the main force for weathering. Concomitantly, the biofilm layer stabilizes the rock surface via coating and protects the weathered front. Our hypothesis contradicts common models, which typically consider biofilms to be agents that promote weathering. In contrast, we propose that the microbial colonization of mineral surfaces acts to mitigate geomorphic processes in hot, arid environments.

1 Introduction

In arid and hyperarid stony deserts, bedrock surfaces are typically barren and free of vegetation or continuous soil mantle. When these surfaces are exposed to atmospheric conditions, they undergo weathering processes that shape the landscape (Smith, 2009). Weathering is an in situ set of processes that includes physical, chemical, and mechanical forces that result in the breakdown and transport of shattered fragments of the parent rock. Weathering can appear in a range of sizes and morphologies (Smith et al., 2005), including gravel shattering (Amit et al., 1996), surface crazing (Smith, 1988), ventifacts (Smith, 1988), micro-rills (Smith, 1988; Sweeting and Lancaster, 1982), and cavernous patterns ,also known as tafoni, honeycomb, or pitting (Mustoe, 1983; Viles, 2005). Weathering is an essential, although often neglected, element in the overall denudation of hot deserts.

Cavernous weathering is one of the most frequently occurring weathering patterns that has been observed in various regions across the globe, including humid, arid, cold, hot, coastal, and inland sites (Bruthans et al., 2018). In the Negev Desert, Israel, cavernous weathering patterns are common in carbonate rocks in arid and hyperarid regions. Upon exposure to the atmosphere, these rocks develop a carbonate coating, termed calcrete or dolocrete (on limestone or dolomite, respectively) by displacive and replacive cementation of calcium or dolomite onto the rock surface (Wright and Wacey, 2004; Alonso-Zarza and Wright, 2010). Following the cementation processes, typical honeycomb features are formed on the exposed parent rock, typified by pits separated by thin walls that are coated by the calcrete or dolocrete. Recent studies suggest that microbial activity also promotes the processes of calcrete and dolocrete formation (Alonso-Zarza et al., 2016; Alonso-Zarza and Wright, 2010).

The accepted conceptual model for the formation of cavernous rock weathering in hot deserts involves the presence of permeable rocks that are subjected to soluble salts and repeated episodes of drying–rewetting cycles (Goudie et al., 2002; Smith, 1988; Smith et al., 2005). The proposed mechanisms assume that cavernous weathering results from physico-chemical processes including salt crystallization (Cooke, 1979; Scherer, 2004), incipient fractures (Amit et al., 1996), exfoliation (Shtober-Zisu et al., 2017), or stress erosion (Bruthans et al., 2014; McArdle and Anderson, 2001). Recently, Bruthans et al. (2018) conclusively demonstrated that moisture flux followed by salt crystallization at the boundary layer govern the case hardening model in temperate climates. In addition, biological mechanisms have been proposed to promote rock weathering through mechanisms such as flaking via colony growth (Viles, 2012), acidification by bacterial extractions (Garcia-Pichel, 2006; Warscheid and Braams, 2000), or alkalization during photosynthesis by cyanobacteria (Büdel et al., 2004). In contrast, it has been proposed that micro- and macro-organism colonization can mitigate weathering in temperate, coastal regions (McIlroy de la Rosa et al., 2014; Mustoe, 2010) via encrustation or protection from direct rain impact. However, it is not clear which of these mechanisms dominates or what the relative contribution of chemical vs. biological processes to weathering in arid environments is.

Microorganisms colonizing rocks form a hardy biofilm known as the biological rock crust (BRC), which is common in most arid and hyperarid regions worldwide (Gorbushina, 2007; Lebre et al., 2017; Pointing and Belnap, 2012). Epilithic communities colonizing rock surfaces are ubiquitous in arid environments, whereas hyperarid rocks, which experience increased radiation and desiccation, are dominated by endolithic communities that colonize internal rock pores (Makhalanyane et al., 2013; Pointing and Belnap, 2012; Viles, 1995). The BRC communities include cyanobacteria and other phototrophic and heterotrophic bacteria, but very low abundances of archaea, fungi, or algae (Lang-Yona et al., 2018). However, the BRC inoculum has not been resolved and has been proposed to originate from settled dust (Viles, 2008) or the surrounding soil (Makhalanyane et al., 2015).

The goal of this study was to illuminate the origin and role of BRCs in cavernous weathering of exposed limestone and dolomite rocks in arid and hyperarid regions. We predicted that the BRC communities on exposed rock surfaces would resemble either the ever-present dust or the surrounding soil, supporting a subset of adapted taxa from both sources. We further hypothesized that the cavernous weathering morphologies of exposed rocks result from salt mobilization by dew, causing crystallization pressure under atmospheric conditions. The developed rock biofilms clog the surface rock pores through secretion of extracellular polymeric substances (EPSs), lowering evaporation and slowing the salt crystallization, but also stabilizing the exfoliated rocks and preventing further weathering. Thus, the presence of a BRC mitigates the geomorphic processes. To test our hypotheses, we applied a holistic approach combining field observations, and geological, geotechnical, and molecular microbiology characterization to elucidate the morphology, origin, and role of BRCs in arid cavernous weathering.

2 Materials and methods

2.1 Study site

We focused on two sites in the Negev Desert, Israel: Sede-Boqer – an arid site, and Uvda Valley – a hyperarid site (Fig. S1, Table S1 in the Supplement). Both sites have rocky terrains underlaid predominantly by carbonate rock slopes consisting of limestone, dolomite, chalk, marl, clay, and chert from the Cretaceous period to the Eocene. Our analyses compared samples from the limey Turonian age Shivta Formation located in the arid region with samples from the dolomitic Turonian age Gerofit Formation located in the hyperarid environment. The Negev Desert, Israel, has maintained arid to hyperarid conditions since the Holocene and has an aridity index (P/PET) of 0.05–0.005 (Amit et al., 2010; Bruins, 2012), which is similar to other arid and hyperarid areas worldwide, e.g. the Namib and Atacama deserts (Azua-Bustos et al., 2012; Viles and Goudie, 2007). The long-term aridity of the Negev Desert makes it a reliable site for testing the relationships between BRCs and geological substrates.

2.2 Field sampling

A total of 24 rock samples were collected along rocky slopes facing northward during November and December 2014; these samples were comprised of 12 limestone samples from the limey Turonian age Shivta Formation at the arid site (30.88^∘ N, 34.78^∘ E; WGS 84 Grid; the samples were named SB 1–12) and 12 dolomite samples from the limey, dolomitic Turonian age Gerofit Formation at the hyperarid site (29.94^∘ N, 34.97^∘ E; WGS 84 Grid; the samples were named UV 1–12). Concomitantly, six soil samples (ca. 500 g each) were collected: three samples from the arid site (referred to as SBSoil 1–3) and three samples from the hyperarid (referred to as UVSoil 1–12) site. Each rock or soil sample was a composite of four sub-samples that were pooled and homogenized in the lab.

We also collected settled dust samples using glass bead traps (Goossens and Rajot, 2008). The traps were placed during a dust storm on December 2013 and collected 3 months later at both the arid (the samples were named SBDust 1–2) and hyperarid (the samples were named UVDust 1–2) sites. In total, two dust samples were collected per site. Each dust sample was a composite of two sub-samples that were pooled and homogenized in the lab.

2.3 Geological analyses

The geological methods used in this study are based on direct field observations and a detailed characterization of the subject lithologies (i.e. limestone and dolomite) which included the morphology (thin sections), mineral components (X-ray powder diffraction – XRD), porosity, and permeability (automatic gas permeameter porosimeter), and the elastic properties (Schmidt hammer). Petrographic thin sections, 30 µm thick, were prepared for each lithology to test the main components in both the BRC and the host rocks, and the sections were then examined under a light microscope (Zeiss, Oberkochen, Germany). The XRD analysis of mineral components (Sandler et al., 2015) was conducted on the BRC and host rocks using three replicates of each. Powdered samples were scanned using an X'Pert³ Powder diffractometer equipped with a PIXcel detector (Malvern Panalytical, Almelo, Netherlands); the scanning range was 3–70^∘ 2θ, with a step size of 0.013^∘, and a speed of 70.1 s per step. The total effective porosity (Φ) and the permeability (k) tests (Scherer, 1999) were performed on 12 rock core cylinder samples that had a 18.5 mm radius and 26.5 mm height using an automatic gas permeameter porosimeter (Core Laboratories, Houston, Texas, USA). Six cores were prepared from each lithology, and from the six cores three were cut parallel and three were cut perpendicular to the bedding plane. Next, 5 cm layer was mechanically removed from the top of the core using a diamond saw (Dremel, Racine, WI, USA) to ensure the complete removal of the BRC. Before testing porosity and permeability, samples were oven dried at 110 ^∘C for 24 h. Schmidt hammer (Lassen, Aarhus, Denmark) tests were applied in the field (Goudie, 2016; Viles et al., 2011). A total of 20 measurements were carried out for each lithology.

2.4 FTIR and stable isotope analysis

Fourier transform infrared spectroscopy (FTIR) analysis was conducted to test for the presence of extracellular polymeric substances (EPSs) on the rock surfaces, and the host rock was used for comparison. The spectra were recorded using a Vertex 70 FTIR spectrometer (Bruker, Billerica, MA, USA) with a 4 cm⁻¹ scan resolution. A total of 1–2 mg of pulverized rock was taken from each sample (n=2), and the spectra were measured twice; the spectra were collected over a wavenumber range from 4000 to 600 cm⁻¹, and a baseline correction was carried out. The spectral absorption bands indicative for EPSs were identified according to published information (Ferrando et al., 2018).

For δ¹³C and δ¹⁸O analyses, 1–2 mg of rock surface powder (i.e. calcite or dolomite) was obtained using a micro drill (Dremel, Racine, WI, USA) along with a cross section of the rock crust and its host rock. Four profile measurements of δ¹³C and δ¹⁸O were performed on samples UVSL 5 and UVSL 6 from the hyperarid site and samples NWSH 1 and NWSH 2 from the arid site. Measurements (in duplicate) of δ¹⁸O-H₂O and δ¹³C-DIC were performed using a gas source isotope ratio mass spectrometer (GS-IRMS; Thermo Fisher Scientific, Waltham, MA, USA) coupled to a Gas Bench II interface (Thermo) after CO₂ equilibration or CO₂ extraction by acidification for δ¹⁸O-H₂O and δ¹³C-DIC, respectively. The samples were calibrated against internal laboratory standards: Vienna Standard Mean Ocean Water (VSMOW) and NBS-19 carbonate standard. δ¹³C values were also referenced against VSMOW and valued for carbonate relative to the Vienna PeeDee Belemnite (VPDB) standard, as previously described (Uemura et al., 2016), with a SD of 0.1 ‰. All values are reported in per mille (‰).

2.5 Desiccation experiment

To test the effect of BRCs on water transport rates in the rock, clogging and desiccation experiments were performed on 16 rock core cylinders from both lithologies (limestone and dolomite). Rock cylinders (with a radius of 37 mm and height of 6.5 cm) were drilled using a rock core drill. Each set of eight rock cores from the two different lithologies included four rock cores that were kept intact, and four rock cores that had their BRCs mechanically removed using a diamond saw (Dremel, Racine, WI, USA) to a depth of 5 cm. Each cylinder was immersed in distilled water for 72 h. Next, the entire surface of the cylinder, except the upper base, was covered with epoxy (Devcon) and then with aluminium foil, leaving only the top of the cylinders (with or without BRC) uncovered to allow evaporation. The cylinders were then weighed (t=0), incubated in an oven (dried) at a temperature of 44 ^∘C for 48 h, and weighed every 2 h during the first 12 h and then every 6 h to determine the residual water content. Second-degree polynomial functions were fitted using the lm{stats} function to determine evaporation rates and then compared using an ANOVA (both analyses were carried out in R).

2.6 DNA extraction PCR amplification and sequencing

For DNA extraction from all rocks, the surface (ca. 100 cm²) was scraped using a rasp (66–67 HRC hardness; Dieter Schmid, Berlin, Germany) that was cleaned with 70 % technical-grade ethanol before each sampling. DNA was then extracted from 0.5 g of the homogenized sample using bead-beating in the presence of a CTAB buffer and phenol, according to a previously published extraction protocol (Angel et al., 2012). A 466-bp fragment of the 16S rRNA gene was amplified using the universal bacterial primers 341F (CCTAYGGGRBGCASCAG) and 806R (GGACTACNNGGGTATCTAAT) flanking the V3 and V4 regions (Klindworth et al., 2012). Library construction and sequencing were performed at the DNA Services Facility at the University of Illinois, Chicago, using an Illumina MiSeq sequencer (Illumina, San Diego, CA, USA) in the 2×250 cycle configuration (V2 regent kit).

2.7 Sequence processing and analysis of bacterial communities

Paired reads generated by the MiSeq platform were quality filtered and clustered into operational taxonomic units (OTUs) using the UPARSE pipeline (Edgar, 2013), with modifications. Contig assembly was carried out using the fastq_mergepairs command. Contigs were then dereplicated with the derep_fullength command, and singleton sequences were removed. Following this, OTU centroids were determined with the cluster_otus command (set at 3 % radius). Abundances of OTUs were determined by mapping the filtered contigs (before dereplication, including singletons) to the OTU centroids using the usearch_global command (set at 0.97 % identity). Following these steps, a total of ca. 1.4 G reads remained. OTU representatives were classified using the mothur software package's implementation of a naïve Bayesian sequence classifier (Schloss et al., 2009; Wang et al., 2007) against the SILVA 119 SSU NR99 database (Quast et al., 2013). All downstream analyses were performed in R V3.4.4 (R Core Team, 2016). Data handling and manipulation were undertaken using the phyloseq package (McMurdie and Holmes, 2013). For α-diversity analysis, all samples were subsampled (rarefied) to the minimum sample size using bootstrap subsampling at 1000 iterations, to account for library size differences, whereas for β-diversity analysis library size normalization was carried out using the geometric mean of pairwise ratios (GMPR; Chen et al., 2018). The abundance-based coverage estimator (ACE) richness estimate (O'Hara, 2005) and the Shannon diversity index (H) were calculated using the estimateR function in the vegan package (Oksanen et al., 2018) and tested using an ANOVA and a Tukey HSD test in the stats package. Variance partitioning and testing were undertaken using PERMANOVA (McArdle and Anderson, 2001) with the vegan-function adonis using Horn–Morisita distances. Differences in phyla composition between the sample types were tested using the non-parametric Scheirer–Ray–Hare test (Mangiafico, 2018; the scheirerRayHare function in the rcompanion package), followed by the post hoc Mann–Whitney test (the wilcox.test function in the stats package), and FDR corrected using the Benjamini–Hochberg method (Ferreira and Zwinderman, 2006; the p.adjust function in the stats package). Detection of differentially abundant OTUs was carried out using ALDEx2 (Fernandes et al., 2014). Plots were generated using the ggplot2 (Wickham, 2016) and ggtern (Hamilton, 2017) packages.

3 Results and discussion

3.1 Field and mineralogical observations

Weathering features were observed in about 30 % of the exposed rocks sampled from both arid and hyperarid sites. Neither the prevalence of weathering nor its morphology seemed to differ between sites despite the different climates and underlying geology. In all cases, the weathering type was classified as either tafoni or honeycomb weathering (Goudie et al., 1997; Groom et al., 2015; Fig. 1a); furthermore, weathering was coupled with the presence of sub-aerial biofilm, which was burrowed underneath the surface and protected by sedimentary deposits (Fig. 1b). The weathering and presence of the crusts were restricted to the atmospherically exposed parts of the rock. The presence of an identical weathering morphology and the prevalence of this morphology in different climates and lithologies challenges the current model, which assumes that surface permeability, moisture, and the presence of salts are primary factors that control weathering rates (Goudie et al., 2002; Smith, 1988; Smith et al., 2005). However, the lack of a correlation between the magnitude of tafoni weathering and climate has already been reported (Brandmeier et al., 2011).

Figure 1 (a) Comparable weathering features in the exposed limestone and dolomite rocks at both sites as noted in field outcrops (hammer for scale, 30 cm long). (b) Visual presence of a rock crust with similar thickness (3–6 mm) in both rock types. The crust's mineralogical composition matched that of the host rock. (c) Thin-section analysis of the rocks showing lamination structure in the BRCs. Dashed lines indicate the interface between BRC host rocks. BRC's mineralogy includes micritic to microsparitic dolomitic or calcitic crystals.

To study the possible differences between these sites, we performed geological characterization of 10 limestone and dolomite rocks collected from the arid and hyperarid sites, respectively, testing for mineral content, porosity, permeability, and elasticity. As expected, our results showed different lithological parameters between the limestone and dolomite rocks (Table S2), although they displayed similar weathering features. Moreover, petrographic thin-section analysis showed that crusts had developed to a similar thickness of 1–6 mm on both rock types, irrespective of climatic conditions including mean annual precipitation (Fig. 1c). However, microclimatic conditions, such as dew or surface temperature may impact local morphologies. Furthermore, the thin sections showed that the crusts are composed of masses of micritic to microsparitic minerals that form a laminated structure (Fig. 1c). Such laminated structures indicate that the crusts are stage four terrestrial calcretes and dolocretes, suggesting a mature crust phase (Alonso-Zarza and Wright, 2010). The calcretes and dolocretes identified on the rocks' surface reject the previously suggested impact of mineralized networks or case hardening (McBride and Picard, 2004). In fact, the detection of mature calcretes could serve as an indication of atmospheric exposure but has also been suggested to result from biogenic activity (Alonso-Zarza and Wright, 2010; Goudie, 1996).

3.2 Composition and chemical characteristics of the rock crusts

To test our hypothesis that the crusts are biogenic and involved in rock weathering processes, we characterized their origin and nature. An XRD analysis of the crust layers and bedrock at the locations showed that the crusts are composed of a similar mineralogy to their respective host rocks, indicating that local weathering, rather than dust deposition, is the source of crust generation (Table S2).

The biogenic nature of the crusts was confirmed using a cross-section analysis of the stable carbon and oxygen isotope ratios in the crust and host rock (Fig. 2a). For both limestone and dolomite, values of δ¹³C increased between the crust and the host rock layers and ranged from −4.1 ‰ in the calcrete to −0.9 ‰ in the limestone bed, and from 0.2 ‰ in the dolocrete to 2.0 ‰ in the dolomite bed. Such values are typical indicators of carbon isotope exchange of primary marine CaCO₃ (abundant in the bedrock) with CO₂ released by microbial respiration (i.e. of carbon originating from photosynthesis) with the subsequent precipitation of pedogenic calcrete (Brlek and Glumac, 2014; Mora et al., 1991). Analysing δ¹³C, along with δ¹⁸O, compositions of pedogenic carbonates are a useful way of reconstructing palaeo-vegetation (e.g. the C3∕C4 plant ratio (Ehleringer et al., 1997; Mora et al., 1991). Our δ¹³C results concur with δ¹³C values collated from speleothems (secondary mineral deposits formed in caves) collected in the central and southern Negev Desert (Vaks et al., 2010), which were also dated to end of the Pliocene (the past 2.5 million years). The low ratio detected here (Fig. 2a) and by Vaks et al. (2010) suggest that the Negev region has only been able to support limited vegetation for at least 2.5 Gyr; if this is the case, the role of the crust in shaping the morphology of the rock surfaces has been considerable. These results support the hypothesis that calcretes and dolocretes are of biogenic origin; therefore, the crust can be referred to as BRC. Moreover, they indicate a similar developmental trajectory for both the calcretes and dolocretes, which is independent of aridity or lithological parameters.

Figure 2 (a) Carbon and oxygen isotope ratio depth profiles of the limestone (top) and dolomite (bottom) BRCs compared with to their host rocks. The blue horizontal line indicates the border between BRC and host rock. (b) Fourier transform infrared (FTIR) analysis of limestone (top) and dolomite (bottom) BRCs indicating the presence of extracellular polymeric substance (EPS) molecules through the distinctive peak ranging between 1020 and 1040 cm⁻¹, which was absent from the host rocks.

In contrast, the trend in values of δ¹⁸O differed between rock types. In the limestone rocks, the ratio ranged from −3.0 ‰ to −6.8 ‰ between the BRC and the host rock, whereas in the dolomite it was higher, ranging from −5.4 ‰ to −0.6 ‰. The decrease in δ¹⁸O in the host limestone rock could be explained by meteoric water substitution (Sandler, 2006). In contrast, the more negative δ¹⁸O values in the dolocrete compared with the host dolomite are attributed to isotopic differentiation of meteoric water due to condensation (Rayleigh distillation) and could result from the large distance from the Mediterranean Sea (which is the primary source of rainfall in the area) compared with closer limestone rocks. In speleothems, similar patterns of δ¹⁸O values have been reported in the central and southern Negev Desert (Vaks et al., 2010). The results suggest that the calcrete and dolocrete studied here have been experiencing arid to hyperarid climates since the Pleistocene, alluding to the possible source of rain. A similar study conducted in the Thar Desert in India also inferred sedimentary rock stable isotope patterns to paleoclimate (Andrews et al., 1998).

To study the potential role of BRCs in the weathering process, their composition was characterized using FTIR, as was previously reported (Sheng et al., 2010). We focused on the functional groups and elemental compositions in the EPSs or microbial aggregates and found a distinct peak in the BRC layers ranging between 1020 and 1040 cm⁻¹ in both limestone and dolomite rocks that was absent in the host rocks (Fig. 2b). This peak is indicative of the presence of EPSs of bacterial origin (Shirshova et al., 2006), pointing to the significant components of asymmetric and symmetric stretching of ${\mathrm{PO}}_{\mathrm{2}}^{-}$ and P(OH)₂ in phosphate as well as vibrations of C–OH and C–C bonds, which are found in polysaccharides and alcohols (Jiang et al., 2004). These results provide strong support for the biogenic nature of the crust, as EPS is a common feature of many if not most biofilms (Drews et al., 2006). The detected EPS could serve several functions in the BRC such as a dust-particle trap to collect the dust and its nutrients, a binding agent to individual members of the biofilm (Davey and O'toole, 2000), or a protective agent that decreases evaporation, retains moisture, and shields from radiation (Or et al., 2007; Roberson and Firestone, 1992).

Based on these findings, we hypothesized that a BRC could in fact act as a mitigator during the weathering process by clogging the pores on the surface of the rock and thereby minimizing capillary rise. Consequently, crystallization of dissolved salts, considered to be the primary mechanism for rock weathering, is mitigated. To test this hypothesis, we performed a drying experiment to estimate water loss from the rock surfaces covered with BRCs. The results suggest that (in both limestone and dolomite rocks) water moves through the rock and is lost to evaporation two or three times faster in the absence of BRC than when it is present (Fig. 3). Considering that salt transport due to hydraulic movement is a dominant weathering mechanism (Huinink et al., 2004), reduced evaporation due to BRC coverage will also inevitably lead to decreased weathering rates. Moreover, our results differ from those performed on temperate sandy stones that showed no significant effect of BRC on decreasing water transport rates (Slavík et al., 2017).

Figure 3 Desiccation of rock cores in the presence and absence of BRCs as a function of time, following full hydration. The curves indicate a second-degree polynomial line fitting (all fitted curves were significantly different from one another in ANOVA tests with P values < 0.01 and R²values>0.95).

3.3 The microbial composition and origin of the BRCs

To elucidate the identity of the bacterial communities on the limestone and dolomite BRCs, we applied multiplexed barcoded amplicon sequencing of the small subunit RNA gene (SSU rRNA). In addition, we compared the BRC communities to those in the samples of the surrounding soil and settled dust in order to deduce the origin of the rock biofilm. As expected, we found simple BRC communities (i.e. low-richness and low-diversity). The communities of the BRC showed an average of 182 observed phylotypes, 354 predicted phylotypes, and a Shannon H value of 3.8 (Fig. 4a; Table S3) for arid limestone, and 129 observed phylotypes, 315 predicted phylotypes, and a Shannon H value of 3.3 for hyperarid dolomite, with no significant difference between the rock types.

The surrounding soil was significantly richer and more diverse (P<0.05) at the arid site (416 and 746 observed and predicted OTUs, respectively, and the Shannon H value was 5.6 on average), and equally rich but slightly more diverse at the hyperarid site (221 and 466 observed and predicted OTUs, respectively, and the Shannon H value was 3.8, on average). The diversity of the dust samples was as poor as that of the BRCs (169 and 107 observed and predicted OTUs, respectively, and the Shannon H values were 3.0 and 1.5, on average) and did not differ between sites (Fig. 4a; Table S3). The number of observed OTUs in the soil and their diversity scores were somewhat lower in this study compared with other reports from similar environments (Barberán et al., 2014; Lang-Yona et al., 2018; Šťovíček et al., 2017), which may be due to sequencing technologies and depth. The lower richness and diversity in hyperarid vs. arid samples and the BRC and dust vs. soil samples is expected and is comparable with trends reported in other studies (Angel and Conrad, 2013; Barberán et al., 2014; Lang-Yona et al., 2018).

Figure 4 Microbial community features of the BRCs, the surrounding soils, and settled dust at the two study sites. (a) Comparison of the richness in the form of the observed number of OTUs (S obs.) and the predicted number of OTUs (the ACE index), in addition to a comparison of α-diversity (Shannon diversity index – H) between the different sample types. Identical lower-case letters indicate no statistical difference between groups in a Tukey HSD test. (b) Clustering of sample types using a PCoA ordination based on the Horn–Morisita distance matrix. Identical lower-case letters indicate no statistical difference between groups in a pairwise PERMANOVA test. Ellipses denote 95 % confidence intervals around the arid and hyperarid samples assuming multivariate normal distribution. (c) Composition of bacterial phyla in the different sample types (see Table S5 for results of statistical tests in the relative abundance of different phyla between sample types). (d) Relative contribution of each bacterial OTU to the community composition of each sample type: top – arid site, bottom – hyperarid site (see Fig. 5 for statistical detection of preferentially abundant OTUs between each sample-type pair).

β-diversity analysis, using variance partitioning, showed statistically significant differences between samples at the OTU level, based on climate, sample type (i.e. rock, soil, or dust), and, to a small extent, also on their interaction. Because previous field studies in dryland regions have showed that soil microbial communities typically exhibit rapid changes after hydration events but little to no seasonal changes (e.g. Angel and Conrad, 2013; Barnard et al., 2013; Šťovíček et al., 2017); therefore, we disregarded the seasonal aspect in our analysis. These variables were found to significantly contribute to the differences in bacterial communities, accounting for 22 %, 40 %, and 3.8 % of the total variance, respectively (Fig. 4b, Table S4). Pairwise comparisons further showed that the two BRCs significantly differed from one another (P<0.01) and also from their surrounding soil and dust samples (P<0.05 in all cases; Table S4). The bacterial community in the samples was typical for dryland areas, mostly dominated by members of the following phyla: Proteobacteria, Actinobacteria, Deinococcus-Thermus, Chloroflexi, Bacteroidetes, Cyanobacteria, Acidobacteria, Firmicutes, and Gemmatimonadetes (Fig. 4c, Table S5). Similar communities have repeatedly been reported for arid and hyperarid soils and rocks (Angel and Conrad, 2013; Barberán et al., 2014; Lang-Yona et al., 2018). While cyanobacteria are typically the main primary producers in the soil and rock communities (Weber et al., 2016), recent studies have shown that other autotrophs may also significantly contribute to the energy balance of these biofilms (Ji et al., 2017).

The BRCs of the two rock types differed in the relative abundance and composition of the major phyla. Most notably, Proteobacteria were significantly more dominant in the hyperarid compared with the arid samples (P=0.02), comprising 21 % and 44 % of the community on average in the limestone and dolomite BRC, respectively. In contrast, the Actinobacteria showed an opposite trend (P=0.03), comprising 42 % and 21 % of the community on average in the limestone and dolomite BRC, respectively. The two BRCs also differed in their composition of Firmicutes, Gemmatimonadetes, and Chloroflexi (P<0.03; Fig. 4c, Table S5).

The soil samples generally showed similar trends on the gross taxonomic level as their respective BRC samples. While none of the phyla differed significantly between the hyperarid BRC and the soil, the Deinococcus-Thermus, Acidobacteria, Firmicutes, and Gemmatimonadetes phyla differed significantly between limestone BRC and the surrounding arid soil (P<0.04; Table S5). Lastly, the arid and hyperarid dust samples were dominated by members of the Proteobacteria, with other phyla comprising only a minor fraction of the community (with the notable exception of Bacteroidetes, which dominated one of the dust samples). However, these differences were not significant, probably due to the small sample size (Table S5).

Despite the general similarities in the community composition between samples on the phylum level, many of the OTUs found in each sample were unique to the BRC, soil, or dust as evident from the ternary diagrams (Fig. 4d). Direct analysis of the differences in the OTUs detected 130 (10 %) differentially abundant OTUs in the dolomite BRC and 74 (6 %) differentially abundant OTUs in the limestone BRC (Fig. 5). Similarly, several differentially abundant OTUs were also detected when comparing the BRCs with their respective soil and dust samples. However, these differentially abundant OTUs were fewer, probably due to the small dust sample size (Fig. 5).

Figure 5 Detection of differentially abundant OTUs between each sample-type pair using ALDEx2. Each circle denotes a single OTU, and its size is its average relative abundance across all samples. The x axis shows the classification of each OTU, whereas the y axis shows the effect size in terms of the log2 fold difference in relative abundance between each sample-type pair. Red circles are OTUs that show significant differential abundance at the P<0.05 level. Numbers next to the arrows indicate the number of significant differentially abundant OTUs that are either more abundant (up arrow) or less abundant (down arrow) in the sample-type appearing in the column compared with the sample-type appearing in the row. Numbers in brackets indicate the total number of OTUs in each pair-wise comparison.

Other BRC bacterial communities have previously been described (e.g. rock varnish, quartz pavement) (Kuhlman et al., 2006; Lang-Yona et al., 2018; Wong et al., 2010a, b) but their origins and roles in geomorphological processes were not considered. Our results suggest that despite the similarity in the morphology and magnitude of rock weathering features in the arid limestone and hyperarid dolomites, the two BRCs we studied harboured distinct microbial communities, differing in over 16 % of the OTUs in addition to their composition at the phylum level. Moreover, despite the spatial proximity and continuous interaction between the limestone and dolomite surfaces and their respective surrounding soils and dust particles, the bacterial communities of the BRCs were distinct. The abilities of bacteria to disperse, settle, and persist in a given location could be an important factor resulting in the biogeographic patterns observed here. The difference between arid and hyperarid soil communities could result from the local contribution of aeolian material, which might affect the loess soil diversity (Crouvi et al., 2008). Alternatively, the hyperarid site experiences slow pedological processes, whereas the arid soil formation was enhanced (Amit et al., 2011) resulting in disparate bacterial communities. The three matrices (BRC, soil, and settled dust) studied here sparsely shared their bacterial communities. In particular the BRC communities had little in common with either the soil or the dust communities (Fig. 4). This demonstrates the ecological filtering effect of the rock surfaces, which imposes unique abiotic challenges on the microbes living on it (Horner-Devine and Bohannan, 2006). This also suggests that the BRCs cannot be regarded as passive deposits of microbial cells originating from the surrounding soil or dust, but should rather be seen as a specific subset of adapted microbes that can persist and form a biofilm under these unique conditions.

3.4 The role of BRC in arid rock weathering – synthesis

Honeycomb weathering patterns are prevalent worldwide and are found in both humid and dry ecosystems. According to contemporary models, this form of weathering is the result of the transport of dissolved salts through the rock and their eventual crystallization in surface pores, which leads to fractures and eventual flaking of rock material (Rodriguez-Navarro et al., 1999). In this study, we found that the weathering patterns and magnitude are similar on rocks from both arid and hyperarid sites, despite the differences in precipitation and lithologies. In arid and hyperarid regions, BRCs have been shown to form once the rock is exposed to the atmosphere (Pointing and Belnap, 2012). A developed crust of biological origin was microscopically and isotopically apparent on all weathered rocks and was shown to be supported by EPSs (Fig. 2). Similar to weathering magnitude, the BRCs showed no observable differences in form or depth despite the different aridity and lithology. Both BRCs comprised bacterial taxa that are typical for xeric environments (Pointing and Belnap, 2012) and included many heterotrophs but also dominant phototrophs or otherwise autotrophic members (Fig. 4). The two BRCs did differ in their bacterial communities at the OTU and higher taxonomical levels, demonstrating a discrepancy between community structure and ecological function (forming BRCs). The BRC communities also differed from their surrounding soil and dust, indicative of the specialization of the colonizing taxa to rock environments. In the absence of mineralized networks or case hardening (i.e. the addition of cementing agent to rock matrix material) we conclude that calcrete and dolocrete were formed through the colonization of microorganisms and the secretion of EPSs, which served as a thin biofilm (Brantley et al., 2011; Weber et al., 2016).

Our results further suggest that this biogenic layer mitigates evaporation and reduces water transport, thereby alleviating salt crystallization pressure in the rock pours (Scherer, 2004). Crystallization of calcium sulphate and sodium chloride solutions, which are abundant in these soils, was shown to build pressure within pores and stress rocks (Scherer, 2004; Sperling and Cooke, 1985). This process is enhanced under low relative humidity and rapid evaporation and compromises the durability of the rocks (Rodriguez-Navarro et al., 2003; Rodriguez-Navarro and Doehne, 1999). Our results suggest that the presence of a BRC decreases evaporation rates (Fig. 3) and thus attenuates the crystallization pressure and reduces damage to the rocks. Moreover, the BRC may also stabilize the rock following exfoliation, preserving the weathered structure.

Arid weathering features, which lead to debris formation result from a dynamic balance between the erosive salt forces and the mitigating effects of the BRC. The role of microbial biofilms in the protection of surfaces from mineral weathering was extensively studied for biomineralization and sedimentation processes (Adams et al., 1992; Dupraz et al., 2009). However, the role of BRC in weathering processes under atmospheric conditions in the desert has not been considered before. We propose that microbial colonization of mineral surfaces protects the rocks from weathering by mitigating salt crystallization and stabilizing the weathered front. Rock weathering processes are typically believed to be controlled at different scales ranging from the climatic scale, down to local conditions at the site, and eventually the micro-scale (Smith, 2009; Sperling and Cooke, 1985; Viles, 2001). The results presented here suggest that in arid environments, micro-scale conditions determine the magnitude of weathering that shape the landscape.