Phylogenetic and Multivariate Analyses To Determine the Effects of Different Tillage and Residue Management Practices on Soil Bacterial Communities

Characterization of the El Batán experimental station and the long-term field experiment.The study described here is part of long-term trial, which began in 1991. The research station is located near the former lake Texcoco, in the semiarid, subtropical highlands of central Mexico. The experimental design consisted of a randomized complete block with two replications, where treatments combined different wheat-maize rotations, tillage/planting methods, and residue management practices, but only four treatments were used in this research, which included the factor of disturbing soil (yes/no), referred to as conventional tillage (CT) and zero tillage (ZT), combined with the factor of residue management, i.e., removal (−R) or retention (+R) of crop residues. All four treatments were under rotation of maize and wheat. The sampling was done in the same season after the maize cropping cycle. Details of this experiment can be found in reference 22.

Soil sampling and characterization.Soil samples were taken from replicated plots at the end of the fallow period. Ten subsamples (0 to 15 cm) were collected at random from each plot with a 2-cm-diameter auger and mixed to yield one composite sample per treatment (7). Soil used for DNA extraction was stored immediately at −70°C to avoid modification of the composition of bacterial communities. Details of the techniques used to measure pH, electrolytic conductivity, total N and C, water-holding capacity, and soil microbial biomass (SMB) can be found in reference 36.

Soil DNA extraction.Total DNA was directly extracted from soil by using a hybrid technique based on methods described earlier (15, 55, 58). These techniques were modified to improve the purification of extracted DNA. Five grams of soil from the composite samples of each treatment was used for this procedure (43). Soil was added to a centrifuge tube containing 25 ml 0.15 M sodium pyrophosphate (24). The tube was shaken on a vortex for 1 min, left to settle for 10 min, and centrifuged at 7,700 × g for 10 min. The supernatant was decanted and discarded. This washing step was repeated three times. Five milliliters of 0.15 M phosphate buffer (0.15 M NaH₂PO₄, pH 8) was added to the tube containing the sediment, shaken on a vortex for 1 min, left to settle for 10 min, and centrifuged at 7,700 × g and 15°C for 10 min. The supernatant was decanted and discarded. This step was repeated once. The sediment was resuspended in 5 ml of lysis solution I (0.15 M NaCl, 0.1 M EDTA [pH 8.0], 10 mg of lysozyme ml⁻¹), mixed, and incubated at 37°C for 1 h. Five milliliters of lysis solution II (0.1 M NaCl, 0.5 M Tris-HCl [pH 8.0], 12% sodium dodecyl sulfate) was added, and the solution was treated with two cycles of freezing at −70°C for 20 min and thawing at 65°C for 20 min. Two milliliters of 0.15 M aluminum ammonium sulfate [AlNH₄(SO₄)₂] was added in the first thawing cycle (6). The sample was centrifuged at 7,700 × g and 15°C for 10 min. The supernatant was transferred to a clean tube and resuspended with 2.7 ml 5 M NaCl and 2.1 ml 10% Triton X-100 in 0.7 M NaCl. The tube was shaken and incubated at 65°C for 10 min. Twelve milliliters of chloroform-isoamyl alcohol (CH₃Cl-C₆H₁₂O, 24:1) was added, mixed, and centrifuged at 3,000 × g for 20 min at 15°C. The solution was transferred to a clean tube, 12 ml 13% polyethylene glycol (MW, 8,000) dissolved in 1.6 M NaCl was added, and the mixture was incubated on ice overnight and centrifuged at 12,000 × g and 4°C for 30 min. The sediment and DNA extract were resuspended in 500 μl of sterile deionized water and transferred to a sterile microcentrifuge tube. One volume of absolute ethanol was added, and the mixture was incubated at 4°C for 30 min. The tube was mixed and centrifuged at 10, 000 × g and 4°C for 30 min. The sediment was washed with 1,000 μl of cool 70% ethanol and dried at 65°C. The sediment was resuspended with 300 μl of sterile deionized water and stored at −20°C.

PCR amplification, cloning, and sequencing of 16S rRNA genes.The universal bacterial primers 46F (5′ GCC TAA CAC ATG CAA GTC 3′) and 1540R (5′ GGT TAC CTT GTT ACG ACT T 3′) were used to amplify highly variable regions V1 to V9 of 16S rRNA genes (14, 57). The PCR products were inserted into the pCR2.1 vector using the TOPO TA cloning kit (Invitrogen, Carlsbad, CA) and then transferred into Escherichia coli DH5α (Invitrogen, Carlsbad, CA). Sequencing reactions were performed by using Big Dye Terminator version 3.1 and an ABI Prism model 3100 DNA sequencer (Applied Biosystems). The specific conditions for these procedures can be found in reference 10.

Phylogenetic analysis.The sizes of the sequences obtained were 1,513 to 1,537 bp. All sequences were compared with reference sequences by BLAST search (2). The presence of chimeric sequences was checked with the RDP Chimera Check program. Multiple cycles of alignment based on phenetic procedures were performed with Clustal X version 1.7 to establish the best alignment conditions (54). The following settings were used: slow/accurate; gap-opening penalty, 15.00; gap extension penalty, 6.66; delay divergent sequences, 30%; DNA transition weight, 0.50; DNA weight matrix, IUB; negative matrix, off.

Phylogenetic relationships between sequences were determined by two methods, distance and maximum parsimony using PAUP* 4.0b10 software (52). The MODELTEST program (39) was used to select the best nucleotide substitution model. Distance matrixes were generated, and branch support was obtained with 100 bootstrap replicates.

For maximum-parsimony analyses, only the parsimonious informative characters were used, i.e., 820 bp for ZT/−R, 850 bp for ZT/+R, 867 bp for CT/−R, and 795 bp for CT/+R. Heuristic tree searches were done using a tree bisection reconnection model and a branch-swapping algorithm with 100 random stepwise replications. One hundred trees were saved for each pseudoreplicate. A rescaled consistency index, derived from trees obtained by unweighted analysis, was used to generate an a posteriori weighted data set. The same heuristic search conditions as used for the unweighted data were used to analyze the weighted data set. Branch support was obtained with 100 bootstrap replicates. Desulfurobacterium thermolithotrophum was used as the external group in every phylogenetic reconstruction.

Taxonomic assignment was obtained by using the Roselló-Mora and Amann prokaryote criteria (45). Phylogenetic groups were conformed at the order level for multivariate analysis.

Community richness and composition analysis.The distance matrixes (generated with the PAUP* 4.0b10 software) were used to obtain the operational taxonomic units (OTUs) for each library. A 3% distance level between sequences was considered the cutoff among different OTUs. Rarefaction curve, richness estimator (abundance-based coverage estimator, bootstrap richness estimator, bias-corrected Chao1, and diversity indices [Shannon's H and Simpson's D]) were determined using DOTUR software version 1.51 (47) for each treatment.

LIBSHUFF analysis was used to compute the differences between the clone libraries and estimated homologous and heterologous coverage within and between libraries from different treatments (48).

AMMI model to analyze the responses of the phylogenetic groups across agricultural treatments.The AMMI model (18, 60) was used to understand the responses of the 19 phylogenetic groups to different agricultural management practices represented by agricultural treatments ZT/+R, ZT/−R, CT/+R, and CT/+R. The AMMI model is usually used in plant breeding for studying genotype response patterns across different environmental conditions (11). In this case, the response variable was the relative proportion of the different phylogenetic groups in the four agricultural treatments and it can be represented as ȳ_ij = μ + τ_i + δ_j + Σ_k^t= 1λ_kα_ikγ_jk + ε̄_ij, where μ is the grand mean over all genotypes and environments, τ_i is the additive effect of the ith genotype, δ_j is the additive effect of the jth environment, the constant λ_k is the singular value of the kth multiplicative component that is ordered λ₁ ≥ λ₂ ≥ … ≥ λ_j, and α_ik are elements of the kth left singular vector of the true interaction and represent the genotypic sensitivity to hypothetical environmental factors represented by the kth right singular vector with elements γ_jk. α_ik and γ_jk satisfy the orthonormalization constraints Σ_iα_ikα_ik′ = Σ_iγ_jkγ_jk′ = 0 for k ≠ k′ and Σ_iα_ik²= Σ_jγ_jk²= 1. This model was called the AMMI model by Zobel et al. (60) and Gauch (18). Gabriel (17) described the least-squares fit and explained how (i) the residual matrix of the genotype-environment interaction (GEI) term, Z = ȳ_ij−ȳ_i.−ȳ_.j−ȳ_.., is subjected to a singular-value decomposition after adjustment for the additive (linear) terms and (ii) the results can be interpreted in a graph called a biplot which reflects, in a reduced-dimensional graph, the most relevant information in the genotype-environment matrix.

Representation of the biplot obtained by AMMI analysis.The biplot is a graphical representation of the first two components of the AMMI model, that is, α₁, α₂, γ₁, and γ₂. While other multivariate analysis methods, such as a conventional principal-component analysis (PCA), are applied to the raw data, the AMMI model calculates the first two components by applying the additive analysis of variance to the raw data to determine the residuals of interaction (GEI term), which are then subjected to a PCA to represent, in a low dimension, the GEI effects of the different agricultural treatments on the phylogenetic groups (18).

The results of the analysis of interactions between bacterial (phylogenetic) groups and agricultural treatments (environments) are represented as a series of vectors in a two-dimensional space drawn from the origin (axis 1 = 0, axis 2 = 0) to the endpoints determined by their scores, i.e., contribution to data variation. The agricultural treatments located toward the center of the biplot (axis 1 = 0, axis 2 = 0) do not discriminate the bacterial groups, which means that the groups have the same response in terms of GEI in those agricultural treatments, whereas agricultural treatments far away from the center (axis 1 = 0, axis 2 = 0) discriminate the phylogenetic groups differently. In the same way, bacterial groups whose presence showed a negligible response to soil management were located toward the graph center (axis 1 = 0, axis 2 = 0). Conversely, bacterial groups located far from the center (axis 1 = 0, axis 2 = 0) had a noticeable response to agricultural practices.

Bacterial groups located in the same direction of a particular agricultural practice were positively affected (increased their proportion) by the corresponding practice, whereas those bacterial groups located in the opposite direction of the vector representing the agricultural practice were negatively affected by the conditions prevailing in that environment. An angle of less than 90° or greater than 270° between a phylogenetic group vector and an agricultural treatment vector indicated that the phylogenetic group had a positive response to the treatment, i.e., an increase in its relative proportion within the population. Conversely, a negative response is indicated if the angle between the vectors is between 90° and 270°.

Nucleotide sequence accession numbers.The 16S rRNA gene sequences obtained in this study have been deposited in the GenBank database under accession numbers EU440544 to EU440729 , FJ889182 to FJ889353 , EU665003 to EU665174 , and EU449555 to EU449738 .

Physical-chemical soil characterization.Initial characterization of the soil at the experimental site showed characteristics similar to those normally found in an agricultural soil of the central highlands of Mexico (Table 1). The SMB, organic C, and total N were greater in soils with retained residues (ZT/+R, CT/+R) than in soils where residues were removed (ZT/−R, CT/−R).

Phylogenetic reconstructions.One hundred seventy-two clones per treatment were used for phylogenetic reconstructions using both the distance and maximum-parsimony methods. Since similar tree topologies were obtained with both methods, only results from the maximum-parsimony analysis will be discussed here.

The results of the phylogenetic analyses (see Fig. S1 to S4 in the supplemental material) show the phylogenetic relationship between the sequences amplified from each soil treatment and 16S rRNA gene GenBank sequences. While we were able to assign numerous sequences at the species level (see Fig. S1 to S4 in the supplemental material), phylogenetic groups were constructed only at the order level to reduce the number of bacterial groups for statistical analysis. The number, description, and relative proportions of groups by treatment are shown in Table 2. The most abundant groups were Acidobacteriales, with relative proportions of 0.43 to 0.54; Sphingomonadales, with relative proportions of 0.14 to 0.18; Myxococcales, with relative proportions of 0.058 to 0.814; and Xanthomonadales, with relative proportions of 0.035 to 0.098. The bacterial groups found in all treatments with lower relative proportions were Actinomycetales, with relative proportions of 0.011 to 0.040; Burkholderiales, with relative proportions of 0.011 to 0.029; Gemmatimonadales, with relative proportions of 0.005 to 0.058; and Oceanospirillales, with relative proportions of 0.005 to 0.017. The Pseudomonadales group, of which the main representatives were related to fluorescent Pseudomonas spp., was found in higher proportions in soils with retained residue and zero or conventional tillage (ZT/+R and CT/+R). Conversely, the highest proportions of Acidobacteriales and Oceanospirillales group members were found in those soils with removed residues and zero or conventional tillage (ZT/−R and CT/−R).

Composition, diversity, and richness analyses.LIBSHUFF analysis of the sequence libraries generated for each treatment indicated that the differences in their bacterial composition were due to the distinct makeup of the communities (P < 0.001). Rarefaction analysis showed that the treatments presented different degrees of diversity (Fig. 1). A decline in the rate of OTUs from rarefaction curves was observed only at 10 and 20% cutoffs (data not shown), indicating that at the family and order levels, the major bacterial groups were detected. Analysis with the DOTUR software also showed that ZT/+R soils had greater richness and diversity than CT/+R soils and those with removed residues (ZT/−R and CT/−R) (Table 3).

• Open in new tab
• Download powerpoint

FIG. 1.

Rarefaction analysis of 16S rRNA gene sequences from the four treatments: ZT/−R, ZT/+R, CT/−R, and CT/+R. The total number of sequences (172 per treatment) is plotted against unique OTUs defined by using a distance level of 3% applying the furthest-neighbor assignment algorithm in DOTUR.

AMMI analysis: interactions between bacterial groups and management practices.Figure 2 shows the relationship between bacterial groups and management practices according to AMMI analysis. The ZT/+R, ZT/−R, and CT/+R treatments had the most important effect on phylogenetic groups, while a discriminable effect was revealed by CT/−R, as shown by its lower vector magnitude. The phylogenetic groups that showed a noticeable response to the soil treatment were Acidobacteriales, Actinomycetales, Burkholderiales, Caulobacteriales, Gemmatimonadales, Pseudomonadales, Rhizobiales, Rhodospirillales, Sphingomonadales, and Xanthomonadales. For example, treatment ZT/+R had a positive effect on the presence of Pseudomonadales, Xanthomonadales, Rhizobiales, Gemmatimonadales, Sphingomonadales, and Rhodospirillales and a negative effect on the presence of Acidobacteriales and Actinomycetales (Fig. 2). Similarly, treatment CT/+R had a positive effect on the presence of Burkholderiales, Caulobacteriales, Pseudomonadales, Sphingomonadales, and Acidobacteriales and a negative effect on the presence of Rhizobiales, Gemmatimonadales, and Rhodospirillales (Fig. 2). Treatments with zero tillage and residue removal (ZT/−R) had a positive effect on the presence of Rhodospirillales, Gemmatimonadales, Rhizobiales, Actinomycetales, and Acidobacteriales and reduced the abundance of Xanthomonadales, Pseudomonadales, Sphingomonadales, Burkholderiales, and Caulobacteriales. Treatment CT/−R had a negative effect on the presence of Acidobacteriales and Actinomycetales. Other phylogenetic groups, i.e., those that were grouped near the center of the graphs, such as Bacillales, Caldilineales, Chromatiales, Clostridiales, Legionellales, Myxococcales, Oceanospirillales, Rubrobacteriales, and unclassified Gammaproteobacteria, were less affected by management practices. These analyses support the conclusions of previous studies that agricultural management practices can affect the bacterial composition of soil.

• Open in new tab
• Download powerpoint

FIG. 2.

Biplot of the two-way interactions of 19 phylogenetic groups (numbered 1 to 19 as in Table 2) and four treatments (ZT/−R, ZT/+R, CT/−R, and CT/+R) associated with the relative proportions of the phylogenetic groups.