An assessment of the hypervariable domains of the 16S rRNA genes for their value in determining microbial community diversity: the paradox of traditional ecological indices

Abstract

Amplicon length heterogeneity PCR (LH-PCR) was investigated for its ability to distinguish between microbial community patterns from the same soil type under different land management practices. Natural sagebrush and irrigated mouldboard-ploughed soils from Idaho were queried as to which hypervariable domains, or combinations of 16S rRNA gene domains, were the best molecular markers. Using standard ecological indices to measure richness, diversity and evenness, the combination of three domains, V1, V3 and V1+V2, or the combined V1 and V3 domains were the markers that could best distinguish the undisturbed natural sagebrush communities from the mouldboard-ploughed microbial communities. Bray–Curtis similarity and multidimensional scaling were found to be better metrics to ordinate and cluster the LH-PCR community profiling data. The use/misuse of traditional ecological indices such as diversity and evenness to study microbial community profiles will remain a major point to consider when performing metagenomic studies.

Introduction

Molecular microbial ecologists often use DNA profiling techniques and the prokaryotic ribosomal genes as phylogenetic markers to assess community structural ‘diversity’ quickly (Torsvik, 1996; Hill, 2003). Among the various molecular profiling techniques (Diez, 2001; Fennell, 2004), amplicon length heterogeneity PCR (LH-PCR) is a method that queries the hypervariable domains of the 16S rRNA genes to identify structural patterns in microbial communities (Suzuki, 1998; Bernhard & Field, 2000; Litchfield & Gillevet, 2002; Mills, 2003; Tiirola, 2003; Bernhard, 2005). In addition, the data are phylogenetically relevant in that the amplicon lengths (i.e. the natural length variation within sequences) generated can be directly associated with specific taxonomic sequences archived in the databases using common sequence alignment and analysis tools (Suzuki, 1998; Crosby & Criddle, 2003).

Once a community pattern is obtained, microbial ecologists have often analysed the data using traditional ecological indices (Watve & Gangal, 1996; Dunbar, 2000; Osborn, 2000; Hill, 2003). These indices are based on the clear definition and ecological description of an individual species, an entity that is often difficult to define in microbiology (Istock, 1996; Watve & Gangal, 1996; Hughes, 2001). Thus, microbial ecologists are often faced with the dilemma of deciding if these traditional indices, designed for discrete macro-community analyses, are appropriate measures to apply to microbial community profiles.

This study investigated soil bacterial communities using LH-PCR molecular profiling. Previous studies have shown that soil type and its inherent resources determine microbial community structure (Zhou, 2002), and that soil degradation or manipulation (e.g. mouldboard ploughing) can negatively impact soil biodiversity (Girvan, 2003). The objectives of this study were: (1) to interrogate three of the nine hypervariable domains of the bacterial small ribosomal subunit (16S rRNA) genes using LH-PCR; (2) to assess which domain or combination of domains provides the highest discrimination; and (3) to compare the traditional ecological indices of richness, diversity and evenness to non-metric multidimensional scaling in the study of microbial communities. To date, LH-PCR studies (Suzuki, 1998; Bernhard & Field, 2000; Mills, 2000; Ritchie, 2000; Tiirola, 2002; Bernhard, 2005) have used only the first two or three hypervariable regions of the 16S rRNA genes. In addition to amplifying two hypervariable domains inclusively (V1+V2), we utilized two other sets of universal primers to amplify the hypervariable domains V1and V3 exclusively. Therefore each combined profile was a concatenation of data obtained using different pairs of primers on the same samples.

Materials and methods

Soil samples

Soil samples were collected at the USDA Agricultural Research Service's Northwest Irrigation and Soils Research Laboratory, Kimberly, ID. All sites are classified as a temperate high, semidesert ecosystem. Soil on all sites has been classified as a coarse-silty, mixed, superactive, mesic durinodic xeric haplocalcid, with 0.1–0.21 grams of clay per gram of soil clay and 0.6–0.75 grams of silt per gram of soil silt, and organic matter of approximately 13 g per kg of soil. The soil had a pH of 7.6–8.0. Plots were 30 m × 30 m and the experiments were arranged in a completely randomized design. Soil samples were taken in triplicate at various depths (0–5 cm, 5–15 cm and 15–30 cm) from three native sagebrush (NSB) sites and three irrigated agricultural cropland sites under mouldboard ploughing systems (IMP) (depth: 0–30 cm). All samples were transported on ice to the laboratory, where they were homogenized, and aliquots frozen at −80°C until processed for molecular analysis.

DNA extraction, quantification and LH-PCR

Five hundred milligrams of homogenized soil samples were extracted in triplicate using the Qbiogene BIO 101 FastPrep^® instrument and FastDNA^® SPIN kit for soil (QBiogene, Irvine, CA) with slight modifications to the manufacturer's protocols previously described (see Mills, 2003). The metagenomic DNA was quantified using a Hoefer DNA Quant fluorometer (Hoefer, San Francisco, CA) and diluted to 10 ng μL⁻¹ working stocks.

For the LH-PCR reactions only the forward primer of each set was labelled with a fluorescent dye. Primers for the V1 domain were P1F_6-FAM 5′ GCG GCG TGC CTA ATA CAT GC 3′; P1R 5′ TTC CCC ACG CGT TAC TCA CC 3′; and for the V3 domain were 338F_HEX 5′ ACT CCT ACG GGA GGC AGC AG 3′ (Cocolin, 2001). The V1+V2 domain primers were 27F_6-FAM 5′ AGA GTT TGA TCM TGG CTC AG 3′; 355R 5′ GCT GCC TCC CGT AGG AGT 3′ (Suzuki, 1998). Final concentrations for all PCR reaction mixes were: 1 × PCR buffer, 2.5 mM MgCl₂, 0.25 mM each dNTP (Promega, Madison, WI), 0.5 μM forward and reverse primers (Invitrogen Life Technologies, Carlsbad, CA) 0.25 U AmpliTaq Gold^® DNA polymerase (Applied Biosystems, Foster City, CA), 0.1% (weight in volume) bovine serum albumin [BSA, fraction V, non-acetylated (Fisher Scientific, Pittsburgh, PA)], 10 ng of genomic DNA and DEPC-treated water (Sigma, St Louis, MO) to a final volume of 20 μL. The MJ DNA Engine PTC-200 programmable thermocycler (MJ Research, Waltham, MA) was used and PCR cycling conditions were as follows: one initial denaturing step of 95°C for 11 min followed by 25 cycles of 95°C for 1 min, annealing at 55°C for 1 min, and extension at 72°C for 1 min. The final extension step was for 72°C for 10 min. These parameters had been previously optimized (see Mills, 2003) so as to minimize any PCR template or kinetic biases that are inherent to PCR and amplification of metagenomic DNA samples (Suzuki & Giovannoni, 1996; Suzuki, 1998).

Electrophoresis

A master mix of a 5 : 2 : 1 ratio of deionized formamide (98%, Sigma), ABI Blue Dextran-EDTA loading dye and ABI GeneScan™ 500 ROX internal standard (Applied Biosystems) was used to denature LH-PCR products before loading on the ABI 377 genetic analyser (Applied Biosystems). The 3 : 2 mixture of loading buffer to PCR product was denatured at 95°C for 4 min, snap-cooled, and held on ice until loading 0.75 μL of product onto a 36-inch well-to-read (WTR) 5% long-range polyacrylamide gel (FMC, Philadelphia, PA). Standard run conditions, using Filter D were used for a 3.5-h run. GeneScan^® and Genotyper^® software (Applied Biosystems) were used to collect and analyse the data.

Binning and analysis

Binning and normalization were performed in order to eliminate errors that are inherent to the collection and analysis software (see Mills, 2003). Raw data for each subplot and domain(s) were averaged and the relative ratios of peak heights (intensities) were calculated by dividing each individual peak height in each electropherogram by the total intensity. The mean relative ratios were used in all subsequent analyses as described in (Dunbar, 1999) and for the multidimensional-scaling (MDS) plots. The traditional ecological indices used were

1. 1

   Shannon information index

   H′=−Σp_i (ln p_i), where p_i is the proportion of the individuals in the ith species.

2. 2

   Evenness index

   E=H′/H_max where H_max=ln (S), where S is the richness or number of amplicons.

anovas were used to compare ecological indices among the four groups (SPSS ver. 13.0, SPSS Inc., Chicago, IL). When significant differences were observed, Bonferonni post-hoc comparisons were performed and all confidence intervals were set at 95%.

The ratio data were imported into Primer 5 (Primer-E, Ltd, http://www.primer-e.com) and square-root-transformed, and Bray–Curtis similarity matrices were calculated using

1. 3

   Bray–Curtis similarity

   Cz=2w/(a+b), where a is the sum of abundances of all amplicons found in a given sample; b is the sum of amplicon abundances in another sample; and w is the sum of the lower of the abundance values for each amplicon common to both samples (Pohle & Thomas, 2001).

To visualize similarity among samples, non-metric MDS was employed using the Bray–Curtis similarity matrix. Analysis of similarity (ANOSIM) was performed for comparisons of similarity between groups, and similarity percentage (SIMPER) analysis was used to determine which amplicons were responsible for discriminating between treatments.

For calculating the concatenated profiles, the hypervariable-domain raw data were realigned, combined and treated as one data set in various domain combinations. Relative ratios and binning were recalculated for each analysis, as described above in order for the total intensity to equal one; the relative ratios for each peak changed because the sum of the total intensities of the combined domains changed depending on the domains queried. Peaks that fell below 0.01 when combined with another data set were eliminated. Thus, with some domain combinations, there was an attenuation of the number of data points used in the analysis.

Results and discussion

Ecological indices and LH-PCR domains

Ecological indices were calculated for each separate domain profile, i.e. V1, V1+V2, and V3 hypervariable domains (Table 1). There was significant variation among groups in all domains (see anova results, Table 1). Using only the V1 domain data, diversity and evenness for the NSB samples at all depths were significantly different from those of the IMP samples (here and subsequently in this section, ‘significant differences’ between treatments refers to P<0.05 for Bonferroni post-hoc comparisons). For V1+V2 data, NSB 5–15 and NSB 15–30 were not significantly different when compared with each other but were distinct when compared with the other samples. When querying the V3 domain, the diversity indices were not significant.

The three data sets, V1, V3 and V1+V2, were combined and the indices recalculated (Table 2). IMP was significantly different from NSB 5–15 and NSB 15–30 with respect to diversity. All NSB sample evenness indices were not significantly different from each other, but IMP was significantly different from NSB 15–30. Indices were then recalculated using only two domains. For the combined V1 and V3 domains, there were no significant differences in richness between samples, but IMP diversity and evenness were significantly different from all NSB samples. The same was true for the V1, V1+V2 data set. When V3 and V1+V2 were combined IMP was not significantly different from NSB 5–15 or NSB 15–30; IMP was, however, significantly different from NSB 0–5. Evenness did not differ between any of the samples.

Multidimensional scaling

Multidimensional scaling tightly clustered the NSB samples apart from IMP with V1 and the three-domain data (Figs 1a and 1e). Similarly, distinct NSB and IMP groups were present with V1 and V3 combined data (Fig. 1d). The V1 and V3 data also indicated the subtle differences related to depth. The V1+V2 or V3 domain data were able to distinguish NSB from IMP samples but the clustering was not as pronounced. All other domain combinations followed the same MDS clustering trend as the V1 and V3 data, with only slight variations in the spacing of the clusters within the three different depth NSB groups (data not shown).

ANOSIM is a measure of the dissimilarity of a priori defined groups. Global R-values near zero indicate no difference among groups, and R=1 or near one indicates that samples within groups are more similar to each other than samples from different groups. In this study, the global R-values for NSB groups compared with IMP for any combination of domains was R≥ 0.80, P<0.001, indicating that NSB always grouped distinctly from IMP. SIMPER analysis of the combined V1 and V3 domains (Table 3) shows the major amplicons that contributed to the dissimilarity among the soil groups. How consistently an amplicon was able to differentiate among the groups were indicated by the ratio of the average dissimilarity to the standard deviation (Dis:SD). A large ratio indicates that the amplicon consistently and substantially contributed to the differences between LH-PCR profiles.

Choice of diversity measures and hypervariable domains

A study by Hill (2003) addressed the issue of the suitability of ecological indices for microbial ecology studies. This report included a comprehensive review of the index selection issue using microbial community clone libraries and concluded that the Shannon index appeared to be more sensitive to overall community change and the loss of rare populations than did the Simpson dominance index (Hill, 2003). Since LH-PCR profiles reflect a combination of dominant and rare data points, the Shannon index was selected to measure diversity.

Unlike large clone libraries, however, microbial community profiles generated by any profiling method, including LH-PCR, only represent the minimum detectable diversity. When any of the (two or three) domain data were combined, NSB samples had higher diversity than IMP, but within NSB, based on depth, in general the indices did not significantly differ. The individual domain data did not, however, support this trend, and the interpretation differed depending on the domain queried. Careful consideration when selecting both which hypervariable domain(s) to query and which indices to use, with respect to their suitability and limitations, is required if diversity measures are to be applied appropriately to LH-PCR community profiles.

Other microbial ecology studies have addressed the issue of diversity using other molecular techniques and these same ecological indices. Dunbar. (2000) concluded in their soil study using terminal restriction fragment length polymorphism (TRFLP) analyses that the use of multiple restriction enzymes to generate community patterns did not necessarily increase the overall resolution of diversity. Instead, multiple digests were used to increase the confidence that similarities between samples were not a result of technical biases that may come from using only one enzyme (Dunbar, 1999). In addition, diversity indices were calculated individually for each enzymatic pattern. They observed that diversity indices varied based on the restriction enzyme used. Likewise, in our study, the choice of domain queried had implications for diversity measures. Unlike Dunbar, we observed an increase in discrimination of diversity using the combined domain data.

Similarity and MDS

Patterning the similarity/dissimilarity of the community proved to be more informative than the traditional ecological indices in discriminating treatment groups. Bray–Curtis similarity is used principally with continuous numerical data and species abundance data similar to those generated by LH-PCR profiling (Gotelli & Ellsion, 2004). Clearly, in this study, the Bray–Curtis similarity and subsequent MDS ordinations strongly reflected between site dissimilarities. The advantage of using Bray–Curtis similarity over other similarity indices based on presence/absence calculations and/or Euclidean distances is that there is less bias introduced by shared absences of amplicons (Gotelli & Ellsion, 2004). The intrinsic importance of both dominant and rare community members within the community can still be reflected using this similarity index and ordination method.

In a recent study by Bernhard. (2005), LH-PCR data patterned the microbial communities across a salinity gradient. Using MDS, these authors were able to discriminate clearly between the bacterial communities associated with fresh water and those in estuarine and marine waters. These different aquatic systems drive both physiological and habitat adaptation, and were reflected in structural changes (i.e. amplicons) in the microbial communities. The present study differs from Bernhard's in that the same soil ecosystem (no physiological gradient per se) was being studied, and disturbance and depth were assumed to be the ecological drivers. It has been shown previously that microbial diversity changes, and is often lower with cropping, land management and use (Ibekwe, 2002) and with depth (Blume, 2002). Even though similar micro organisms were no doubt associated within a given Idaho soil type, LH-PCR and MDS analysis were better able to pattern the subtle community differences associated with disturbance and somewhat with depth. MDS strongly supported the clustering of the samples based on tillage. However, the discriminatory pattern produced was also a function of which domains were queried. Even though the diversity indices seemed to need only two domains (i.e. V1 combined with V3) to discriminate adequately between samples, the more robust scaling output was reflected with either the V1 or the three-domain data. However, V1 and V3 data were able to differentiate NSB from IMP samples and also to identify differences associated with depth in the NSB soil groups. Herein lies the paradox of choosing which domains and measures better represent community structural differences.

Conclusions

Univariate measures such as Shannon's information index and its associated evenness index appear to be more appropriate measures for community clone libraries since individual bacterial identification can be ascertained through sequencing (Wintzingerode, 1997; Hill, 2002). Multivariate and ordination measures such as Bray–Curtis similarity and MDS that use iterative procedures to map the similarity (or dissimilarity) produced by profiling were better measures to apply to this study. ANOSIM can test for differences between multivariate groups, and SIMPER analyses can rank the amplicons contributing the most to the dissimilarities. Because of the inherently lower resolution provided by any molecular profiling method, dynamic patterning of microbial communities can best be analysed using these nonparametric, multivariate methods.

The evolution of profiling techniques over the last few years has moved from publishing only the raw data (electropherograms or gels) to applying traditional ecological analyses to profiling data and now to using multivariate analyses. Perhaps the greatest challenge facing microbial ecologists will be to break away from the traditional ecological paradigms of diversity measures and resolve the paradox by developing new algorithms or unique metrics that will better analyse the complexities still hidden within microbial communities.

Acknowledgements

This study was funded in part by a NSF ADVANCE Fellowship (Award no. 0340695) to D.K.M. and in part by the USDA Agricultural Research Service, Northwest Irrigation and Soils Research Laboratory, Kimberly, ID to J.A.E. Lilliana Moreno, Laurie Richardson, Todd Crandall and the rest of the Mathee and Richardson laboratory members are thanked for their critical reviews and helpful comments.

References