Spatial and Temporal Variation of Archaeal, Bacterial and Fungal Communities in Agricultural Soils

Background Soil microbial communities are in constant change at many different temporal and spatial scales. However, the importance of these changes to the turnover of the soil microbial communities has been rarely studied simultaneously in space and time. Methodology/Principal Findings In this study, we explored the temporal and spatial responses of soil bacterial, archaeal and fungal β-diversities to abiotic parameters. Taking into account data from a 3-year sampling period, we analyzed the abundances and community structures of Archaea, Bacteria and Fungi along with key soil chemical parameters. We questioned how these abiotic variables influence the turnover of bacterial, archaeal and fungal communities and how they impact the long-term patterns of changes of the aforementioned soil communities. Interestingly, we found that the bacterial and fungal β-diversities are quite stable over time, whereas archaeal diversity showed significantly higher fluctuations. These fluctuations were reflected in temporal turnover caused by soil management through addition of N-fertilizers. Conclusions Our study showed that management practices applied to agricultural soils might not significantly affect the bacterial and fungal communities, but cause slow and long-term changes in the abundance and structure of the archaeal community. Moreover, the results suggest that, to different extents, abiotic and biotic factors determine the community assembly of archaeal, bacterial and fungal communities.


Introduction
Understanding temporal and spatial patterns in the abundance and distribution of communities has been a fundamental quest in ecology. Such an understanding is crucial to allow an anticipation of responses of ecosystems such as soil to global changes [1]. Because local conditions are never constant, small disturbances that affect the soil microbial communities might occur [2][3] at different temporal and spatial scales. The assessment of microbial communities at a particular locality may result in patterns that vary greatly both within and between years, and these communities may be subjected to changes over longer time scales as a result of processes such as succession and evolutionary change [4]. One approach to investigate temporal (and spatial) variability in complex systems is to explore patterns of b-diversity. Whereas alpha (a-) diversity represents a measure of the total diversity of a given site, b-diversity is the variation of species composition (turnover) across space or time between paired sites. High bdiversity indicates large differences in community composition among different sites. Such high diversity can result from local as well as regional factors, e.g. changes in the local environmental conditions or limitation of dispersal between sites [5].
Temporal variation of conditions is a very common feature of ecosystems. Ecologists have long been interested in how such variation structures natural communities [6,7]. It can presumably affect the rate of microbial turnover, as microorganisms can process resources and adapt to changes in natural environments on a much faster time scale than macroorganisms [8]. Moreover, many functional microbial groups can show dramatic seasonal changes in soils [9].
The number of studies employing the concept of b-diversity to understand how microbial communities respond to biotic and abiotic parameters has increased substantially in soil ecology. Martiny and co-workers [10] studied the mechanisms driving ammonia-oxidizing bacterial (AOB) communities in salt marsh sediments. They found no evolutionary diversification when comparing the AOB community composition between three continents; although a negative relationship was observed between geographic distance and community similarity. Furthermore, in an attempt to determine to which extent a bacterial metacommunity that consisted of 17 rock pools was structured by different assembly mechanisms [11], the authors studied changes in b-diversity across different environmental gradients over time, including phosphorus concentration, temperature and salinity. They found that there were temporal differences in how the communities responded to abiotic factors. b-diversity allows not only the understanding of temporal but of spatial variations as well. For instance, in a survey of bacterial communities across more than 1000 soil cores in Great Britain [12], no spatial patterns were observed, but instead variations in b-diversity according to soil pH were found, which revealed that b-diversity (between sample variance in a-diversity) was higher in acidic soils (pH 4-5) than in more alkaline soils (pH 7-9) [12]. In the former soils, environmental heterogeneity was highest, calculated as the variance in environmental conditions [12]. In another study, different patterns of bacterial bdiversity were observed between different layers in sediment cores, which could be attributed to historical variation and geochemical stratification [13].
Of the soil microbial groups, bacteria have been mostly studied, as they exhibit an estimated species diversity of about 10 3 to up to 10 6 per g soil [14][15][16]. However, archaea and fungi are also important microorganisms found in soil. Previous studies have shown the ubiquity of archaea in soil, especially the crenarchaeota [17][18][19]. Fungal abundances in the order of 10 4 fungal propagules per g of dry soil were observed in Antarctic soils [20] and 10 7 per g of soil in soil crusts [21]. Fundamental differences in the physiology and ecology of members of such communities would suggest that their patterns of spatial and temporal variation are controlled by distinct edaphic factors.
In this study, we explored the temporal and spatial fluctuations of soil microbial communities and their relation to local environmental conditions. In order to do so, we investigated the spatiotemporal dynamics of the soil microbiota by analyzing the patterns of aand b-diversity of archaea, bacteria and fungi in eight agricultural soils across the Netherlands. We sampled the soils eleven times, from 2009 to 2011. Furthermore, to complement the analyses, we applied TLA (time-lag analysis) [22], a distance-based approach to study the temporal dynamics of communities by measuring community dissimilarity over increasing time lags. TLA provides measures of model fit and statistical significance, allowing the quantification of the strength of temporal community change in a numerical framework [23]. We thus interrogated how the relationship between microbial abundance, species composition and the surrounding environment varies in space and time and how this relates to long-term compositional changes.

Study area and field sampling
The eight soil sites sampled are located in the Netherlands. Their characteristics and geographical coordinates are found in Table 1 and in Table S1. Sampling points were selected to reflect temporal differences in abiotic parameters. For each soil four replicates were taken. Each replicate consisted of 10 subsamples (15-20 cm deep) collected between plots with a spade, away from the plant roots. Soil samples were collected four times over an annual cycle in 2009 (April, June, September and November), three times in 2010 (April, June and October), and four times in 2011 (February, April, July and September). Each sample was placed in a plastic bag and thoroughly homogenized before analysis. A 100-g subsample was kept at 4uC and used for chemical analyses, whereas the remaining soil was kept at 220uC for subsequent DNA extraction and molecular analysis of community composition and total abundance (see below).

Soil chemical analysis
The environmental variables measured included pH, concentrations of nitrate (N-NO 3 2 in mg/kg of soil), ammonium (N-NH 4 + in mg/kg of soil), organic matter (OM in %) and clay content (in %). The pH was measured in CaCl 2 suspension 1:4.5 (g/v) (Hanna Instruments BV, IJsselstein, The Netherlands). Organic matter (OM) content is calculated after 4 hours at 550uC. Nitrate (N-NO 3 2 ) and ammonium (N-NH 4 + ) were determined with a colorimetric method using the commercial kits Nanocolor Nitrat50 (detection limit, 0.3 mg N kg-1 dry weight, Macherey-Nagel, Germany) and Ammonium3 (detection limit, 0.04 mg N kg-1 dry weight; Macherey-Nagel, Germany) according to manufacturer's protocol.

Nucleic acid extraction
DNA was extracted from 0.5 g of soil using Power Soil MoBio kit (Mo Bio Laboratories Inc., NY), according to the manufacturer's instructions, after the addition of glass beads (diameter 0.1 mm; 0.25 g) to the soil slurries. The cells were disrupted by bead beating (mini-bead beater; BioSpec Products, United States) three times for 60 s. Following extraction, the DNA preparations were electrophoresed over agarose gels in order to assess DNA purity, quality (average size) and quantity. The quantity of extracted DNA was estimated on gel by comparison to a 1-kb DNA ladder (Promega, Leiden, Netherlands) and quality was determined based on the degree of DNA shearing (average molecular size) as well as the amounts of coextracted compounds.

Real-time PCR quantification (qPCR)
Absolute quantification was carried out in four replicates on the ABI Prism 7300 Cycler (Applied Biosystems, Germany). The 16S rRNA gene was amplified by qPCR using diluted extracted DNA as template. Specific primers for archaea (group 1 crenarchaeota) 771F/957R [24] and for V5-V6 region of bacteria 16SFP/16SRP [25] were used. We have chosen to focus on Crenarchaeota, as this group is often more common in soil environments than Euryarchaeota [26]. For Fungal community primers 5,8S/ITS1f [27] were chosen. Cycling programs and primer sequences are detailed in Table S2. The specificity of the amplification products was confirmed by melting-curve analysis and on 1.5% agarose gels. Standard curves were obtained using serial dilutions of plasmid DNA containing the cloned 16S rRNA gene obtained from Burkholderia terrae BS001 or ITS region of Rhizoctonia solani AG3. Dilutions ranged from 10 7 to 10 2 gene copy numbers/ml. The archaeal standard curve was obtained by serial dilution of PCR product generated from Cenarchaeum symbiosum with the aforementioned archaeal specific primers [24].

Data analyses
The diversity of each of the soil bacterial, archaeal and fungal communities was determined on the basis of the PCR-DGGE profiles. Total diversity (a) of the dominant community members was estimated from these data using the Shannon index, as recommended by Hill et al. [37], as well as the number of DGGE bands (species richness). We calculated the temporal b-diversity of archaeal, bacterial and fungal communities as the mean of all pairwise Bray-Curtis dissimilarities based on the relative abundance of DGGE bands, as previously described [38,39,11]. To support results from the calculated b-diversity and to test the statistical significance and strength of community dynamics we used time-lag analysis (TLA) [22] by plotting Hellinger-transformed [40] distance values against the square root of the time lag for all lags. The time-lag analytical approach can produce a number of general theoretical patterns with time-series data [22]. The square root transformation reduces the probability that a smaller number of points at larger time lags will bias the analysis [41]. The Bray-Curtis matrices as well as Hellinger-transformed distances were determined in PRIMER-E (version 6, PRIMER-E Ltd, Plymouth, UK; [42]).
To test how a-diversity, b-diversity and microbial abundance varied in relation to environmental variables, parametric Pearson correlation coefficients were calculated between a and b diversities, soil pH, organic matter, nitrate, ammonium, clay content and soil moisture, as well as between total abundances and TLA slopes using SPSS v18.0.3 (SPSS Inc., Chicago, IL, USA). All variables except pH were transformed (Log(x+1)) prior to all analyses. Moreover, we applied variance partitioning to evaluate the relative contribution of the drivers of the microbial assemblages. Forward selection was used on CCA (Canonical Correspondence Analysis) to select a combination of environmental variables that explained most of the variation observed in the species matrices. For that, a series of constrained CCA permutations was performed in Canoco (version 4.0 for Windows, PRI Wageningen, The Netherlands) to determine which variables best explained the assemblage variation, using automatic forward selection and Monte Carlo permutation tests (permutations = 999). The length of the corresponding arrows indicated the relative importance of the chemical factor explaining variation in the microbial communities.

Variability of environmental parameters
Soil pH, nitrate, ammonium and organic matter levels were determined in triplicate across all soil samples. Soil pH was significantly higher (P,0.05) in soils K, G and L (7.  (Table S1).
Considering each soil individually, they had characteristically different values, with higher levels of nitrate and ammonium found in soils B, V, D and S than in soils W, K, G and L ( Temporal variations in the abundance of archaeal, bacterial and fungal communities and their responses to abiotic variables We studied the variations in the abundances of archaeal, bacterial and fungal communities over time and across all samples in three years. The total bacterial abundance showed significant temporal variation during the whole period, ranging between  ), and significantly higher abundance was observed in the sandy soils depending on the sampling time (Fig. 1). Overall and on average, the abundance of bacteria was higher than that of the fungi, except in September 2009 and during 2011. We used Pearson's correlation to examine how soil parameters influenced the abundances of bacterial, archaeal and fungal communities. Whereas the archaeal abundances were positively correlated with soil pH (r = +0.883, P,0.001), they were negatively influenced by nitrate (r = 20.764, P,0.05). A positive relationship was observed between fungal abundance and soil organic matter (r = +0.722, P,0.05), and a negative relationship was observed between fungal abundance and archaeal abundance (r = 20.484, P,0.05). Relationships between the abundance of bacteria and fungi, as well as between bacteria and archaea, were not significant. Interestingly, none of the soil parameters measured influenced bacterial abundance significantly.

Patterns of a-diversity and response to abiotic variables
Understanding how species are distributed in space and time may yield a first avenue towards their assembly rules [43]. We used two ecological measures, i.e. the Shannon index (H9) and species richness, as proxies to study the variations in the adiversities of the archaeal, bacterial and fungal communities. Differential patterns of archaeal, bacterial and fungal a-diversities were observed, as measured by H9 (Fig. 2). The H9 values of the archaeal communities ranged from 1.  (Fig. 2). In general, the differences observed between sandy and clayey soils for the bacterial and fungal diversities (Shannon index) were time point-dependent. For archaea, a higher Shannon index was noticed in the sandy soils compared to the clayey ones in 2009 and 2010 but not in 2011.
Concerning correlations with edaphic factors, a positive effect of OM content was observed on the archaeal a-diversity (r = +0.691, P,0.05) ( Table 2). When using the number of DGGE bands as a measure of a-diversity (species richness), a significant and strong positive correlation was found between archaeal a-diversity and nitrate levels (r = +0.962 , P,0.001) ( Table 2). None of the soil parameters measured correlated significantly with bacterial or fungal a-diversity.

Patterns of temporal b-diversity and responses to abiotic variables
The patterns of temporal b-diversity of the archaeal, bacterial and fungal communities (taking into account the variations in community composition of each microbial group in individual soils over time) showed small but significant variations across soils (Fig. 3A). Bacterial b-diversities were in general higher than fungal ones across soils, except for soil V. There were slight but significant differences (P,0.05) between sandy and clayey soils regarding the temporal b-diversity of archaeal and bacterial but not of fungal communities (Fig. 3B).
Although chemical parameters might show variability over time, significant correlations could still be observed. The patterns in the archaeal temporal b-diversities observed were mainly due to positive correlations with nitrate (r = +0.874, P,0.05) ( Table 2). None of the soil parameters measured were correlated with bacterial and fungal temporal b-diversities. Canonical correspondence analysis was used to test the significance of the influence of soil parameters on the community parameters. We used variance partitioning to control for the effect of each individual parameter, while all others are defined as covariables in the constrained analyses [44]. Considering the whole data set, soil parameters explained 45%, 6.6% and 6.9% of the temporal variability in archaeal, bacterial and fungal community structures, respectively. The archaeal community was mostly affected by changes in OM (11.9%) and nitrogen specimens (nitrate + ammonium; 7.8% each), whereas the bacterial and fungal community variations were mostly related to ammonium (2.1% and 2.2% for bacteria and fungi, respectively) (Fig. S1).
Significant relationships were observed between variation in bdiversity and H9 for bacterial (r = +0.602, P,0.05) and fungal communities (r = 20.481, P,0.05) but not for archaeal communities.  Table 2. Pearson's correlation coefficient between soil chemical parameters, biotic parameters (total abundance, alpha diversity, beta diversity and slopes from TLA analysis), calculated from the eight soils over time. Quantifying temporal changes of archaeal, bacterial and fungal communities The temporal changes of the microbial guilds were quantified and statistically tested via TLA. TLA analyses were performed separately per year and also considering all three years. For both analyses, the results and conclusions were similar. Therefore, we decided to include only the whole three year dataset. A statistically significant regression line (P,0.05) was observed for the archaeal community, with an overall slope of 1.835 (Fig. 4A and Table 3). Moreover, all eight soils showed indications of directional changes in community composition, yielding regression lines that were statistically different from zero (P,0.05) with the exception of the G soil (Table 3). Although the slopes were small (Table 3), they were mainly reflected in the positive Pearson correlations with nitrate levels (r = +0.814, P,0.05). The bacterial communities showed a similar trend as observed for the archaeal ones, with significant regression lines and a slope of 0.785 (Fig. 4B and Table 3). Although analyses of the fungal communities in the eight soils showed that only three soils were undergoing directional changes (B, V an L soils), the overall result based on the simultaneous analysis of all soils yielded statistically significant regression lines (slope of 0.638, P,0.05, Fig. 4C). None of the soil parameters measured had significant effects on the rates of change of bacterial and fungal communities. Significant and contrasting relationships were observed between the TLA slopes and the H9 values of archaeal (r = +0.629, P,0.05) and bacterial (r = 20.523, P,0.05) communities but not of fungal communities.

Temporal variation in the abundance of soil microbial communities
In our study, population sizes of archaea, bacteria and fungi, estimated using quantitative PCR, were found to be within the range observed in other soil systems [24,45]. Quantitative PCR of soil DNA, as any PCR based approach, has its inherent limitations, whether it be the biases of soil DNA extraction, PCR, or the core genes targeted. However, the method is highly reproducible and sensitive, enabling the quantification of microbial abundance changes across temporal and spatial scales. Moreover, this study performed multiple qPCR runs in order to ensure results were statistically significant. In our calculations, we also took into account the efficiency and amount of extracted DNA from the soil samples. Therefore, we argue that our results are representative of the fluctuations observed between different times rather than pure noise.
A high abundance of crenarchaeota in soils has been previously observed [24,46], possibly indicating a crucial functional role for such organisms in agricultural soils. Furthermore, the bacterial abundance was often higher across soils than the fungal abundance (except at the end of 2009 and the end of 2011), supporting the finding that bacterial:fungal (B:F) ratios are quite high in agricultural or grassland soils as compared to, say, forest soils [45,[47][48][49]. These comparisons are very important in the context of whether soils are thought of a being fungal (more ''natural'') or bacterial (more highly cultivated) dominated. Indeed, such elevated B:F ratios may also reflect anthropogenic disturbances due to agricultural practices.
The variations in microbial abundances could be explained by several parameters, depending on the target group. Soil pH and nitrate explained more than 75% of the variation in archaeal abundance. Previous studies have reported negative effects of pH on group 1.1c Crenarchaeota [50] in acid forest soils and negative relationships between nitrate and archaeal abundance [21]. The positive correlations between archaeal abundance and soil pH observed here suggest that our soils might be dominated by crenarchaeal species that are adapted to conditions of higher soil pH (7.0-7.5) [51], which may be linked to the long agricultural history of the plots studied here.
Interestingly, the bacterial abundances didn't respond to soil pH or any other measured abiotic parameter, although several studies have reported pH as the main determinant of bacterial community composition [52][53][54][55][56][57]. It has been shown that some specific bacterial taxa abundances decrease or increase with a changing pH, for instance members of the Acidobacteria and Actinobacteria [56]. Although pH may have driven changes in the relative abundance of some bacterial classes, the abundance of total bacteria remained quite constant in the different pH ranges, indicating that the carrying capacity of the soil was not strongly affected by pH. Fungal abundance was also not affected by pH, but this was expected since the pH range in our soils was within the (wide) pH optimum for this group, often covering 5-9 pH units without significant inhibition of growth [58,59]. We also observed that when conditions apparently favored increases in fungal abundance, archaeal abundance decreased, suggesting that fungi and archaea might compete for similar niches. Nonetheless, fungal abundance was positively affected by OM content, which is consistent with the saprophytic status of most fungi [60].

Temporal variation in a-diversity
None of the soil parameters measured in this study were able to explain the patterns of a-diversity observed for bacteria and fungi. It might be that the taxonomic scale was too broad and a deeper analysis would allow a better understanding of the observed patterns, as observed by Rasche et al. [61]. With PCR-DGGE, only the most abundant taxa, comprising no more than 0.1-1% of the community, can be detected. In other words, only the most abundant organisms are detected with PCR-DGGE. Because of these caveats, the parameters calculated from PCR-DGGE fingerprints and correlations based thereon should be interpreted as indications and not as absolute conclusions.
Archaeal a-diversity, on the other hand, was shown to respond to nitrate and OM levels. Nitrate had opposite effects on archaeal richness and diversity, depicting a community that responds to increasing nitrate with an increase of richness but a great decrease of evenness. This most likely indicates the outgrowth of previously undetectable OTUs. Thus, in addition to the strong negative effect on archaeal abundance observed by qPCR, nitrate availability seems to be a crucial factor determining archaeal community structure. The positive correlation between archaeal diversity and soil OM content indicate that the OM provides a substantial fraction of carbon to the local archaeal communities. Recently, genomic analyses of Crenarchaeum symbiosum and Nitrosopumilus maritimus suggested that these organisms are capable of mixotrophy [62,63] and that group 1.1c Crenarchaeota are able to grown on methanol and methane [64]. This suggests that archaea might not be solely sustained by ammonia oxidation [65,66].

Variation in b-diversity over time (species turnover)
To assess how dynamic each soil microbial group was over time, we calculated the temporal b-diversity (average Bray Curtis dissimilarity) for each soil over time and for each microbial group. We observed a higher temporal b-diversity for archaea than for bacteria and fungi across all soils. This indicates that the archaeal communities are much more dynamic than the bacterial or fungal ones along a time gradient. These differences are probably due to the differential physiologies and sensitivities to environmental perturbations of these microorganisms. It has been shown that changes in temperature, moisture [61,67] and resource availability due to seasonal variation [61] can affect soil archaeal as well as bacterial communities. Moreover, a clear pattern was observed for bacterial b-diversity in the metacommunity of 17 rock pools, with higher variations during summer and lower during autumn [11]. The temporal variations of archaeal and bacterial communities were also higher in the clayey soils than in sandy ones, suggesting that the latter harbors more dynamic communities.
One main finding of this study is that, although the b-diversity patterns of the three microbial domains investigated are related with the same set of abiotic factors, the total percentage of variation able to explain those patterns was much higher for archaeal (45.0%) than for bacterial (6.6%) or fungal (6.9%) communities. This suggests that the archaeal communities might be much more sensitive to environmental changes than the bacterial or fungal ones. Based on these results, we hypothesize that the archaeal communities of agricultural soils with a long history of N-fertilization are more sensitive to disturbances than the corresponding bacterial or fungal communities.

Quantification of b-diversity
To be able to quantify community dynamics, allowing comparisons and providing a general overview of long-term trends in the complex soil system, we used the slopes obtained from TLA. TLA has been intensively used to identify directional changes and to quantify temporal dynamics of macroorganisms [69][70][71][72][73]but very few studies have focused on microorganisms. Although the TLA slopes for archaea and bacteria were small, they were significantly different from neutral. Clearly, even small changes can be part of a long-term trend. On the contrary, changes in the fungal communities were non-significant, suggesting stochastic species dynamics.
Changes in environmental variables within soil sites determine how time affects turnover (b-diversity), as different microbial interactions are favored if prevailing conditions change [74]. Only archaeal communities responded to changes in environmental parameters, being strongly correlated with nitrogen availability and with the degree of temporal variation quantified by TLA. This might suggest that at some level strongly deterministic processes are acting on the archaeal but not on the bacterial and fungal communities in these soils. Another explanation is that archaea are much more limited in their ecoversatility, whereas bacteria and fungi are highly functionally redundant. The observed relation between bacterial richness and TLA slopes, e.g. high turnover at low richness, was also noticed in a study of the distribution of British birds [75]. The authors discuss that low species richness areas tend to have relatively more random mixtures of species than high species richness areas. The relation observed between archaeal turnover and species richness suggests a less random distribution of species, caused mainly by nitrate contents.
In this study we demonstrate that changes in the community composition of bacteria and fungi could be linked to both environmental and biotic factors (e.g. species-species interactions), as their (a-) diversity co-varied significantly with their b-diversity over the time period of the study. Conversely, archaea showed no significant correlation between aand b-diversity, and the community shifts were mainly driven by the surrounding environment, mostly by the effects of soil pH and nitrate concentrations. This might indicate that changes in archaeal community are mostly driven by environmental factors, as previously observed by Zinger et al. [68] in a study on the patterns of archaeal, bacterial and fungal communities in an alpine landscape. Furthermore, we propose that different environmental and biological mechanisms act on each microbial niche. A more comprehensive understanding of the rules governing these important soil microorganisms will require additional field work as well as microcosm experiments to identify the key environmental and biotic factors driving the assemblage of these communities.
Ethic statement. No specific permits were required for the described field studies. The locations are not protected. The field studies did not involve endangered or protected species. Figure S1 Biplots of canonical correspondence analysis (CCA) of Archaeal, Bacterial and Fungal similarity matrices and vector fitting of the environmental variables. Similarity matrices from DGGE data were obtained from eight soils over three years (2009, 2010 and 2011). Physicochemical data, soil moisture (Humidity), soil nitrate (NO3), soil ammonium (NH4), organic matter (OM), clay content (clay + silt) and soil pH (pH) are presented with black arrows. (DOCX)