Characterizing the fecal microbiota of infants with botulism

Background Infant botulism is the most prevalent form of botulism in the USA, representing 68.5 % of cases reported from 2001–2012. Infant botulism results when botulinum toxin-producing clostridia (BTPC) colonize the infant gut with concomitant in vivo production of the highly potent botulinum neurotoxin (BoNT). The gut microbiota of infants with botulism is largely uncharacterized; therefore, it remains unclear whether the microbiota profile of these patients are distinct in composition, abundance, or diversity. To address this uncertainty, we employed 16S rRNA gene profiling to characterize the fecal microbiota in 14 stool samples among laboratory-confirmed and non-confirmed infant botulism cases. Results Seven bacterial phyla were identified among all 14 infant stool samples examined. Compared to samples from non-confirmed cases, the fecal microbiota of infant botulism patients displayed significantly higher Proteobacteria abundance. Of the 20 bacterial families identified, Enterobacteriaceae was significantly more abundant in samples from infants with botulism. Firmicutes abundance and the abundance ratio of Firmicutes/Proteobacteria was significantly lower in samples from infants with botulism. Lactobacillus spp. abundance was notably reduced in 12 of the 14 samples. Clostridium botulinum and Clostridium baratii were identified in low relative abundances in confirmed and non-confirmed samples based on their 16S rRNA gene profiles, although their toxigenicity remained undetermined. No significant differences were observed in the number of operational taxonomic units (OTUs) observed or in fecal microbiota diversity between laboratory-confirmed and non-confirmed samples. Correlations between individual phylum abundances and infant age were variable, and no significant differences were shown in number of OTUs observed or in fecal microbiota diversity between samples delineated by overall mean age. Conclusions Significant differences in Proteobacteria, Firmicutes, and Enterobacteriaceae abundances were identified in the fecal microbiota of infants with botulism when compared to samples from non-confirmed cases. Fecal microbiota diversity was not significantly altered in infants with botulism, and a limited presence of BTPC was shown. It could not be determined whether the fecal microbiota profiles shown here were comparable prior to patient illness, or whether they were the direct result of infant botulism. The results of this study do, however, provide a detailed and descriptive observation into the infant gut microbiota after intestinal colonization by BTPC.


Background
From 2001 to 2012, the National Botulism Surveillance System reported an average of 142 cases per year of botulism in the USA [1]. Infant botulism is responsible for 68.5 % of these cases. Each year, the National Botulism Laboratory Team (NBLT) at the Centers for Disease Control and Prevention (CDC) investigates approximately 20 suspected cases of infant botulism submitted by states across the USA. First described in 1976 [2], infant botulism occurs in children less than 1 year old, and infant susceptibility to gut colonization of botulinum neurotoxin (BoNT) producing clostridia (BTPC) may be associated with perturbations of the developing infant gut microbiota [3]. Strains of Clostridium botulinum, Clostridium baratii, and Clostridium butyricum have each been implicated as the causative agents in infant botulism cases. Positive confirmation of infant botulism is established when BoNT and/or isolates of BTPC are detected in the stool.
Compared to adults, infants harbor an intestinal microbiota (i.e., the assemblage of microorganisms which occupy the gut) that is transient and minimally complex [4,5]. A number of factors are known to influence the microbial composition within the infant gut including age [6][7][8], mode of delivery [9,10], diet [11][12][13], antibiotic use [14][15][16], and disease [17,18], making it difficult to delineate their respective effects. An important function of the microbiota beyond digestion and metabolism is to prevent the propagation of pathogenic microorganisms [19][20][21]. It has been suggested that alterations of the developing infant gut microbiota may suppress BTPC colonization resistance, thus increasing susceptibility to botulism in certain individuals. Although studies using mouse models appear to support this claim, the etiology of this process remains unclear [22,23].
Some risk factors identified for infant botulism include honey consumption, breastfeeding status, and living in a rural area [24][25][26][27][28]. Many of these risk factors have been shown to independently influence the microbiota in both healthy and diseased infants. However, delineating the degree of influence of botulism on the infant gut microbiota from disease-independent factors is a difficult challenge and not the focus of this study. Our aim was to provide a detailed descriptive analysis of the fecal microbiota of infants with botulism and to identify any significant alterations in bacterial abundance, composition, and/or diversity among samples from confirmed and non-confirmed infant botulism cases. To characterize the prokaryotic members of the infant gut microbiota, we targeted the hypervariable V3 region of the 16S rRNA gene for high-throughput sequencing and performed extensive downstream taxonomic analysis. This methodological approach allowed for a thorough and descriptive taxonomic characterization of the patient fecal microbiota and provided a means to evaluate disparities between the gut microbiota of infants with botulism and the gut microbiota of infants from nonconfirmed botulism cases.

Evaluation of amplification and sequencing bias
Each of the five phyla from the mock community sample were successfully identified from sequence read analysis (Fig. 1). Overall mean absolute percent error (MAPE) for abundance quantification of the microbial mock community was 14.4 % (SD ± 33.1). Overall error was heavily skewed (MAPE = 74 %; SD ± 10) toward two phyla (Bacteroidetes and Deinococcus-Thermus), which represented the bottom 10 % of total phylum abundance. Abundance quantification of the top 90 % of phyla resulted in a MAPE of 7.6 % (SD ± 0.27). Results from mock community sequencing indicated that the study methodology employed was largely effective for phylum- Fig. 1 Histograms of mock community DNA-sequencing results. Expected phylum abundance percentages were calculated by classifying each of the 20 known bacterial species into their respective phylum. Observed phylum abundance percentages were calculated by assigning 16S rRNA gene sequence reads to representative taxonomy by alignment with Greengenes ribosomal database using an OTU definition of 97 % sequence homology or greater with a minimum cluster size of 3 level taxonomic analysis; however, the percent error may increase when quantifying underrepresented phyla.

Sample and patient characteristics
We analyzed 14 infant stool samples from 14 patients representing a broad range of ages. Eight samples were associated with laboratory-confirmed botulism cases representing BoNT serotypes A, B, and F (designated as "confirmed"), and six samples were from cases where no laboratory confirmation of botulism could be made (designated as "non-confirmed"). Mean patient age across all 14 samples was 87.9 days (SD ± 98.3). Mean patient age in confirmed samples was 68.0 days (SD ± 88.7), and mean patient age in non-confirmed samples was 114.3 days (SD ± 120.8). There was no significant difference in patient age between confirmed and nonconfirmed samples (p = 0.22; alpha = 0.05). Coded sample names, BoNT type detected, quantitative PCR results, and patient ages at the time of sample collection are listed in Table 1. No additional patient or sample metadata was available for these cases.

Taxonomic and statistical analysis of fecal microflora
Sequence read statistics and alpha diversity metrics from 16S rRNA gene sequencing are listed in Table 2. The number of operational taxonomic units observed (OTU; defined as a subset of 3 or more sequence reads that share 97 % or greater sequence homology) was lowest for CDC68053 and highest for CDC68083 with an average number of OTUs per sample of 902.4 (SD ± 502.0). There was no significant difference in number of OTUs observed between confirmed and non-confirmed samples (p = 0.61; alpha = 0.05). Shannon diversity index, a measure of bacterial richness and evenness, revealed that the fecal microbiota from sample CDC68083 was highest in bacterial diversity while sample CDC68116 displayed lowest diversity. Fecal microbiota diversity was not significantly different between confirmed and nonconfirmed samples (p = 0.52; alpha = 0.05).
Family-level taxonomic analysis of the fecal microflora in the infant stool samples revealed the presence of 20 bacterial families (Table 3). Sequence reads that could not be assigned to a bacterial family, and families with <1.0 % abundance in each of the 14 samples, were designated as "Other". The Gammaproteobacteria family Enterobacteriaceae displayed significantly higher abundance (p = 0.02, alpha = 0.05) in confirmed samples (x ̅ = 44.0 %; SD ± 33.5) compared to non-confirmed samples (x ̅ = 8.4 %; SD ± 15.7) (Fig. 3). No other bacterial family displayed a significant difference in abundance between confirmed and non-confirmed samples.
Based on 16S rRNA gene taxonomy, C. botulinum was identified in 6 of the 8 confirmed samples and 2 of the 6 non-confirmed samples, with relative abundances ranging from <0.001 to 0.01 % of the overall bacterial community where detected (Table 4). C. baratii was identified in 5 of the 8 confirmed samples and 2 of the 6 non-confirmed samples with relative abundances ranging from <0.001 to 0.003 % of the overall bacterial community where detected (Table 4).   Samples were re-grouped relative to mean patient age of all 14 samples to examine the association between age and phylum-level abundance. Nine samples were collected from patients younger than the overall mean patient age (collectively referred to as the "younger" group), and 5 samples were collected from patients older than the overall mean patient age (collectively referred to as the "older" group). Proteobacteria abundance was significantly higher (p = 0.04; alpha = 0.05) among the younger group (x ̅ = 50.3 %; SD ± 28.5) compared to the older group (x ̅ = 8.8 %; SD ± 8.2) (Fig. 5). No other phylum displayed a significant difference in abundance between the younger and older group. The average ratio of Firmicutes/Proteobacteria was 170 times lower (p = 0.13; alpha = 0.05) in the younger group (x ̅ = 4.1; SD ± 9.3) compared to the older group (x ̅ = 698.3; SD ± 1526.2), and the average ratio of Firmicutes/Bacteroidetes was 6.0 times lower (p = 0.57; alpha = 0.05) in the younger group (x ̅ = 13671.8; SD ± 20643) compared to the older group (x ̅ = 2287.1; SD ± 3028.2) (data not shown). There was no significant difference (p = 0.28; alpha = 0.05) in number of OTUs observed between the younger infant group and the older infant group. Fecal microbiota diversity was not significantly different between the younger infant group and the older infant group (p = 0.20; alpha = 0.05).
Redundancy analysis (RDA) was used to explore associations among infant microbiota composition and abundance, and to examine the influence of infant age on sample variance (Fig. 6) [29]. Taxonomically unclassified data were excluded from RDA analysis, as was data representing the phylum Fusobacteria which was only detected in trace abundance (<0.002 %) in one sample. RDA explained 83.7 % of all variation among samples. Overall sample variance was produced by the combined effect of bacterial abundance, bacterial composition, and infant age (F = 3.5; p = 0.048; alpha = 0.05). RDA confirmed a negative correlation between Proteobacteria abundance and infant age, and positive correlations between infant age and Tenericutes, Actinobacteria. Bacteroidetes, and Verrucomicrobia abundance. RDA revealed that Proteobacteria abundance was positively correlated with 7 samples (CDC68116, CDC68058, CDC68068, CDC68129, CDC68053, CDC68083, and CDC69114). Three samples (CDC68126, CDC69078, and CDC68090) were positively correlated to Tenericutes, Actinobacteria, Verrumicrobia, and to a lesser extent Bacteroidetes abundance. The phylum Firmicutes displayed a negative correlation to all other phyla by RDA, and displayed a positive correlation with samples CDC69062, CDC68122, CDC68112, and CDC68192. The influence of age was weighted most heavily on

Discussion
Characterizing the fecal microbiota in infants with botulism provides another step toward understanding the link between the infant gut microbiota and this disease. Although enteric diseases can produce discernable, if not predictable, alterations to the bacterial gut microbiota [30], it can also be influenced by a wide variety of disease-independent perturbations [31]. The findings of our study revealed the fecal microbiota of infants with botulism was characterized by an over-enrichment of Proteobacteria with low Firmicutes abundance. This finding highlighted a significant disparity between confirmed and non-confirmed samples with a clear deviation from a more traditional healthy-infant microbiota profile [16]. Previous studies have reported reduced Firmicutes abundance with a concomitant increase in Proteobacteria abundance in response to disease-mediated alterations of the gut microflora [32,33]. Firmicutes are often highly represented in the gut microbiota of healthy individuals and can be reduced in illness [34], while a significant increase in Proteobacteria abundance can lead to gastrointestinal inflammation in response to environmental and genetic factors [35].
A notable disparity was identified among the relative proportions of Lactobacillus spp. identified among samples. The genus Lactobacillus contains over 180 species of bacteria, many of which are commonly found propagating the microbiota of healthy infants. Lactobacillus spp. abundance was markedly reduced (mean relative abundance = 0.2 %; SD ± 0.004) in the fecal microbiota of all 8 confirmed samples, and 4 of the 6 non-confirmed samples (data not shown). Similar reductions have been reported in previous studies investigating the fecal microbiota of infants [5,[36][37][38]. Among these reports, it was suggested that a variety of factors (e.g. caesarian delivery, aggressive antibiotic treatment, long-term incubation) may combine to suppress Lactobacillus colonization, significantly delay the establishment of beneficial bacteria, and increase susceptibility to enteric pathogen colonization.     Lactobacilli have been shown to provide some inhibition activity against several common food pathogens [39]. Although it remains unclear whether the reduction of beneficial bacteria such as Lactobacillus spp. can increase susceptibility to BTPC colonization in infants, this finding illustrates a distinct feature of the gut microbiota profile in infants once BTPC colonization has been established. Mariat et al. [40] illustrated the utility of using phylum-level abundance ratios to profile the infant gut microbiota. In that study, it was reported that the ratio of Firmicutes/Bacteroidetes in healthy infants was significantly lower than that of adults. Although we did not observe a significant difference in Firmicutes/Bacteroidetes abundance ratios between our sample groups, we do report a significant decrease in the abundance ratio of Firmicutes/Proteobacteria in confirmed infant botulism samples; a ratio which was 515-times lower than that observed in samples from non-confirmed cases. Whether this reduction preceded BTPC colonization, or whether it was a consequence of botulism could not be concluded. It should be noted that phylum abundances could have been influenced by gaps in sample collection (the time between onset of illness and sample collection) and sample storage times prior to DNA extractions. Nevertheless, taxonomic analysis indicated that the primary disparity between confirmed and non-confirmed samples resulted from significant differences in the relative abundances of Firmicutes and Proteobacteria within the fecal microbiota.
Only a handful of studies have reported their efforts to characterize the fecal microflora in infants with botulism. Long et al. (1985) [27] examined stool samples from seven infant botulism patients and found a high proportion of enterobacteria in the fecal microflora. Similarly, we report significantly higher Enterobacteriaceae abundance in the fecal microflora of infants with botulism compared to samples from non-confirmed cases. A comparison of family-level abundances revealed Enterobacteriaceae to be the only bacterial family which displayed a significant disparity between confirmed and non-confirmed samples. Enterobacteriaceae was the  In another study, Wilcke et al. [26] quantified C. botulinum from stool samples of four infants with laboratoryconfirmed infant botulism and found that C. botulinum abundance ranged from 0.01 to 3.3 % of culturable fecal microbiota. In our study, neither C. botulinum nor C. baratii ever exceeded 0.01 % of the fecal microbiota. Although it is likely that BTPC intestinal abundance fluctuates temporally over the course of the illness, these findings do support the widely-held assumption that botulism can manifest from a relatively limited presence of BTPC in the infant gut. Unfortunately, there is no available data to indicate whether there is a minimum threshold for BTPC colonization needed to produce an amount of toxin required for disease onset. Given the extremely high potency of BoNT these organisms are capable of producing, coupled with the low weight of infants relative to adults, it was not surprising to find no more than 0.01 % BTPC abundance in stool samples from infant botulism patients.
The identification of C. botulinum and C. baratii 16S rRNA gene sequences in both botulism-confirmed and non-confirmed samples highlight the current limitations of using a sequencing-based approach for botulism diagnostics. Infant botulism is confirmed by the laboratory detection of BoNT and/or the detection of BTPC in infant stool. Although C. baratii can produce BoNT subtype F, non-toxigenic strains are more commonly isolated from stools. Non-toxigenic C. baratii do not harbor the serotype F botulinum neurotoxin gene (bont/F) within their genome, and as a result, 16S rRNA gene sequencing does not provide the necessary resolution to distinguish between toxigenic and non-toxigenic isolates of this organism. Similarly, amplicon sequencing cannot resolve C. botulinum from non-toxigenic strains of clostridia (such as C. sporogenes) that share >99 % 16S rRNA gene identity.
The temporal development of the infant gut microbiota throughout the first year of life is not consistent across all individuals [7,8]. This yearlong progression from near sterility to a fully developed microbiome is guided by a myriad of external factors and can be interrupted by illness [6,10,44]. With the exception of Proteobacteria abundance, the microbiota in the younger infant group did not significantly differ from that of the older infant group in composition, abundance, or diversity. RDA indicates a negative correlation between infant age and samples with higher Proteobacteria abundance, and a positive correlation between infant age and samples with higher abundances of Actinobacteria, Tenericutes, Verrumicrobia, and Bacteroidetes. Excluding age as a co-variable in RDA (data not shown) did not alter the cluster patterns shown in Fig. 6, nor the correlations approximated between samples and phylum composition and abundance. Together, these findings indicate that differences in fecal microbiota abundance and composition between botulism-confirmed and non-confirmed samples were not due solely to patient age disparities.

Conclusions
The infant gut microbiota undergoes a dramatic transition throughout the first year of life, and a variety of factors guide its development. Perturbations of the infant gut microbiota have been shown to increase susceptibility to opportunistic pathogens, yet it remains unclear whether similar disturbances can be associated with increased susceptibility to BTPC colonization. Researchers have speculated that infants display an increased susceptibility to botulism due to their transient and minimally complex gut microbiota; however, there is little published data to support this claim. In an effort to resolve some of the uncertainty surrounding this speculation, we conducted the first detailed analysis into the gut microbiota of infants with botulism using high-throughput, nextgeneration sequencing technology.
We have identified several distinct features among BTPC-colonized microbiota which are comparable to those previously identified in infants suffering from botulism and other enteric diseases. Our findings revealed that, compared to that of non-confirmed samples, the fecal microbiota of infants with botulism contained significantly higher abundances of Proteobacteria and Enterobacteriaceae, while Firmicutes abundance remained significantly lower. Additionally, our findings indicated that botulinum toxin-producing clostridia were present in very low abundances, indicating that botulism can manifest from a relatively limited presence of these organisms in the infant gut.
These findings can be used to guide future infant botulism research efforts by directing focus toward several distinct features of the infant gut microbiota which have been identified here. The data generated from this study has greatly expanded our insight into the gut microbiota profiles in infant botulism patients. Ultimately, our understanding of infant botulism will continue to increase as sequencing technologies and computational analyses are further developed. These developments should greatly increase our capabilities as investigators interested in gaining a deeper understanding of botulism and its influence on children impacted by this disease.

Sample properties and patient characteristics
In compliance with a human subjects exemption protocol (#4991.0) approved by the CDC Human Research Protection Office, coded stool samples were selected for investigation from 14 infants suspected to have botulism with patient ages ranging from 14 to 350 days old. The samples were selected to represent a broad array of sample properties and patient characteristics. Of the 14 samples sent for botulism testing, NBLT confirmed 8 samples for infant botulism (confirmed), while 6 could not be confirmed (non-confirmed). Refer to Table 1 for complete sample properties and patient characteristics.

Genomic DNA extraction and sequencing
Genomic DNA was extracted from each stool sample using the PowerFood ® Microbial DNA Isolation Kit (MO BIO Laboratories, Inc., Carlsbad, CA). Libraries were prepared for Ion Torrent PGM sequencing following the Fusion Method-Ion Amplicon Library Preparation for bidirectional sequencing (Ion Amplicon Library Preparation [Fusion Method] User Guide; Publication number 4468326, Revision C). Extracted gDNA was purified using the DNA Clean and Concentrator-5 kit (Zymo Research, Irvine, CA) and quantified on a NanoDrop 2000 UV-vis Spectrophotometer (Thermo Scientific, Wilmington, DE). For bidirectional sequencing, two direct PCR amplifications targeting the hypervariable V3 region of the 16S rRNA gene were performed on each sample gDNA. PCR target primers were previously described by Zhang et al. [45] and synthesized with Ion Torrent compatible adaptors and key sequences to generate approximately 280 bp PCR amplicons (Table 5). Each 50 μL PCR reaction consisted of 45 μL of Platinum PCR SuperMix High Fidelity buffer (Life Technologies, Grand Island, NY), 1 μL each of 10 μM primers and 50 ng of purified gDNA. PCR thermocycling conditions was set up with an initial denaturation at 94°C for 2 min followed by 35 cycles of 94°C for 30 s, 58°C for 30 s, and 68°C for 30 s with a final extension of 68°C for 2 min. Following amplification, PCR products were purified using the DNA Clean and Concentrator-5 kit (Zymo Research), and size selection on the purified amplicons was performed using the E-Gel Agarose Gel Electrophoresis System (Life Technologies). Each sizeselected library was quantified using the Qubit dsDNA HS Assay Kit (Life Technologies), diluted to 18 pM, and pooled into an equimolar solution. Each DNA fragment library was templated onto Ion Sphere Particles (ISPs) via emulsion PCR using the Ion OT2 instrument, and the resulting template-positive ISPs were quality checked to assess templating efficiency. Template-positive ISPs were enriched on the Ion ES instrument, loaded onto Ion 314 v2 chips, and sequenced using the Ion PGM.
Genomic DNA from Microbial Mock Community B (Even, High Concentration), v. 5.1H (bei Resources) was used to evaluate the effects of amplification and sequencing bias. The mock community sample of gDNA consisted of a single pool of equimolar RNA operon counts from 20 bacterial strains representing 5 bacterial phyla. Library preparation, sequencing, and taxonomic analysis of the mock community DNA was performed in the same manner as stool sample gDNA. A phylum-level comparison of actual abundance (known phylum percent distribution of the mock community) versus observed abundance (PGM sequencing results) was performed. The formula below was used to calculate mean absolute percentage error (MAPE) in abundance quantification (F t = observed abundance, A t = actual abundance, n = 5): Raw sequence data and the associated metadata for all 14 samples have been deposited in MG-RAST metagenomics analysis server (http://metagenomics.anl.gov). All sample data is freely available and openly accessible from the MG-RAST website under project number 13930, or by navigating directly to http://metagenomics.anl.gov/ metagenomics.cgi?page=MetagenomeProject&project= 13930.

Sequence analysis and taxonomic characterization
Sequence reads were analyzed with the open source software package Quantitative Insights into Microbial Ecology v1.8.0 (QIIME) [46]. Prior to analysis, reads with a read length of less than 150 bp or greater than 350 bp, reads with missing quality scores or mean quality scores below 25, reads with primer mismatches, and reads with greater than 6 ambiguous bases were removed from analysis. The remaining reads were aligned using Python Nearest Alignment Space Determination (PyNAST), and clustered into operational taxonomic units (OTUs). OTU clustering is a metric of binning sequence reads where a single OTU is defined as a subset of reads (minimum cluster size of 3) which share 97 % or greater sequence homology. The QIIME workflow used for this analysis can be accessed via https://github.com/tbshirey/QIIME.
The representative taxonomic identities were assigned by alignment with the Greengenes ribosomal database [47], and stacked histograms of the relative taxonomic abundances were generated showing phylum-level taxonomic representatives for each sample. Family-level percent abundances, Enterobacteriaceae abundance, and abundance of two expected BTPC (C. botulinum and C. baratii) were also recorded for each sample. Phylumlevel ratios of Firmicutes/Proteobacteria and Firmicutes/ Bacteroidetes were calculated based on average abundance percentages for each sample. Bacteroidetes was not detected in samples CDC68129 and CDC68192 which prevented calculation of the Firmicutes/Bacteroidetes ratio for these samples.

Statistical analyses
Significance testing was performed using the Mann-Whitney U rank-sum test in Statistical Package for the Social Sciences (SPSS; IBM Corporation, Armonk, NY). A non-parametric 0.95 confidence interval (α = 0.05) accompanied these estimates. Linear regression (reported as the adjusted R 2 value; R ̅ 2 ) was performed in SPSS to model the correlation between taxonomic abundance and infant age. Shannon diversity index illustrating inter-sample fecal microflora taxa diversity was calculated in QIIME by rarifying all samples. This index is a commonly used measure of diversity which uses an algorithm to account for both species richness (i.e., the number of species present) and species evenness (i.e., the distribution of species) within a sample. The Shannon diversity index estimates diversity on a scale beginning at 0 (no diversity observed) and increasing where greater diversity is estimated.
RDA was performed using the software package Canoco 4.5 (Microcomputer Power, Ithaca, New York) to examine the relationships between sample microbiota community structure, disease confirmation, and infant age. Bacterial phyla (dotted blue vectors) served as response variables, and patient age at the time of sample collection (solid red vector) was used as a predictor co-variable. Sample symbols were generated to indicate disease confirmation (crosses = botulism confirmed samples; circles = nonconfirmed samples). Vector length is proportional to abundance (phylum vectors) or age (age vector), with each variable increasing in the direction of the vector. Correlations between samples and vectors are approximated based on sample proximity to each vector (indicating correlation between the sample and phylum composition and/or age), and the location of the sample point along the length of each vector (indicating correlation between the sample and phylum abundance, and/or age). Prior to ordination, Monte Carlo permutation tests (499 permutations under reduced model) were performed on nontransformed data. RDA was visualized using CanoDraw™ (Microcomputer Power). All other data visualized in tables, histograms, and bar charts were managed in Microsoft Excel (Microsoft, Redmond, WA), and box plots were generated using R (http://cran.r-project.org/).