Microbial secondary succession in soil microcosms of a desert oasis in the Cuatro Cienegas Basin, Mexico

Ecological succession is one of the most important concepts in ecology. However for microbial community succession, there is a lack of a solid theoretical framework regarding succession in microorganisms. This is in part due to microbial community complexity and plasticity but also because little is known about temporal patterns of microbial community shifts in different kinds of ecosystems, including arid soils. The Cuatro Cienegas Basin (CCB) in Coahuila, Mexico, is an arid zone with high diversity and endemisms that has recently been threatened by aquifer overexploitation. The gypsum-based soil system of the CCB is one of the most oligotrophic places in the world. We undertook a comparative 16S rRNA 454 pyrosequencing study to evaluate microbial community succession and recovery over a year after disturbance at two sites. Results were related to concurrent measurements of humidity, organic matter and total C and N content. While each site differed in both biogeochemistry and biodiversity, both present similar pattern of change at the beginning of the succession that diverged in later stages. After one year, experimentally disturbed soil was not similar to established and undisturbed adjacent soil communities indicating recovery and succession in disturbed soils is a long process.


Introduction
In ecological theory, succession is defined as the predictable manner by which communities change over time during the colonization of a new environment or following a disturbance (primary and secondary succession respectively) (Begon et al. 2006).
Bacterial succession has been explored in a variety of environments and over different timescales (Fierer et al. 2010). However, there are many theoretical and methodological questions that remain unansweredsolved. For example the order of occurrence and species turnover inaccording to the successional stages can differ and be explained by stochastic or functional factors depending on the use of ecological or functional classifications instead of 16S rRNA taxonomic classifications (Burke et al. 2011a;Burke et al. 2011b). Notwithstanding, there are theoretical models for microbial community succession and assembly, most of them derived from those already well defined for macroorganisms (Prosser et al. 2007). Due to the complexity of the patterns observed in the succession of microbial communities, alternative explanations started to appear based on carbon (C) inputs to the systems, suggesting that species richness and the biomass of specific ecological groups change along successional stages (Fierer et al. 2010).
In recent years, the role of deterministic and stochastic processes in community assembly has been tested. It has been shown that bacterial community assemblages can be regulated by the local environment (species sorting) (Langenheder & Székely 2011), massive immigration that can prevent competitive exclusion of species (mass effect) (Lindström & Langenheder 2012), and neutral process, which assumes that all species are similar in their competitive ability and in dispersal. In the neutral process scenario, this scenario changes in community composition result from neutral drift over PrePrints PrePrints time. (Woodcock et al. 2007;Ofiteru et al. 2010). These mechanisms can co-occur, resulting in microbial communities being structured by more than one process (Langenheder & Székely 2011). However, t The development of such theoretical models are based mostly on work done in laboratory environments, and little is known of what happens in natural systems such as soils (Caruso et al. 2011;Navarro et al. 2009;Nemergut et al. 2007;Schmidt et al. 2007). Data on how microbial communities change through time in natural environments are needed in order to refine these theoretical models.
Soils are among the most diverse microbial environments analyzed to date (Youssef & Elshahed 2009;Daniel 2005), which makes the identification of differences in community diversity patterns between stable and disturbed soils challenging. Arid soils are thought to have lower productivity and diversity than other soil habitats, and therefore may provide a better opportunity to validate molecular ecosystems approaches to examine the genetic and functional organization of native microbial consortia. Earth's arid regions represent today nearly one third of total continental ecosystems (Collins et al. 2008) and are thought to be more vulnerable under most scenarios of global climate change (Intergovernmental Panel on Climate Change, 2007 report).
The Cuatro Cienegas Basin (CCB) in Coahuila, Mexico, is an arid zone that has recently been threatened by aquifer overexploitation. The gypsum-based soil system at CCB is one of the most oligotrophic environments in the world despite the high diversity in comparison with others arid soils (López-Lozano et al. 2012). Nevertheless, due to recent and ongoing overexploitation of the deep aquifer in this oasis, we decided to evaluate the community sensitivity after disturbance. In this study we experimentally In January of 2007, 20 kg of soil were collected from each site. The soil was mixed and used to construct microcosms of 1 kg of soil in permeable bags of nylon mesh (25 x 25 cm). The microcosms were sterilized by autoclaving followed by a dose of 25 kGy of gamma rays (at the Instituto Nacional de Investigaciones Nucleares, La Marquesa, México). One microcosm sample from each site was immediately used for DNA extraction as a control.
In February 2007, plots of 8 x 8 m were established at each site. The plots were divided in 64 quadrants of 1 m 2 , and 40 microcosm bags were placed randomly within the plots as replicates. Every three months, three microcosm bags were collected randomly during a period of one year (3, 6, 9 and 12 months), 12 in total for each site. The remaining microcosms were deposited in the site in order to guarantee enough samples in each sampling period but not analyzed for this experiment. In addition, three samples from 15 cm 3 of undisturbed soil next to the microcosms in each plot were randomly collected at the beginning and end (12 months) of the experiment to describe the community that potentially could colonize the microcosms. 50 g of each sample were stored in liquid nitrogen until DNA extraction. The remaining soil was stored in black plastic bags at 4ºC during one month until processing in the laboratory for biogeochemical analysis. using a Bran-Luebbe Auto Analyzer III (Norderstedt, Germany). Total and inorganic C were determined by dry combustion and coulometric detection (Huffman 1977). Organic C was calculated as the difference between total and inorganic C. N and P concentrations were determined following acid digestion; N was determined using a modified Kjeldahl method (Bremmer & Mulvaney 1982), and P was determined with the molybdate colorimetric method following ascorbic acid reduction (Murphy & Riley, 1962).
Microbial C (micC) and microbial N (micN) concentrations were determined from field-moist samples by a chloroform fumigation extraction method (Vance et al. 1987).
Inorganic N (NH4+ and NO3-) was extracted from fresh sub-samples with 2 M KCl, followed by filtration through a Whatman No. 1 paper filter (Robertson et al. 1999), and determined colorimetrically by the phenol-hypochlorite method. Dissolved organic C, N and P were extracted in two grams of fresh material with deionized water after shaking for 1 h, and filtering through a Whatman# 42 filter and a 0.45µm nitrocellulose membrane (Jones y Willett, 2006). Dissolved organic C (DOC) was determined by combustion and coulometric detection (Huffman 1977). Dissolved organic N (DON) and dissolved organic P (DOP) were determined after acid digestion. DON was calculated as the difference between digested soluble N and NH4+ in deionized water extracts. The DOP was calculated as the difference between digested dissolved P and inorganic P (as orthophosphate).

PrePrints PrePrints
(RMANOVA) with one between-subject factor (site: Dry Lagoon and River) and one within factor (sampling date: 3, 6, 9 and 12 months). In order to compare biogeochemical soil parameters between sites, and between disturbed and undisturbed soil in the last sampling date (12 months), we carried out a two-way analysis of variance (ANOVA) (factor 1 levels: sites Dry Lagoon and River; factor 2 levels: undisturbed soil and microcosm).
The relationship between microbial community composition in terms of the most abundant bacterial families, soil characteristics and samples (date and site) was analyzed by canonical correspondence analysis (CCA). In this analysis, tThe samples are represented by a centroid. Its position is indicative of the relationship between a specific sample and either of the ordination axes. Soil characteristics are represented by vectors. Vectors of greater magnitude and forming smaller angles with an ordination axis are more strongly correlated with that ordination axis. High scores of absolute value for a given family or a given site on a CCA axis indicate that it is highly related to the axis and to the environmental variable exhibiting high correlation to the axis. All soil characteristics were tested for significant contribution to the explanation of the variation in bacterial family community composition with an ANOVA like permutation test to assess the significance of constraints. Only variables that were significant by the permutation test at the P≤0.05 level were included in the CCA biplot. The CCA was performed using the package Vegan in R (http://www.r-project.org/). In the supplementary material, we added a DCA analysis without the restrictive variables for comparison ( figure S2).

PrePrints PrePrints
(EPICENTRE Biotechnology) according to the manufacturer's instructions, with an additional step of bacterial isolation using the fractionation centrifugation technique described in (Holben et al. 1988). This step was performed on frozen 50 g soil samples before DNA extraction as a way of reducing the remaining concentrations of salts, polysaccharides and secondary compounds. DNA was stored at -20 °C.
In order to confirm sterilization in the control microcosms, the DNA region coding for 16S rRNA was PCR amplified using universal primers (Tables S1 and S2). Since no band was recovered in a PCR gel run using the sterilized soil, the sterilization was assumed to be effective.
Pyrosequencing of 16S rRNA tags. 16S rRNA genes were amplified from each sampling date (pooled DNA from the three subsamples was used as template) using the 939F (TTGACGGGGGCCCGCAC) and 1492R (TACCTTGTTACGACTT) paired primers.
Sequencing was undertaken using the standard Roche 454 Titanium LIB-A kit with multiplex identifier sequence (MID) tags (Table S3)  Lagoon are: before sterilization SRS346170, at 3 months SRS346171, at 6 months SRS346172, at 9 months SRS346173 and at 12 months SRS346174. The samples ID for microcosms in the River are: before sterilization SRS346176, at 3 months SRS346177, at 6 months SRS346178, at 9 months SRS346179 and at 12 months SRS346180. The samples ID for the undisturbed adjacent soil without disturb at 12 months are SRS346175 for the Dry Lagoon, and SRS346181 for the River.

Results
Nutrient fluctuations during the secondary succession. In order to understand the relationships between the changes in community composition and nutrient transformations that occurred during the secondary succession process in the Churince soils, the soil samples were biogeochemically characterized and compared. We found that the CCB soil at the Dry Lagoon site was more oligotrophic (C:N:P ratio of 125:5:1) than the River site (C:N:P ratio 300:16:1) and in comparison with a general "average" soil C:N:P of 186:13:1 (Cleveland & Liptzin 2007). RMANOVA analysis revealed consistent differences in the soil nutrients between sites across time ( Table 1). The Total Organic C was higher at the site adjacent to the River throughout the sampling, even though it showed a fluctuating pattern during the secondary succession experiment, probably due to both seasonal effects and succession stages ( Table 1) (Table 1).
Values increased during the first three and six months after sterilization, then micC decreased at 9 months and increased again in the last sampling date. Ammonium showed contrasting patterns between sites. At the Dry Lagoon, ammonium increased showing higher values at 6 months, and then decreased gradually. At the River, ammonium increased gradually to drop in the last sampling. The microbial N (micN) increased in samples near the River, but presented fluctuations in the Dry Lagoon. Only the nitrate was higher in the Dry Lagoon microcosms but there were not significant differences by time. ANOVA analysis of soil nutrients in microcosms samples and undisturbed adjacent soil at 12 months showed significant differences between treatments in Total Organic C, dissolved organic N (DON), dissolved organic P (DOP) and ammonium (Table 2). However, undisturbed soil samples were similar to the perturbed samples in their total N and P content (TN, TP and DIP) ( Table 2). These data suggests that, in general, the sampling sites are different and as such, respond differently to perturbation.   PrePrints PrePrints any eukaryotic specie which mitochondria and chloroplast would be amplified by these techniques. Suggesting that while wind and water could help the colonization processes; nematodes and other invertebrates did not play an important role in the migration of bacteria between disturbed and undisturbed sites.
Despite the differences observed in nutrient content between sites, the 16S rRNA sequence analysis showed similar variations in the diversity and community composition of the sites, even though the overall diversity was higher in the River during the entire sampling period (Table 4, Figure 1). At both sites, during the first six months (3 and 6 months samples) there was an increase in diversity, then the diversity decreased by nine months and increased again after one year (12 months) following a similar fluctuating pattern to the biomass. Contrasting the diversity in the last sampling date with the diversity of the adjacent undisturbed soil, we observed that the undisturbed soil had higher diversity than microcosms samples (Table 4).
Prior to disturbance, both Dry Lagoon and River communities were dominated by Acidobacteria (44% and 60%, respectively) ( Table 3). The Dry Lagoon had a large proportion of bacteria that were unclassified at the 97% identity threshold (around 12 to 33 % of the sequences). Proteobacteria were also abundant, comprising 9% of the Dry Lagoon sequences and 35% of the River sequences. Overall there was a gradual increase in the abundance of the Bacteroidetes during the succession process, which were not abundant in the first sampling (less than 0.24%) but increased with time (around 7-9%). Autotrophic phyla such as Clorobi, Chloroflexi and Cyanobacteria, increased in abundances the first nine months and then decrease their abundance in the last sampling (12 months) (Table 3).

PrePrints PrePrints
A more detailed analysis at the family level showed that the families present before the sterilization in low abundance were not found after 3 months, but many of them were recovered by 6 months and onward (Table S4). The initial colonizers at 3 months were different between sites, but in both cases, the majority of the identified families were related to known opportunistic heterotrophs. In the Dry Lagoon the most abundant families were Moracellacea, Streptococcaceae, Alcaligenaceae, Pseudomonadaceae, Micrococcaceae and Bacillaceae. At the River site at 3 months, Holophagaceae, Acidobacteriaceae, GIF3 (Chloroflexi), Enterobacteriaceae, Streptococcaceae, Listeriaceae and Bacillaceae were found. By 6 months the abundance of unclassified taxa increased, but also the abundance of families that include previously described members with other metabolic capabilities besides heterotrophy (Table S4). The relative abundances of some bacterial groups at different taxonomic levels correlated significantly with Total Organic C, Total N, humidity and NO3 - (Table 5), and these trends were corroborated with the CCA (Figure 2).
Similarity between sites and samples through the time, measured with the Bray-Curtis distance and Jaccard coefficient using OTUs at 97% of identity, achieved a peak three months after sterilization, while it was lower before sterilization and six, and nine months after sterilization. After a year the microcosms communities were more similar to the communities from undisturbed soil at each site ( Figure 3).

Discussion
Microbial diversity and secondary succession in a gypsum based soil. In this study we Even with the higher sequencing effort afforded by 454 pyrosequencing, we still only captured 52% to 81% of the microbial communities in our samples (Figure 1 and Table   4). This was surprising, since in more fertile and wet soils analyzed with similar sequencing effort, the maximum number of OTUs at 97% of identity seldom exceeds  Roesch et al. 2007), while in our more diverse sample, we estimated ~ 9,000 OTUs at 97% identity. However, the sites analyzed in this study differed in diversity, with soil adjacent to the River being more diverse than the Dry Lagoon if we compare the OTUs at 97% (Table 4). In the previous clone library study the same pattern was observed (López-Lozano et al. 2012). We also evaluated the community composition in terms of "ecological groups", considering groups in which abundance correlated with high or low C mineralization rates (copiotrophs or oligotrophs, respectively) in general soil surveys (Fierer et al. 2007(Fierer et al. , 2010, and widely known autotrophs (phototrophs such as Cyanobacteria and Chloroflexi). Also we determined some metabolic capabilities by comparison of the sequences with the closer cultivated organisms. Despite the broad ecological classification based only in 16S rRNA gene sequences, our results showed interesting trends in groups considered as oligotrophs, copiotrophs (organism that need more nutrients) as well as heterotrophs and autotrophs.
The 16S rRNA libraries of both sites are dominated by Acidobacteria in all samples. This is one of the most common phyla found in soil libraries worldwide (Janssen 2006). The Acidobacteria are in general oligotrophic (Eichorst et al. 2007;Fierer et al. 2007) and have been shown to comprise 50% of clone libraries in arid soils (Dunbar et al. 1999;Kuske et al. 1997), but they are less abundant in more nutrient-rich agricultural soils (Nagy et al. 2005;Roesch et al. 2007). Hence, it is not surprising that they represent ~30-60% of bacteria in all our libraries in CCB sampling sites. In contrast, Bacteroidetes have been classified as copiotrophs in general (Fierer et al. 2007). In our libraries Bacteroidetes become much more abundant at the 6 months sampling, and it is possible that these phyla appear when the accumulation of organic material is enough to sustain the mineralization rate of this group. In fact, between 6 and 9 months there was a drop PrePrints PrePrints in the Dissolved Organic C while the microbial C increase in the soil microcosms.
Autotrophic groups, such as Cyanobacteria are abundant at both sites. This phylum is common in most environments including soil and biological soil crust of arid zones (Gundlapally & Garcia-Pichel 2006;Nagy et al. 2005). The abundance of yet another autotrophic group, the Chloroflexi, increased during the initial stages of succession (3-9 months) but decreased in the last sampling date (12 months). With the exception of high Chloroflexi in the River sample during the first samplings, the majority of the colonizers are opportunistic heterotrophs (a greater part members of Firmicuttes,

Betaproteobacteria and Gammaproteobacteria). It is possible that especially in the Dry
Lagoon site, the initial colonizers were benefited by the Dissolved Organic C, released during the sterilization process. The nutrient pulse could have facilitated colonization by heterotrophs that benefited from available resources in the disrupted microcosm. The finding of preferential growth of opportunistic heterotrophs during early succession agrees with the findings from another study using outdoor sterile microcosms seeded by rainwater bacteria (Langenheder & Székely 2011). These authors found that neutral and species sorting processes interacted during the assembly of bacterial communities, and the importance of each depended on how many generalists and specialists were present in the community. We suggest that in our microcosms, the initial faster growing community (generalists) depleted the nutrients, and then were out competed by groups that have alternative energy sources (specialists). At 6 months the abundance of groups with more specialized metabolic capabilities increased, and the taxa with low abundances before the sterilization, "the rare biota" that were not present after three months returned to the community from this sampling date. The observed nitrate accumulation in the Dry Lagoon site could be due to an increase in nitrification, an alternative mechanism of obtaining energy that could persist under those conditions PrePrints PrePrints (Montaño et al. 2007). Also, the higher DOP in the Dry Lagoon site suggests that P mineralization is not occurring as efficiently as in the River site, perhaps because of limited energy for exo-enzyme production necessary for acquisition (F. Garcia-Oliva, personal communication). Hence, the Dry Lagoon site might be supporting a more oligotrophic microbiota than the River site. Conversely, changes in DOC at the River site correlated with changes in microbial C. Additionally, ammonium correlated with microbial N, suggesting that the initial surge of nutrients could have facilitated a faster migration of the neighboring heterotrophic biomass in this site than in the Dry Lagoon site.
Photosynthetic autotrophs (such as Chloroflexi) seem to be sufficient for nurturing the community in the River site providing C sources, since nitrates and DOP are scarce, indicating low nitrification rates.
Regarding the diversity shifts, both sites showed similar patterns. There is an initial increase in the Shannon index, followed by a decrease at 9 months and an increase at 12 months. We do not discard the possibility that these convergent patterns only in the diversity index shifts can be related to climatic factors. At nine months a decrease in diversity was observed, which coincided with an unusually heavy precipitation event (Table S5).
The Bray-Curtis similarity index and Jaccard similarity coefficient analyses, suggest that even though the initial communities are similar in composition, they differ in later stages.
In addition, the analysis by phylogenetic and ecological groups, suggest that the communities respond in a similar trajectory of initial colonization first by heterotropic generalists and later specialists. However the River site was characterized by higher nutrients and diversity, together with the continuous presence of Cloroflexi in high PrePrints PrePrints abundance suggest an early food web based on primary production. This pattern is similar to the early stages of the autotrophic succession suggested by Fierer et al.
(2010). Our experimental design cannot be used to test this autotrophic model, but our data show that the colonizers were a mix of heterotrophs and autotrophs.
In this study we found community composition of the soil microcosms after perturbation did not recover to similar communities found in the undisturbed soil community after one year. We are not assuming that there is only one "climax" community in terms of composition. What we suggest is that a community would be "recovered" when it has similar performing conditions to the neighboring undisturbed community. The performance of the community can be inferred by the soil nutrient content and characteristics. Due to the significant differences in physicochemical parameters between microcosms and control sites at 12 months, we can conclude that despite the diversity found at the Churince soils, the small patches (1 kg mesh bags) did not recover to resemble undisturbed soil in either site, after a year of migration and succession.
Nevertheless, it seems they followed a parallel paths for such recovery at the beginning of the succession and diverged in later stages.

Conclusions
This is the first study on bacterial soil dynamics conducted at the CCB. Our results provide important insight into microbial community dynamics in response to a disturbance (secondary succession). In this study we found changes in community composition across time that were indicative of the successional process. Our data showed evidence of similar initial colonizers in both sites, but the divergence in later successional stages reflects stochastic factors suggesting a species sorting and neutral Notwithstanding, it is interesting that succession of small, perturbed spots is very slow; this reveals the temporal scale is important for this community in terms of resilience, but general long-term monitoring is necessary to better understand the temporal patterns and natural variability of this area. Changes in microbial communities due to disturbance may directly affect ecosystem processes, which are vital in a protected natural area, threatened by over-exploitation of aquifers such as is occurring in the CCB.

Figure 1
sampling effect Rarefaction curves of the A) Dry Lagoon and B) River in all sampling dates (3, 6, 9 and 12 months).
OTUs were determined at 97% sequence identity.

Changes in biogeochemical parameters
Mean ± SD of biogeochemical parameters of microcosms and undisturbed soil samples of each site for the last sampling date, 12 months (n=3). A two way ANOVA was used for statistical comparison and results are presented as F values with significance levels: * = P < 0.05, ** = P < 0.01, *** = P < 0.0001, ns = not significant.   Table 5. Spearman's rank correlations between the relative abundances of the most abundant bacterial phyla, proteobacterial classes and bacterial families, and the soil properties in CCB. Bold numbers: P < 0.05; Bold and underlined numbers P < 0.001.