Analysis of the Diversity of Substrate Utilisation of Soil Bacteria Exposed to Cd and Earthworm Activity Using Generalised Additive Models

Biolog EcoPlates™ can be used to measure the carbon substrate utilisation patterns of microbial communities. This method results in a community-level physiological profile (CLPP), which yields a very large amount of data that may be difficult to interpret. In this work, we explore a combination of statistical techniques (particularly the use of generalised additive models [GAMs]) to improve the exploitation of CLPP data. The strength of GAMs lies in their ability to address highly non-linear relationships between the response and the set of explanatory variables. We studied the impact of earthworms (Aporrectodea caliginosa Savigny 1826) and cadmium (Cd) on the CLPP of soil bacteria. The results indicated that both Cd and earthworms modified the CLPP. GAMs were used to assess time-course changes in the diversity of substrate utilisation (DSU) using the Shannon-Wiener index. GAMs revealed significant differences for all treatments (compared to control -S-). The Cd exposed microbial community presented very high metabolic capacities on a few substrata, resulting in an initial acute decrease of DSU (i.e. intense utilization of a few carbon substrata). After 54 h, and over the next 43 h the increase of the DSU suggest that other taxa, less dominant, reached high numbers in the wells containing sources that are less suitable for the Cd-tolerant taxa. Earthworms were a much more determining factor in explaining time course changes in DSU than Cd. Accordingly, Ew and EwCd soils presented similar trends, regardless the presence of Cd. Moreover, both treatments presented similar number of bacteria and higher than Cd-treated soils. This experimental approach, based on the use of DSU and GAMs allowed for a global and statistically relevant interpretation of the changes in carbon source utilisation, highlighting the key role of earthworms on the protection of microbial communities against the Cd.


Biolog Ecoplates
Biolog MicroPlates were developed in the late 1980s to assist in the identification of bacterial strains [1]. These 96-well plates contained carbon sources and a tetrazolium violet redox dye that turned purple if inoculated microorganisms utilised these sources. By comparing the obtained carbon substrate utilisation patterns with databases, it was possible to establish a probable identification [2]. Later, microbial ecologists used Biolog plates to investigate patterns at the community level. Therefore, a new plate specifically designed for community analysis and microbial ecological studies was created; this new plate was referred to as the EcoPlate [3]. The EcoPlate contains 31 of the most useful carbon sources for soil community analysis [3], allowing for community-level physiological profiling (CLPP) of heterotrophic bacterial assemblages. This technique has also been widely used to assess the toxicological impacts of different pollutants [4], including different heavy metals [5][6][7][8][9].
The use of EcoPlates results in a CLPP, yielding a very large amount of data that may be difficult to interpret. These difficulties are related to the alteration of the original microbial communities due to soil sampling and pre-treatment as well as bacterial extraction [4,10]. Moreover, the Biolog procedure can be considered a culture method in which the originally inoculated community is altered. The contribution of a certain species population to the colour profile will depend on its culturability in the Biolog wells and on its interactions with other species (differential growth and competition). The tetrazolium dye also introduces some bias in the profile because not all bacteria are able to reduce it [4]. Thus, the Biolog method is more useful for comparing soil microbial communities than for community characterisation [4]. CLPPs provide little insight regarding the function of the community in situ unless they are combined with other microbial methods that do not rely on the culturing of the soil microflora [11,12].
The pattern of positive and negative responses as well as substrate oxidation rate and extent are highly reproducible for simple microbial communities, particularly if the inoculum densities are similar [10]. Nonetheless, data corrections are required if inoculum density is not controlled. A kinetic approach must be used to extract data from the plates to avoid problems related to incubation time. However, there is no consensus regarding which statistics are most appropriate, and the interpretation of their meaning remains unclear [4].

Earthworms and cadmium (Cd) on soils
Earthworms are key members of the soil macrofauna in temperate soils [13] that directly or indirectly modulate resource availability (quality, quantity, and distribution) for other organisms. Earthworms are considered ecosystem engineers because they modify, maintain, or create habitats [14] as they build soil biostructures, which consist of aggregates and macropores [15][16][17]. The soil zone influenced by earthworm burrowing and casting was first termed the drilosphere by Bouché, 1975 [18,19], where microbial biomass is greater than in the surrounding soil [20]. In general, the mechanical and biological activities of earthworms favour organic matter mineralisation and humification and stimulate microbial activity [13].
Cd levels in the environment vary widely, and the average natural abundance of this element in the earth's crust varies from 0.1 to 0.5 mg kg 21 . However, human activities, such as mining, metal smelting, urban and industrial emissions, waste incineration, coal combustion, traffic dust, and particularly the use of phosphate fertilisers and sewage sludges, are quoted as the primary reason for the increase in soil Cd content over the last several years in Europe [21]. As heavy metals cannot be degraded, they tend to accumulate in soils. Their disappearance indicates that they have leached into deep soil layers and transferred to subterranean waters or transported to other locations by runoff [22]. The fate of Cd in soils will depend on a balance between adsorption, leaching, and biological uptake. In most soils, more than 99% of Cd is found in the solid phase, and less than 1% is present in the soil solution [21].
Earthworm species belonging to the three main ecological groups, i.e. epigeics, endogeics and anecics, have been shown to increase the mobility and bioavailability of some heavy metals in soils [23,24]. As a result, soil ecosystem function can be altered by heavy metal pollution through changes in microbial communities and earthworms in the drilosphere [25]. These changes in microbial community composition may lead to increased heavy metal tolerance in the bacterial community [26,27]. However, a decrease in microbial diversity also leads to less soil resilience to natural or human-induced perturbations [28]. Thus, studies on the joint effect of earthworm engineering and soil pollution on the bacterial community may help to prevent and remediate human impacts at the ecosystem level.
Experimental approach: the use of generalised additive models (GAMs) Our goal is to demonstrate that the joint action of earthworm activity and Cd toxicity induce changes in the bacterial community that are measurable by changes in the CLPP. To this end, we have assessed and measured these changes using a combination of descriptive and multivariate statistics together with GAMs [29,30]. GAMs are semiparametric extensions of generalised linear models that allow for the flexible modelling of covariate effects, avoiding the limitations of parametric models.
This methodology is appropriate when the relationship between the variables is expected to be complex and not easily fitted by standard linear or non-linear models, which is the case in most biological studies. One of the most frequent uses of GAMs in ecology is the modelling of species distribution and species abundance [31], related to environmental and geographical variables [32], being plant ecology and fisheries the two fields where the use of GAMs have been generalized in the last decade. Studies of air pollution and health are another domain where these models have been widely used as standard methods since its flexibility in modelling spatial and temporal processes [33]. One of the great advantages of these models is its ability to determine the nature of the relationship between the response and explanatory variables rather than assuming a fixed relationship, allowing for fitting complex relationships in a straightforward manner.

Soil and earthworm sampling
The earthworm Aporrectodea caliginosa is commonly used in soil ecotoxicological studies [34]. This earthworm is present in the natural soils surrounding our facilities, is easy to culture and sensitive enough against different toxicants, and being a key actor of the matter cycling and soil formation, the results obtained with it are environmentally relevant.
A total of 102 individuals were hand-sorted in an experimental field near the Pyrenean Institute of Ecology (Montañ ana, Zaragoza) in June 2011. Experimental soils were collected in an experimental field near Jaca (northern Spain). No specific permissions were required for these activities since both experimental fields from which the soil and the earthworms were collected belong to the Spanish National Research Council. The activities neither involve vertebrates nor endangered or protected species.
The soil pH was determined from extracts obtained with a 1:2.5 ratio (g soil:g distilled H 2 O [pH w ] or in 0.1 M KCl [pH KCl ]). Organic matter (OM) content, cation exchange capacity (CEC), and soil texture were measured by calcination, atomic absorption spectrometry, and aerometry, respectively [35]. Field capacity (0.3 bar) and wilting point (15 bar) were determined in two pressure chambers [36]. The pH values at the earthworm sampling site were 7.82 (water) and 7.43 (KCl). The prevailing soil conditions (pH, OM content, and CEC) are provided in Table 1 All earthworms remained in a terrarium with the same soil from which they were collected for two weeks (maintenance period). Then, they were moved to another terrarium that contained the experimental soil sieved at , 2 mm. The duration of this preconditioning period was 16 days, and its purpose was to allow the earthworms to empty their guts of the previously ingested soil.

Design and duration of the experiment
We tested the effect of Cd and earthworms on the ability of bacteria to degrade the 31 different sources of carbon included in Biolog EcoPlates. Accordingly, a bifactorial experiment with four treatments was designed (S: soil control; Cd: soil irrigated with Cd; Ew: soil engineered by earthworms; EwCd: soil engineered by earthworms and irrigated with Cd). A detailed schema of the design is provided in Figure 1. Six replicates (pots) per treatment were filled with 1.5 L of sieved soil. The soil surface exposed to air in each pot was approximately 180 cm 2 . The earthworms were exposed to a natural light/dark cycle from June to August. The air humidity as well as air and soil temperatures were controlled. The soil moisture was measured by a Tektronix 1502C TDR Cable Tester reflectometer (Bracknell, Berkshire, United Kingdom) [37]. After the preconditioning period, specimens of a similar size and weight were selected and introduced into each pot (three individuals). Each pot was subsequently watered with distilled water to reach field capacity (i.e., 30% w/w).
All pots were watered with 200 mL of distilled water from Days 1-16. From Day 17 onward, the pots were watered with 200 mL of distilled water or Cd solution depending on the treatment ( Figure 1). Cd was added to the soil as CdCl 2 + 2.5 H 2 O (Panreac), and the pots were watered with a 343 mg L 21 solution. The irrigation was designed to achieve soil concentrations below the median lethal concentration (LC 50 ) to A. caliginosa, which has been estimated in 540 mg Cd g-1 for a 56-day exposure [38]. The percolated water was collected daily 24 h after the watering that generated it and before watering the pots again. The soils were exposed to both Cd and earthworms for 16 days (see Figure 1).  (1), soil bacteria were exposed to earthworms and irrigated with clean water. Later, soils were irrigated with Cd or clean water depending on the treatment (2). At the end of the experimental period, bacteria were extracted by centrifugation and suspended in liquid media (3). This media was inoculated in Biolog EcoPlates (4). The results of colour development over 167 h were recorded and normalised (5). The results were fitted to a sigmoidal model (6), and their parameters were calculated (slope, time to reach 50% of the maximum colour intensity, and maximum intensity reached). Parameters were compared to assess the effects of the different treatments (7), and the carbon sources that were more informative (accumulating more discriminating power) were selected for multivariate analysis (8). Finally, GAMs were fitted to assess the effect of the different treatments on the use of diverse carbon sources over time, determined by the Shannon-Weaver index (9). doi:10.1371/journal.pone.0085057.g001

Measurement of heavy metal
Following Hobbelen et al. [39], three soil sample types were obtained on Day 32. First, pore water was collected by placing 40 g samples of fresh soil into 50 mL Falcon tubes and centrifuging for 40 min (2,000 g). The supernatants were filtered with a Büchner funnel and cellulose filters (0.45 mm pore size; Sartorius Stedim Biotech) and fixed with 69% HNO 3 (Hiperpur; Panreac) at a final concentration of 4% (v/v). Next, the soil was air dried. To obtain CaCl 2 extracts, a 5:1 mixture (5 g of 0.01 M CaCl 2 :1 g of soil) was prepared and shaken for 2 h at 200 rpm. The extracts were filtered and fixed. Finally, the soil samples were acid-digested in a microwave (Milestone digester) by adding 4 mL of HNO 3 (PA-ISO 69%) and 1 mL of H 2 O 2 (33%) to each 0.1 g of dried soil (20 min at 220uC, final temperature). All samples were stored at 4uC until further analysis. The Cd concentration in the soil and percolated water was determined by inductively coupled plasma-optical emission spectrometry (ICP-OES).

Bacterial community-level physiological profiling (CLPP)
The microbial community population structure and changes in carbon source usage were assessed by live/dead staining and Biolog tests, respectively. Bacterial cells were extracted from soils following Katayama [40] and modified as described in Ikeda [41]. Briefly, 10 g of fresh soil was diluted in 95 mL of Milli-Q water. The samples were homogenised for 15 min, and 10 mL of the soil suspension was pipetted into a 50 mL Falcon tube. The sample was sonicated for 1 min and subsequently centrifuged at 1,000 g for 10 min to obtain 9.5 mL of sample. The remaining 0.5 mL was resuspended in 10 mL and submitted to a new dispersioncentrifuge cycle to obtain additional microbial inocula. This cycle was repeated five times until a 47.5 mL sample was obtained, which was refrigerated at 4uC until analysis.
Bacterial extracts were later filtered using 60 mm cellulose filters to reduce interference during OD measurements (for Biolog). Viable and dead cells (live/dead staining) were stained following the LIVE/DEAD protocol [42]. The cells in 10 randomly selected fields per sample were counted using NIS-ELEMENTS software (Nikon ß). The number of viable cells was calculated as totaldead.
For CLPP, integrated soil samples were used (sampling was performed at different depths). One Biolog TM EcoPlate was inoculated per treatment and incubated at 25uC in a dark chamber. Each plate contained three sets (replicates) of carbon sources, and each plate was filled with samples from one pot. An Anthos 2010 microplate reader (Biochrom) was used to read the absorbance in each well at 590 nm (OD 590 ) at even intervals from 0 to 167 h.

Diversity of substrate utilisation (DSU)
DSU was determined by the Shannon-Weaver index (H), which is one of the most commonly used diversity measures [43,44]. This index considers both the number of species and their abundance. In this case, the H value describes the ability of the bacterial community to degrade more or fewer types of carbon sources, thus being an index of physiological diversity of the bacterial community. Microbial communities that are able to degrade more substrates or/and to degrade them with similar efficiency would have higher values of H. The index was calculated as follows: where p i is the normalised value of OD 590 (see point 2.6). This equation is commonly used to process Biolog data [45,46].

Statistical analysis
The carbon degradation data from every well were normalised by being divided by the initial OD 590 at time 0 after averaging the three wells containing the same carbon source. This procedure was performed to reduce the variability due to differences in inocula densities [47,48].
The results of colour development showed sigmoidal relationship between the OD 590 and time for all carbon sources. Therefore, the curves were fitted to a non-linear model [47,49] with the four-parameter log-logistic function given by the formula: Where parameter b is the slope of the curve (S), c is the lower asymptote of the sigmoidal curve, d is the upper asymptote of the curve representing the maximal degradation increase (M), and parameter e is the time to reach 50% of maximum degradation (TM 50 or ED 50 ). T-tests were performed to evaluate which ratios between parameters from two different treatments were significantly (p,0.05) different from 1; p-values were adjusted using Bonferroni correction for multiple test.
For the CLPP analysis, the carbon sources shown to be more informative regarding the effects of heavy metals on bacterial communities were selected. Other sources (such as D-xylose and D-cellobiose) identified in previous studies were also considered [50][51][52]. Canonical discriminant analysis (CDA) was performed on the selected carbon sources (b-Methyl-D-Glucoside, i-Erythritol, D-Mannitol, a-Cyclodextrin, N-Acetyl-D-Glucosamine, a-D-Lactose, and the two previously mentioned).
Thin-plate regression splines [53] were used as smoothers, with optimal effective degrees of freedom chosen automatically using the generalised cross-validation criteria [30].
All analyses were performed with the free statistical software R [54]. The drm and compParm functions of the drc package [49] were used to estimate and compare kinetic parameters, the candisc package [55] was used for CDA, and the gam function of the mgcv library [30] was used for GAM modelling.

Kinetic parameters
The results of colour development revealed sigmoidal relationships between the OD 590 and time for all carbon sources (see point 6 in Figure 1); thus, the results were adjusted to the four parameters log-logistic model [47]. Three parameters of this model were selected for their usefulness in comparing treatments: the maximum value attained (M), the time at which the colour development is half of its maximum value (TM 50 ), and the slope (S). Those parameters presenting significant differences between treatments are shown in Table 2.
This process is the first step of the experimental approach proposed in this study; the carbon sources exhibiting significant differences between treatments (i.e., discriminant sources) would be potential candidates for the subsequent data analysis (see Figure  1) with the aim of increasing the discriminant power of the procedure.
Treatments provoked differences in the kinetic parameters of 23 of the 31 carbon sources. Specifically, the largest differences in CLPP were promoted by the activity of earthworms in Cd-treated soils, which accounted for the differences in the kinetics of 13 carbon sources. On the other hand, the addition of Cd to earthworm-treated soils only resulted in differences in two carbon sources (Table 2).    To assess the effect of the different treatments in the degradation of the different carbon sources, the parameters from the 4 parameters log-logistic models, were divided and the resulting ratio compared to 1 using t-tests. The letter indicates the parameter compared (T = TM 50 , S = slope, M = maximum). Values between brackets are 95% CI. Table shows only those comparisons with significant p-values (adjusted using Bonferroni correction for multiple test), meaning that these ratios significantly (p,0.05) different from 1 and thus the parameters compared were significantly different. As example, the first value in the column S/Ew indicates that the slope of the log-logistic model of degradation for A3 in treatment S was four times lower than in treatment Ew. doi:10.1371/journal.pone.0085057.t002 functions is 97%. The first axis accounts for 87% of the variance and is strongly influenced by D-cellobiose. The second function is strongly influenced by D-cellobiose and is also influenced by Nacetyl-D-glucosamine. Figure 3 presents the CDA based on M values. The cumulative percentage of the first two functions is very high (approximately 99%). In this case, the sources with the most influence are also N-acetyl-D-glucosamine and D-cellobiose for both axes. Figure 4 presents the CDA based on slopes. Although the cumulative percentage is also very high (94%), in this case, the standardised coefficients indicate that the weight of the different sources is more similar than that in the previous CDA.

Canonical discriminant analysis (CDA)
Assessing the time course changes in diversity of substrate utilization (DSU) based on the Shannon-Wiener index using generalised additive models (GAMs) Results from GAM model shows that treatments exhibit significantly different time course changes in DSU, with a explained deviance of 87%. The effective degrees of freedom, chi  Figure 5, the red dotted line represents a hypothetical community with a DSU that was constant as a function of time; in other words, a community exhibiting such a time-course trend would degrade all carbon sources with similar affinity (S) and capacity (M). Thus, this red line can aid in our interpretation of the results. All treatments caused a decrease in DSU during the first 54 h. This decrease was more acute in the bacterial community from the Cd-treated soils. Thereafter, differences were evident between treatments. The communities from soils containing earthworms (Ew and EwCd) presented similar values until the end of the experiment (143 h) regardless the presence of Cd. On the other hand, the community from the Cd-treated soil exhibited an increase in DSU during the next 43 h and did not change thereafter. The community exhibiting the most diverse behaviour was that obtained from control soils (S). In this community, the decline in DSU was constant over time. The model validation did not show any violation of variance homogeneity (i.e. no patterns in residuals) or normality.

Kinetic parameter analysis
Cd exposure resulted in the enhanced degradation of eight carbon sources, as shown by increases in S, M, and/or TM 50 (see column 2 in Table 2). These results indicate the Cd tolerance of soil microbial communities. It has been shown that long-term exposure to Cd (months to years) may cultivate a more tolerant community [26,27], exhibiting enhanced enzyme activities [56]. However, our data indicate that the acclimatisation of the communities to Cd may be achieved after only a few days, which has also been shown in communities exposed to mercury [6]. The low percentage of dead cells in Cd-treated soils also supports this observation.
On the other hand, earthworm activity had a modulating effect on Cd impact. Although earthworm activity allowed for the survival of a higher number of bacterial cells, it also caused reductions in bacterial catabolic activity. As shown in column Cd/ EwCd in Table 2, earthworm activity resulted in lower M values for 12 carbon sources. Because EwCd contained three times more active bacteria than Cd, this result indicates a marked reduction in their enzymatic activity. The reduction of bacterial growth and activity by earthworms has also been recently described in the absence of toxicants [57].

Canonical discriminant analysis (CDA)
In all three CDAs (based on TM 50 , M and S), treatments resulted in different CLPPs, as shown by the centroids. The centroids are the weighted position of each treatment in the multidimensional space created by the eight carbon sources used in the analysis (Figs. 2, 3, and 4); arrows indicate the effect of each source in the discriminatory space. The clear separation between centroids means that these sources gathered a sufficient amount of information to differentiate the effect of each treatment on the CLPP. The CDA illustrates how the CLPP is affected by both Cd and earthworms. Compared with Cd, earthworms had a greater effect on the physiological profile of bacterial communities. The three parameters used (TM 50 , M, and S) allow for a clear discrimination between the effects of the four treatments; however, these parameters have different sensitivities for each experimental factor.
CDA on TM 50 (Figure 2). This parameter allowed for discrimination of the control soils and can be used to detect the interactive effects of Cd and earthworms. The EwCd treatment is characterised by faster degradation of D-cellobiose (G1), acyclodextrin (E1), and b-methyl-D-glucoside (A2), whereas Cd is characterised by faster degradation of N-acetyl-D-glucosamine (E2) and D-xylose (B2).
CDA on S (Figure 4). Slopes also allow for a clear discrimination of Cd impacts. Thus, the Cd treatment appeared clearly separated from the other treatments. The S values of a-Dlactose (H1) and D-cellobiose (G1) were enhanced by Cd.
Sources indicative of specific impacts. The increase in TM 50 for the degradation of D-cellobiose (code G1 in Biolog plates and in Figure 2) appears to be a promising indicator for earthworm activity. Changes in the slope of a-D-lactose are also indicative of earthworm activity regardless of the presence of Cd. In contrast, an increase in S for 2-hydroxybenzoic acid may be indicative of Cd impacts on bacterial-earthworm consortia.

Time course changes in diversity of substrate utilization (DSU) by GAM analysis
GAMs revealed significant differences for all treatments (compared to control -S-) in the time trend of DSU. However, all treatments shown a common initial decline (the first 54 h) in the DSU; that is generally attributed to changes in the community inoculated [11], which becomes less diverse. The CLPP would be dominated by fast growing species that are positively selected by high nutrient levels [4]. Other mechanisms related to metal tolerance and earthworm's effects would explain the differences shown beginning at 54 h and continuing onward.
The acute initial DSU decrease shown by Cd treated soils, suggests the dominance of a few Cd-tolerant microbial taxa [26,27]. That community presented very high metabolic capacities on a few substrata (as indicated by the results shown in Table 2), and that may explain the acute decrease of DSU (i.e. intense utilization of a few carbon substrata). After 54 h, and over the next 43 h the increase of the DSU suggest that other taxa, less dominant, may reach high numbers in the wells containing sources that are less suitable for the Cd-tolerant taxa.
Earthworms were a much more determining factor in explaining time course changes in DSU than Cd. Accordingly, Ew and EwCd soils presented similar trends and values, regardless the presence of Cd. Moreover, both treatments presented similar number of bacteria (between 3 and 3.6e7 cells g 21 DW) and higher than Cd-treated soils (1.1e7). Altogether indicates that earthworm activity provided some protection to the bacterial community against Cd. The bacterial community promoted by earthworms activity during the first 16 days of the experiment, maintained its degrading diversity even after the exposure to a new stressor (Cd). Differently to other studies [58], Cd did not exerted an additional selection pressure. This experimental approach, based on the use of DSU and GAMs allowed for a global and statistically relevant interpretation of the changes in carbon source utilisation, highlighting the key role of earthworms on the protection of microbial communities against the cadmium.