Climate controls prokaryotic community composition in desert soils of the southwestern United States

Abstract

Aridisols are the dominant soil type in drylands, which occupy one-third of Earth's terrestrial surface. We examined controls on biogeographical patterns of Aridisol prokaryotic (bacterial and archaeal) communities at a regional scale by comparing communities from 100 Aridisols throughout the southwestern United States using high-throughput sequencing of the 16S rRNA gene. We found that microbial communities differed among global biomes and deserts of the Southwest. Differences among biomes were driven by differences in taxonomic identities, whereas differences among deserts of the Southwest were driven by differences in relative sequence abundance. Desert communities were dominated by Actinobacteria, Proteobacteria and Crenarchaeota, supporting the notion of a core set of abundant taxa in desert soils. Our findings contrast with studies showing little taxonomic overlap at the OTU level (97% sequence similarity) across large spatial scales, as we found ∼90% of taxa in at least two of the three deserts. Geographic distance structured prokaryotic communities indirectly through the influence of climate and soil properties. Structural equation modeling suggests that climate exerts a stronger influence than soil properties in shaping the composition of Aridisol microbial communities, with annual heat moisture index (an aridity metric) being the strongest climate driver. Annual heat moisture index was associated with decreased microbial diversity and richness. If the Desert Southwest becomes hotter and drier as predicted, these findings suggest that prokaryotic diversity and richness in Aridisols will decline.

Aridisols, biodiversity, biogeography, 16S rRNA gene, structural equation model, climate

INTRODUCTION

In arid environments, extreme temperature fluctuations, soil moisture content and plant traits influence soil microorganisms (Collins et al.2008; Bell et al.2013; van Gestel, Reischke and Bååth 2013). Aridisols (from the Latin aridus, meaning ‘dry’), the third most abundant soil type in the world (Walther 2014), are characterized by water deficiency (i.e. evapotranspiration far exceeds precipitation for a majority of the year; Brady and Weil 2008), low organic matter content, accumulations of inorganic minerals (e.g. salt, gypsum, and carbonates), and short-statured vegetation that often includes bunchgrasses and shrubs. While there is a large body of work on a subset of desert soil microorganisms—for example, the autotrophic component of biological soil crusts (e.g. Belnap and Gardner 1993; Garcia-Pichel and Belnap 1996; Belnap 2003; Yeager et al.2004)—less is known about Aridisol microbial communities, particularly regarding the factors controlling their biogeographic distributions at large spatial scales (but see Fierer et al.2012; Wang et al.2012a; Maestre et al.2015). The limited number of molecular examinations of prokaryotic communities in desert soils indicates they are taxonomically and functionally distinct from those found in other biomes (Drenovsky et al.2010; Fierer et al.2012).

Soil bacterial communities across ecosystems are often strongly influenced by pH, which explains a large proportion of the variance in soil bacterial diversity at local, regional and continental scales and across ecosystem types (Fierer and Jackson 2006; Lauber et al.2009; Griffiths et al.2011). Most Aridisols, however, have a neutral to alkaline pH; with little variation in the pH of these arid soils, it is likely that other environmental parameters influence patterns in microbial community composition. Although some evidence suggests that plant communities shape microbial communities at the local scale in non-desert ecosystems (Zak et al.2003; Wardle et al.2004; Gao et al.2013), this is not always the case (Ramirez et al.2010; Bastida et al.2013; McHugh and Schwartz 2015), even at larger, continental scales (Fierer and Jackson 2006).

Aridity has been invoked as an important variable shaping soil microbial communities in desert soils for bacteria (Ben-David et al.2011; Köberl et al.2011; Maestre et al.2015) and fungi (Tedersoo et al.2014; Maestre et al.2015) due to the typically negative relationship between aridity and resource availability (e.g. water and nutrients). A recent study (Maestre et al.2015) used a structural equation modeling (SEM) approach on high-throughput sequencing data to separate the direct and indirect effects of aridity on microbial diversity and total abundance in arid systems globally; however, metrics of diversity and abundance do not address whether communities are compositionally distinct along climate and soil gradients.

The purpose of this study was to assess the spatial structuring of soil microbial communities across the arid Southwest of the United States, and to determine what factors drive these patterns. We formulated five a priori hypotheses. First, we hypothesized that (1a) microbial communities are distinct among global biomes and (1b) among deserts of the Southwest. Additionally, we predicted that microbial community separation would be more pronounced among biomes than deserts. To test this set of hypotheses, we integrated our sequencing dataset with publicly available sequences from a cross-biome prokaryotic community comparison (Fierer et al.2012) to evaluate how the composition of desert prokaryotic communities compared to those found in other biomes. Second, we hypothesized that (2) geographic distance would structure soil microbial communities of the Desert Southwest. Specifically, we predicted that microbial communities in geographic proximity would be compositionally more similar. We made this prediction because this distance-decay pattern has been observed in all domains of life (e.g. Nekola and White 1999; Green et al.2004; Horner-Devine et al.2004). Finally, we hypothesized that (3a) climate would be the most important variable structuring soil microbial communities of the Desert Southwest because it has been identified as a common driver of soil microbial communities at regional and continental scales (Ettema and Wardle 2002; Fierer and Jackson 2006; Lauber et al.2008, 2009; Fierer et al.2009). Furthermore, because our study examined microbial communities of Aridisols, we hypothesized that (3b) among climate variables, aridity would be the most important predictor variable. We addressed this set of hypotheses using SEM to determine whether the global pattern that aridity structures microbial communities in drylands was also observed at a regional scale. By accounting for other variables associated with aridity, we assessed whether aridity was a causal factor influencing microbial communities and evaluated the relative importance of aridity in comparison to other factors. By considering soil microbial communities across a large geographic extent in extreme Aridisols, where the controls on microbial distributions are likely to be different from temperate ecosystems, this work adds substantially to our understanding of microbial biogeography and provides new context for identifying the mechanisms controlling distributions of microbial communities (see Fierer et al.2012 and Maestre et al.2015 for other examples of large-scale studies in arid systems).

MATERIALS AND METHODS

Study area

The southwestern United States encompasses extensive drylands containing portions of the Mojave, Sonoran and Great Basin Deserts. These three deserts are distinct in terms of climate variables, such as temperature, aridity, and the amount and timing of precipitation (Huxman et al.2004), and also in terms of vegetation (Shreve 1942) and soil characteristics (Table 1). Soils were sampled in the southwestern United States, across an area of 460 000 km² (Fig. S1a, Supporting Information). The median pairwise distance between any two sites was 350 km linea recta (minimum = 0.4 km; maximum = 888 km; Fig. S1b). The study area spanned three of the four US deserts (Mojave Desert, Sonoran Desert and Great Basin Desert), encompassing 100 Aridisol sites selected using soil survey maps from the Natural Resources Conservation Service (U.S. Department of Agriculture). Across these sites, soil texture varied from sandy clay loam to loamy sand.

Table 1.

Climate and soil chemical characteristics of each desert in our study area.

	Climate				Soil characteristics
Desert	MAP (mm)	MAT (°C)	Elevation (m)	AHM	pH	SOM (%)	NO3− (μg g−1)	P (μg g−1)	K (μg g−1)	Ca (meq/L)	Mg (meq/L)	Na (meq/L)
Great Basin	268 ± 9.9	11.6 ± 0.3	1599 ± 34	88 ± 4	7.8 ± 0.1	1.3 ± 0.1	3.1 ± 0.8	8.8 ± 1.4	341 ± 47	6.4 ± 1.2	10.5 ± 9.0	6.3 ± 5.4
Mojave	145 ± 11	20.0 ± 0.5	646 ± 73	230 ± 19	7.6 ± 0.1	0.5 ± 0.1	3.2 ± 0.5	9.1 ± 0.8	184 ± 21	3.0 ± 0.2	0.7 ± 0.09	0.7 ± 0.5
Sonoran	248 ± 20	20.6 ± 0.4	634 ± 86	161 ± 22	7.7 ± 0.1	0.8 ± 0.1	37 ± 18	7.3 ± 1.5	479 ± 75	9.5 ± 4.3	4.0 ± 2.1	11.4 ± 7.5

Data are averaged across sampling sites in the Great Basin Desert (n = 57 sites), Mojave Desert (n = 20 sites) and Sonoran Desert (n = 23 sites). Climate data for MAP, MAT, elevation and AHM, calculated as (MAT + 10)/(MAP/1000), were obtained from www.climatewna.com for the 1981–2010 period using the coordinates of each sampling location. Values represent means ± standard error.

Climate and vegetation data

Climate data for individual sites were obtained from an open-data platform hosted by the Center for Forest Conservation Genetics at the University of British Colombia (www.climatewna.com) for the period from 1981 to 2010. Climate data were generated via PRISM high-resolution spatial climate data for the conterminous United States (Wang et al.2012b). Sites ranged in elevation from 80 to 2197 m, mean annual temperature (MAT) from 6.1°C to 23.0°C, and mean annual precipitation (MAP) from 77 to 420 mm. Great Basin Desert samples were collected at an elevation ∼2.5 times greater than that of the other two deserts. The Great Basin Desert is the coolest and wettest desert in our study, whereas the Mojave Desert is the driest, but just as hot as the Sonoran Desert (Table 1). The annual heat-moisture index (AHM) for our sites ranged from 39.3 to 422. AHM is a biologically relevant indicator of aridity that is calculated as the ratio of heat and precipitation (Wang et al.2012b): (MAT + 10)/(MAP/1000). Soils in areas of high AHM have hot and dry soil environments that experience large daily fluctuations in temperature, compared to soils under cooler and wetter conditions (Nobel and Geller 1987). In essence, AHM is an indicator that captures how much soils can store and lose heat, making it a relevant environmental variable for soil microbial communities (van Gestel et al.2011).

Data for vegetation type and cover were collected from LANDFIRE (http://www.landfire.gov/). The most common existing vegetation types were mixed desert scrub (e.g. with cacti, palo verde, creosotebush), desert grassland, shrubland and piñon-juniper woodland. Vegetative cover ranged from 0% to 59% (mean = 31%).

Soil sample collection

Soil samples were collected from 16 to 20 June 2014 from areas between plants. Desert plants can have extensive root systems (Rundel and Nobel 1991), so plants can still influence soil in interplant spaces. Sampling occurred during the dry period preceding monsoon rains, and no precipitation was observed across the study region during the period of collection. To ensure that we captured the variability in soil at each site, we collected five replicate samples, with a minimum distance of 10 m between each replicate (n = 100 sites, 500 samples). Climate and vegetation data for these sites (described above) were taken at GPS coordinates corresponding to the center of the sampling area for each site. Samples were collected after removing visible desert varnishes or biological soil crusts when present. Samples were taken from the top ∼5 cm of soil using sterile scoops and transferred into sterile Falcon tubes. Samples were stored on dry ice and transported to the laboratory, where DNA was extracted.

Physical and chemical soil properties

The five replicates at each site were combined to create a composite soil sample. The composite samples were sent to Colorado State University's Soil, Water, and Plant Testing Laboratory to determine soil physical and chemical properties. Soil chemistry variables were selected based on known or predicted biological influences they have on soil microbes. Soil pH, electrical conductivity (mmhos cm⁻¹), organic matter (%) and sodium adsorption ratio were measured. Nitrate (NO₃⁻-N), P, K, Zn, Fe, Mn and Cu were assessed using the NH₄HCO₃-diethylenetriaminepentaacetic acid (AB-DTPA) method (with values provided in ppm). Additionally, the exchangeable bases Ca, Mg and Na were extracted using ammonium acetate (values provided in meq L⁻¹). All macro- and micronutrients were analyzed by inductively coupled plasma-atomic emission spectroscopy.

DNA extraction and sequencing of 16S rRNA genes

Because soil microbial communities are spatially more variable than soil chemical properties (Ettema and Wardle 2002), we sequenced all replicates at a site to better capture the true diversity within these soils. We then summarized compositional data across the replicates to obtain representative site-level information. Total genomic DNA was extracted using the PowerLyzer PowerSoil DNA extraction kit (MoBio Laboratories, Carlsbad, CA). Extracted DNA was quantified by PicoGreen (Molecular Probes, Eugene, OR) fluorescence and normalized to 1 ng μL⁻¹ prior to amplification.

Bacterial communities were prepared for high-throughput sequencing on an Illumina MiSeq using a two-stage PCR approach (Berry et al.2011). Each sample was first amplified using primers 515f and 806r, which target the hypervariable v4 region (Bates et al.2011). This was done in triplicate 8 μL PCR amplifications containing 1 mM of each primer, 0.01 U μL⁻¹ Phusion HotStart II Polymerase (Thermo Fisher Scientific, Waltham, MA), 1X Phusion HF buffer (Thermo Fisher Scientific), 3.0 mM MgCl₂, 6% glycerol and 200 μM dNTPs. PCR conditions were 95°C for 2 min; 15 cycles of 95°C for 30 s, 55°C for 30 s and 60°C for 4 min. Initial PCR products were checked on a 1% agarose gel. Triplicates were then pooled, diluted 10-fold and used as template in the subsequent tailing reaction with region-specific primers that included the Illumina flowcell adapter sequences and a 12 nucleotide Golay barcode (15 cycles identical to initial amplification conditions). Products of the tailing reaction were purified with carboxylated SeraMag Speed Beads (Sigma-Aldrich, St. Louis, MO) at a 1:1 (volume to volume) ratio (Rohland and Reich 2012) and quantified by PicoGreen fluorescence. Equal quantities of the reaction products were then pooled, and the library was bead-purified once again (1:1 ratio) and quantified by qPCR using the Library Quantification Kit for Illumina (Kapa Biosciences, Woburn, MA). The amplicon library was denatured and loaded at 11 pM (including a 30% PhiX control) onto an Illumina MiSeq instrument (San Diego, CA) using 2 × 150 paired-end read chemistry at Northern Arizona University's Environmental Genetics and Genomics Laboratory. A total of three amplicon pools were sequenced, returning 14.24, 15.6 and 14.75 million reads passing filter for sequencing runs 1, 2 and 3, respectively. All sequences were submitted to MG-RAST, project ID 14035.

Sequence processing and community comparisons

Sequence data were analyzed using Quantitative Insights into Microbial Ecology v 1.7 (QIIME, Caporaso et al.2010b). Read 1 of each library was demultiplexed, and the reads and downloaded libraries were combined for subsequent analysis. For quality filtering, the default score was changed from 25 to 30. Open reference OTU picking was performed at 97% identity using uclust (Edgar 2010). The most abundant sequence for each OTU was aligned with PyNAST (Caporaso et al.2010a) against the Greengenes v13_5 database (DeSantis et al.2006), and taxonomy was assigned using Ribosomal Data Project classifier (Wang et al.2007). A phylogenetic tree was built using FastTree (Price, Dehal and Arkin 2010). Additional quality filtering was applied, discarding any OTUs that accounted for less than 0.005% of the total sequences (Bokulich et al.2013). After quality filtering, 19 266 450 16S rRNA gene sequences remained, representing 2 744 unique OTUs. Of the initial 100 geographic locations, 12 were removed due to low-sequencing depth in the samples. To remove heterogeneity in sampling depth, the bacterial libraries were rarefied to a depth of 5 000 sequences prior to calculations of alpha (Chao1 and Shannon) and beta diversity. Our quantitative analyses utilized site-level community data, based on the five replicates per site.

Crossbiome prokaryotic community composition

To compare hot desert prokaryotic community data to community data from other biomes, additional 16S rRNA gene libraries from the study ‘Cross-biome metagenomic analyses of soil microbial communities and their functional attributes’ (Fierer et al.2012) were accessed (MG-RAST project ID 10307). These libraries were from soils collected at 16 sites: three hot deserts (Mojave and Chihuahuan), six Antarctic cold deserts, two tropical and temperate forests, and one prairie, tundra and boreal forest. Our Aridisol dataset was collected using the same sampling methodology, primers, sequencing approach and principal coordinates analyses as those utilized by Fierer et al. (2012), with the exception that we did not sample soils during the peak of the plant growing season. However, desert soil microbial community composition has been shown to be temporally stable (Armstrong et al.2016), and hence, the timing of sampling would have little effect. To further enable the best comparison between the data sets, OTU picking and taxonomic assignments were performed on the combined data set.

Data analysis

Statistical analyses were conducted in R (3.0.2, R Core Team 2012) unless otherwise noted. For the prokaryotic cross-biome analyses, principal coordinates analyses based on weighted Unifrac distances were conducted in PRIMER v. 6 using pairwise distances between prokaryotic communities from different biomes (Fierer et al.2012 dataset) and hot deserts (our own study plus desert communities from Fierer et al.2012). To test if hot desert communities differed significantly from one another, permutational multivariate analysis of variance (PerMANOVA) and pairwise t-tests were performed in PRIMER. To determine whether site-level climate and soil characteristics differed among the Great Basin, Mojave and Sonoran Deserts, we obtained P-values from randomization tests (5 000 bootstrap iterations). Bonferroni corrections were applied to the P-values because we used three pairwise comparisons per variable. Chao1 and Shannon diversity indices were calculated using the Vegan package in R (Oksanen et al.2013). To assess the degree of endemism, we determined presence and absence of taxa at the OTU level for each desert. A Venn diagram was then constructed to visualize the shared taxa (present in more than one desert) and unique taxa (found only in one desert).

We created an interpolated map of community similarity by first determining how similar the microbial communities were at the OTU level. The dissimilarity matrix was based on Bray-Curtis distances (Legendre and Gallagher 2001) among sites using relative abundances from all replicates per site. We then performed non-metric dimensional scaling (NMDS) analysis (‘ecodist’ R package; Goslee and Urban 2007) on the dissimilarity matrix to summarize the data into one dimension using 1 000 iterations, and used the scores of the iteration that produced the lowest stress. The NMDS scores were interpolated across the study area using inverse distance weighting (‘gstat’ package; Pebesma 2004), and visualized using color (‘ggmap’ R package; Kahle and Wickham 2013), such that sites with similar microbial communities had similar color on the map. We also created an interpolated map using the original data of the climate variable that was the best predictor of community similarity (i.e. AHM), as determined from the SEM approach.

For the SEM, we used two NMDS axes to describe the soil microbial community because this was the least number of axes that resulted in the greatest decrease in stress (final stress = 0.29, r² = 0.71). Because of the constraints of SEM related to our limited sample size, we used partial-Mantel tests (‘ecodist’ package) to screen for predictor variables that went into our a priori SEM. We did this in two ways. First, we tested whether geographic distance controlled microbial community composition (either directly or through environmental variables). Second, we tested whether microbial community composition was controlled by environmental, climatic or plant community variables, with geographic distance held constant. Additionally, to determine the relative importance of the three types of variables tested (i.e. climate, soil properties and plant communities), we conducted separate analyses with each variable type grouped into a composite variable. Partial-Mantel tests allowed us to examine associations between variables while accounting for the influence of all other variables. We tested whether dissimilarity matrices (Euclidean distance) for geographic distance, climate variables, soil property variables, vegetative cover and vegetation type correlated with the microbial community dissimilarity matrix (Bray-Curtis distance). The type of distance matrix used was based on long-standing recommendations (McCune, Grace and Urban 2002). For all partial-Mantel tests, dissimilarity matrices were compared and P-values were calculated using distributions estimated from 999 permutations (Jackson and Somers 1989). Partial-Mantel correlation P-values were corrected for multiple comparisons using the False Discovery Rate technique in R (Benjamini and Hochberg 1995; Benjamini 2010).

The SEM was conducted (IBM® SPSS® Amos 22.0.0; Arbuckle 2006) because it combines factor analyses with multiple regression (Bollen 1989), thereby allowing us to determine if our causal inferences about what controls the soil microbial community were supported by our data. We generated an a priori model that showed our hypothesized interactions among soil properties, climate and the microbial community. Vegetation characteristics and geographic distances were excluded because partial Mantel tests showed these variables did not correlate well with the microbial community dissimilarity matrix. We tested our proposed causal links by confronting our model with observations, which allowed us to falsify or support our proposed model and evaluate the fit between our model and the observational data (Grace and Pugesek 1998). Next, we developed a measurement model (Bollen 1989; Bowker et al.2005) describing the hypothesized relationships among our measured variables and underlying, unmeasured factors (i.e. latent variables). This measurement model was then tested iteratively using the bootstrap goodness-of-fit test until a satisfactory fit with the data was achieved (Bollen and Stine 1992). Our final model (see summarized representation in Fig. 5) included MAT, MAP, elevation and AHM, which were highly intercorrelated (all r ≥ 0.60−0.97) and therefore suitable as indicators for the latent ‘climate’ variable, and NO₃⁻, pH and soil organic matter (SOM) as soil property variables.

RESULTS

Climate and soil chemical characteristics

The Sonoran and Mojave Deserts had similar MATs, which were at least 8.4°C higher than the Great Basin Desert (Table 1). The Mojave Desert was the driest, receiving up to 58% of the MAP of the other two deserts (Table 1). Thus, the aridity index, AHM, was lowest in the cold and wet Great Basin Desert and greatest for the hot and dry Mojave Desert (Table 1).

The Great Basin Desert had significantly higher SOM (P < 0.001) than the other two deserts (Table 1). Despite higher SOM, soil pH was similar for all three deserts (mean = 7.7). Soil NO₃⁻ was highest for the Sonoran Desert (Table 1), but this was due to three sites (of 23) having exceptionally high concentrations (≥190 μg g⁻¹) that greatly influenced the mean. Soil P, Ca, Mg and Na concentrations were similar among the three deserts, exhibiting high within-desert variability for the latter two cations (Table 1). The only nutrient to differ among deserts was K, with lower values in the Mojave Desert (Table 1).

Cross-biome and cross-desert microbial community composition

In support of hypothesis 1a, prokaryotic communities from our Aridisols study were compositionally distinct from communities found in other biomes, but similar to desert soil microbial communities previously characterized within and near our study region (Fierer et al.2012; Pseudo-F_94,2 = 8.7, P = 0.001; Fig. 1a). Compared to non-desert biomes, our Aridisol communities had lower abundances of Acidobacteria and Verrucomicrobia. Instead, Aridisol microbial communities from our study were dominated by Actinobacteria (32.9%), Proteobacteria (22.5%) and Crenarchaeota (9%; Fig. S2a, Supporting Information), and all other phyla were relatively rare (<9%). In support of hypothesis 1b, there were also significant differences among microbial communities of the three deserts (Pseudo-F_2,382 = 51.4, P = 0.001), with the hottest and driest Mojave Desert having communities of distinct composition. Compared to the other two deserts, communities from the Mojave Desert had the lowest relative abundances of Actinobacteria, Crenarcheaota, Acidobacteria and Firmicutes, but highest relative abundances of Proteobacteria and Cyanobacteria (Fig. S2a). The higher relative abundances of Proteobacteria in the Mojave Desert were driven by higher abundances of Alphaproteobacteria (Fig. S2b), which was the most dominant class of Proteobacteria across the three deserts. The most notable difference in Alphaproteobacteria was that the order Rhizobiales (members fix atmospheric nitrogen and have beneficial relationships with plants) comprised 64% of that class in the Mojave Desert and accounted for 16.8% of the total prokaryotic community. The dominant order in the Great Basin and Sonoran Deserts, making up about 10% of the total prokaryotic community, was Nitrososphaerales within the phylum Crenarchaeota. Actinobacteria was the second most dominant phylum in the Mojave Desert, and the dominant phylum in the Sonoran and Great Basin Deserts. The most notable compositional difference within this phylum was that the order Actinomycetales comprised half of the Actinobacteria in the Mojave Desert, whereas it comprised only 28% of the Actinobacteria in the other two deserts (Fig. S2c). Other dominant orders in the Actinobacteria were Rubrobacterales (23−28%) and Solirubrobacterales (11−19%; Fig. S2c).

Figure 1.

Ordination plots derived from principal coordinates analysis of weighted Unifrac distances. (A) Aridisol (hot desert) prokaryotic communities (open symbols) are contrasted with communities from other biomes (closed symbols) using sequencing data acquired from Fierer et al. (2012). (B) Aridisol communities are isolated and assigned to the four United States deserts for comparison of community composition among deserts in our study (Great Basin, Mojave and Sonoran; open symbols) and those from Fierer et al. (2012) (Chihuahuan and Mojave; closed symbols).

Diversity and degree of endemism

Community differences were also apparent in terms of diversity and richness, with the Great Basin Desert having the highest phylotype diversity (both Chao1 and Shannon; Fig. 2a) of the three deserts. In contrast, the diversity in the Sonoran Desert was similar to the Mojave Desert, but because of large within-desert variability in diversity, this resulted in lower Chao1 diversity, but not Shannon diversity, compared to the Great Basin Desert. Observed OTU richness was also highest for the Great Basin Desert (3 659 OTUs) compared to the other two deserts (2 369 OTUs for the Mojave Desert and 3 485 OTUs for the Sonoran Desert; Fig. 2b). Eight per cent (288 OTUs) of the phylotypes found in the Great Basin Desert were unique to that desert, whereas only 0.04% (1 OTU) was unique to the Mojave Desert (Fig. 2b). Furthermore, we found little support for high endemism because community differences were related to differences in the relative abundance of OTUs rather than by presence or absence of particular OTUs, as a high percentage of OTUs (87% of phylotypes) were found in at least two deserts (Fig. 2b) and, of these, 70% (i.e. 61% of all OTUs) were found in all three deserts. Shared OTUs among the deserts spanned all phyla, and they exhibited similar relative abundances at the phylum level compared to whole communities (shared and unique) in each desert (i.e. compared to Fig. S2a; data not shown).

Figure 2.

Differences in microbial communities in Great Basin, Mojave and Sonoran Deserts in terms of (A) diversity using Shannon and Chao1 indices based on relative abundances (±standard deviation), and (B) the number of shared OTUs among deserts and the number of OTUs unique to each desert based on presence/absence.

Drivers of Aridisol community composition

In support of hypothesis 2, geographically similar sites in the southwestern United States supported more similar microbial communities (Mantel r = 0.36; P < 0.001). However, geographic distance did not directly influence the composition of soil microbial communities in the Desert Southwest, as this correlation was not significant when the partial correlations of other variables were considered (partial Mantel r = 0.0052; P = 0.45); rather, the influence of geographic distance was indirect, through associations with climate (partial Mantel r = 0.30; P < 0.001) and soil properties (partial Mantel r = 0.49; P < 0.001). Furthermore, the influences of climate and soil properties on microbial communities were much stronger than the influence of vegetation (partial Mantel r = 0.027; P = 0.35). When we examined the effects of climate and soil properties on vegetation (using composite variables), we found that vegetation cover was weakly associated with soil properties (partial Mantel r = 0.069, P = 0.094), but not climate (partial Mantel r = –0.0078, P = 0.42). Vegetation type, however, was associated with neither soil properties (partial Mantel r = 0.018, P = 0.37) nor climate (partial Mantel r = 0.035, P = 0.37). Because geographic distance and vegetation were not predictors of microbial community composition, we simplified our SEM by removing them. Thus, our SEM examined the relative influence of climate versus soil properties on prokaryotic community composition (Fig. 3).

Figure 3.

Simplified diagram of structural equation modeling results illustrating that climate (latent variable, circle), indirectly and directly, has a larger influence on soil microbial community composition than soil properties (composite variable, gray box) directly. Results also show how climate affects the individual properties of soil (NO₃⁻, pH and SOM), with SOM being the strongest affected by climate. Highly intercorrelated indicators of the latent variable ‘climate’ were MAT, MAP, elevation and AHM. The soil microbial community dissimilarity was determined using NMDS of the Bray-Curtis distance matrix of OTU-level community composition of two axes (NMDS 1 and NMDS 2). The high P-value indicates a good model fit. Standardized regression weights are shown for each path, with arrow widths indicating the strength of the relationships; dashed arrows indicate non-significant relationships.

In the simplified SEM, we allowed climate to directly and indirectly (through soil properties) influence microbial community composition. Our model explained 50% of the variance of NMDS axis 1 and 25% of the variance of NMDS axis 2. In support of hypothesis 3a, we found stronger direct and indirect effects of climate on community composition in comparison to soil properties (Fig. 3). There was only a weak direct influence of soil properties on microbial community composition once we partitioned out the strong direct influence of climate. Soil properties affected NMDS 1 (β_stand = –0.26), but not NMDS 2 (β_stand = –0.07; Fig. 4). The influence of soil properties on NMDS 1 was mostly driven by SOM (β_stand = 0.22) and NO₃⁻ (β_stand = –0.16), and not pH (β_stand = –0.08). Soil properties themselves were influenced strongly by climate, most notably for SOM and pH (Fig. 3). In support of hypothesis 3b, the strongest climate predictor of SOM was AHM (β_stand = –0.97): high AHM (i.e. extreme hot and dry conditions) coincided with low SOM. The strongest predictor of pH was MAP (β_stand = –0.78), with MAP negatively associated with pH.

Figure 4.

Interpolated map of community composition. The map is based on the inverse-distance-weighted interpolation of the NMDS scores onto one axis using the Bray-Curtis distance matrix of OTU-level community composition. Similar colors indicate similar prokaryotic communities. Sampling sites are indicated by black symbols.

Climate showed the greatest direct and indirect (via soil properties) influences on community composition, affecting both NMDS axes of the community ordination (with absolute values of direct paths with a standardized regression coefficient (β_stand) > 0.5 on NMDS 1 and NMDS 2). The most important climate variables directly influencing community composition were AHM (β_stand of –0.65 on NMDS 1) and MAT (β_stand of –0.55 on NMDS 2). When visualized spatially, geographic regions with similar AHM also had similar soil community composition (Fig. 4). Notably, the hot and dry Mojave Desert, with high AHM, had soil microbial community composition scores that were distinct from the other two deserts (Fig. 4).

Associations between relative abundances of microbial taxa and AHM

Because our SEM showed AHM strongly influenced microbial composition, we regressed relative abundances of the dominant phyla against AHM in order to understand the responses of specific taxa (Fig. 5). Relative abundances of Cyanobacteria increased with increasing AHM, associating with hotter and drier sites, whereas the relative abundances of several phyla (e.g. Crenarcheaota, Acidobacteria, Gemmatimonadetes, Planctomycetes, Verrucomicrobia and Nitrospirae) decreased with AHM, associating with moister and cooler conditions (Fig. 5). Of these, Acidobacteria, Planctomycetes and Verrucomicrobia showed the clearest, non-linear, declines with increasing hotter and drier conditions. Four phyla had a unimodal relationship with AHM (i.e. highest relative abundances at an intermediate AHM): Proteobacteria (class Alphaproteobacteria), Chloroflexi, Armatimonadetes and FBP. The remaining phyla showed no clear relationship with AHM: Actinobacteria, Proteobacteria (classes Beta-, Delta-and Gammaproteobacteria), Firmicutes, Bacteroidetes, Thermi, BRC1 and Euryarchaeota.

Figure 5.

Relative abundances of phyla with AHM across the three deserts. Panels are ordered in decreasing dominance (i.e. average relative abundances across samples), except for the second most dominant group, the Proteobacteria, which were subdivided into Alpha-, Beta-, Delta- and Gammaproteobacteria and grouped together following the most dominant phylum, the Actinobacteria. The line represents a locally weighted polynomial regression using weighted least squares, a LOESS curve fitting procedure, with the shaded area representing the 95% confidence interval.

DISCUSSION

Cross-biome and cross-desert comparison

Our study indicates that Aridisols have distinct microbial communities from those found in soils of other biomes (supporting hypothesis 1a). Non-desert biomes contained higher relative abundances of Acidobacteria and Verrucomicrobia compared to desert biomes in this study and other desert surveys (Andrew et al.2012; Fierer et al.2012; Maestre et al.2015). Consistent with other studies (Andrew et al.2012; Maestre et al.2015), our desert communities were dominated by Actinobacteria and Proteobacteria. Together, these phyla comprised over half of the total taxa we detected. Generally, and even in hyperarid soils (Neilson et al.2012), Actinobacteria dominate soil communities.

Specialized survival strategies are necessary for prokaryotes to withstand the physiological stresses imposed by low water potential and nutrient concentrations, intense heat and solar radiation, and soil instability, which are characteristic of desert soils. Actinobacteria can grow in low humidity environments (Doroshenko et al.2005) and are known to tolerate drought well (Barnard, Osborne and Firestone 2013; McHugh, Koch and Schwartz 2014). This is in part because of the filamentous morphology of some Actinobacteria, enabling them to exploit soils at low soil water potential with very thin and discontinuous water films (Stark and Firestone 1995) or with water in distant soil pores (Torsvik and Øvreås 2008). Additionally, their filamentous nature aids in the formation of stable soil aggregates (Torsvik and Øvreås 2008), thereby increasing soil stability. Other drought adaptations include spores that may persist for a long time in the absence of water (Chen and Alexander 1973). The Firmicutes, which can form drought-tolerant endospores, were also well represented in Aridisols from our study (Bueche et al.2013). Besides drought, prokaryotes in hot deserts must withstand large temperature fluctuations, with average diurnal temperatures ranging from 12°C to 45°C, compared to 4°C−8°C in coastal and temperate regions (Dai, Trenberth and Karl 1999; van Gestel, Reischke and Bååth 2013). Phyla including Crenarchaeota and Chloroflexi are able to withstand substantial temperature ranges (Cavicchiolo 2006; Hatzenpichler et al.2008) and extreme heat (Gladden et al.2011). Collectively, these traits likely contribute to the establishment of a core set of abundant taxa in desert soils (Andrew et al.2012), and enable microbes to withstand extreme environmental stress.

Prokaryotic communities in the Mojave Desert were the most distinct of the three deserts (supporting hypothesis 1b). Surprisingly, of the deserts in our study, the hottest and driest Mojave Desert was the sole desert dominated by Proteobacteria and not by Actinobacteria, even though Actinobacteria are typically the most abundant phylum in desert soils (Andrew et al.2012; Neilson et al.2012; Maestre et al.2015). Not only were community differences in the Mojave Desert reflected by this difference in the dominant phylum, but also in the dominant order; Rhizobiales (phylum Proteobacteria) dominated in the Mojave Desert, whereas Nitrososphaerales (phylum Crenarchaeota) was the dominant order in the other two deserts. Members of the Rhizobiales are plant symbionts commonly occurring in semiarid ecosystems, even in bare soil (Hortal et al.2013), and they supply usable nitrogen to their hosts by fixing it from the atmosphere. Members of the Nitrososphaerales, which have only been described since 2005 (Könneke et al.2005; Treusch et al.2005), are ammonia oxidizers that can fix atmospheric CO₂ and are hence important in both carbon and nitrogen cycling (Kerou et al.2016). Interestingly, the second most dominant order, Actinomycetales (phylum Actinobacteria), was present in equal relative abundances across all deserts. Members of this order can degrade lignin and other recalcitrant litter (Heuer et al.1997) and are thus important decomposers in deserts. The Mojave Desert was further compositionally distinct from the other two deserts by having lower relative abundances of Acidobacteria and Firmicutes, but higher relative abundances of Cyanobacteria.

The Acidobacteria, Planctomycetes and Verrucomicrobia showed the strongest, non-linear, declines with increasing AHM (i.e. they were associated with moister and cooler conditions). These phyla are more common in non-desert biomes, such as forests and tundra (Fierer et al.2012). Only Cyanobacteria increased in relative abundance with increasing AHM. Cyanobacteria are inhabitants of all desert soils (Hagemann et al.2015), where they are important producers. They are also key constituents of biological soil crusts (Belnap 2003; Hagemann et al.2015), contributing to soil stabilization and nitrogen fixation. Frost heaving that leads to tall, pinnacled crusts makes biological soil crusts more pronounced in the cooler Utah desert soils, whereas they are less conspicuous in hotter and drier climates (Belnap 2003). Perhaps intentionally avoiding visible soil crusts in our sampling approach biased against Cyanobacteria in Utah. This could explain an apparent positive relationship of Cyanobacteria with AHM. The four phyla that had unimodal relationships with AHM (i.e. highest relative abundances at intermediate AHM) were Proteobacteria (class Alphaproteobacteria), Chloroflexi, Armatimonadetes and FBP. Other phyla were found in similar relative abundances across the range of AHM in our study. These ‘insensitive’ phyla included the dominant, drought-tolerant filamentous Actinobacteria, and less dominant phyla (i.e. typically at relative abundance <5%), such as Proteobacteria (classes Beta-, Delta- and Gammaproteobacteria), Firmicutes, Bacteroidetes, Thermi, BRC1 and Euryarchaeota. These findings suggest that these prokaryotic groups cope equally well across environmental gradients and therefore have a wider environmental niche.

Low endemism

Our work contrasts with studies showing little taxonomic overlap at the OTU level across large spatial scales (Fulthorpe et al.2008; Griffiths et al.2011). Instead, we determined that nearly 90% of the OTUs were found in at least two deserts, with 70% of these common to all deserts, thereby demonstrating low prokaryotic endemism. Even in the Great Basin Desert, which had the highest OTU level diversity, a mere 8% of the phylotypes found therein were unique to that desert. Our findings at the OTU level provide support for a core desert microbiome common to the Southwest desert region, an idea previously proposed for Sonoran Desert soils (Andrew et al.2012). This suggests that deserts exert a strong selection pressure that resulted in this set of core desert phyla that are adapted to extreme desert environments. In our study, this core set spanned all phyla. Because differences in Aridisol prokaryotic community composition (Fig. S2) arose mainly from differences in relative abundances rather than differences in taxonomic identity, our findings further support the notion of ubiquitous dispersal within this region (van der Gast 2015). One mechanism explaining the large taxonomic overlap is the dispersal of bacteria on dust particles throughout desert ecosystems. Deserts are a large source of these airborne particles that accommodate transport for bacteria over large distances (Prospero et al.2005; Barberán et al.2015). Alphaproteobacteria, which comprised 76% of the Proteobacteria in our study, can be abundant on dust particles in the continental United States (Barberán et al.2015). Additionally, it is possible that desert bacteria are dispersed across deserts via rain clouds (Delort et al.2010). For example, species from the two most abundant phyla in our study (Actinobacteria and Proteobacteria) have been cultivated from cloud water (Amato et al.2005, 2007), suggesting that monsoon rains might be an underappreciated dispersal mechanism in the southwestern United States.

Drivers of Aridisol prokaryotic community composition

A growing body of evidence indicates that microbial communities follow a distance-decay pattern (Horner-Devine et al.2004; Martiny et al.2011; Monroy et al.2012), whereby communities closer in space are more similar in composition than communities that are farther apart. This is observed in many other ecological systems (Nekola and White 1999). We found support for the distance-decay pattern in our Aridisol study (hypothesis 2), but this pattern disappeared once other variables, namely climate and soil properties, were factored out with partial-Mantel tests and SEMs. Our results are consistent with other studies that show community similarity decay with geographic distance was weak (Martiny et al.2011) or non-existent when other environmental factors (e.g. climate, soil properties) were accounted for (Horner-Devine et al.2004; Fierer et al.2007), demonstrating that an apparent relationship between geographic distance and community similarity is merely indirect. Because there was no residual effect of geographic distance once climate and soil properties were partitioned out, we suggest our study captured the main factors driving the divergence of prokaryotic communities across geographic distance.

Among the factors we measured, climate and soil properties best explained Aridisol microbial community differentiation, with climate having a much stronger influence on microbial communities than soil properties (supporting hypothesis 3a). The SEM allows causal interpretations, but does not preclude other interactions. Nevertheless, consistent with our findings, climate does tend to be an important factor controlling microbial communities at regional and continental scales, as well as other factors such as topography and soil pH (Ettema and Wardle 2002; Fierer and Jackson 2006; Lauber et al.2008, 2009; Fierer et al.2009); vegetation type, land use, soil nutrient status and the quality and quantity of SOM can also be important (Högberg, Högberg and Myrold 2007; Jangid et al.2008; Lauber et al.2008; Fierer et al.2009). We found no indication of vegetation influencing the prokaryotic communities, perhaps because the vegetation data were too coarse at the spatial scale and lacked sufficient species detail to effectively link vegetation to prokaryotic composition. Of the climate variables tested in our SEM, AHM, the ratio between MAT and MAP, had the strongest influence on prokaryotic community composition (supporting hypothesis 3b). The higher AHM of the Mojave Desert likely lead to extremely high daily temperature fluctuations, as dry soils heat up and cool down much faster and to a greater degree than wet soils (Nobel and Geller 1987). High daily soil temperature fluctuations have been shown to be an important stressor for microbial communities (Maestre et al.2015; van Gestel et al.2016) and could therefore exert a strong selective pressure on soil microbial communities. Of the soil properties, organic matter content was the most influential for microbial community similarity. The importance of SOM in arid ecosystems is consistent with a field study of global drylands that suggests microbial communities are limited by carbon (Maestre et al.2015). Soil pH had the weakest effect on microbial community composition in our study, which contrasts with many studies in terrestrial (Fierer and Jackson 2006; Jangid et al.2008; Lauber et al.2008, 2009; Fierer et al.2009; Rousk et al.2010; Rousk, Brookes and Bååth 2010) and aquatic (Fierer et al.2007) ecosystems. However, we only observed small variation in pH across our soil samples (Table 1), which is common for deserts (Maestre et al.2015).

Conclusion

Results from our study suggest that climate plays a larger role than do soil properties in structuring microbial communities at regional and continental scales. However, while past studies in temperate ecosystems have found that precipitation and temperature are the major climatic variables driving soil microbial communities (Castro et al.2010; de Vries 2012), our study showed the combination of these two (i.e. AHM) is the most important climate factor in Aridisols of the Southwest. The extremes in temperature and precipitation, and not merely the means, are likely more important to Aridisol microbial communities than to microbial communities in mesic ecosystems. Our results suggest that prokaryotic communities in the Desert Southwest could thus be disproportionately more susceptible to climate change, as this region is predicted to become hotter and drier, with more intense but less frequent precipitation events (Bernstein et al.2007). Soil microbes influence element cycling, plant community dynamics and the evolutionary responses of ecosystems to global change (Bardgett and van der Putten 2014). Accordingly, loss of soil biodiversity and simplification of community structure are expected to result in loss of ecosystem function and stability (Wagg et al.2014). Because the ecological strategies of microbes related to the moisture regime are likely linked to traits that influence their functional potential (Evans and Wallenstein 2014), microbial communities in the arid Southwest could contribute to shifts in ecosystem function as wetting and drying cycles fluctuate with climate change. This could be particularly true for ecosystem processes that are governed by a narrow range of phylogenetic groups (Graham et al.2016).

SUPPLEMENTARY DATA

Supplementary data are available at FEMSEC online.

Acknowledgements

We thank Matthew A. Bowker for his help and guidance with the structural equation model.

FUNDING

This work was supported by a National Science Foundation Grant, No. EF-0747397. TAM was supported by a National Science Foundation Postdoctoral Research Fellowship in Biology under Grant No. 1402451. NvG was supported by the Department of Energy's Biological and Environmental program, Office of Science, Grant No. DE-SC0010632.

Conflicts of interest. None declared.