New Threshold and Confidence Estimates for Terminal Restriction Fragment Length Polymorphism Analysis of Complex Bacterial Communities^†

Samples and DNA extraction.

Fifteen soil and two sediment samples, collected between June 1997 and November 2003, were frozen at −20°C within 24 h of sample collection, after which total community DNA was extracted using the protocol of O'Farrell and Janssen (27). The samples used were: A, marine sediment from Williamstown Hatt Reserve seagrass bed; B, marine sediment from below the Williamstown Anglers Club pier; C, rhizosphere soil from a eucalypt in Wombat State Forest; D, rhizosphere soil from a pine tree on Moonee Ponds Creek; E, soil from a vineyard at Bunyip; F, rhizosphere soil from a camellia bush in Ascot Vale; G, soil from a heavily manured vegetable garden in Bunyip; H, soil from a pasture (0 to 10 cm below ground) at Ellinbank; HL, soil from a pasture (3 to 10 cm below ground) at Ellinbank; HR, roots and associated soil from pasture grasses at Ellinbank; I, soil from a petrol- and oil-contaminated grassed area in Ascot Vale; J, soil from a native grassland in Sunbury; KB, soil from a forage brassica plot at Ginninderra; KW, soil from a wheat plot at Ginninderra; L0, rhizosphere soil from a eucalypt on Moonee Ponds Creek; L2, soil under grass (2 m away from sample L0) on Moonee Ponds Creek; and L4, soil under grass (4 m away from sample L0) on Moonee Ponds Creek. All but the two Ginninderra sites, which were located in the Australian Capital Territory, were in Victoria, Australia.

Generation of T-RFLP profiles and data sets.

T-RFLP profiles were generated using the primers FAM27f (5′-GAGTTTGATCMTGGCTCAG-3′), labeled at the 5′ terminus with 6-carboxyfluorescein, and 519r (5′-GWATTACCGCGGCKGCTG-3′) and otherwise followed the protocol of Sait et al. (35). Unincorporated primers and reaction components were removed using a Wizard SV Gel and PCR Clean-Up System column (Promega, Annandale, New South Wales, Australia) according to the manufacturer's instructions. Approximately 100 ng of purified PCR product was digested overnight, at the specified temperature, with 5 U of one of the following restriction endonucleases: BstUI, HaeIII, HhaI, HinfI, MspI, or Sau96I (New England Biolabs, Inc., Beverly, Mass.). Digested PCR products were then precipitated and analyzed on a model 377 DNA sequencer (Applied Biosystems, Foster City, Calif.) at the Australian Genome Research Facility, Parkville, Victoria, Australia (35) to generate profiles of fragments up to 530 nucleotides (nt) long.

A separate profile was generated for each sample with each of the six restriction endonucleases. Each profile consisted of T-RFs that each had a reported fragment length, in nucleotides, and a reported peak area of fluorescently labeled product, in fluorescence units (FU). The raw data set of fragment lengths and corresponding peak areas obtained from the Genotyper software (Applied Biosystems) was first compiled so that all T-RFs in each profile of the data set were aligned with the T-RFs of the same inferred length (rounded to the nearest nucleotide) in every other T-RFLP profile generated by that restriction endonuclease.

All of the profiles generated with any one RE constituted a data set. The fluorescence integrated under any one peak is referred to as the area of that peak, and the total area for any one profile is the sum of the areas of all of the peaks excluding those generated by fragments of less than 30 nt or greater than 500 nt. To standardize data sets, small peaks that may have been detected in a profile as a result of loading large amounts of DNA were removed by application of a threshold area. Peaks with an area smaller than this threshold area were removed from the data set, using three different standardization methods: the constant percentage threshold, the constant baseline threshold, and the variable percentage threshold.

Constant percentage threshold calculations.

Sait et al. (35) determined the threshold area for a data set by applying increasing area thresholds, as percentages of the total area on each trace, until the minimum percentage that resulted in independence of the number of peaks remaining in each profile and the total area of each profile before threshold application was found. This standardization method is referred to as the constant percentage threshold calculation. All peaks that contributed less than this constant percentage threshold to any profile within a data set were removed before analyzing relationships between the profiles in the data set. A constant percentage threshold value was determined for each data set (referred to as constant percentage [different]) and also by combining all the profiles in all six data sets into a global data set to calculate a global constant percentage (referred to as constant percentage [global]).

Constant baseline threshold calculations.

The second method for standardization was the constant baseline threshold calculation of Dunbar et al. (10), applied to peak areas. The total areas in all of the profiles in a data set were normalized to the same value as that of the profile having the smallest total area, and then all of the peaks in each trace were reduced proportionally by the factor required to yield that normalized total area. The first constant baseline threshold was set as the smallest peak area detected in the unmanipulated data sets (rounded up to 50 fluorescence units), and peaks with an area equal to or smaller than this were removed from profiles after normalization (referred to as constant baseline [50 FU]), before analyzing relationships between traces in the data set. A second analysis was performed on data sets where the constant baseline threshold was set at 100 fluorescence units (referred to as constant baseline [100 FU]).

Variable percentage threshold calculations.

The third method for standardization was termed the variable percentage threshold method. The total area of each profile was divided by different values (divisors) to yield numbers that were used as percentage thresholds. For each divisor, all peaks that contributed less than the percentage threshold calculated for that profile were removed. Then, for each divisor, the remaining number of peaks was plotted against the total area, so that each profile contributed one point on that plot. We have found that, sometimes, if there is very little variation in the total area of the profiles, there is no detectable relationship between the number of peaks and the total area, and no threshold needs to be applied (M. Sait and P. H. Janssen, unpublished data). If standardization is required, a useful divisor to start with is 1,000 times the mean total area for the data set. Different divisors were then tested, and the divisor that resulted in the weakest relationship between the number of peaks remaining and the initial total area was considered to be the optimal divisor (Fig. 1). The unique percentage threshold value for each profile was calculated by dividing the total area of that profile by the optimal divisor. Peaks that contributed less than that percentage threshold were removed from that profile before analyzing relationships between traces in the data set.

Comparison of profiles.

Matrices of Sorensen-Dice pairwise similarity coefficients (8, 40) were calculated for all possible comparisons of profiles within a data set using a program written in the C++ programming language. These similarities were converted to distances, where distance = 1 − similarity, to produce distance matrices. These matrices were represented graphically as dendrograms using the Fitch-Margoliash least-squares (14) and neighbor-joining (36) algorithms in PHYLIP (13). Consensus dendrograms were generated from the six data sets using the CONSENSE program in PHYLIP, and the consensus values reported are the numbers of REs that generated coherent clusters of samples. Shannon-Wiener diversity indices (41) were calculated using the peak areas, expressed as a proportion of the total area, as measures of abundance.

Multivariate statistical analysis.

Matrices of Sorensen-Dice similarity coefficients were generated for each T-RFLP data set from each RE and for each threshold calculation procedure. A first-stage nonmetric dimensional scaling (MDS) plot was then derived from each of the similarity matrices. The second-stage MDS (2sMDS) analysis was then used to examine similarities between each of the first-stage MDS plots. For the second-stage analysis, a Spearman rank correlation was computed between the corresponding elements of each pair of first-stage similarity matrices. The resulting Spearman rank correlations then became the elements of a second similarity matrix, which was used to generate the 2sMDS plot. Each point on the 2sMDS represents a first-stage MDS plot and therefore essentially provides a visual interpretation of the similarity between all of the original MDS plots.

We also calculated a relative dispersion index for each data manipulation procedure (5, 39). The dispersion index compares the average rank similarity across a given group of samples. In our case, we defined each of the data manipulations or REs as the sample groupings. For interpretation, the lower the dispersion index, the tighter is the grouping of samples within the given treatment.

All multivariate statistical routines were carried out with the PRIMER 5 software package (Primer-E Ltd., Plymouth, United Kingdom).

Experimental design.

In this study, we used only one method of determining the size of peaks (by area according to Genotyper software outputs), one method of comparing T-RFLP profiles (Sorensen-Dice similarity coefficients), and two methods of representing the manipulated data sets (dendrograms and 2sMDS). We tested how different ways of calculating the threshold area and how the choice and number of REs used can affect the conclusions drawn about the similarity of samples. Peaks with areas smaller than the threshold area are removed from the data set to produce a manipulated data set. This is required to standardize profiles so that differences due to sample loading, resulting in differences in sensitivity for different samples, do not affect comparisons based on the presence or absence of peaks.

DNA was extracted from the 17 samples, and PCR products were generated from all DNA samples. The pooled products from two PCRs from each DNA template were separated into six aliquots and digested using six different REs. To create T-RFLP profiles that were similar to some of these, additional PCRs were performed using three of the DNA samples (D, F, and G), and these were also aliquoted and digested with the same six REs. In this way, a total of 120 T-RFLP profiles were generated. We chose REs that spanned the range of the number of T-RFs of 50 to 500 nt predicted for a set of about 4,600 16S rRNA gene sequences and therefore should produce profiles that vary in the number of sequences theoretically contributing to any one peak (11).

To ensure that T-RFLP profiles more accurately represented the bacterial community of a sample and to reduce method-associated noise, it has been suggested that the products from multiple PCRs from the same template be pooled to minimize the effects of amplification bias in individual PCRs (6) and that consensus data sets be generated from multiple electrophoretic separations of the digestion products from one sample and RE to minimize variations due to profile generation (10). However, for this study we relied on variations in PCRs and electropherograms to generate less-than-ideal data sets. We reasoned that such data sets would allow us to test the effects of thresholds and choices of REs, since we knew that the pairs of replicates (D and D′, F and F′, and G and G′) were derived from the same communities. Analyses that resulted in these replicates being deemed to be similar can therefore be considered more useful than analyses that failed to recognize replicates as being similar. This avoided a priori assumptions of sample similarity to test analysis methods, which may not be valid (15).

Similarities between raw T-RFLP profiles.

The twenty profiles (17 samples and 3 extra replicates) generated with each RE were compiled into one data set per RE. The similarity of each pair of profiles within each data set was expressed as a Sorensen-Dice pairwise similarity coefficient. These values were then used to generate distances, and dendrograms were constructed from these distance values to depict the similarities between different profiles. Profiles that were similar to each other grouped closely in these dendrograms. The samples and replicates grouped differently in the analyses of different raw data sets, so that few groupings were consistently found in analyses of all six data sets, i.e., the use of different REs generally resulted in different pairs of profiles being deemed to be similar. For example, comparison of the data sets generated using BstUI and MspI reveals that few groupings are shared between the two (Fig. 2A and B). The two sediment samples (A and B), the two Ginninderra agricultural soil samples (KB and KW), and the two parkland soil samples separated by 2 m (L2 and L4) did group together in the analyses of both data sets. In fact, these three pairings were the only ones found consistently in the analyses with all six REs (Fig. 3A).

Two of the replicate PCR pairs, D and D′ and F and F′, were paired only by four of the six REs, while the third replicate pair, G and G′, paired in only one data set (Fig. 3A). The only RE to group all three replicates in analyses of the raw data sets was HaeIII. D and D′ were not grouped by BstUI and HhaI, but F and F′ were not grouped by BstUI and Sau96I. The different groupings of replicates meant that it was difficult to be sure whether some of the groupings formed by nonreplicate samples were meaningful, and, even with HaeIII, it could not be ruled out that some of the pairings of nonreplicates detected might be fortuitous rather than reflective of real similarities between the original bacterial communities.

Effects of different threshold calculations.

We attempted to standardize the data sets and so produce greater similarity between the replicates. Three different methods were used to calculate a threshold area value to remove minor peaks that may have been detected purely as a result of the amount of DNA applied to the separation gels. A constant percentage (different) threshold (35) was applied to each data set from each individual RE. In addition, this method was also used to calculate a constant percentage (global) threshold area for all of the data sets combined into one data set, so that one threshold area was applied to all six data sets. Sait et al. (35) expressed the area under each peak as a percentage of the total area for all of the peaks in that profile and determined the minimum threshold as a percentage for the data set, so the profiles have comparable numbers of peaks while still retaining enough peaks to analyze the microbial community. Constant baseline thresholds, based on the method of Dunbar et al. (10), were applied to each data set from each RE, using the threshold areas of 50 and 100 FU. Dunbar et al. (10) actually compared peak heights, rather than area, but their method of proportionally reducing the heights of each peak in larger profiles and removing peaks that then fell below the threshold can also be used on peak areas. The third method combined elements of both of these published methods, and we term this the variable percentage threshold method. A unique threshold area, as a percentage of the total area, was calculated for each trace in a way that is dependent on the total peak area for each individual profile. A greater total peak area in any profile will increase the threshold area value of the baseline for that profile. This compensates for the increased sensitivity due to more labeled product being present in that profile. The constant percentage and variable percentage methods require enough samples to allow a reasonable estimate of the relationship between total peak area and the number of peaks in each profile to be determined (Fig. 1).

The variable percentage threshold method appeared to be the best for determining threshold values with these samples because it grouped all three pairs of PCR replicates (D and D′, F and F′, and G and G′) in all six data sets generated with different REs when the Sorensen-Dice similarity coefficient was represented by the Fitch-Margoliash least-squares treeing method (Fig. 3B and Table 1). The other methods of threshold determination resulted in the replicates grouping together less often. When the groupings were analyzed using the neighbor-joining method of depicting the outcomes as trees, the constant percentage (global) method resulted in the most consistent grouping of replicates, but some of the other methods also appeared satisfactory (Table 1). The variable percentage threshold method also resulted in the largest number of REs generating similar groupings between profiles generated from different samples of all the methods used (Table 1). Groupings uncovered in analyses generated using a number of different REs are more likely to consist of truly similar samples than groupings found using only one RE, indicating that this threshold method may be better for standardizing data sets and subsequently detecting true similarities between complex communities.

MDS is a robust approach to examining community data that has been used extensively to analyze data sets based on macrofaunal community structures (5) and, more recently, to analyze microbial community structures (32, 43). Here we used an extended MDS approach, termed 2sMDS, as an additional method to examine T-RFLP data for each RE and for each threshold calculation procedure. To date, 2sMDS analysis has been applied only to the study of marine benthic invertebrates (28, 29, 38) but has proven to be very useful in allowing metascale pattern analysis. Applying this technique to our data allowed us to place our results within a robust statistical framework. Each data point on second-stage ordination essentially represents an MDS ordination plot for one of the REs and one of the standardization procedures, and the position of each point is determined by the similarity of that MDS plot to all other plots. The most consistent grouping is indicated in such plots by the smallest spread of points within a given treatment.

The variable percentage threshold method gave rise to the most consistent grouping of the samples (Fig. 4A) and therefore represents the preferred approach to analyzing data. That is, the relationship of the samples to each other as determined by the six REs was most similar when the variable percentage threshold was applied to standardize the data. The relative dispersion index, a measure of the closeness of grouping, was smallest for the variable percentage threshold method (Table 1). At the other extreme, the constant percentage (different) threshold method (where each data set had its own unique threshold value determined) gave rise to the least consistent grouping of the samples (Fig. 4A) and had the largest relative dispersion index (Table 1). These results are generally consistent with the recovery of groupings of pairs of samples in consensus dendrograms (Table 1). The constant percentage (different) method resulted in a large loss of T-RFs (Table 1) (see the supplemental material), and this loss of information, and consequent simplification of the profiles, seems to have resulted in less-accurate detection of similarities between communities because the analyses are made on a small number of characters (T-RFs). Since this method was nearly as poor as the use of the raw data with no threshold, this suggests that overmanipulation of the data must be avoided. The three methods that resulted in the smallest dispersion all resulted in intermediate numbers of peaks (mean of 54 or 55 peaks) remaining after application of thresholds.

To test the effect of suboptimal threshold values on the data, and so the effects of under- and overmanipulating the data sets, we applied different thresholds to the data sets and compared the effects on the relative dispersion indices in 2sMDS. We applied threshold values at 0.5, 1.0, 1.5, and 2.0 times the optimal value that had been calculated using the variable percentage threshold method. As expected, the number of peaks decreased more rapidly than remaining total area as the threshold value increased (Fig. 5A and B), since application of the threshold eliminated the smallest peaks. The dispersion index decreased as the threshold value increased (Fig. 5C), but the consensus of agreement between the six different REs was greatest at the optimal threshold. In this instance, the dispersion indices are showing one outcome of steadily removing peaks from a given data set. As rare or small peaks are steadily removed from the data sets, the resultant data sets are becoming more similar; hence, there is less scatter with the ordinations, seen here as smaller dispersion indices. The optimum threshold identified by the variable percentage threshold method therefore results in the highest level of agreement between the different REs. Since we can expect that each RE will tend to produce similar patterns when the starting community is similar but also that experimentally induced noise will tend to obscure the true pattern of relationships between samples, we can assume that the highest consensus values occur when the optimum threshold is applied and a balance between noise elimination and information retention is achieved.

These outcomes suggest that the use of any method to calculate a threshold area allows the subsequent analysis to detect real groupings by eliminating noise. Overall, the variable percentage method appeared to be the most useful of the methods tested, probably because of its ability to set an optimal threshold area value.

Choice of restriction endonuclease.

The profiles generated by the different REs contained similar numbers of T-RFs, and measures of diversity also suggested that the amounts of information in the profiles were not greatly different when they were generated using different REs (Table 2) (see the supplemental material). However, the 2sMDS ordinations plotted by RE showed that individual enzymes can generate different community patterns (Fig. 4B), as evident from the scatter of the enzyme clustering throughout the 2sMDS (Table 2). It is clear that data sets generated using some enzymes, for example BstUI, are more strongly affected by the choice of threshold method than others, such as those generated using HaeIII (Fig. 4B). It is important to note that REs that resulted in more consistent groupings, like HaeIII and HinfI, did not generate profiles that looked more similar to each other (Table 2); instead, the relationships between samples were more consistent when these REs were used, regardless of which standardization method was used (Fig. 4B). This suggests that data sets generated using some enzymes are less susceptible to noise. In contrast, the choice of standardization method has a large impact on the outcomes of analyses of data generated using other REs. REs that generated profiles in which many small peaks were common to most profiles and large peaks were mainly unique should result in more differentiation of profiles when increasing thresholds were applied, while REs that generated many large peaks that were common to most profiles should be less sensitive to increasing thresholds. However, there were no obvious differences in the distribution of peak sizes in profiles generated using different REs (see the supplemental material). Instead, the distribution of the number of larger peaks relative to smaller peaks was almost identical for profiles generated with each of the six REs. When the six REs are ranked by relative dispersion index across all threshold methods (Table 2), there is no correlation with the characteristics listed for the same REs by Engebretson and Moyer (11). Other studies report sometimes-contradictory conclusions regarding the best REs in experimental comparisons (3, 6, 10, 20, 22, 34). We suggest that the analysis of data generated using multiple REs reduces the likelihood of a lack of resolving power due to RE selection.

Generation of confidence values using multiple REs.

Computer simulations of T-RFLP analyses have shown that the choice and number of REs used affects the profile's representation of the community, and the use of multiple REs was suggested to allow a more accurate assignment of peaks to 16S rRNA gene sequences (11). We have extended this by generating multiple data sets and constructing a consensus tree depicting relationships from the results of the individual data sets for six REs (Fig. 3). This allows an estimate of the confidence of different groupings to be made. Many studies employing T-RFLP use only one RE. Conclusions about the similarity of different samples could be attempted based on data sets generated from just one RE. For example, profiles of samples C and HL generated using BstUI appeared similar (Fig. 2C), allowing a conclusion to be made that their bacterial communities were similar. However, the consensus of the six separate data sets indicates that this is likely to be an incorrect conclusion (Fig. 3B), since these samples were found to be similar in only two of the six data sets (generated with BstUI and HinfI). This agrees with results found in a previous study (9), where the combined data from REs were not consistent, and indicates that conclusions about complex communities based on T-RFLPs generated from just one RE may be erroneous.

Each peak in a T-RFLP profile generated from a complex bacterial community such as those found in soils (11, 33) is likely to represent more than one species. The use of different REs will group these species together in different combinations (into different peaks), so that each RE produces a different simplified representation of the community. DNA from different soils may generate the same-sized T-RF that originates from different sets of 16S rRNA genes. These sets can be expected to be distributed in different peaks if a different RE is used. In our analyses, samples that were grouped together by all six REs are therefore more likely to be truly similar, while the confidence in groupings found in data sets generated by fewer REs is lower. The use of multiple REs generates a consensus value, which estimates a confidence level of the grouping. This is different from using a bootstrap analysis to generate confidence values for groups. During a bootstrap analysis, a new data set is generated by sampling characters randomly, with replacement, so that the resulting bootstrapped data set is of the same size as the original, but with some characters removed and others duplicated (12). The random variation of the results from analyzing these bootstrapped data sets is considered to be typical of the variation that might arise when collecting new data sets. As such, it tests the size and robustness of the data set. Our confidence values express the actual consensus of multiple real data sets. This represents a useful means of assessing the significance of any groupings observed.

Although some studies applying T-RFLP to complex microbial communities have used a number of different REs to generate more information for sample comparisons (2-4, 6, 9, 10, 26, 34), this is the first study to generate consensus values from the data generated using multiple REs. Ayala-del-Río et al. (1) initially tested a subset of samples with six different REs, but carried out community analysis on profiles from two REs, which were selected as the best. Within this study we have been unable to determine criteria for selecting the best REs for comparing complex microbial communities, and we believe the information from all six REs is valuable in our community analyses. Using six REs will also allow the detection of some problems with T-RFLP, such as the generation of artifacts that may occur during PCR (31).

Conclusions.

The factors governing the effect of RE choice are poorly understood, and the outcomes of computer-generated simulations (11) did not correlate well with the outcomes of our studies on real samples. We conclude that the use of multiple REs, employed individually, overcomes possible effects of RE choice on generating useful T-RFLP information for complex bacterial communities. Our results suggest that our variable percentage threshold method is a useful addition to the limited range of standardization methods available, because it allows determination of an optimal threshold for each profile, and so minimizes information loss. We also suggest that generating multiple profiles using different REs for each sample, and so giving confidence values for sample groupings, may be a means of determining the significance of relationships detected. We used six REs, but this could be varied. In conjunction with pooling PCRs (6) and generating consensus profiles from multiple separations (10), more confident interpretation of T-RFLP analyses can be made. The consensus approach could also be used with other gene profiling techniques (21), where multiple profiles could be generated using different primer sets, and consensus values generated as indicators of confidence.