OUP user menu

Soil microbial diversity patterns of a lowland spring environment

Sotirios Vasileiadis, Edoardo Puglisi, Maria Arena, Fabrizio Cappa, Johannes A. van Veen, Pier S. Cocconcelli, Marco Trevisan
DOI: http://dx.doi.org/10.1111/1574-6941.12150 172-184 First published online: 1 November 2013


The Po river plain lowland springs represent unique paradigms of managed environments. Their current locations used to be swamps that were drained 6–7 centuries ago, and they have been in constant use ever since. Our aims were to identify the effects of land use on the microbial communities of these soils, look for associated diversity drivers, and assess the applicability of ecology theories with respect to identified patterns. We screened the microbial diversity across a land use transect via high-throughput sequencing of partial 16S rrRNA gene amplicons. Land use had a major effect on soil properties and microbial community structures. Total organic carbon and pH were major diversity drivers for Bacteria, and pH was important for Archaea. We identified the potential contribution of soil amendments to the indigenous microbial communities, and also gained insights into potential roles of taxa in the organic carbon turnover. Verrucomicrobia coincided with the higher values of the recalcitrant organic carbon. Actinobacteria and Acidobacteria correlated with the more labile organic carbon. Finally, the higher diversity found in the soils less enzymatically active and relatively poorer in nutrients, may be explained to an extent by niche-based theories such as the resource heterogeneity hypothesis and Connell's intermediate disturbance hypothesis.

  • disturbance
  • land use
  • microbial community
  • next generation sequencing
  • soil


Soil is a highly complex and important matrix encompassing immense bio-diversity and a large number of biological processes (Barrios, 2007). Many of these processes contribute to ecosystem services (e.g. nutrient cycling, soil erosion control and biological pest control) supporting human activity in both natural and managed environments. However, human input and related perturbations, from attempts to enhance ecosystem services, alter soil traits close to or even beyond points at which natural attenuation mechanisms may not longer lead to functional restoration (Shennan, 2008).

A characteristic example of human intervention in natural ecosystems is the drainage of wetlands for agricultural use. Particular environments where this strategy was applied in the past are the lowland springs and their surroundings in the northern area of the Po river plains (Kløve et al., 2011), in Italy. Such drastic change in the physico-chemical circumstances may have a strong impact on the structure and functioning of organisms living in the system. This holds in particular for microbes, which, in parallel, are the primary agents responsible for key ecosystem processes. Although the effects of land use and management on the soil microbial community structure have long been acknowledged (Guo & Gifford, 2002; Kelliher et al., 2005; Acosta-Martínez et al., 2008; Wallenius et al., 2011), relevant detailed information is scarce (Garbeva et al., 2004). The vast numbers of microbial inhabitants of complex soil environments (Schloss & Handelsman, 2006) hampered efforts for detailed diversity screening during previous decades. Recent high-throughput sequencing technologies have overcome some methodological problems (Sogin et al., 2006; Roesch et al., 2007; Lauber et al., 2009; Gloor et al., 2010; Bartram et al., 2011; Bates et al., 2011; Delmont et al., 2012; Uroz et al., 2013), therefore increasing the depth of investigation and filling associated gaps. However, more issues need to be addressed to uncover the total existing microbial diversity in targeted, polymerase chain reaction (PCR)-based approaches, e.g. with selected nucleic acid extraction strategies, the primer choices, the PCR setup and sequencing design strategies, rendering the achievement of such a goal still expensive in terms of time and effort (Berry et al., 2011; Delmont et al., 2011; Lombard et al., 2011; Vasileiadis et al., 2012). Therefore, studies using approaches based on a single primer set and/or DNA extraction method, like the presented study, should restrict their claims to the part of the existing diversity that has been uncovered. To this extent, when using the term ‘total’ in conjunction with microbial community-associated terms hereafter, we refer to values and estimates that have been uncovered.

Our aim was to study the microbial communities of soils along a land use transect in the Gaverina (Bergamo, Italy) lowland spring environment (Supporting Information, Fig. S1). Combining the screening depth provided by novel high-throughput sequencing technologies and appropriate data analysis approaches, we explored 16S rRNA gene-based microbial diversity and associations with land use traits, soil chemical and biological properties. Further, we assessed the applicability of existing ecological hypotheses in explaining the patterns of shaping microbial communities in these soils.

Materials and methods

Soil sample physical-chemical-biochemical properties

Historically, all soil environments studied here originated from the bed of a swamp owing its existence to underground water tension. Drainage, which took place 6–7 centuries ago, led to a blossoming of agricultural activity in the area located next to the river Po and generated the land use and management gradients of the soils studied here (Kassen & Rainey, 2004) (Fig. S1: Fonti della Gaverina 45°27′55.69″N, 9°38′20.05″E, elevation 97 m, Italy – sampling carried out in May 2009). The gradients include (1) the low-land spring banks (saturated with water during the high precipitation season), (2) maintained meadow zones around the springs (plant biomass is a mixture of grass, legumes and native plants, with the above-ground biomass being removed twice a year and the plant composition being supported by no tillage sowing of legume seeds when their production decreases – every 3–5 years); and (3) the surrounding fields, cultivated during the sampling period with maize (part of a rotation program with legumes, grasses and maize, receiving cow slurry every 2 years). Top-soil samples (the upper 5–10 cm) were collected in triplicate with an auger from these three environments and soil from the selected depth range was sieved through a 2-mm mesh (soil properties and management traits are provided in Table 1996). Particle-size analysis was carried out using the pipette method (Day, 1965). Three carbon fraction measurements (labile, moderately labile and recalcitrant) were based on the Walkley–Black method as modified elsewhere (Chan et al., 2001) for obtaining the different fractions. To investigate potential links between the carbon fractions and the energy flow associated with the soil microbial communities we measured widespread associated biological activities: β-glucosidase (EC and acid phosphatase (EC activities were determined by the p-nitrophenol method of Eivazi & Tabatabai (1988) and Margesin & Schinner (1994), with p-nitrophenyl-β-d-glucoside and p-nitrophenylphosphate as substrate, respectively. Further, nitrate reductase, being genetically widespread (Philippot, 2002), was determined using KNO3 as a substrate according to Fu & Tabatabai (1989).

View this table:

Soil properties and qualitative land use and management traits

Management traitsMaizeMeadowRiparian
Slurry – TillagePlant biomass removalSeasonal soil saturation
Soil perturbation frequencyFrequentNot frequent
TextureLoam – Clay loamLoamLoam – Clay loam
Chemical propertiesAVGSDAVGSDAVGSDF statistic
pH***(c) 6.4± 0.10(b) 6.9± 0.00(a) 7.16± 0.1258.4
CEC (mEq 100 g−1)1.07± 0.181.79± 0.251.1± 0.712.5
Humidity (%)14.57± 2.2618.06± 2.2824.31± 9.982
N (%)0.17± 0.010.32± 0.030.2± 0.240.9
TOC (%)***(b) 1.43± 0.06(a) 2.98± 0.26(c) 0.74± 0.04158
TOC/N8.25± 0.089.34± 0.057.97± 5.57N/A
Labile OC proportion***(a) 0.74± 0.10(a) 0.63± 0.01(b) 0.37± 0.0328.3
Moderately labile OC proportion0.05± 0.020.11± 0.070.12± 0.032.1
Recalcitrant OC proportion**(b) 0.21± 0.08(b) 0.26± 0.08(a) 0.51± 0.0415.6
Enzymatic activities
Nitrate reductase (μg N g−1 24 h−1)59.43± 35.68120.07± 156.35159.2± 72.090.7
β-glucosidase (μg PNP g−1 h−1)***(a) 55.61± 10.35(a) 66.99± 6.80(b) 13.70± 5.5938.4
Phosphatase (μg PNP g−1 h−1)***(b) 61.29± 0.73(a) 170.07± 43.24(c) 6.22± 2.9833.3
  • Post-hoc pairwise comparisons: Tukey HSD (α < 0.05). Significant pairwise differences are denoted with different letters in parenthesis.

  • N/A: conditions were not met and non-parametric test was not significant.

  • anova significance: *** 0.001, ** 0.01, * 0.05.

DNA isolation, PCR conditions, multiplexing and sequencing

DNA was extracted from 500 mg of each soil sample using the FastDNA® SPIN kit for soil with a FastPrep® 24 instrument (MP Biomedicals, LLC, Solon, OH) according to the manufacturer's instructions. Extracts were quantified using the Quant-iT™ HS dsDNA Assay kit (Invitrogen, Paisley, UK) and the Qubit™ fluorometer (Invitrogen). DNA purity and shearing were tested using a Biophotometer (Eppendorf, Hamburg, Germany) and 0.8% agarose gel, respectively. DNA 2 ng extracts were used for the bacterial 16S rRNA gene PCR amplifications; 20 ng was used for the archaeal ones. The targeted V5 region of the 16S rRNA gene was selected on the basis of the available sequencing technology read length screening abilities (100-bp paired-end reads, related information is provided in the end of the present section) and also conservation of candidate V region flanking loci, as indicated in a previous database search and simulation study (Vasileiadis et al., 2012). PCR 50 μL reactions were performed according to the following program: 94 °C for 5 min, 35 cycles X [94 °C for 30 s of denaturation; 50 °C (for bacterial primer-sets) and 60 °C for 30 s of primer annealing for the bacterial and the archaeal 16S rRNA gene targeting primer-sets, respectively (Table S1); 72 °C for 30 s of elongation] and 72 °C for 10 min. PCR 50 μL reactions were performed using the following mixtures: 1 × PCR buffer, 2.5 mM MgCl2, 2.5 U AmpliTaq® Taq polymerase (Applied Biosystems, Foster City, CA), 0.4 mM of each dNTP, 0.5 μM forward primer and 0.5 μM of reverse primer, 2 and 20 ng respectively of template DNA when bacterial and archaeal 16S rRNA gene targeting primers were used. DNA extracts or PCR product were stored at −20 °C until further use.

Polymerase chain reaction (PCR) products were labelled with 6-bp indices (one for each soil sample, Table S2) as previously described by Meyer et al., (2008) and were concomitantly pooled in a single sample along with other indexed samples not related to the analysed sequences. The sequencing process of the PCR products pool was carried out by Fasteris SA (Geneva, Switzerland), using the TruSeq™ SBS v5 kit (Illumina Inc., San Diego, CA) for the sequencing library generation. The sequencing reaction was performed with a HiSeq 2000 Illumina instrument (Illumina Inc.). The samples were further multiplexed with bulk DNA samples in the same lane using the Illumina barcoded adapters, to avoid Illumina technology-related amplicon sequencing issues (Bartram et al., 2011).

Dataset preparation

The software used for the base-calling and demultiplexing (according to Illumina barcodes) processes was hiseq control soft v1.1.37, rta v1.7.45 and casava v1.7. Sequences were further demultiplexed according to the sample indexes provided in Table S2 and quality-controlled using the Fastx-toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). Good quality forward and reverse amplicon reads that could unambiguously re-construct the amplicons of their origin were used in follow-up analyses. Considering that most V5 region amplicons are shorter than 150 bp (Vasileiadis et al., 2012), and also the method used for sample indexing, the 100-bp length of each of the two obtained reads per amplicon was sufficient for re-constructing the vast amplicon majority. This task was performed with the ‘pandaseq’ script, previously made available by Bartram et al., (2011), with the restrictions of at least 10 bp of overlap between read pairs and no allowed mismatches, resulting in 3 969 194 assembled reads. Further quality control steps resulted in 2 989 302 high quality reads (average and minimum PhredQ values of > 28 and ≥ 20, respectively) with a mean length of c. 162 bp (including index and primer sequences) for all samples (Archaea and Bacteria). These reads were further demultiplexed according to soil sample of origin and targeted microbial group, and their unassembled versions were submitted to the Sequence Read Archive (SRA) database of the National Center for Biotechnology Information (NCBI, http://www.ncbi.nlm.nih.gov/Traces/sra) and are available under accession numbers SRS396685, SRS396687 and SRS396690-705.

Data analysis

Qualified sequence reads were used for operational taxonomic unit (OTU) and taxonomy-supervised analyses using mothur software v1.29 (Schloss et al., 2009). Sequences were pre-processed for analysis by: (1) sequence chimera control with Uchime (Edgar et al., 2011) with the exception of archaeal sequences (due to dependence of the method in existing databases which are still poor for Archaea); (2) removing non specific amplicons that were classified with the Bayesian classifier (Wang et al., 2007) for 50% bootstrap cutoff (Claesson et al., 2009) in non-targeted kingdoms, according to the available mothur version of the ribosomal database project (RDP) database; and (3) down-sampling each sample sequence set to the sample represented by the lowest sequence number per taxonomical kingdom in order to achieve equal sequence representation of samples. This choice was previously shown to reduce sequencing artifact effects on the data analysis (Schloss et al., 2011) and also biases associated with the calculation of several diversity indices (Gihring et al., 2011). All referred steps resulted in analysis of a total of 781 560 (86 840 per sample) sequences for Bacteria and 35 775 (3975 per sample) sequences for Archaea for all samples. In the OTU-based analysis, sequences were aligned with the nast algorithm (DeSantis et al., 2006; Schloss, 2009, 2010) and a kmer search approach vs. kingdom-specific SILVA pre-aligned reference databases as amended by the mothur software developers (Schloss et al., 2009; Quast et al., 2013). Pairwise distances of the aligned sequences and the average linkage algorithm (Sun et al., 2009; Schloss, 2010; Schloss & Westcott, 2011) were used for sequence clustering into OTUs encompassing sequence reads with a distance equal or < 3%. A taxonomy-supervised approach was also followed to evaluate relationships among samples and environmental variables with identified taxa. In this approach, the sequences sampled for equal representation, as described above, were classified using the Bayesian classifier (Wang et al., 2007) with 50% bootstrap cutoff value (Claesson et al., 2009). The greengenes (McDonald et al., 2012) and the silva databases (Quast et al., 2013) were used during classification of the bacterial and archaeal datasets, respectively, due to their coverage of the respective taxa.

Calculated indices related to α-diversity were the Good's diversity coverage estimate (Good, 1953), the observed richness (S), the non-parametric Shannon (npH′) (Chao & Shen, 2003) and Shannon-based Equitability [EH = H′/Hmax = H′/ln(S)] (Sheldon, 1969). Analysis of variance of means (anova) and the Tukey's honestly significant difference (HSD) pairwise comparison test (α < 0.05) were performed to assess significant differences between land use and management with respect to total microbial community diversity indices. The anova assumptions were tested using the Shapiro normality test and the Levene's test of homogeneity of variance. Where anova was not applicable, the Kruskal–Wallis non-parametric test for significant differences estimation and the Nemenyi–Damico–Wolfe–Dunn joint ranking test (for confidence intervals of 99%) were performed, along with a Tukey post-hoc analysis for pairwise comparisons. Total community correlation with the environmental variables was performed with Mantel tests, using the Bray–Curtis distance metric for microbial data and the Spearman r coefficient. Correlations between individual OTUs and between the diversity indices with the measured environmental variables were assessed with pairwise correlations using the Spearman r coefficient. The microbial data matrices were subjected to Hellinger transformation (Legendre & Gallagher, 2001) prior to the performance of the multivariate analyses. To assess total community tendencies we performed unconstrained multivariate analyses such as hierarchical clustering using the average linkage algorithm, principal components analysis, principal coordinates analysis (PCoA) on the Bray–Curtis sample distance matrix and analysis of similarity (anosim) of the Bray–Curtis distances. The effect of environmental variables was assessed with distance-based redundancy analysis (db-RDA; Legendre & Anderson, 1999) on the Bray–Curtis distance estimations of the microbial abundance matrices, and permutation tests (9999 permutations) were carried out to assess significance of results. Prior to beginning db-RDA analysis, the measured environmental variables were normalized using the Box–Cox power law transformations (Osborne, 2010) and the output data values were expressed as deviations from the mean (z-scores) to eliminate the effect of unit and distribution range differences. Significant constraints were selected by the forward selection method (Blanchet et al., 2008). Performance of the db-RDA was preferred over canonical correspondence analysis when the length of the first axis of the performed detrended correspondence analysis was lower than 2, indicating linear responses of microbes in the observed environmental gradients (Lepš & Šmilauer, 2001).

The α-diversity-associated indices were calculated with modules of mothur software v1.28.0 (Schloss et al., 2009), and the rest of the statistical analyses were performed in r v2.14.2 (R Development Core Team, 2011) using the packages Vegan (Dixon, 2003), Hmisc (http://biostat.mc.vanderbilt.edu/wiki/Main/Hmisc) and Scales (http://cran.r-project.org/web/packages/scales/index.html).


Soil properties

Soil pH, TOC and its fractions were the chemical properties that showed statistically significant differences among land use types (Table ). Higher average values for most properties were derived from meadow (highest TOC, CEC, N content) or riparian samples (highest recalcitrant OC and pH). The exception was the labile OC proportion, where the highest value was observed for the maize samples. Statistically significant differences were shown for the measured biochemical activities of β-glucosidase and acid phosphatase, showing similar patterns between environments. Both of them had higher values in meadow, followed by maize, with β-glucosidase being significantly different between the riparian and the other two environments, whereas for acid phosphatase the differences were significant between all three tested soil environments. Nitrate reductase, on the other hand, did not show significant differences between environments.

Taxonomy-supervised analysis

About 80% of the bacterial sequences were classified to family level, whereas only slightly over 40% of the archaeal sequences was classified to class level for a 50% bootstrap cutoff with the Bayesian classifier (Fig. S2). The dominant bacterial sequences were Acidobacteria,Actinobacteria,Firmicutes,Proteobacteria and Verrucomicrobia (Fig. a). Db-RDA performance with land use as constraining factor showed significantly different communities among sample groups. Highly correlated taxa within each sample group were fermenters such as Bacilli and Clostridia for the maize group, Rhizobiales for meadow samples, and Planctomycetes,Bacteroidetes and Rhodocyclales for the riparian group. TOC and pH were significantly correlated environmental variables, with the shifts observed for the bacterial communities. Dominant archaeal taxa were Thermoplasmata,Methanobacteria,Methanomicrobia,Thaumarchaeota,Halobacteria and several unclassified Euryarchaeota (Fig. b). There was a significant distinction between land use types for Archaea, with the maize soils being clearly separated from the other two soil groups. Thermoplasmata were more abundant in the maize samples, while sequences in unclassified Euryarchaeota were more abundant in meadow and riparian samples. pH was the sole environmental variable significantly associated with community shifts according to db-RDA results for Archaea.

Relative participation of taxa per sample (upper), db-RDA performance with the land use as constraining factor (middle) and db-RDA performance with variable constraints according to the forward selection approach (lower), for the bacterial (a) and archaeal (b) taxonomical classification datasets. Significance of the results was based on permutation tests (9999 permutations); the percentage of variance explained by each db-RDA model and the significance is provided in the top right corner of each plot; it is analyzed among variables in the low corner; axes percentages provide the proportions of the plotted constrained variance; dash-dotted lines show the standard deviations per soil group.

Correlation of OTU diversity with land use type

The achieved coverage of the identified total diversity according to the dedicated sequencing effort in the data analysis, was over 95% and 80% for the bacterial and archaeal datasets, respectively (Fig. ) for the 3% sequence–distance defined OTUs (referred to as OTUs hereafter). α-Diversity screening showed that bacterial communities were more diverse, OTU-rich and even for the riparian and the maize soils compared with the meadow soils. However, statistically significant differences were observed only between the riparian and meadow soils (Fig. a). An inverse relationship between soil groups and indices was observed for Archaea. Riparian soils were still more OTU-rich, diverse and even, but with maize soils having the lowest values (Fig. b). β-Diversity analysis of the sampled soils using unconstrained ordination showed clear differences between the microbial communities of the different land use soil groups (Fig. S3). The OTU datasets demonstrated a lower level of association of the archaeal samples from the same land use groups than was found for the bacterial samples. In both cases the meadow samples had higher within group distances compared with the other soil groups.

Boxplots of the coverage estimate, the observed richness (S), the npH′, and the EH′, for the bacterial (a) and archaeal (b) analysed datasets. Significant differences among soil groups (ma: maize, me: meadow, ri: riparian) are indicated by different letters according to the performed anova and Tukey (HSD) test (P≤ 0.05).

OTU-based assessment of diversity drivers

Sample distances according to their OTU composition were compared with differences of the measured soil property values. Of the measured environmental variables, Mantel tests showed significant correlations between the total bacterial community shifts and TOC, together with its labile and moderately labile fractions, as well as pH (Table ). Furthermore, shifts in the phosphatase and β-glucosidase activities were significantly correlated with the total bacterial community shifts. For the archaeal community matrices, the correlated variables were pH and β-glucosidase. Spearman r correlation coefficient tests between the total community diversity indices and measured variables, showed significantly inverse ranked relationships between the npH′ index and the Shannon Equitability (EH′) with TOC and its labile fraction for Bacteria (Fig. S4). Three measured variables were significantly correlated with three archaeal diversity indices (Fig. S5). The pH and the recalcitrant OC fraction were positively correlated with the S, the npH′ and the EH′, whereas the β-glucosidase activity was negatively correlated with these three indices.

View this table:

Mantel r-values and test significance (999 permutations) for generated OTU relative abundance matrices

Bacteria Archaea
r Significance r Significance
pH0.33 * 0.50 *
TOC0.36 ** 0.25.
OC-labile0.32 * 0.20
OC-moderately-labile0.23 * 0.11
N0.29 * 0.27.
Nitrate reductase0.060.03
Phosphatase0.34 ** 0.18
β-glucosidase0.33 * 0.47 *
  • Alpha value: 0.1, 0.05 (*), 0.01 (**), 0.001 (***).

Individual OTU associations with the environmental variables were examined via the Spearman r coefficient rank correlation tests. For P cutoff values of 0.01, the r coefficient values were above 0.7 or below –0.7 for all correlated OTUs (Fig. S6). Sequences of these OTUs were extracted from each dataset and assigned to taxa (Fig. ). Overall, higher proportions of the archaeal genotypes showed correlations with chemical and enzymatic variables (up to c. 52.6%) compared with bacterial genotypes (c. 3.8%). However, correlations yielding high sequence numbers were observed in a narrower set of the considered variables for Archaea (e.g. pH and negatively correlated OTUs with moderately labile OC) than for Bacteria. Although the taxonomy levels used in this analysis were quite broad, some taxa showed, for the majority of correlated sequences, distinct responses to the measured variables.

Taxonomical assignments of sequences residing in the highly correlated OTUs (P≤ 0.01) with measured variables according to the Spearman test, for the bacterial (a) and the archaeal (b) datasets. The sizes of the stacked bars are proportional to the total sequence reads used per dataset.


Several rhizobia were positively correlated with total N. Actinobacteria were correlated with soils having less labile carbon and higher measured phosphatase activity. Acidobacteria were more associated with lower pH values (although there were genotypes with opposite responses) and increased TOC and labile OC. Firmicutes were correlated with the lower humidity, lower pH and less labile OC levels. Verrucomicrobia were positively associated with phosphatase and all carbon fractions except for the moderately labile carbon. Furthermore, several Proteobacteria were correlated with higher pH and CEC levels, whereas the lower values of these variables were dominated by other taxa (e.g. Firmicutes, Acidobacteria, Actinobacteria). The increased phosphatase and β-glucosidase activities were mostly correlated with Acidobacteria,Actinobacteria,Verrucomicrobia and some proteobacterial taxa. High levels of nitrate reductase activity coincided with some Bacteroidetes and proteobacterial taxa.


Methanobacteria,Methanomicrobia and Thermoplasmata were abundant in soils with a lower pH. Methanobacteria were well correlated with total and labile OC and also phosphatase and β-glucosidase enzymatic activities. Thaumarchaeota were positively associated with the measured nitrate reductase activity: a higher unclassified euryarchaeotal sequence incidence was found in soil samples with higher pH and less labile OC. Sequences classified as Halobacteria showed a preference for soils with higher pH values.


The land use differences between the three studied soil groups had profound effects on the biogeochemical properties of soil and its microbial community structures. We have used both the OTU and the taxonomy-supervised analysis in an attempt to fill methodological gaps (Schloss & Westcott, 2011; Sul et al., 2011) and exploit the information potentials of each of these approaches with respect to our objectives.

The three land use types led to the development of distinct microbial communities for Bacteria and Archaea, as supported by both the taxonomy and the OTU-based approaches (Fig. 2008; Rousk et al., 2009; Kuramae et al., 2011; Lopez-Sangil et al., 2011). Similarly, factors significantly correlated with total bacterial or archaeal community shifts were identified in these variables, but with differences concerning taxon responses. Total Bacteria were correlated with TOC and pH gradients according to the taxonomy-supervised approach (Fig. 2011). Previously, the participation of archaeal members was estimated to be between 0% and 10% of that of the prokaryotic members (Bates et al., 2011), whereas pH was indicated to affect this balance drastically by causing reduction of archaeal member abundance with lower values of pH (Bengtson et al., 2012). The latter might also contribute to the increased difference of the archaeal maize communities from the other two soil environments compared with the difference between meadow and riparian soils. The corresponding pH differences were ≥ 0.5 between maize and meadow or riparian areas, whereas the meadow–riparian difference was 0.22. In the total bacterial community, results did not demonstrate equal dependence on the pH changes, showing distinct communities between all three soil groups. Furthermore, the previously identified abundance ratios (Bates et al., 2011) can explain the lower percentage of bacterial sequences that participated in the responsive OTUs, compared with Archaea (see Fig. ) and also the total community distance-based correlation with variables (Fig. , Table ).

An interesting finding of this study is the high abundance of halobacterial assignments (24.5% of total) in the archaeal sequence dataset (Fig. 2004; Oren & Ventosa, 2009; Oxley et al., 2010).

In an attempt to identify the diversity drivers in the studied soils, we performed modeling of the microbial responses to environmental variables with db-RDA and correlation tests of the variables with individual OTUs. Results showed that maize soils were rich in microbial taxa commonly found in animal excrement (Wang et al., 2010) such as Clostridia and Bacilli and the euryarchaeotal Thermoplasmata (Figs 2003; Snell-Castro et al., 2005; Mao et al., 2011). The distinction of these taxa is also associated to the depth of our strategy, as previous works using increased depth approaches of their time, such as PCR-single-strand-conformation-polymorphism analysis (Peu et al., 2006), failed to detect the associated taxa in soil shortly after amendment with manure. A high rhizobial incidence in meadow is consistent with the presence of leguminous plants among the meadow plants (Welbaum et al., 2004), whereas an increased presence of Planctomycetes and Rhodobacterales, the Alphaproteobacteria known to thrive in freshwater environments, was observed in riparian soils (Crump et al., 1999; Schlesner et al., 2004; Imhoff, 2005; Buesing et al., 2009; Fuerst & Sagulenko, 2011). The significantly increased presence of these taxa along with the periodic exposure to water flow in the riparian soils, indicate more open and dynamic ecosystems compared with the other soils. With respect to bacterial diversity drivers, potential roles have been suggested by the results of the individual OTU–environmental variable correlation tests (Fig. 2009; Rui et al., 2009; Deangelis et al., 2011; Martinez-Garcia et al., 2012) have shown associations with different OC fractions. This could suggest a niche sharing concerning organic matter decomposition between these taxa, with Acidobacteria and Actinobacteria dominating the less recalcitrant OC decomposition stages, whereas the more versatile taxon of Verrucomicrobia, some of whose members are able to decompose more complex organic matter (Martinez-Garcia et al., 2012) and participate in methylotrophy (Chistoserdova et al., 2009), are positively correlated with the recalcitrant OC and (in a lesser degree) with the labile OC.

Another interesting finding was the reduced bacterial richness, diversity and evenness found in meadow soils as compared with the other soil environments. Such effects, when identified during fallow periods in crop rotation, have been interpreted as an outcome of resources depletion (Welbaum et al., 2004). However, measured enzymatic activities related to energy flow in soil systems (β-glucosidase, acid phosphatase), had higher values in meadow samples, whereas close relationships were found in meadow soils with root-associated (Rhizobiales) or complex organic molecule decomposing taxa (Verrucomicrobia,Actinobacteria) (Figs 2006). Therefore, rhizodeposits may favor a restricted number of microbes more efficient in metabolizing these compounds. Such a phenomenon is compatible with the previously proposed hump-shaped relationship of microbial diversity in response to the increase of dominant resources according to the resource heterogeneity hypothesis (Lynch et al., 2004). Thus, a reduced variety in resources may lead to a corresponding reduced number of niches, and therefore to reduced microbial diversity, as also identified in the meadow samples.

Another potential explanation of the identified diversity patterns is related to underlying disturbance phenomena. Disturbance effects are difficult to assess in field experiments because they are usually confounded with other factors contributing to measured variables, and a more disturbance-focused study strategy may be necessary before making associated claims. However, observed diversity indices (Fig. 1978) proposing that perturbation events, when not deleterious for the majority of existing populations, reduce dominance and promote diversity. Originally, the theory described the increase in under-canopy plant diversity after an intermediate disturbance, e.g. fire incidents with reduced effects on plant propagative material, due to post-disturbance lack of competition for resources. The proof-of-principle concerning the applicability of the theory in the microbial world was provided by a laboratory microcosm setup (Buckling et al., 2000). In that study, a unimodal response of pseudomonad population diversity according to colony phenotypes was demonstrated when static broth cultures were subjected to different disturbance frequencies. Similarly, soil bacteria with more autochthonous lifestyles may become more competitive when disturbance events affect the relative abundances of dominant populations, therefore contributing to a diversity increase as in the reported cases of the maize and riparian soils when compared with the less perturbed meadow soils.

Concluding remarks

In the present study we attempted an in-depth per-sample analysis of the microbial communities under the selective pressures associated with different land uses. The chemical and microbiological patterns identified, strongly suggest that land use and management decisions can have a profound effect on chemical and biological soil properties. In the soils studied, bacterial primary producers were more affected by related interventions in soils than were their archaeal counterparts. pH was a principal diversity driver for both Bacteria and Archaea, the latter being more dependent on it. Finally, we speculate that the increased diversity found in the energy-wise poorer and the more perturbed environments may suggest that soil microbial diversity is to a significant part regulated in line with niche-based theories, such as the resource heterogeneity hypothesis (Lynch et al., 2004) or Connell's intermediate disturbance hypothesis (1978).

Supporting Information

Additional Supporting Information may be found in the online version of the article:

Table S1. Primer sequences used in the present study. References correspond to Materials section reference list.

Table S2. Adapters used in the present study for sample indexing of PCR products provided in reference Meyer et al., (2008).

Fig. S1. Satellite photograph of the Gaverina area lowland spring (Source: ‘Gaverina’. 45°2757.32″-69″ North, 9°3824.21″-27.20″ East, Google Earth, March 30, 2012).

Fig. S2. Percentages of classified sequences per taxon level with the Bayesian classifier for 50% bootstrap cutoff, of the bacterial (A) and archaeal (B) sequence data.

Fig. S3. Unconstrained multivariate analyses performed on the bacterial (A) and archaeal (B) OTU data matrices.

Fig. S4. Significant Spearman r correlations (P ≤ 0.05) between bacterial dataset diversity indices and measured environmental variables.

Fig. S5. Significant Spearman r correlations (P ≤ 0.05) between archaeal dataset diversity indices and measured environmental variables.

Fig. S6. Spearman r values for significant correlations (P ≤ 0.01) between individual OTUs and measured environmental variables for the bacterial (A) and archaeal (B) datasets.


This study was financially supported by the European Community 7th Framework Project GENESIS (226536) on groundwater systems; the Cariplo Foundation under project No. 2011-1088 ‘SNAC- Synthetic and Natural Agrochemical Compounds: ecological impacts on the soil system and effects on plant production’; and the Doctoral School on the Agro-Food System (Agrisystem) of the Università Cattolica del Sacro Cuore (Italy). The authors are thankful for the critical evaluation of the manuscript by Dr Ioannes Stergiopoulos (Department of Plant Pathology of the University of California, Davis, CA) and Dr Eiko E. Kuramae (Microbial Ecology group of the Netherlands Institute of Ecology – NIOO-KNAW, Wageningnen, Netherlands). Moreover we would like to acknowledge the assistance during the study setup by Dr Stefano Longo and the technical support of Giulia Maria Parisi and Pierluisa Fantini.


  • Editor: Cindy Nakatsu


View Abstract