The role of ecosystem engineers in shaping the diversity and function of arid soil bacterial communities

Ecosystem engineers (EEs) are present in every environment and are known to strongly influence ecological processes and thus shape the distribution of species and resources. In this study, we assessed the direct and indirect effect of two EEs (perennial shrubs and ant nests), individually and combined, on the composition and function of arid soil bacterial communities. To that end, topsoil samples were collected in the Negev desert highlands during the dry season from four patch types: (1) barren soil; (2) under shrubs; (3) near ant nests; or (4) near ant nests situated under shrubs. The bacterial community composition and potential functionality were evaluated in the soil samples (14 replicates per patch type) using 16S rRNA gene amplicon sequencing together with physico-chemical measures of the soil. We have found that the EEs affected the community composition differently. Barren patches supported a soil microbiome, dominated by Rubrobacter and Proteobacteria, while in EE patches Deinococcus-Thermus dominated. The presence of the EEs similarly enhanced the abundance of phototrophic, nitrogen cycle, and stress-related genes. In addition, the soil characteristics were altered only when both EEs were combined. Our results suggest that arid landscapes foster unique communities selected by patches created by each EE(s), solo or in combination. Although the communities' composition differs, they support similar potential functions that may have a role in surviving the harsh arid conditions. The combined effect of the EEs on soil microbial communities is a good example of the hard-to-predict non-additive features of arid ecosystems that merit further research.

1 Introduction

Hot desert environments are characterized by long droughts interspersed by intermittent and unpredictable rain events. Water and nutrients in hot desert environments are scarce and unevenly distributed across the land, resulting in patches of contrasting productivities. High-productivity patches, also called resource islands, are defined by large concentrations of organic matter and nutrients (Bachar et al., 2012; Ben-David et al., 2011; Schlesinger et al., 1996; West, 1981). These resource islands can be formed through the redistribution of nutrients and water by ecosystem engineers (EEs), such as perennial plants or invertebrates (Wilby et al., 2001; Wright et al., 2006). EEs are also known for impacting many components of a given environment, such as soil features, annuals distribution, or community composition of microorganisms (De Graaff et al., 2015; Oren et al., 2007).

An EE is an organism that, directly or indirectly, modifies the availability of resources to other organisms by transforming the physical state of abiotic and/or biotic components of the ecosystem sensu Jones et al. (1994). The impacts of EEs range from physical, through the creation of biogenic structures (e.g. tunnels) (Lavelle, 2002), to chemical, through the production of compounds that have physiological effects (e.g. root exudates) (Lavelle et al., 1992), to biological, through organism behaviour (e.g. seed dispersal) (Lavelle et al., 2006). In drylands, resources, such as nutrients or water, are often concentrated around EEs, boosting the development of diverse populations of annual plants and invertebrates (Wright and Upadhyaya, 1996) as well as microbial communities (Bachar et al., 2012; Ginzburg et al., 2008; Saul-Tcherkas and Steinberger, 2011). This taxonomical response to changes in the physico-chemical conditions is linked to the potential function of the community (Narayan et al., 2020). This implies that the variation in taxonomy by the presence of an EE could potentially be associated with changes in functionality.

In desert ecosystems, ants are a notable example of an EE (Ginzburg et al., 2008). They redistribute resources by tilling the soil, bringing soil from the deep layers to the upper layers (bioturbation), and by gathering, storing, and ejecting food items, such as plant material or dead invertebrates, in and around the nest (Filser et al., 2016; Folgarait, 1998; MacMahon et al., 2000). EEs in arid environments also include perennial shrubs (Callaway, 1995; Schlesinger and Pilmanis, 1998; Segoli et al., 2012; Shachak et al., 2008; Walker et al., 2001). Their root systems create a soil mound that traps litter and seeds, allowing for higher water infiltration. The root exudates increase the content of organic matter and the shrub canopies decrease evaporation, prolonging water availability following a rain event (Bachar et al., 2012). In addition, the presence of shrubs alters the course of water run-off (Oren et al., 2007), which impacts the locations of available water for soil microbial communities. In addition, root systems have their own microbiome, which interacts with the soil microbial community (Steven et al., 2014).

The roles of both ants and perennial shrubs as EEs were reported in various ecosystems (Facelli and Temby, 2002; Farji-Brener and Werenkraut, 2017; Frouz et al., 2003; Gosselin et al., 2016; Pariente, 2002; Schlesinger et al., 1996). However, we know little about their joint effect on arid ecosystems. We hypothesized that each EE would shape a unique soil bacterial community via changes in the soil physico-chemical properties. We further predicted that since shrub canopies and ant nests may affect soil properties differently, their combined effect on the microbial community is non-additive and thus cannot be predicted by the contributing components. To test our hypotheses, we explored arid soil bacterial microbiomes and soil chemical features during the dry season of 2015. We sampled four different patches: under Hammada scoparia shrubs; near the nest openings of the harvester ant Messor ebeninus; in combined patches of the ants' nests under shrubs; and in barren soil.

2 Materials and methods

2.1 Sampling

The study was conducted in a long-term ecological research (LTER) site in the central Negev, Israel (Zin Plateau, 34^∘80^′ E, 30^∘86^′ N). It is characterized by 90 mm annual rainfall and average monthly temperatures fluctuating from 13 ^∘C (January) to 35 ^∘C (August). Vegetation is scarce and dominated by the perennial shrubs Hammada scoparia and Atriplex halimus (Gilad et al., 2004).

Sampling was conducted as previously described (Baubin et al., 2019) with slight modifications, such as the inclusion of Shrub & Nest samples. To summarize, we sampled four distinct patch types: (1) barren soil (Barren); (2) under the canopy of H. scoparia (Shrub); (3) 20–30 cm from the main opening of the nest of M. ebeninus (Nest); and (4) 20–30 cm from an ant nest's opening that was situated under a shrub canopy (Shrub & Nest). Samples were collected in October 2015, after an 8-month drought.

We sampled 14 random experimental blocks from each of the four patches (4 patch types ×14 blocks = 56 samples). The samples were collected using a scoop that was sterilized between each sampling using 70 % technical ethanol. Soil was collected from the top 5 cm after removal of the crust and debris. Three subsamples of ∼100 g were collected from each block and pooled together. In the lab, samples from two adjacent blocks were composite and homogenized using a 2 mm sieve. The samples were then separated for consecutive analyses: 15 g of each soil sample was stored in −80 ^∘C for bacterial analysis, 25 g was used to determine the water content in the soil, and the rest was used for the measurements of physico-chemical properties.

2.2 DNA extraction, amplification, and sequencing

Total nucleic acids were extracted from 0.5 g of soil as previously described (Angel, 2012), purified with the ExgeneTM Soil SV kit (GeneAll, Seoul, South Korea) according to the manufacturer's instructions. The 16S rRNA encoding genes V3–V4 region was amplified using 341F and 806R primer (Klindworth et al., 2013). The PCR consisted of 2.5 µL 10× standard buffer, 10 µM primers, 0.8 mM dNTPs, 0.4 µL DreamTaq DNA polymerase, 4 µL template, 1 mM bovine serum albumin (Takara, Kusatsu, Japan), and 12.6 µL Milli-Q water. Triplicate PCR reactions (95 ^∘C for 30 s; 28 cycles of 95 ^∘C for 15 s, 50 ^∘C for 30 s, 68 ^∘C for 30 s; 68 ^∘C for 5 min) were pooled and amplicon concentration and purity were measured by electrophoresis (Nanodrop ND-1000, Thermo Fisher Scientific, Waltham, MA, USA). The amplicon libraries were constructed and sequenced on the Illumina MiSeq platform (2×250, pair-end) at the Research Resources Centre at the University of Illinois.

2.3 Soil physico-chemical analysis

The physico-chemical parameters of the soil samples were assessed following the standard methods (SSSA, 1996). Water content was measured by gravimetry. Other parameters were measured as follows by the Gilat Hasade Services Laboratory (Moshav Gilat, Israel). The pH was measured in saturated soil extract (SSE). Phosphorus (P) was extracted by the Olsen method using a 0.5 M sodium bicarbonate solution (NaHCO₃) and the absorbance of the final solution was measured at 880 nm using a spectrophotometer. Nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) and ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) were extracted with a 2 N potassium chloride (KCl) solution and measured at 520 and 660 nm, respectively. Organic matter (OM) content was determined by the Walkley–Black method using a dichromate oxidation (${\mathrm{Cr}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{7}}^{-\mathrm{2}}$) and the amount of oxidizable OM was measured at 600 nm.

2.4 Bioinformatic analysis

The reads were quality-checked with MultiQC and trimmed using TrimGalore. Briefly, all reads with a quality of less than 20 and shorter than 150 bp were removed, and the rest were analysed further. The reads were then gathered into amplicon sequence variants (ASVs) (99 % identity cutoff) and merged using Dada2 (Callahan et al., 2016) in QIIME2 (Bolyen et al., 2018) following the NeatSeq-Flow pipeline (Sklarz et al., 2018). ASV counts were normalized to equal sampling depth (9100 reads). The taxonomic assignment was done using Silva (version 132) (Quast et al., 2013) through QIIME2 and all non-bacterial data have been characterized as unclassified.

2.5 Statistical analysis

The statistical analysis was done using R (R Core Team, 2016). To visualize the differences between patch types, a non-metric multidimensional scaling (NMDS) plot was created using the Bray–Curtis dissimilarity and the significance of these differences was analysed using a non-parametric analysis of similarity (ANOSIM) (“vegan” package, Oksanen et al., 2014). The envfit function (“vegan” package, Oksanen et al., 2014) was applied on the NMDS data to evaluate the effect of soil parameters on the bacterial community. The NMDS was plotted using the “ggplot2” package (Wickham, 2016) and the arrows representing the effect of each soil parameter as well as the centroids for each patch type, calculated using envfit, were added to the plot. The bacterial data were analysed using the “phyloseq” package (McMurdie et al., 2017). The relative abundance, whenever higher than 0.05 %, of each phylum, class, and order was calculated and then plotted using a stacked bar plot (“ggplot2” package, Wickham, 2016). The significance of difference between patch types was assessed using a non-parametric test: a Kruskal–Wallis test and a post hoc Dunn test (Dinno, 2017; Dunn, 1964; Kruskal and Wallis, 1952). All sequences retrieved in this study were uploaded to BioProject (https://www.ncbi.nlm.nih.gov/bioproject, last access: 30 September 2018) under the submission number PRJNA484096.

2.6 Functional prediction

The prediction of function of the 16S amplicons was done with Piphillin using the KEGG database (October 2018). Piphillin generates a genome abundance table that is normalized to the 16S rRNA copy number for each genome (Iwai et al., 2016; Narayan et al., 2020). To analyse the arid soil microbial functionality, we selected metabolisms and respective genes related to arid soil using groups and genes from the KEGG database (Kaneshisa and Goto, 2000). We selected steps in metabolic pathways for different methods of harvesting energy (organotrophy, lithotrophy, and phototrophy) (Cordero et al., 2019; Greening et al., 2016; León-Sobrino et al., 2019; Tveit et al., 2019), for parts of the nitrogen cycle (Galloway et al., 2004), and for the survival of the individual during a drought (DNA conservation and repair, sporulation, and reactive oxygen species (ROS)-damage prevention) (Borisov et al., 2013; Hansen et al., 2007; Henrikus et al., 2018; Preiss, 1984; Preiss and Sivak, 1999; Rajeev et al., 2013; Repar et al., 2012; Slade and Radman, 2011). Then, we looked for each step in the KEGG database and picked out genes of interest to build our own database. The assignment of function to the KEGG numbers was done in R. The significance of the differences between patch types in predicted functionalities was evaluated using a non-parametric test, a Kruskal–Wallis test and a post hoc Dunn test (Dinno, 2017; Dunn, 1964; Kruskal and Wallis, 1952), and boxplots were created in R.

3 Results

3.1 Soil physico-chemical characteristics

Table 1 depicts the differences in the soil characteristics (full list of values in Table A1) between the patches (barren, nest, shrub, and Shrub & Nest). Shrub & Nest patches have higher concentrations of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and P (30 and 54 mg kg⁻¹, respectively) than the average of the other patches combined (4.7 and 22 mg kg⁻¹, respectively). When verifying with a Kruskal–Wallis test and a Dunn test on the values of these soil variables (Table A2), we see that the differences between patch types are significant (Shrub & Nest vs. all other patches, p<0.05). Patches with two EEs also have a significantly higher concentration of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (9.72 mg kg⁻¹) and OM (8.21 %) compared to all other patches (${\mathrm{NH}}_{\mathrm{4}}^{+}$ mean: 5.62 mg kg⁻¹, p value < 0.05; OM mean: 5.51 %, p≤0.05). However, the water content and pH did not show significant differences between patches (Table A2).

Table 1 Soil parameters presented as mean ± standard deviation (${\mathrm{NO}}_{\mathrm{3}}^{-}$: nitrate, P: phosphorus, ${\mathrm{NH}}_{\mathrm{4}}^{+}$: ammonium, OM: organic matter content, water: water content).

3.2 Beta diversity

The summary of the sequence analysis can be found in Table A4. Dada2 analysis yielded 2318 ASVs and the ANOSIM results (Fig. 1, Table A3) suggest that there are significant differences in the microbial community between patch types (R=0.28; p=0.001). The envfit function shows that most soil parameters correlated with the barren patches but not with the other three patch types.

Figure 1 Non-metric multidimensional scaling (NMDS) of the soil 16S microbial communities in the dry season under different patch types. The centroid for each patch type is represented by a dashed circle. The arrow vectors represent the effect of each soil physico-chemical characteristic on the bacterial community calculated with the envfit function. ${\mathrm{NO}}_{\mathrm{3}}^{-}$: nitrate, ${\mathrm{NH}}_{\mathrm{4}}^{+}$: ammonium, OM: organic matter content, P: phosphorus, water: water content. The patch types are significantly different from each other (ANOSIM, R=0.28247; p value = 0.001). P, OM, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and ${\mathrm{NH}}_{\mathrm{4}}^{+}$.

3.3 Community composition

The community was mostly composed of Actinobacteria, Proteobacteria, Deinococcus-Thermus, Bacteroidetes, and Firmicutes (Fig. 2). The relative abundance for each phylum is detailed in Table A5. We focused on the results of the three main phyla: Actinobacteria, Deinococcus-Thermus, and Proteobacteria. Using pair-wise comparisons, we saw that shrub patches and nest patches had similar communities (no significant differences, p>0.05); therefore, we considered them to be single EE patches. For these patches, an average relative abundance of nest and shrub patches was used for statistical data. For the Actinobacteria phylum, patches with one EE had significantly lower relative abundance than barren patches (one EE: 9 % vs. Barren patch: 35 % p<0.005) or patches with two EEs (17 %, p value: 0.02). For the Deinococcus-Thermus phylum, barren patches had significantly lower relative abundance than patches with one or two EEs (Barren: 3 % vs. one EE: 25 % vs. two EEs: 9 %, p<0.05). A similar pattern was detected in the Proteobacteria phylum (Barren: 38 % vs. one EE: 44 % vs. two EEs: 39 %, p<0.05). Additionally, we looked at the next three most abundant phyla: Firmicutes, Bacteroidetes, and Chloroflexi. For Firmicutes, the relative abundance of this phylum was significantly higher in the Shrub & Nest patch than in the Barren and Shrub patches. For Bacteroidetes, the Nest patch had a significantly lower relative abundance than the other patches. For Chloroflexi, there was a significant decrease in relative abundance in the Shrub, Nest, and Shrub & Nest patches compared to the Barren patch. All p values can be found in Table A6. The class and order plots show differences between patch types. However, the resolution is not high enough to enable us to draw significant conclusions (Figs. B1 and B2).

Figure 2 Bar plot of the relative abundance (%) of the most abundant phyla in the soil microbial community in the dry season under different patch types (phyla with a relative abundance >0.05 %). The Proteobacteria have been separated into their classes (represented here in shades of green). The relative abundance of Deinococcus-Thermus increases when one EE is present, while the population of Actinobacteria decreases.

3.4 Functional prediction

The abundance of each gene group has been normalized to the 16S rRNA copy number for each genome. The functional prediction results focus on eight distinct gene groups: phototrophy, lithotrophy, organotrophy, DNA conservation, DNA repair, nitrogen cycle, sporulation, and ROS-damage prevention (listed in Table A7). Figure 3 shows the pattern of the obtained functions. It shows higher abundances of the gene group encoding for DNA conservation, DNA repair, nitrogen metabolism, ROS-damage prevention, sporulation, and phototrophy in patches associated with at least one EE compared to the barren patches (Table A8). Therefore, we analysed the results as barren vs. average of the other three patch types that were not significantly different from one another (Table A9), and significant differences (p<0.04) between barren and EE patches were detected. The genes related to lithotrophy only differed between patches with one EE and the barren patches (p<0.03), but patches with two EEs were similar to the barren plots. Finally, for genes related to the organotrophy, there were no significant differences between the patches (p>0.05).

Figure 3 Boxplots of the functional prediction of the 16S sequences. Each panel (boxplot) represents a different group of genes associated with a certain functionality. The full list of genes can be found in Table A7. The patch types are represented by distinct colours and patterns. The y axis is the abundance in copy number (CN) normalized to the 16S rRNA copy number for each genome.

4 Discussion

In desert environments, during the dry season, a large portion of the microbial community is dormant or shows reduced metabolic activity (Bay et al., 2018; Cordero et al., 2019; Lennon and Jones, 2011; Schulze-Makuch et al., 2018). However, the presence of EEs enhances the potential for functions related to metabolism and to survival functions (Fig. 3). EEs create havens of resources and water, which can be affiliated with the concept of resource islands (Schlesinger and Pilmanis, 1998). However, their individual, and combined, effects do not always lead to significant changes in the composition of the soil microbial community (Fig. 2). While the soil parameters might be modified by the presence of both EEs, the microbial community might take a longer time to change due to their slow turnover in the dry season. However, these communities experience more habitable conditions due to the modulating effects of the EEs on the environmental conditions. The increase in the activity of gene groups can be explained by an increase in nutrients in the joint EE patches (Table 1).

Both Actinobacteria and Deinococcus-Thermus were abundant in all patches, but their relative abundances were negatively correlated. Each phylum featured a dominant genus that is well adapted to stressful conditions: Rubrobacter dominated the barren soil, while Deinococcus dominated the EE patches (Fig. 2 and Table A5). Rubrobacter are specialized in surviving strong desiccation and low nutrients (Bull, 2011; Ferreira et al., 1999), showing high relative abundance in arid barren soils of the Negev desert highlands (Meier et al., 2021). Deinococcus are versatile organisms, highly adapted to a wide range of extremes, such as radiations, temperatures, and xerification (Chanal et al., 2006; Prieur, 2007; Slade and Radman, 2011). This versatility allows them to thrive in EE patches as they can better adapt to perturbations compared to Rubrobacter.

Only the combination of EEs resulted in significant changes (p values: Table A2) in ${\mathrm{NO}}_{\mathrm{3}}^{-}$, P, and, to a lesser extent, NH${}_{\mathrm{4}}^{+}$, pH, and OM (values: Table A1). When located under a shrub, ants can increase their seed consumption, which enhances the amount of leftovers around the nest (Wagner, 1997) and increases the concentrations of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and P. These macronutrients are important drivers of the biological processes, as they are often the limiting factors of microbial growth and activity in the terrestrial environments (FAO et al., 2020). The physico-chemical measures, including soil water content, OM, nitrogen, P, and pH, did not match the changes observed in bacterial composition or function (Tables A1, A2, and A9 and Fig. 1), as was previously reported (Angel et al., 2010; Bachar et al., 2012; Vonshak et al., 2018). Indeed, there was no significant link between the changes in the bacterial communities and the measured soil parameters (Table A10).

The EE patches analysed in this study share the same habitat and resources, but their impacts are distinct (Passarelli et al., 2014), and, thus, their joint impact is non-additive. The behaviour of each EE is important as it becomes a feature of the combined impact of both EEs (Alba-Lynn and Detling, 2008). However, the effect of both EEs together cannot be inferred from their individual environmental impact or from their mutual interaction (Gilad et al., 2004). Here, we investigated a sessile organism with a passive and slow impact (the perennial shrub) and compared it to a motile organism (the ants) with an active and transient impact. Ants may have both a short-term impact, through the seasonal accumulation of seeds and organic matter, and a lasting impact, due to the alternation of the nest mound which remains in the same place for decades (Wagner and Jones, 2004). We have previously proposed that the observed differences in communities could be mediated by microclimatic characteristics under shrub patches (Bachar et al., 2012). It has been reported that the desert dwarf shrubs affect the physical features of their immediate soil patch. Shrubs were shown to divert water flow and reduce evapotranspiration rates following rain events (Segoli et al., 2008; Whitford and Duval, 2002) and reduce temperature and radiation year round (Kidron, 2009). Likewise, ants aerate the soil, thus increasing infiltration during rain events (Berg and Steinberger, 2008), and mix the layers through bioturbation (Folgarait, 1998). Therefore, the prolonged water availability and altered physical conditions from the wet season may have lasting effects on the communities' structure (Baubin et al., 2019), shaping the composition and functions observed here (Figs. 2 and 3).

5 Conclusions

In this study, we focused on two EEs only, but there are many EEs in one ecosystem, and knowing their joint impact would help explain the nutrient turnover and the bacterial communities in this ecosystem. The main stress-resistant phyla (Actinobacteria and Deinococcus-Thermus) react differently to the presence of EEs. The presences of these EEs also lead to a higher potential activity in the microbial communities. However, even though they have similar impacts, when together, EEs have non-additive effects.

Appendix A

Table A1 Soil characteristics data. NH${}_{\mathrm{4}}^{+}$ and P show the highest discrepancy between Shrub & Nest patches and the other three types.

Table A2 P values of the Dunn test between patch types on the soil characteristics variables. Bold numbers are significant (<0.05).

Table A3 Results of the pair-wise adonis test between patch types. Bold numbers are significant (<0.05).

Table A4 Number of reads before and after the trimming stage and during the Dada2 stage.

Table A5 Relative abundance (%) of the taxonomic community per patch type.

Table A6 P values of the Dunn tests between patch types on the relative abundance of the five most abundant phyla. Bold numbers are significant (<0.05).

Table A7 List of the genes used for function prediction ordered by groups and subgroups.

Table A8 Abundance (in copy number (CN)) of each patch type within each group of gene.

Table A9 Chi-square values and p values of the Dunn tests between patches done on the functional prediction results. Bold numbers are significant (<0.05).

Table A10 Results of the adonis analysis of the impact of soil parameters on the bacterial community.

Appendix B

Figure B1 Bar plot of the relative abundance (%) of the most abundant classes in the soil microbial community in the dry season under different patch types (classes with a relative abundance >0.05 %). The resolution is too low to draw significant conclusions.

Figure B2 Bar plot of the relative abundance (%) of the most abundant orders in the soil microbial community in the dry season under different patch types (orders with a relative abundance >0.05 %). The resolution is too low to draw significant conclusions.