Supplementary Methods
Description of the BC LTSP study sites
Detailed information about the Long-Term Soil Productivity installations in British Columbia can be found in the establishment report at http://www.for.gov.bc.ca/hfd/library/FIA/2009/FSP_Y092042b.pdf. Here we provide a synopsis of the most important characteristics. As shown in Figure 1, the six BC LTSP sites are located in two different biogeoclimatic zones, i.e. the Sub-Boreal Spruce (SBS) and Interior Douglas-Fir (IDF) zone respectively. The SBS is characterized by snowy, cold winters, and short, warm, and moist summers with mean annual temperature ranges from 1.7 to 5 °C. Mean annual precipitation can range from 415 to 1650 mm, with snow accounting for approximately 25-50% of the total precipitation. Climax tree species in the SBS are hybrid white spruce (Picea engelmannii x glauca) and subalpine fir (Abies lasiocarpa), whereas lodgepole pine (Pinus contorta) is common in maturing forests in the drier and more southern parts of the zone. Other seral tree species include Douglas-fir (Pseudotsuga mensiesii), trembling aspen (Populus tremuloides), and common paper birch (Betula papyrifera). The IDF is characterized by warm, dry summers and cool winters with mean annual temperatures range from 1.6 to 9.5 °C. Mean annual precipitation is between 300 and 750 mm but can exceed 1000 mm in the wettest subzones, and 20-50% of the precipitation falls as snow. Douglas-fir is the most common climax tree species in the IDF, while ponderosa pine (Pinus ponderosa) or hybrid white spruce and western redcedar (Thuja plicata) can be found on drier and wetter sites, respectively. Seral species include lodgepole pine, trembling aspen, and paper birch.
Figure 1. Locations of the six LTSP experimental sites examined in this study, i.e., SBS1-3 and IDF1-3. The Sub-Boreal Spruce (SBS) and Interior Douglas-Fir biogeoclimatic zones are indicated in blue and red, respectively.
Each LTSP site of the respective biogeoclimatic zone represents an experimental block in a randomized complete block design. For the SBS experiment, three sites were selected in zonal ecosystems (maturing seral to maturing climax stands) in three different subzones with an average distance of 285 km between sites. For the IDF experiment, three sites with predominantly zonal ecosystems (maturing seral to maturing climax stands) were selected in one subzone with an average distance of 9 km between sites. The average distance between the six sites in both biogeoclimatic zones is 395 km. Uniformity of stand structure as well as soil and site properties critical to plant growth were important criteria in the site selection process (Table 1, Hope 2006, Kranabetter and Sanborn 2003, Sanborn et al. 2000). All sites have deep (> 1 m), medium-textured soils and mor humus forms. Within each site, understory vegetation was uniform prior to treatment application.
PCR amplification and sequencing
The DNA concentration of samples was adjusted to 10 ng DNA µl-1 with ddH2O containing bovine serum albumin (BSA, New England Biolabs, Ipswich, MA) to yield a final concentration of 1 µg BSA µl-1. These samples were heated for 5 min at 90°C to bind PCR inhibiting substances like humic acids to the BSA (Rådström et al. 2004). PCR was performed using 50 ng soil DNA in a total volume of 50 µl containing 1× PCR-buffer (Qiagen Inc., Mississauga, ON), 2 mM MgCl2, 0.2 µM of each primer, 0.4 mM dNTP (Fermentas, Burlington, ON), 15 µg BSA and 2 U of HotStar Taq polymerase (Qiagen). Amplification was performed with initial denaturation for 15 min at 95 oC followed by 30 cycles with denaturation for 40 s at 94 oC, annealing for 40 s at 60 oC (bacterial 16S) or 58 oC (eukaryal ITS) respectively, and extension for 1 min at 72 oC, followed by a final extension for 10 min at 72 oC. The hypervariable region V1-V3 of the bacterial 16S rRNA gene was amplified using primers 27F (A-AGAGTTTGATCMTGGCTCAG, Lane 1991) and 519R (B-GWATTACCGCGGCKGCTG, Amann et al. 1995) where A and B are the two standard FLX primers (Roche 454 Life Sciences). The internal transcribed spacer 2 (ITS2) of the fungal ribosomal operon including parts of the 5.8S and large subunit rRNA gene was amplified using primers ITS3 (B-GCATCGATGAAGAACGCAGC, White et al. 1990) and ITS4 (A-TCCTCCGCTTATTGATATGC, White et al. 1990). Each sample was amplified in triplicate and subsequently pooled prior to purification. Each pooled PCR product was purified using Illustra GFX™ PCR purification columns (GE Healthcare, Piscataway, NJ). Purified PCR products were checked for quality and quantity with a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies, Wilmington, DE). Concentrations were adjusted to 40 ng µl-1 with Tris-HCl buffer. Barcodes and 454 GS FLX Titanium adapters were retrofitted to the amplicons at the Genome Québec Innovation Centre (Montréal, Canada) prior to sequencing. Barcodes were added to FLX primer A and tags were uni-directionally sequenced from that end.
Target extraction and verification
After quality filtering of the raw sequences using mothur (Schloss et al. 2009), the bacterial 16SV1-V2 and fungal ITS2 regions were verified and extracted using v-xtractor (Hartmann et al. 2010) and its equivalent tool for fungal ITS (Nilsson et al. 2010), respectively. The ITS primers used (White et al. 1990) amplify the ITS region for several other groups of eukaryotes in addition to fungi, but the extraction software tool only targeted and extracted fungal ITS (Nilsson et al. 2010). As for 16S, the V1-V2 region was selected as the segment of choice within the V1-V3 amplicon. In an initial assessment, the V1-V2 region provided a better trade-off between number of reads and phylogenetic information when compared to shorter (V1, higher number of reads but lower success of taxonomic assignment) or to longer segments (V1-V3, higher success of taxonomic assignment but lower number of reads) of the amplicon. Importantly, these software tools extract gene segments that are comparable for community analysis and at the same time confirm basic authenticity of the target (Hartmann et al. 2010).
OTU Clustering with Crunchclust
Sequences were clustered into OTUs using crunchclust, a novel scalable open-source software tool written in C++. crunchclust uses a greedy incremental clustering algorithm optimized for large amplicon data sets derived from 454 pyrosequencing and accounts for frequent homopolymer detection errors considered critical when using this technology (Behnke et al. 2011, Huse et al. 2010, Quince et al. 2009). Since pyrotag data sets comprise a very large number of reads and encompass mostly variable domains, clustering based on multiple sequence alignments are frequently inaccurate (Liu et al. 2010), difficult to compute as data sets grow (Liu et al. 2010), and artificially inflate OTU richness (Sun et al. 2009, Huse et al. 2010). Aligning pyrotags against well-curated reference alignments can improve accuracy and computational practicability (Schloss 2009), but does not work well for new deeply branched lineages, domains with high divergence rates, or sequences with homopolymer errors. Therefore, obtaining reasonable reference alignments for many hypervariable genetic regions such as the ribosomal internal transcribed spacers is unlikely. As a result, clustering algorithms independent of multiple sequence alignments have been developed (Edgar 2010, Huang et al. 2010), which however do not solve the problems associated with the inherent noise. For this matter, dedicated denoising pipelines have been released that try to eliminate erroneous reads (Huse et al. 2010, Quince et al. 2009, Quince et al. 2011, Reeder and Knight 2010). Unfortunately, these pipelines are often slow and complex to run.
In one single step, crunchclust de-replicates raw reads sorting them according to their abundance and clusters the de-replicated reads into OTUs using a global exact Needlemann-Wunsch (Needleman and Wunsch 1970) pairwise alignment to calculate the exact number of differences between two sequences while ignoring dissimilarities due to homopolymer lengths and terminal gaps when the distal primers are not reached. Although real variations in homopolymer length will go undetected using the optional homopolymer-filter, this strategy has proven to be an efficient means to reduce artificial inflation of OTU richness caused by sequencing errors, which have been reported to account for a major part of errors in such datasets (Behnke et al. 2011). These homopolymer errors by far exceed the true biological variation in homopolymer length among sequences. crunchclust does not require high performance computing resources (although it is MPI compatible) and runs efficiently on laptop computers using a single CPU core. Memory requirements are negligible and the software is supported on Linux, Windows, and Mac operating systems. crunchclust is available at http://code.google.com/p/crunchclust/. We benchmarked crunchclust against related “denoising” pipelines (Huse et al. 2010, Quince et al. 2009, Quince et al. 2011, Reeder and Knight 2010) and usearch (Edgar 2010) using 454 data sets of mock communities equivalent to one titanium run. Where other pipelines required days to process this information, both crunchclust and usearch required only minutes, the former using an exact method and the latter a fast heuristic one. However, usearch does not account for homopolymer errors, whereas crunchclust does. In conclusion, crunchclust represents a rapid algorithm dedicated to accurately clustering massively parallel pyrotag sequence data sets to produce denoised OTUs for use in downstream analytic pipelines such as mothur.
Multivariate analysis of community structures and diversity
Data were standardized by dividing the number of reads in each taxonomic unit by the total number of reads in each sample. Standardized data were square root transformed to downweight the contribution of quantitatively abundant OTUs to the similarities calculated between samples. A resemblance matrix was calculated using Bray-Curtis similarities (Bray and Curtis 1957) based on the standardized square root transformed read abundance data. The Bray-Curtis similarity coefficient is one of the most widely applied resemblance measures and meets many criteria desired for application to ecological data (Clarke et al. 2006 #10473). Principal Coordinate Analysis (PCO, Gower 2005) implemented in the primer6+ package was used to display similarities in microbial community structures among all samples. PCO is an unconstrained metric multidimensional scaling ordination that extracts major variance components from the multivariate data set in order to reduce dimensionality of the data cloud by minimizing the residual variation in the space of any chosen resemblance measure.
Tests of the multivariate null hypotheses of no differences among a priori defined groups were examined using analysis of similarities (ANOSIM, Clarke 1993), permutational multivariate analysis of variance (PERMANOVA, Anderson 2001), and canonical analysis of principal coordinates (CAP, Anderson and Willis 2003) implemented in primer6+. These non-parametric discriminative methods based on permutation tests do not rely on assumptions that are commonly too stringent for most ecological data sets (Anderson 2001) and base the analysis on multivariate distance measures that are reasonable for these data (McArdle and Anderson 2001). ANOSIM is a univariate non-parametric analogue of analysis of variance (ANOVA) and provides the test statistic R that equals 0 when there are no differences among groups and 1 for the maximum difference among groups (Clarke 1993). PERMANOVA is analogous to multivariate analysis of variance (MANOVA), which allows partitioning the variability of the data according to a complex design or model, and provides F-ratios that are analogous to Fisher's F-ratio in MANOVA (Anderson 2001). CAP uses PCO followed by canonical discriminant analysis to provide a constrained ordination that maximizes the differences among a priori groups and reveals patterns that can be cryptic in unconstrained ordinations (Anderson and Willis 2003). The canonical correlation (δ2) of each CAP axis is an estimate of the amount of shared variance between the respective canonical variates of dependent and independent variables. CAP also implements the ‘leave-one-out allocation’ approach to determine the misclassification error or alternatively the classification success provided by the model (Lachenbruch and Mickey 1968). The ratio between source (known affiliation) and successfully classified (predicted affiliation) data provides an quantitative estimate of the degree of discrimination among the groups achieved by the canonical axes or, more general, a goodness of fit of the model (Anderson and Willis 2003). Ordination and classification success have to be examined in combination to draw conclusions about separation of a prior groups. Finally, CAP tests the null hypothesis of no difference among the groups by calculation of a trace statistic traceQ_m'HQ_m and a P value based on permutations (Anderson and Willis 2003). Permutational analysis of multivariate dispersions (PERMDISP, Anderson 2006) was used to test for heterogeneity of community structure in a priori groups. ANOSIM, PERMANOVA, CAP, and PERMDISP analyses were performed with 9,999 permutations and run as routines included in the primer6+ software package. In the PERMANOVA routine, residuals of raw data were permuted under a reduced model and partitioning was performed using Type I sums of squares, which is recommended when covariates are defined (Clarke and Gorley 2006). Location factors such as horizon and site were defined as covariates in order to subtract the variance explained by these factors. Since Type I is sequential and the input order of the variables matters, the analyses was re-run after swapping the order of the treatment factors. The number of PCO axes included in the CAP analysis was automatically selected by the diagnostics routine in primer6+, choosing the number of axis showing the lowest misclassification error. Pearson correlations coefficients between species abundance and CAP coordinates were calculated in primer6+. All ordinations were plotted with statistica 8.0 (StatSoft, Tulsa, OK).
mothur was used to calculate Good’s coverage (Good 1953) as well as the Simpson’s index for evenness (E1/D, Simpson 1949). Simpson’s indices are among the measures that perform the best on ecological data (Magurran 2004). Importantly, E1/D is largely independent from species richness and sampling size (Magurran 2004), making it well suited for pyrotag data where the accuracy of richness estimates is questionable. In order to account for variability among different experimental sites, Simpson’s evenness was normalized relative to the reference plots within each site. The effect of harvesting was examined by only including samples at compaction level C0 and using one-way ANOVA followed by Fisher’s Least Significant Difference (LSD) post-hoc test. The effects of the different treatment levels within the factorial design was examined by excluding the reference samples and using factorial ANOVA followed by a desirability approach (Derringer and Suich 1980) to find the combination of OM removal and soil compaction where evenness and diversity are maximized.