Spatial heterogeneity of physicochemical properties explains differences in microbial composition in arid soils from Cuatro Cienegas, Mexico

Study area

The study site is locally known as “Churince system.” It is located in the western part of the CCB (26°50′N, 102°08′W; Fig. 1) at 740 m a.s.l. The system consists of a spring, an intermediate lagoon, and a dry desiccation lagoon connected by short shallow creeks. The annual precipitation in the area is less than 250 mm, occurring mainly from May to October. Temperatures fluctuate from 0 °C in January to 45 °C in July, with a mean annual temperature of 21.4 °C (CCB weather station). The vegetation is mainly halophile and gypsophile grasslands (Challenger, 1998). The area is also dominated by physical and biological soil crusts. The soil is predominantly basic, rich in calcium and sulfates, but very poor in nutrients, and belongs to Gypsisol type (IUSS Working Group WRB, 2015).

Sampling design

The sampling plot was approximately at 50 m from the dry desiccation lagoon. The gypsophile grass Distichlis spicata covered 10% of the plot (Fig. 1) and it was the only species present in this area. Physical and biological crusts occupied the open areas between plants. In order to explore the spatial relationship between soil chemical distribution and microbial community structure at local scale, we used a plot of 8 m × 8 m that consisted of a nested system of four quadrats (A–D quadrats) of 16 m², which were divided in eight 1 m² “replicates,” following a checkerboard pattern (Fig. 1). Vegetation cover for each sampling “replicate” or site (1 m²) was registered qualitatively in order to have further ecological context for the results. In August 2007, we collected soil samples (500 g) from the first 10 cm at each site, to a total of 32 samples (eight samples for each 4 m² quadrat), under the SEMARNAT collection permits 06590/06 and 06855/07. Soil samples were homogenized in the field and divided in two subsamples, which were stored at −20 °C (for molecular analyses) and at 4 °C (for chemical analyses), respectively. Analyses were performed upon arrival to the laboratory.

Physicochemical and biochemical analyses

Soil samples were air dried and sieved through a 2 mm mesh prior to physicochemical and biochemical determinations, which were performed twice for each sample. Total carbon (TC) was determined by dry combustion and coulometric detection (Huffman, 1977) using a Total Carbon Analyzer (TOC, UIC Mod. CM5012; UIC, Inc., Chicago, USA). For total nitrogen (TN) and phosphorus (TP), the samples were acid digested and determined colorimetrically using a Bran-Luebbe Auto-analyzer III (Germany), according to Bremner (Bremmer & Mulvaney, 1982) and Murphy & Riley (Murphy & Riley, 1962), respectively. Inorganic N forms (NH${}_4^{+}$ and NO${}_3^{-}$) were extracted with 2 M KCl after shaking for 30 min, followed by filtration through a Whatman #1 filter, and measured colorimetrically by the phenol–hypochlorite method. Inorganic P (Pi) was extracted with sodium bicarbonate, and determined colorimetrically by the molybdate-ascorbic acid method (Murphy & Riley, 1962). Dissolved organic C (DOC), N (DON) and P (DOP) were extracted with deionized water after shaking for 1 h and then filtered through a Whatman #42 filter. DOC was determined with a TOC module for liquids (UIC-Coulometrics), while DON and DOP were acid digested and measured colorimetrically.

Electrical conductivity and pH were determined in soil with deionized water (soil solution ratio 1:2). To quantify water-soluble cations (Ca²⁺, Mg²⁺, K⁺, Na⁺) and anions (HCO${}_3^{-}$, Cl⁻, SO${}_4^{2-}$), soil samples were shaken with deionized water for 19 h, centrifuged at 2,500 rpm and filtered through a Whatman #42 filter. Ca²⁺ and Mg²⁺ were analyzed by atomic absorption spectrophotometry with an air/acetylene flame (Varian SpectrAA 110), while Na⁺ and K⁺ by flamometry (flame photometer Corning PFP7) (Bower, Reitemeier & Fireman, 1972). Anions were determined by liquid chromatography (Waters Mod. 1525) with a mobile phase of borate sodium glucanate (Bower, Reitemeier & Fireman, 1972).

Molecular analyses

Microbial community structure was characterized using terminal restriction fragment length polymorphisms (T-RFLPs) of bacterial 16S rRNA gene.

Genomic DNA was extracted from the soil samples using the Soil Master DNA Extraction Kit (Epicentre Biotechnology), with an additional previous step based on the fractionation centrifugation technique in order to reduce the high salt concentration (Holben et al., 1988). After extraction, genomic DNA was cleaned with Microcon columns (Fisher Scientific) with the purpose of removing any substance that could inhibit PCR amplification. These protocol modifications gave the best results from various methodologies tested; however, we were only able to amplify the 16S rRNA gene from 21 out of the 32 soil sampling sites (amplicons obtained in each quadrat: A = 3; B = 3; C = 7; D = 8). The low yield in the DNA amplification could be attributed to molecular applications inhibitors of unknown nature (López-Lozano et al., 2012).

Amplification of the bacterial 16S rRNA genes was carried out in a final volume of 50 µL containing: 0.2 µM of each fluorescently labeled domain-specific primers (VIC-27F and FAM-1492R) (Lane, 1991), 0.2 mM of each dNTP, 1 U of Taq Platinum DNA polymerase (Invitrogen), 2.5 µL DMSO, 2.5 µL BSA, 1 mM MgCl₂, 1 mM buffer, and 20 ng of DNA. Five independent PCR reactions were performed for each sample with the following program: 5 min at 94 °C; 30 cycles at 94 °C for 1 min, 52 °C for 2 min, 72 °C for 3 min; and 72 °C for 10 min. PCR products were pooled and purified from 2% agarose gel (Gel extraction kit, Qiagen Inc.). The amplicons were restricted with AluI enzyme (Promega) at 37 °C for 3 h and 65 °C for 20 min. Three independent readings of the size and abundance of fluorescently labeled terminal restriction fragments (TRFs) were performed for each sample using an ABI 3100-Avant Prism Genetic Analyzer (Applied Biosystems), as described previously (Coolen et al., 2005). For each profile of TRFs, we established a baseline and only those TRFs with peak heights ≥50 fluorescent units were used in subsequent analyses (Blackwood et al., 2003). Each unique TRF was considered to be an operational taxonomic unit (OTU). Estimations of diversity were derived from matrices constructed based on the presence and abundance of TRFs using relative peak area as an estimate of abundance calculated as: ${\displaystyle \mathrm{Ap}=\left(n_i∕N\right)\times 100}$ in which n_i represents the peak area of one distinct TRF and N is the sum of all peak areas in a given T-RFLP pattern (Lukow, Dunfield & Liesack, 2000).

Statistical analyses

Statistical and diversity analyses were performed in R (R Development Core Team, 2011), mainly with vegan (Oksanen et al., 2012), ggplots (Warnes, 2012) and BiodiversityR (Kindt & Coe, 2005) packages.

All soil properties data were expressed on a dry-weight basis. Non-normal data were log-transformed to normalize the distribution of the residuals. We analyzed the data using both multivariate (MANOVA to detect patterns: whether there were significant effects of the four quadrats on overall soil variables) and univariate (to detect significant differences in individual variables) methods. These analyses were followed by multiple pairwise tests, using Tukey’s honestly significant difference (HSD), at the 5% level of significance, to identify possible differences in the soil variables between quadrats. The correlations between each pair of variables were calculated using Pearson’s correlation coefficient. Soil properties were then standardized and ordered by principal components analysis (PCA), and the sampling points from the four quadrats were visualized with the two first principal components.

Alpha diversity indices (Shannon, Simpson, and Berger-Parker) and richness estimates were calculated for each quadrat using the T-RFLPs profiles. Microbial diversity indices were also analyzed using ANOVA type III for unbalanced data and evaluated using Renyi’s entropy profiles for eight scales (α = 0, 0.25, 0.5, 1, 2, 4, 8, infinite) (Rényi, 1961; Chao et al., 2014) with the BiodiversityR package. These profiles provide a comprehensive analysis of the diversity, giving a parametric measure of the uncertainty of predicting the OTUs richness, as well as the relative abundance of OTUs, at different scales between the four quadrats. To evaluate the sampling effort, rarefaction curves were constructed for each quadrat using EstimateS v.9.1.0 (Colwell, 2005). Microbial community structure between quadrats was examined through Venn diagrams (with the matrix of OTUs presence) and ordination analyses (with the matrix of OTUs abundance). To visualize communities’ structure, Bray-Curtis dissimilarity distances were calculated with the relative abundance of T-RFLPs profiles. Similar communities were then clustered using the Ward’s hierarchical clustering algorithm, which tries to minimize variances in agglomeration. A heatmap of the relative abundance of OTUs was constructed with dual hierarchical clustering.

Community structure was also investigated for correlations with chemical parameters following a multivariate analysis. For this, the relative abundance of T-RFLPs profiles were ordered by Detrended Correspondence Analysis (DCA) with Hellinger transformation (Blackwood et al., 2003), and correlations between the ordination axes and soil properties were calculated. This eigenvector-based ordination technique uses a chi-square distance measure and assumes that TRFs have a unimodal distribution along ecological gradients (Legendre & Legendre, 1998), which is a more appropriate assumption than linearity for ecological analysis of T-RFLPs data (Culman et al., 2008). Permutation tests under reduced model were used to identify significant explanatory soil variables (p < 0.05), which were plotted onto the ordination map as vectors. Only the soil variables corresponding to the same sampling sites as the T-RFLPs data were used for this analysis.

To further evaluate which soil abiotic properties could be used as predictors of the distribution of OTUs diversity, a multiple linear regression analysis was performed. Significant explanatory variables were chosen as the best predictors of OTU abundance by stepwise multiple regressions and by minimizing the Akaike Information Criterion (AIC). Collinearity between explaining variables was also tested using the Variance Inflation Factor (VIF).

Spatial heterogeneity in physicochemical and biochemical parameters

The total plant cover in the experimental plot was only 10%. However, quadrats C and D were more densely and homogeneously covered than A and B (Fig. 1). Overall there was a high presence of soil crusts and biocrusts, particularly in the A quadrat.

Soil chemical properties varied significantly between the four quadrats (Wilks’ lambda = 0.000, F = 7.896, p < 0.05). Soil samples were alkaline (pH between 8.6–8.8) due to the high presence of salt in this arid ecosystem (Table 1). Cations (Ca²⁺, K⁺) and anions (HCO${}_3^{-}$, Cl⁻, SO${}_4^{2-}$) were the most variable parameters in this small plot showing significant differences between quadrats, which means that ions significantly contribute to soil heterogeneity in this system. C and D quadrats had the greatest concentration of cations (0.64 and 0.65 cmol Ca²⁺ kg⁻¹, 1.27 and 1.33 cmol K⁺ kg⁻¹, respectively), except for Mg²⁺ and Na⁺, while A and B quadrats had the highest concentration of anions (2.8 and 2.3 cmol HCO${}_3^{-}$ kg⁻¹, 2.8 and 2.5 cmol Cl⁻ kg⁻¹, 15.6 and 16.4 cmol SO${}_4^{2-}$ kg⁻¹, respectively). The high concentration of Na⁺ found in these soils (mean value of 147 cmol kg⁻¹) indicates salinity stress. Total forms and nutrients content were very low in this arid soil (TC: 2.4–2.8 mg g⁻¹; TN: 0.48–0.6 mg g⁻¹; TP: 0.03–0.04 mg g⁻¹; NH${}_4^{+}$: 3.6–4.2 µg g⁻¹; NO${}_3^{-}$: 1.5–2.3 µg g⁻¹), as expected, and they did not show significant differences among quadrats. The total C/N ratio (a quality index for soil organic matter) was also very low (from 4.6 to 6.1) and did not show significant differences between the four quadrats. On the other hand, C availability (DOC) was higher in D quadrat (124 µg g⁻¹), which means greater substrate availability for microbial metabolism in this quadrat. As expected, the pH was positively correlated with Mg²⁺ and Na⁺ (Table S1). The TC was only positively correlated with Ca²⁺, while TP was positively correlated with pH and cations, as well as negatively with DOP. The TN was positively correlated with DON and negatively with C:N, NH${}_4^{+}$, NO${}_3^{-}$ and HCO${}_3^{-}$. Finally, N inorganic forms were also positively correlated between them and the C:N ratio.

The complex chemical spatial heterogeneity among these four quadrats was explored using a PCA (Fig. 2). The first component (PC1) explained 54.3%, while the second component (PC2) explained 34.3% of the total variation in the soil parameters among quadrats (Table S2). The variables associated with the PC1 were cations and anions, as well as pH, TP, DOC and DOP. The PC2 was mainly related to soil nutrients (TN, C:N, NH${}_4^{+}$, NO${}_3^{-}$, DON, DOC:DON). A clear separation between quadrats was observed along the PC1 axis, mainly explained by the spatial heterogeneity distribution of ions.

Spatial heterogeneity of microbial diversity

A total of 184 different OTUs were obtained in the four quadrats, of which 121 OTUs had less than 1% of the total maximum relative abundance. Unfortunately, the number of available samples was unbalanced for the microbial diversity study in this plot (amplicons in each quadrat: A = 3; B = 3; C = 7; D = 8) potentially due to the presence of inhibitors of unknown nature, which hampered the 16S rRNA amplification from all samples sites. Despite this constraint, rarefaction curves showed a good community sampling for quadrats A, C and D, with evident subsampling for quadrat B, which is one of the two quadrats for which only three out of eight samples could be analyzed in terms of microbial diversity (Fig. S1). Significant variation in alpha diversity indices among quadrats was detected. Although the A quadrat was also limited in the number of analyzed samples (3), it was the most diverse (H: 3.31) and with the highest evenness (1∕D: 0.944; BP: 0.153), followed by C, D, and B (Table 2). It was also evident the high variability in microbial diversity among replicates (with the exception of A), which reflects the spatial heterogeneity at small local scale of this arid soil. A summary of diversity indices was obtained by calculating Renyi’s community profiles. Ranking based on these profiles is preferred to ranking based on single indices because rank order may change when different indices are used (Kindt & Coe, 2005). These profiles showed the same pattern of diversity, both in terms of richness and evenness, where the highest diversity was found for quadrat A and the lowest for quadrat B, while C and D presented intermediate diversity (Fig. 3).

Despite the high heterogeneity in microbial diversity in such a small plot, Venn diagram revealed a considerable overlap of OTUs among the four quadrats: 18% of OTUs were shared by all quadrats (Fig. 4), which represent 44.2% of the total abundance of the microbial community recovered from this plot with the T-RFLP technique. The C quadrat had the most “unique” OTUs (12%, representing 0.5% of total abundance), followed by quadrat D (11.4% representing 1.7% of total abundance), B (10.3%, representing 3.9% of total abundance) and A (5.4%, representing 0.3% of total abundance) in decreasing order. Interestingly, C and D quadrats shared 84 of the 184 recovered OTUs (46%), suggesting that both quadrats had more similar community composition than A and B quadrats. Moreover, the heatmap (Fig. S2) and the cluster dendrogram (Fig. S3) of OTUs abundance showed the same pattern of grouping as the PCA analysis for soil chemical properties, separating quadrats in two groups: A–B and C–D. As in the PCA analysis, there was a sample in the C quadrat that clearly deviated from the other samples.

Multivariate analyses of microbial community structure

To explore the association between community structure and soil parameters, we performed a DCA analysis. The original 19 soil parameters were reduced to seven non-redundant explanatory variables (TN, DON, Ca²⁺, K⁺, HCO${}_3^{-}$, Cl⁻, and SO${}_4^{2-}$; Table S3), which were the factors that contributed significantly to differences in community composition among the four quadrats in the ordination analysis. The DCA showed a clear separation of the quadrats in two groups, mainly explained by soil salinity: anions (HCO${}_3^{-}$, Cl⁻, and SO${}_4^{2-}$) significantly correlated with OTUs from the A and B quadrats, while TN, DON, Ca²⁺ and K⁺ significantly correlated with OTUs from the C and D quadrats (Fig. 5). Thus, the grouping pattern of these microbial communities showed in the above analyses was also confirmed by the DCA analysis.

To further investigate the relation among OTUs diversity distribution and soil abiotic properties, a full linear model was then tested and reduced for a model where all the coefficients were statistically significant and the best AIC achieved. The best regression model indicated that DON, Ca²⁺, K⁺, and SO${}_4^{2-}$ were statistically significant to explain the distribution of OTUs diversity in this soil (Table S4; R² = 0.59, p < 0.01).

Chemical heterogeneity at local spatial scale is mainly due to ions concentration variability

The values of TC, TN and TP in the experimental plot were lower than the values reported for other deserts (Thompson et al., 2006; Strauss, Day & Garcia-Pichel, 2012), as well as for soils in the CCB (López-Lozano et al., 2012; Tapia-Torres et al., 2015). The very low C:N ratio also suggests deficiency of soil organic C, therefore a low nutrient availability to soil microbes and vegetation, limiting the N cycle due to the lack of C availability. The Redfield ratio in this soil was 71:17:1, which suggests that in these quadrats the C is the limiting nutrient in comparison with a general “average” soil C:N:P of 186:13:1 (Cleveland & Liptzin, 2007). On the other hand, this result also differs from the Redfield ratio of 104:5:1 reported for the same soil system (López-Lozano et al., 2012). These differences could be attributable to the great heterogeneity of this arid environment and the different time of soil sampling in both studies: February 2007 (dry cold season with low evapotranspiration) in López-Lozano et al. (2012) and August 2007 (rainy hot season with high evapotranspiration) in this study.

All soil samples in this study had high alkalinity produced by the elevated concentrations of ions, which is a general pattern in desert soils (Titus, Nowak & Smith, 2002). The high pH decreases P availability, which is very scarce in these soils and it is bonded to Ca²⁺ and Mg²⁺ (Cross & Schlesinger, 2001; Perroni et al., 2014). It is worth to mention that ions varied spatially in identity in this small plot: quadrats A and B were significantly high in anions, while quadrats C and D were significantly high in cations. The huge concentration of Na⁺ in the four quadrats is an indicator of the extremely high salinity in these soils, which negatively affects the soil aggregates stability, as well as nutrients and water availability for plants, favoring the development of soil crusts, which are typical in arid and semiarid soils (Belnap, 2003; Zhang et al., 2007). In particular, salt crusts are abundant in the CCB area. They consist of layers at the soil surface mainly formed by soluble salt crystallizing soil particles at shallow saline groundwater level regions (Zhang et al., 2013).

The high concentration of ions can be attributed to the gypsum-rich nature of the CCB soils, where groundwater rises to the surface by soil capillarity action and water evaporation promotes salt accumulation. This situation results in rivers with a steep salinity gradient (Cerritos et al., 2011) and pools surrounded by saline soils rich in sulfates and extremely poor in nutrients (López-Lozano et al., 2012). Therefore, it is not surprising to find that the soil properties variation in this small plot was mainly explained by ions concentration, grouping the four quadrats in two broad clusters: A–B and C–D. These clusters had a qualitative pattern associated with the vegetation cover, being quadrats C and D more densely and homogeneously covered by vegetation than A and B. Although the present research analyzes soil communities, a previous study of microbial communities of the water system associated with the studied plot showed a correlation of microbial composition and water conductivity gradients (Cerritos et al., 2011). Thus, the spatial variation in these physicochemical properties among the four quadrats may be a consequence of differences in moisture content due to the proximity to a subterranean water flow, indirectly evidenced by the marked patchy distribution of the vegetation cover and the “open” areas occupied by soil crusts (López-Lozano et al., 2012).

Heterogeneity in microbial diversity at local spatial scale is explained by physicochemical factors, not by vegetation cover nor nutrient content

Despite recent important advances in our knowledge of the structure, composition and physiology of biotic components in arid soils (Belnap et al., 2005; Caruso et al., 2011; Maestre et al., 2015; Makhalanyane et al., 2015), little is known about the spatial variability of microbial diversity at local scales and its interactions with chemical heterogeneity in these ecosystems (Housman et al., 2007; Castillo-Monroy et al., 2011; Andrew et al., 2012). T-RFLPs fingerprinting was used in this study to assess the relationship between microbial structure and the small-spatial heterogeneity of soil chemical properties. We are aware that this technique cannot recognize taxonomic groups and accounts mainly for relatively abundant microbial groups. Nevertheless, given the aims of this study of characterizing microbial communities structure and its relationship with abiotic or physichochemical parameters, a fingerprint approach such as T-RFLPs is adequate to provide replicable, valid and sufficient data (Angel et al., 2013), as it has done for many other studies looking at patterns of correlation between microbial diversity/composition and environmental factors (Fierer & Jackson, 2006).

The heterogeneity in OTUs diversity among these quadrats is evident, being the A quadrat the most different with respect to the other quadrats. Despite the fact that the A quadrat had scarce plant cover and similar nutrients and ions concentrations to the B quadrat, it showed the greatest microbial diversity, which could be related to the high presence of biocrusts that may incorporate more resources to the soil, potentially increasing organic C and biomass and, in turn, diversity (Geyer et al., 2013). On the other hand, the B quadrat had the lowest microbial diversity, which could be related to the lowest values of DOC found in this quadrat. Labile organic matter fractions, such as DOC, are the primary energy source for soil microorganisms and are characterized by rapid turnover (Bolan et al., 2011). It has been reported that even in disturbed sites, DOC is the main source of C influencing the composition of the microbial community (Churchland, Grayston & Bengtson, 2013). Then, changes of soil microbial community could be regulated by C availability through labile soil organic matter pools, as have recently been shown to happen in McMurdo Dry Valleys arid soils in Antarctica (Geyer et al., 2013).

Regarding similarity in microbial composition among the four quadrats, cluster dendrogram and multivariate analyses showed two clear groups, which were A–B and C–D, corresponding to the same clustering of quadrats based on soil chemical parameters. A common explanation for the soil microbial composition patterns is related to the presence of plants controlling levels of microbial diversity and driving community assembly (Singh et al., 2007; Berg & Smalla, 2009; Ben-David et al., 2011). However, in our study the observed spatial pattern of microbial diversity distribution at such local scale does not seem to be associated with vegetation cover. For example, the A quadrat is the most diverse in microbial community and the less vegetated, suggesting that microbial diversity in this arid soil could be more related to the presence of “open” areas occupied by biocrusts. On the other hand, abiotic factors, such as ionic content, are statistically explanatory variables in the spatial ordering (DCA) of the microbial communities analyzed. In addition, the multiple regression model selected shows the importance of DON, Ca²⁺, K⁺, and SO${}_4^{2-}$ in explaining the distribution of OTUs abundance in this arid soil. Abiotic drivers of microbial diversity in arid soils has been also reported for the Sonoran desert (Andrew et al., 2012), where location, pH, cation exchange capacity and soil organic C were highly correlated with microbial composition. Therefore, we showed that microbial community diversity and distribution responds to and/or influences local soil physicochemical characteristics at a small spatial scale in this arid ecosystem.