Roots and associated microbes generate acid-forming CO2 and organic acids and accelerate mineral weathering deep within Earth’s critical zone (CZ). At the Calhoun CZ Observatory in the USA’s Southern Piedmont, we tested the hypothesis that deforestation-induced deep root losses reduce root- and microbially-mediated weathering agents well below maximum root density (to 5 m), and impart land-use legacies even after ~70 y of forest regeneration. In forested plots, root density declined with depth to 200 cm; in cultivated plots, roots approached zero at depths >70 cm. Below 70 cm, root densities in old-growth forests averaged 2.1 times those in regenerating forests. Modeled root distributions suggest declines in density with depth were steepest in agricultural plots, and least severe in old-growth forests. Root densities influenced biogeochemical environments in multiple ways. Microbial community composition varied with land use from surface horizons to 500 cm; relative abundance of root-associated bacteria was greater in old-growth soils than in regenerating forests, particularly at 100–150 cm. At 500 cm in old-growth forests, salt-extractable organic C (EOC), an organic acid proxy, was 8.8 and 12.5 times that in regenerating forest and agricultural soils, respectively. The proportion of soil organic carbon comprised of EOC was greater in old-growth forests (20.0 ± 2.6%) compared to regenerating forests (2.1 ± 1.1) and agricultural soils (1.9 ± 0.9%). Between 20 and 500 cm, [EOC] increased more with root density in old-growth relative to regenerating forests. At 300 cm, in situ growing season [CO2] was significantly greater in old-growth forests relative to regenerating forests and cultivated plots; at 300 and 500 cm, cultivated soil [CO2] was significantly lower than in forests. Microbially-respired δ13C-CO2 suggests that microbes may rely partially on crop residue even after ~70 y of forest regeneration. We assert that forest conversion to frequently disturbed ecosystems limits deep roots and reduces biotic generation of downward-propagating weathering agents.
In 2003 a hypothesis was circulated within the scientific community positing that human beings began influencing Earth’s climate with the advent of agriculture ~8K ybp due to effluxes of CH4 and CO2 (Ruddiman, 2003; Ruddiman, 2005). The “Ruddiman hypothesis” has prompted new ways of interpreting fluctuations in atmospheric constituents and stimulated debate among Earth system scientists (e.g. Elsig, et al. 2009; Kaplan, et al. 2010). In contrast to the exploration of human influence on Earth’s atmosphere, relatively little attention has been paid to the deep-soil consequences of long-term human modification of ecosystems at Earth’s surface. We assert that a scientific understanding of Anthropocene signals will be enhanced by circulating and testing a belowground analog to Ruddiman’s hypothesis: that conversion of naturally deeply-rooted ecosystems to cultivation agriculture imparts biogeochemical signatures deep within soil profiles that result from alterations in soil water, gas and solute fluxes, and that soil formation itself has been altered from what would be present sans agriculture even deep within the profile, particularly where cropland has persisted for centuries or millennia. Given the development of science investigating the critical zone (CZ) in recent years – the Earth system extending from the top of canopies to deep within the subsurface (Brantley et al. 2008; Brantley et al. 2017; Wymore et al. 2017) – an appropriate moniker for this hypothesis might be the “Ruddiman-CZ” hypothesis. Ruddiman-CZ designates an embracing of whole-ecosystem (that is, above- and belowground components) responses to land conversion.
This hypothesis is not easy to test. Unlike the original Ruddiman (2003) hypothesis which relies on examination of relatively well-mixed atmospheric concentrations of greenhouse gases trapped in glacial ice over centuries and millennia, the belowground analog requires characterization of paired agricultural and non-agricultural counterparts, with pairs accounting for site-, management-, and plant-specific changes after land conversion to agriculture. Finding multiple locations where these variables can be controlled for any global-scale study is a challenge. However, a first-order step towards probing the Ruddiman-CZ hypothesis can be accomplished by testing a hypothesis narrower in scope and thus more tractable. This becomes feasible if we consider two features common to most agricultural systems growing where forest is the natural land cover: the relatively high frequency of disturbances through cultivation, and the multi-decadal to centuries timescale over which deeply-rooted forests regenerate after agricultural abandonment. Because roots and the microbes dependent on them promote a biological, acid-producing weathering reaction front (Drever, 1993; Maeght et al., 2013; Thorley et al., 2015) that advances downward through ecosystem development, it is reasonable to posit that persistent maintenance of shallowly-rooted ecosystems (i.e. cultivated) where more deeply rooted ecosystems would naturally develop over time will influence the generation of biotic acids and associated soil weathering reactions.
The direct and indirect (via root-dependent microbes) influence of roots on soil weathering stems largely from carbonic acid and organic acids (OAs). Carbonic acid is produced by root- and microbially-derived CO2, and generates acidity with water according to the reaction:
Though a weak acid, carbonic acid can be a potent contributor to soil acidification (Reuss and Johnson, 1986; Amundson and Davidson, 1990; Richter and Markewitz, 1995; Richter and Markewitz, 2001) as soil CO2 concentrations can regularly exceed 5% in deep soil horizons (Richter and Markewitz, 2001; Johnson et al., 2013; Hasenmueller et al., 2015; Lathuillière et al., 2017), driving this reaction forward. Organic acids are also generated by microbes and roots and additionally contribute to soil acidification, as demonstrated by the commonly occurring OA oxalic acid:
Though OAs are typically considered most abundant in surface horizons (Richter and Markewitz, 2001), deep roots and their mycorrhizal symbionts generate OAs that are important agents of mineral nutrient release (Clarkson, 1985; Aoki et al., 2012; Wang et al., 2015). Because rhizospheres exist far below the zone of maximum rooting density (Richter et al., 2006a), the formation of these biotically-derived acids can occur deep within the soil profile (Oh and Richter, 2005; Thorley et al., 2015).
The downward propagation of weathering reaction fronts that are promoted by these root- and microbially-derived acids likely is diminished or even halted by land uses following conversion of naturally deeply-rooted ecosystems. Severe disturbance often kills roots en masse (Hicke et al., 2012), and re-establishment of mature, deep root networks may take longer than the regeneration of aboveground biomass. Plant investments in deep, dense root systems take decades to emerge (Billings, 1936; Zangaro et al., 2008; Knops and Bradley, 2009; Devine et al., 2011; Vickers et al., 2012; Yuan and Chen, 2012) and full regeneration may require many decades to centuries (Sun et al., 2015). Combined with the knowledge that roots and rhizosphere microbes are key soil weathering agents (Maeght et al., 2013) these phenomena suggest that when aboveground disturbances kill roots and limit their regeneration, the production of carbonic acid and OAs by deep roots and their microbial dependents slow or even cease. These effects may be alleviated only via recolonization of deep rooting zones over many decades or centuries (Figure 1). If so, the biogeochemical signatures of aboveground disturbance regimes may be evident far deeper and exhibit far greater persistence than has been assumed.
As a first step towards probing the belowground complement to Ruddiman’s (2003) hypothesis, we hypothesized that frequent cultivation disturbance maintains rooting systems in a persistently shallow state, altering root- and microbially-mediated biogeochemical pools and fluxes important for soil development even well below the zone of maximum root density. Though the presence or absence of roots has important implications for soil physical characteristics such as porosity and the transport of water, microbial substrates, and gases (Scholl et al., 2014), we focus on the propensity for soils harboring root systems of varied ages to generate biologically-derived CO2 and a proxy for biotically-derived OAs, extractable organic C (EOC). EOC is a useful metric of soil weathering agents because it is partially comprised of OAs (Herbert and Bertsch, 1995) and is one metric of mineralizable C substrate (Buscot et al., 2005). Given the long timescales over which forest rooting networks develop (Billings, 1936; Sun et al., 2015), we also hypothesize that effects of agricultural disturbance on these weathering agents and associated parameters will be evident even after many decades of forest regeneration. Unlike studies exploring the influence of fertilizer and pesticide applications in deep soil and groundwater (Kim et al., 2015; Liu et al., 2015) or the influence of land use change on more natural but surficial soil processes (Clements, 1916; Brokaw and Busing, 2000; Richter and Markewitz, 2001; Fornara and Tilman et al., 2009; Kardol and Wardle, 2010; Richter and Billings 2015; Castle et al., 2016; Pulsford et al., 2016), we examine the effects of reductions in deep root density imposed by frequent, long-term disturbance on biotically-derived agents of soil weathering deep in the subsoil. Given the extent of anthropogenic conversion of natural terrestrial ecosystems to annual agricultural crops (up to about 15% of Earth’s terrestrial land surface; Ramankutty et al., 1998; Hooke et al., 2012), deep soil biogeochemical changes due to deep root loss likely has a widespread influence on the biogeochemical environment of weathering profiles observed today.
We leverage replicated plots representing contrasting land use histories and closely related soil types at the Calhoun Critical Zone Observatory (Calhoun CZO) in South Carolina, USA. The work thus invokes space-for-time substitutions across land use histories, while controlling for key ecosystem properties, employing agricultural plots that have supported row-crop agriculture for much of ~200 years and continuously for more than the last 80 years, regenerating pine forests growing for ~70 y on former agricultural land, and old-growth hardwood plots. Though the term “old-growth” can be applied to forests or grasslands (Veldman et al., 2015), here it is used to refer to forests that have had well over a century of root system development. We use this platform to test linkages between root depth distributions and root- and microbially-mediated biogeochemical pools and fluxes relevant to soil development. The work specifically focuses on the availability of soil weathering agents (i.e. CO2 and a proxy for OAs) and multiple indices of soils’ microbial generation of those weathering agents (i.e. C substrate availability and its δ13C values, decay rates, and C mineralization rates). We also investigate whether contemporary soil microbial community composition provides a signal of land use history and the distinct biogeochemical environments it may impose deep within soil profiles. If supported, the hypotheses suggest that widespread, repeated disturbance regimes such as those imposed by annual row crop systems can significantly influence top-down, biotic processes driving evolution of Earth’s critical zone and that these signals persist across multiple decades.
Field sites and soil sampling
The Calhoun CZO (CCZO) is comprised of multiple research areas within several kilometers of 34°36′25″ N, 81°43′13″ W, where contemporary annual precipitation averages 1170 mm y–1 and mean annual temperature (MAT) is 17°C. Research plots at the CCZO represent three historical trajectories of land use change. For two of these land uses, we employed replicated, paired plots within the larger forest complex. We used replicated, “old-growth” (>100 y old) hardwood plots (0.1 ha) to represent old forests growing on soils that have never been plowed. Average tree age at diameter at breast height (DBH) in hardwood plots is 129 ± 21 (standard deviation) y; actual tree age is up to 15 y older. The oldest trees in these plots date from the late 18th century to mid-19th century (W. Cook, pers. comm.). All hardwood stands had intact or nearly intact deciduous canopies visible in 1933 aerial photos taken by the U.S. Forest Service (Brecheisen et al., 2015). Combined with ground-level ecological and geomorphological observations, these images suggest that these stands have experienced limited human disturbance since the late 18th century (Brecheisen and Richter, in prep.). Species composition of these stands is dominated by Quercus alba, Q. rubra, multiple Carya spp., Liquidambar styraciflua, and Liriodendron tulipifera. The currently regenerating pine forests (dominated by Pinus taeda) grow on land originally cleared and cultivated for European-style agriculture and open pastures. After 100 to 150 y of farming and subsequent abandonment, pine stands regenerated via natural regeneration and by planting from the 1930s to the 1950s. Average tree age in the pine plots is 69 ± 4 y, based on known dates of field abandonment and pine planting. Old-growth hardwoods and regenerating pine forests are paired at the same landscape positions (five pairs on upland interfluves and two pairs at mid-slopes). Though comparisons across old-growth hardwood and pine forests may reflect differences in litter and soil chemistry and rooting depths that can occur even between similarly aged hardwood and pine trees, we used these plots as realistic representatives of forest succession in this region.
We also investigated plots cultivated with annual crops (typically sunflower (Helianthus annuus), corn, winter wheat, or sorghum) and exhibiting pronounced Ap horizons. To the best of our knowledge, these plots have been used for cultivated crops for up to 200 y, and aerial photos (Brecheisen et al., 2015) provide robust evidence of these plots being cultivated for more than 80 y. Historically, agricultural crops in this region were dominated by cotton, tobacco, and (to a lesser extent) corn and wheat, and beginning in the 1880s many cropped fields in the region received N, P, and lime inputs. It is impossible to find well-replicated cultivated fields suitable for sampling in this part of South Carolina in the contemporary landscape, so we sampled one complex of fields that covers an area of 8.6 ha. Two of these agricultural plots are located on an interfluve and one is at a midslope position.
Soils were collected during three sampling campaigns in late June 2015, early July 2016, and October to December 2016 (Supplementary Information Table 1). In two sampling campaigns, soils were sampled from 40–50, 100–150, 150–200, 250–300, and 400–500 cm via hand-augering in three paired old-growth forests and pine stands (two interfluves and one midslope), and in the three cultivated plots. During the third sampling campaign, a backhoe excavated soil pits to >200 cm to provide greater replication across the landscape and to develop a quantifiable understanding of rooting depth distributions. These pits were dug in three paired old-growth and pine forests on interfluves, one pair of old-growth and pine forests on a midslope, and in one midslope cultivated field. Within each pit we sampled soils from each genetic horizon, and took high-resolution photographs of a cleaned face for subsequent quantification of root presence (see below). In forested plots, all soil sampling locations were located several meters from trees. All samples were placed in coolers and shipped to the University of Kansas for processing. Upon arrival, samples were stored at 17°C (the mean annual air temperature at the Calhoun CZO), not the more typical 4°C, to ensure that samples were not influenced by temperature outside their normal in situ range. Within hours of arrival, samples were subjected to removal of the generally sparsely occurring roots, hand homogenization, and each biogeochemical assay.
The biota – root density and microbial communities
Photo-derived estimates of root densities within the soil pits were used to quantify differences in rooting system extent and density across the three land uses. On each high-resolution photo of the soil pit faces, we overlaid a grid using Image-J software (Schneider et al., 2012). Each grid cell was 1 cm × 1 cm, sized according to measuring tapes present in each 40 cm wide photo. For each gridded photo, a spreadsheet was generated in which each spreadsheet cell was considered representative of each 1 cm × 1 cm photo grid cell. Presence or absence of a root within each photo grid cell was recorded in each corresponding spreadsheet cell. We counted a root as present in a cell even if it extended from an adjacent cell in which it also was counted. Root density data is expressed as the fraction of each 1 cm row that contains roots.
Microbial genomic data help us understand linkages between root abundances and soil biogeochemical environments. Following manufacturers’ protocols, DNA was extracted from 10 g of each well-mixed soil sample using MoBio PowerMax® Soil DNA Isolation Kit (Carlsbad, CA). Final volumes were concentrated to 100 µL and purified via Agencourt AMPure XP (Beckman Coulter, Brea, CA) beads. The purified DNA was stored at –20°C until further use. Bacterial amplicon PCR was performed on an Eppendorf MasterCycler Gradient thermocycler (Eppendorf, Westbury, NY) using Phusion® High Fidelity DNA Polymerase (New England BioLabs, Ipswich, MA) and primers 341F and 805R, which target the hypervariable regions of the 16S rRNA gene from V1-V3 (Albertsen et al., 2015). Dual indexed libraries were prepared from the amplicons using the Illumina Nextera® XT index kit (Illumina, San Diego, CA). Concentrations of PCR products were measured on Qubit® 2.0 Fluorometer (Life Technologies, Grand Island, NY). Following library preparation and pooling, a paired-end 300 cycle sequence run was carried out on the libraries on an Illumina MiSeq sequencer (Illumina, San Diego, CA) at the Integrated Genomic Facility at Kansas State University, Manhattan, KS. All DNA sequences are publically available through National Center for Biotechnology Information BioProject PRJNA395599 (https://www.ncbi.nlm.nih.gov/bioproject/).
Sequence processing was accomplished via the Quantitative Insights into Microbial Ecology (Qiime v1.9.1) pipeline (Caporaso et al., 2010). First, data were demultiplexed to trim the barcode and primers, to remove low quality sequence reads (Phred score < 20), to truncate end sequences, and to remove all sequences containing ambiguous bases. Using open-reference operational taxonomic unit (OTU) picking in Qiime, singleton OTUs were removed to further de-noise the sequences. Sequences were then clustered into OTUs using uclust (Edgar, 2010) with a 97% sequence identity threshold. PyNAST (Caporaso et al., 2010) was used to align representative sequences from each OTU, and taxonomic identity was assigned using the Greengenes 13_8 database (McDonald et al., 2012). These data permitted examination of multiple measures of α and β diversity (described below).
Solid and extractable soil C pools
Concentrations of soil organic C (SOC) were determined for all soils, as well as the δ13C signature of SOC to characterize the degree to which regenerating forests have transformed C4-influenced pools of SOC in former agricultural plots into C3-dominated SOC pools. Interpretation of δ13C signatures assumes two end members: C3 and C4 vegetation inputs, with approximated δ13C values of –28 and –12‰, respectively. The δ13C signatures are interpreted only semi-quantitatively because though C4 corn and sorghum were and still are periodically planted in Calhoun’s agricultural fields, the historic C4 biomass input to these profiles is unknown. Concentrations of SOC and δ13C of SOC were obtained from soils oven dried at 60°C for >48 h and analyzed on a ThermoFinnigan MAT 253 Continuous Flow System interfaced with a Costech 4010 elemental analyzer at the University of Kansas Keck Paleoenvironmental and Environmental Stable Isotope Laboratory. Instrument precision is greater than ± 0.22‰ for δ13C.
For all soils except one pair of hardwood and pine soil pits, we quantified concentrations of salt-extractable organic C (EOC). EOC is a useful metric of soil weathering agents because it is partially comprised of OAs (Herbert and Bertsch, 1995) and is one metric of mineralizable C substrate (Buscot et al., 2005). For these analyses, five g (fresh weight) of each soil sample were extracted with 25 mL 0.5 M K2SO4 to reveal soil solution OC plus exchangeable, mineral-bound OC (Jones and Willett, 2006; Toosi et al., 2012). Samples were filtered through borosilicate filters (1.5 mm, Sartorius, Germany) and acidified (20 ml 85% H3PO4). At the Marine Stable Isotope Lab at North Carolina State University, extracts were introduced to a total organic C (TOC) analyzer (OI Analytical Aurora 1030D, College Station, USA) interfaced with a Thermo Finnigan Conflo III and a Thermo Electron Delta V Plus isotope ratio mass spectrometer (IRMS).
In situ CO2 concentrations
In spring and early summer of 2015, gas wells were installed by hand-augering to four depths (0.5, 1.5, 3, and 5 m) in each of three paired old-growth and regenerating forest plots and three agricultural plots (SI Table 1). These wells were sampled for in situ [CO2] 39 times, approximately every three weeks from July 31, 2015 to July 26, 2017. We focus on the latter half of the growing season in 2016 and the beginning of the growing season in 2017, months during which no data gaps are present. Gas samples from buried, 750 mL (21 × 6 cm radius) PVC reservoirs were pumped to the surface via minimally diffusive 6.35 mm Bev-a-line tubing (US Plastics, Lima, Ohio). During sampling events, Swagelok-style compression fittings were attached to the tubing and a Vaisala GMP221/GMM221 NDIR CO2 analyzer. Approximately 1.5 L of gas was drawn out to purge stagnant air from the tubing and reservoir and to saturate the gas sampler/analyzer with fresh soil atmosphere. Samples were then drawn into the analyzer.
Ex situ biogeochemical fluxes
Soil samples were subjected to aerobic incubations (17°), during which activities of four extracellular enzymes were quantified following standard protocols. All exo-enzymes assayed are commonly used by soil microbes to decompose soil organic matter prior to substrate uptake and mineralization to CO2: β-1,4-glucosidase (BGase), β-1-4-N-acetylglucosaminidase (NAGase), phosphate-monoester phosphohydrolase (APase), and phenol oxidase (POase). These methods are described in detail elsewhere (DeForest, 2009; Lehmeier et al., 2013; Min et al., 2014). Extracellular enzyme activities are reported as pmol g–1 dry soil h–1.
Microbial efflux rates of CO2 during these aerobic incubations were also quantified. Soil samples (100 g fresh weight) were weighed into Mason jars (0.95 L) and brought to 70% water holding capacity at 17°C. After ~1 h, jars were sealed and we drew initial and 48 h headspace samples (14 mL) and injected them into pre-evacuated glass vials. Headspace samples were injected into a gas chromatograph (Varian CP3800, Agilent) equipped with a thermal conductivity detector. We analyzed replicate gas samples for δ13C of incubation-derived CO2 using a Small Sample Isotope Module connected to a Picarro G2201-i gas analyzer (both Picarro, USA), and corrected the δ13C measurements to account for this system’s CO2 concentration dependence of δ13C.
The hypotheses described above require tests for differences in response variables across depth among land use histories (Supplementary Information Table 2). Because pools of SOC tend to be higher at downslope positions relative to upslope (Papanicolaou et al., 2015; Fissore et al., 2017), it is also important to test for the influence of landscape position where possible. Multiple approaches are required to accomplish these goals (SI Table 2). Except for some measures of microbial community composition, all data sets were first subjected to testing of differences across forest types; separate tests were then performed to assess differences between the pseudo-replicated agricultural soils and the forests (see below).
First, differences in root densities between forest types were assessed, at depths where only forest roots were found – below 70 cm. Root data were binned into 13, 10-cm depth intervals between 70 and 200 cm. Mixed effects ANOVA was used to test for the effect of forest type, depth interval and their interaction on root density, with depth interval considered a fixed effect and forest type as a random effect. To test for the effect of forest type on the shape of rooting depth distributions, root data were fit to a model defined by Arrouays and Pelissier (1994) useful for describing depth distributions of pools of soil organic matter as they decline with depth (Klopfenstein et al., 2015), hereafter referred to as the A-P model:
where b represents the rate of change of root density with depth, d represents depth, z1 and z2 represent the surface and deepest depth, respectively, and C1 and C2 represent root density at those depths, respectively. ANOVA was used to test for any influence of land use on the four modelled parameter estimates generated for each soil profile.
Though a wealth of information is contained in environmental molecular data (Cruz-Martinez et al., 2012; Evans and Wallenstein, 2014), we focused on using genomic data as indices of months if not years (Carini et al., 2016) of dominant soil biogeochemical environments, and whether they exhibited variation with depth and land use history. First, Good’s coverage and rarefaction analysis were used to assess our DNA extraction and sequencing procedures. Four indices of α diversity were calculated for each sample: the Shannon diversity index (H), phylogenetic diversity (PD), Chao1 species richness, and the number of observed OTUs. Variation in these diversity measures with forest type, depth and their interaction was assessed via a mixed effects ANOVA with depth as a fixed effect and forest type as a random effect. Community dissimilarity across all land use histories (both forest types and agricultural soils) was determined using unweighted UniFrac distance. The distances between communities were visualized by Nonmetric Multidimensional Scaling (NMDS) using PRIMER 7 (7.0.5) and the PERMANOVA+ add-on (1.0.3). Canonical Analysis of Principal (CAP) Coordinates (Anderson and Willis, 2003; Jyrkankallio-Mikkola et al., 2016) permitted visualization of the distance matrix, and assessment of statistical differences among microbial community composition with land use history. These analyses also were conducted via PRIMER 7 and the PERMANOVA+ add-on. We explored the relative abundance of a relatively novel phylum, AD3 (Zhou et al. 2003), which is observed in close association with roots (Vik et al., 2013) and thus was of interest as a signal of root-influenced soil volumes. For these analyses, variation in relative abundances of multiple AD3 classes with forest type, depth and their interaction was assessed via a mixed effects ANOVA, with depth as a fixed effect and forest type as a random effect. For those response variables exhibiting a significant interaction, post-hoc protected Welch’s t-tests were conducted. We also examined relationships of AD3 classes’ abundances with EOC concentrations, and assessed whether those relationships varied with land use history.
Because SOC and EOC data reflected a combination of absolute soil depths and genetic soil horizon depths, we generated interpolated, smoothed spline curves following Bishop et al. (1999). These curves permitted contrasts of values within a given depth interval across multiple soil profiles sampled. After visually ensuring that spline curves matched the data, interpolated values were binned into 0 to 25, 25 to 50, 50 to 100, 100 to 300, and 300 to 500 cm depth intervals; a mixed effects ANOVA was used to test for the effect of forest type, landscape position, depth interval, and all possible interactions. Landscape position and depth were considered fixed effects and forest type a random effect. Significant interactions prompted a protected Welch’s t-test to perform post hoc analyses and to assess if the driver of the interaction effect was a contrast of biogeochemical interest. Significant effects of forest type on concentrations of either form of organic C drove assessment of the significance of the relationship between organic C and root densities via regression analysis. ANCOVA was then used to determine if forest type played a significant role in driving the response of organic C to root density, using rooting abundance as a covariate and forest type as a main factor.
In situ CO2 concentrations also were tested for differences between forest types and soil depth. Extensive time series analyses of in situ CO2 and O2 concentration variation with depth and land use history are presented in a forthcoming paper (Z. Brecheisen and D. Richter, pers. comm.). Here, we display heat maps of [CO2] trends observed across depths for the three land use histories for one calendar year, and then focus our analyses on growing season data (18 observations for each site at each depth) at two depths (300 and 500 cm), where [CO2] values commonly reached >4% in the forested plots. A linear mixed effects model (R package lme4) was used to test the influence of forest type and sampling date on [CO2] at 300 and 500 cm, assigning forest type and sampling date as fixed effects and research site and plot as random effects. For each depth, the model was run with and without a forest type*sampling date interaction. We performed a likelihood ratio test using ANOVA to contrast these two models’ output, and used the resulting χ2 statistic and its associated P value to assess the significance of the interaction and to guide model selection. A null model was then constructed, removing forest type as a fixed effect, and a likelihood ratio test performed using ANOVA, this time contrasting our null and working models. Finally, the χ2 statistic from this test and its associated P value were used to infer whether the null model could be rejected and that forest type was a meaningful driver of in situ [CO2].
Variation in the four exo-enzymes, three gas fluxes and δ13C of respired CO2 with forest type, depth and their interaction were tested using a mixed effects ANOVA. For these mixed effects ANOVAs, depth was a fixed effect and forest type was a random effect. For those response variables exhibiting a significant or near significant interaction, we explored the driver(s) of that result via a post-hoc protected Welch’s t-test.
With the exception of β-diversity analyses of the microbial data, for all of the above data sets we assessed whether the relatively less well-replicated agricultural soil data were meaningfully different from old-growth and regenerating forest soil data via t-tests at each depth. For root densities, data were binned into 10-cm intervals down to 70 cm (below 70 cm there were no roots observed in agricultural plots). For SOC and EOC and their δ13C values, data were binned by depth interval as described above for analyses contrasting the two forest types. For root, SOC, and EOC data sets, we also estimated function parameters for agricultural soils using the A-P model, and compared these parameters with those derived for old-growth and regenerating forest soils via t-tests.
With the exception of the NMDS and CAP analyses, all data analyses were conducted in R version 3.4.1 (R Core Team 2017). Except where otherwise indicated, estimates of error are presented as standard errors.
Old-growth forests exhibited root densities averaging approximately twice those of the regenerating forests when considering the entire profile sampled (Figure 2). Below 70 cm, where only forested soils harbored roots, old-growth forests exhibited greater root density (P = 0.053). Root density was also influenced by depth, declining to 200 cm (P < 0.0001). There was no interaction between forest type and depth (P = 0.167).
In contrast to both types of forested plots, we did not observe roots below 70 cm in the agricultural fields. Comparisons using t-tests suggest that agricultural plot root densities were lower than in old-growth forest plots at 10–20, 20–30, 30–40 and 40–50 cm (P = 0.034, 0.032, 0.089, and 0.051, respectively). Agricultural plot root densities were lower than in regenerating forest plots at 10–20 and 20–30 cm (P = 0.019 and 0.035).
A-P modeled estimates of the rate constant describing the rate of change in root density with depth (b) differed between agricultural and old-growth forests (0.043 vs. 0.022 cm–1; P = 0.017), such that the decline in root density with depth was steeper in agricultural soils; in contrast, agricultural and regenerating forest (0.027) soil values of b were similar. Old-growth forest soil values of b exhibited a trend suggesting a less steep decline in root density with depth when contrasted with pine forest b values (P = 0.120).
Microbial community composition
An average Good’s coverage value of 96 ± 2 (a unitless index) indicated that 16S rRNA sequences effectively represented bacterial communities within each sample. Most indices of α diversity, a metric of how diverse a community in an environmental sample is, suggested a decline with depth: across all soils, the Chao 1 index, PD, and the number of OTUs all were significantly influenced by depth (P = 0.031, P < 0.0001, and P = 0.046, respectively). Our mixed effects ANOVA did not indicate any influence of forest type or a forest type*depth interaction, nor did any land use comparison emerge as significantly different when we incorporated agricultural sites into our analyses. However, measures of β diversity, a metric of differentiation in microbial community composition among samples, illustrated that both depth (NMDS analysis, P = 0.117; Figure 3a) and land use (CAP; P < 0.005, Figure 3b) drove differences in microbial community composition.
The AD3 phylum comprised a large fraction of all OTUs in these samples, and in some samples was the largest fraction of all OTUs. Beneath the surface horizon the relative abundance of AD3 class JG37-AG-4 was modestly influenced by forest type (P = 0.111), significantly influenced by depth (P = 0.001), and significantly influenced by the interaction between forest type and depth (P = 0.008; Figure 4a). Driving the interaction between forest type and depth, old-growth forests harbored greater relative abundances of this class than regenerating pine forests from 100 to 150 cm (P = 0.040), and similar trends were observed at two of the other depths assayed. We also observed a significant, positive relationship between EOC concentrations and the relative abundance of AD3 class ABS-6, a class closely related to JG37-AG-4 (Figure 4b), for agricultural plots (linear, P = 0.012, adjusted R2 = 0.57) and for regenerating forests (logarithmic, P = 0.003, adjusted R2 = 0.043). No such relationship was observed for old-growth forests.
Soil organic C and extractable organic C concentrations and δ13C
Concentrations of SOC were relatively low in all mineral soil profiles, which we attribute to these soils providing little protection for soil organic matter from attack by microbially produced exo-enzymes (Richter and Markewitz 2001). Concentrations of SOC ranged from undetectable to 67 mg g–1 and declined with depth for all land use histories and landscape positions (P < 0.0001, Figure 5a). There was no influence of forest type, landscape position, or any interaction on SOC, and we observed no significant effect of land use history on SOC when incorporating agricultural soils into the analyses.
Concentrations of EOC ranged from undetectable to 277 mg g–1 and declined with depth (P < 0.0001, Figure 5b). We observed trends hinting at an influence of forest type (P = 0.110) and landscape position (P = 0.090) on EOC concentrations, likely reflecting [EOC] at the deepest depths where old-growth forest soils exhibited 8.8 times the [EOC] of regenerating forest soils. The landscape position effect suggests that midslope sampling locations possessed greater [EOC] than those on interfluves (64.9 ± 13.3 vs. 43.8 ± 6.26 mg g–1). When we incorporated agricultural plot data into our analyses, we observed lower [EOC] in agricultural soils than in old-growth soils (P = 0.046), and a similar trend for regenerating forest soils (P = 0.117), both at 25–50 cm. At the deepest depth, old-growth forest soils exhibited 8.8 times the [EOC] of regenerating forest soils and 12.5 times the [EOC] of agricultural soils.
We also examined the fraction of SOC comprised of EOC (EOC/SOC). We employed this estimate of the ratio of OA or other labile C pools to total SOC to investigate potential differences across land use histories in belowground C allocation. We observed a significant interaction effect of forest type and depth interval on EOC/SOC (Figure 5c; P < 0.0001), driven by higher EOC/SOC in old-growth forests (20.0 ± 2.6%) compared to regenerating forests (2.1 ± 1.1) between 400 and 500 cm (P = 0.006). Analyses including agricultural soils revealed EOC/SOC ratios in agricultural soils (1.9 ± 0.9%) significantly lower than those in regenerating (P = 0.002) or old-growth (P = 0.0002) forests. We assessed the relationship between root density and [EOC]; all land uses revealed a significant, positive relationship between these parameters. More specifically, ANCOVA revealed a significant interaction between forest type and root density at depths greater than 20 cm (P = 0.047), such that [EOC] increased with root density to a greater extent in old-growth hardwood plots than in regenerating forest plots (Figure 6).
Values of δ13C of SOC ranged from –20.1 to –29.8‰ (Figure 7a) and exhibited no differences between forest types nor an influence of landscape position or any interaction, and only an influence of depth (P < 0.0001). However, when we incorporated agricultural soils into our analyses, we observed a significant 13C-enrichment of agricultural SOC relative to old-growth (P = 0.002) and regenerating (P = 0.0004) forests from 0 to 25 cm, and a similar effect from 25 to 50 cm (P = 0.001 for both old-growth and regenerating forests) and 50 to 100 cm (P = 0.051 and 0.084 for old-growth and regenerating forests, respectively). Thus, in the agricultural fields SOC δ13C values reflect the centuries-old history and contemporary planting of C4 crops. In regenerating forests, SOC δ13C values reflect the dilution of historic C4 crops with contemporary, C3 forest organic matter.
The δ13C of EOC ranged from –20.8 to –30.0‰ (Figure 7b) and exhibited an interaction effect between forest type and landscape position (P = 0.055) driven by old-growth forest δ13C-EOC values being 13C-deplete relative to those in regenerating forests (P = 0.053) by approximately 0.7‰ on interfluves. When we incorporated agricultural soils into our analyses, we observed that agricultural soil δ13C-EOC values were significantly more 13C-enriched than old-growth (P = 0.004) and regenerating (P = 0.004) forest δ13C-EOC values from 0 to 25 cm, and 13C-enriched relative to old-growth forests from 25–50 cm (P = 0.060). From 50 to 100 cm, agricultural δ13C-EOC values were significantly 13C-enriched relative to old-growth forests (P = 0.007), and old-growth forests exhibited a trend of 13C-depletion relative to regenerating forests (P = 0.110). All of these differences were between 2 and 4‰ in magnitude.
CO2 concentrations in situ
Heat maps of in situ soil CO2 illustrate variation across months, land use, and depth: values generally increased with depth and reached peak values late in the growing season (Figure 8a). Values ranged from just above aboveground atmospheric concentrations to more than 7%. At 300 cm, we observed a significant effect of forest type on in situ [CO2]: Old-growth forest values were greater than regenerating forest values (5.49 ± 0.24 vs. 4.07 ± 0.12%, P = 0.042; Figure 8b). At 500 cm, we observed no strong effect of forest type (P = 0.137). Averaged across all sampling dates, t-tests suggested that agricultural soil [CO2] was lower than old-growth forest [CO2] at 300 (2.1 ± 0.1 vs. 5.5 ± 0.8%; P = 0.009) and 500 (1.9 ± 0.1 vs. 5.5 ± .8%; P = 0.014) cm, and lower than regenerating forest values at those depths as well (2.1 ± 0.1 vs. 4.1 ± 0.4% at 300 cm and 1.9 ± 0.1 vs. 4.0 ± 0.4% at 500 cm; P = 0.062 for both depths).
Rates of microbial CO2 production and organic matter decay
Microbial efflux rates of CO2 exhibited a significant influence of depth (P < 0.0001, Figure 9a). Rates ranged from 1.3 to 131.6 ng CO2-C g–1 h–1. No effect of forest type on microbial CO2 efflux rates was observed. δ13C of respired CO2 ranged from –30.4 to –14.1‰; though rates of CO2 generation were low, the CO2 produced exhibited an influence of forest type (P = 0.020) and depth interval (P = 0.100; Figure 9b). The effect of forest type reflects the lower δ13C of respired CO2 from old-growth forest soils (–27.0 ± 0.8‰) compared to regenerating forest soils (–23.5 ± 0.7‰). Because δ13C-CO2 data were unavailable from several agricultural samples, we were unable to assess whether δ13C of microbially respired CO2 in agricultural soils was significantly different from either forest type.
During the aerobic incubations, APase exhibited decreasing activity rates with depth (P = 0.005) and no influence of forest type or a forest type*depth interaction (Figure 10). In spite of the apparent shift in APase activity between agricultural soils and the two forest types, limited replication of agricultural plots limits our abilities to reject the null hypothesis that there is a land use effect of APase activities in these soils. BGase (P = 0.004) and POase (P = 0.067) activities also declined with depth, but neither exhibited an influence of any other factor; activities of NAGase exhibited no influence of any factor tested (data not shown).
We have posited a belowground complement to Ruddiman’s (2003) hypothesis, the Ruddiman-CZ hypothesis: that the advent of agriculture has modified natural biogeochemical signals and associated weathering status more deeply within many of Earth’s soil profiles than is typically appreciated, especially where agriculture has persisted for many years. As a first step, we tested if persistent, disturbance-induced absence of deep roots alters biogeochemical signals that are 1) linked to generation of soil acidity (e.g. in situ CO2 concentrations, ex situ CO2 production rates, organic C substrates and a proxy for OAs) and 2) revealing of dominant biogeochemical conditions (e.g. exo-enzyme activities and microbial community composition), even well below the zone of relatively high rooting density (~30–40 cm; Figure 1). We also tested if these effects persist even after decades of forest regeneration.
Well-developed rooting systems pump photosynthate deep into soil profiles
Depth distributions of Calhoun CZO root systems demonstrate how disturbance regimes and associated ecosystem age can influence soil profiles (Figure 2). The great depths to which tree roots extend (Canadell et al., 1996; Jackson et al., 1996; Maeght et al., 2013; Freycon et al., 2014; Fan et al., 2017) compared to those in most annual croplands (Acuna and Wade, 2013; Ahmadi et al., 2011; Gan et al., 2011; Neykova et al., 2011; Trachsel et al., 2013), and the centennial timescale over which forest roots continue to expand (Zangaro et al., 2008; Knops and Bradley, 2009; Devine et al., 2011; Yuan & Chen, 2012; Sun et al., 2015) are widely observed phenomena. However, to our knowledge the current study is the first to quantify root densities to such depths in replicated plots representing distinct disturbance histories, and to determine a meaningful influence of forest age on rooting depth distributions. We further demonstrate multiple biogeochemical consequences of the greatly reduced root activity below 70 cm in agricultural compared to forested soils, and of the expanded influence of old-growth hardwood vs. regenerating forest root systems. It is well-documented that agricultural management can influence profile chemistry down to groundwater (Kim et al., 2015; Liu et al., 2015), and centuries of agricultural management and associated fertilizer application can leave a long-lasting legacy of nutrient enhancement in soils (e.g. Richter et al., 2006b). In contrast, we demonstrate changes to more naturally-occurring processes as a result of diminished root densities below the zone of maximum rooting density.
The microbial community composition data prompt three observations (note that a more detailed treatment of these data is presented in a forthcoming manuscript). First, we demonstrate a decline in α-diversity of microbial communities with depth similar to multiple studies (Fierer et al., 2003; Eilers et al., 2012; Gabor et al., 2014). Similarities in α-diversity trends with depth across studies give us confidence that such observations may be generalized across diverse landscapes. Second, compositional differences reflected in β-diversity measures (Figure 3) suggest distinct biogeochemical environments in soils subjected to distinct land use histories, even to 500 cm. If these differences were driven by contemporary land cover only, we might expect genomic data from shallower horizons to be more divergent in ordination space than deeper horizons. We observed no such trend, and instead observed similar divergence between land use histories regardless of depth. Periodic fertilizer application likely influences microbial composition in the agricultural plots (Mitchell et al., 2016; Gupta et al., 2016), but microbial compositional distinctions between old-growth and regenerating forests at depth suggest that distinctions in root-influenced environments in deep soils are persistent even after ~70 y of forest regrowth. These results hint that even the deepest microbial communities assayed are reflective of long-term, centuries- and decades-long land management choices, and serve as a reminder of the great depths to which microbes can serve as signals of access to substrates and CZ environmental conditions (Akob and Küsel 2011).
Third, microbial data hint at root densities in these soils far more deeply (500 cm) than the depth to which we directly quantified root densities (200 cm). The AD3 phylum occurs within plant roots (Vik et al., 2013). Consistent with this observation, AD3 comprised more than 70% of bacteria in surface soils, and both classes of AD3 declined with depth in a manner qualitatively similar to root declines with depth (Figure 4a depicts this for class JG37-AG-4). A trend of greater relative abundance throughout soil profiles of the AD3 class found in roots (JG37-AG-4, Vik et al. 2013) in old-growth compared to regenerating forest soils implies that root density varied with forest type down to 500 cm – far deeper than our root counts – and that the presence of roots and their exudates promotes biogeochemical environments distinct from those where roots are relatively scarce. We suggest that deep within soil profiles, where quantifying root biomass is an especially difficult task (Maeght et al. 2013), the abundance of some classes within the AD3 phylum may be a useful index of root density. Our observation that AD3 class ABS-6 increases with soil [EOC] in agricultural and regenerating forest soils, but not in old-growth forest soils (Figure 4b), highlights further distinctions between old-growth systems and those subjected to disturbance. This observation is consistent with the idea that abundances of some AD3 classes are limited by EOC to a greater degree in soils experiencing land use conversion from deeply-rooted systems to agriculture, even after ~70 y of forest regeneration. Given the time-integrated nature of soil-derived DNA (Carini et al., 2016), we infer that these observations reflect robust (i.e. not snapshot) distinctions in dominant biogeochemical environments in these soils.
Biogeochemical consequences of altered root densities with land use were also evident via concentrations and isotopic signatures of soil C pools. As in multiple other studies (Jobbagy and Jackson, 2000; Chantigny, 2003), SOC and EOC concentrations declined with depth in all soils (Figure 5a, b). δ13C of SOC indicated 13C-depletion in both types of forests relative to agricultural plots (Figure 7a), providing evidence of pine organic matter incorporation into soil profiles masking previous inputs of 13C-enriched, C4 biomass. These data are consistent with previous work indicating rapid turnover of SOC in these soils (Richter et al., 1999). The ratio of EOC to SOC (Figure 5c) also indicates a key difference between agricultural and forested profiles: the proportion of EOC comprising SOC increases as root density changes from complete absence in the agricultural plots below ~70 cm to a small but persistent presence through 2 m in the forested plots.
Consistent with the idea that degree of root colonization matters for deep soil C cycling (Maeght et al., 2013), our study also reveals distinctions across forest type in the amount and type of photosynthate pumped to depth. At the deepest depth assayed, [EOC] in old-growth forests was greater than in regenerating forests and agricultural plots, and the greatest EOC/SOC in old-growth forest soils: [EOC] and EOC/SOC in old-growth forests were ~9 and ~10 times those in regenerating forests, respectively, at 400 to 500 cm (Figure 5b, c). This result suggests that EOC pools are sensitive signals of root density and are more responsive to root density in old-growth systems (Figure 6). Though beyond the scope of our study, these results could also reflect changing hydrologic flow paths with altered root density, and highlight a need for examining the soil structural consequences of frequent and persistent disturbance regimes that prevent deep root development. Indeed, our observation of higher [EOC] at midslope positions compared to interfluves emphasizes the sensitivity of soil solutes to flow paths.
Enhanced [EOC] and EOC/SOC ratios in old-growth forest soils at depth are congruent with roots being a key source of sugars and organic acids (Dakora and Phillips, 2002; Aoki et al., 2012; Hunter et al., 2014; Karst et al., 2017), and highlights a potentially important mechanism by which roots can drive soil biogeochemical environments and, ultimately, soil development via acid generation. We do not know what comprises EOC in these soils. However, root sugar exudates can be 13C-deplete relative to biomass (Mary et al., 1992). As a result, the generally 13C-deplete EOC in old-growth soils relative to regenerating forests and the isotopic excursion of EOC between 100 and 200 cm in old-growth forest soils to generally more negative values (Figure 7b) is consistent with old-growth roots exuding sugars to a greater extent than in regenerating forests. This, in turn, may signal a lower investment in more 13C-enriched exudates such as OAs in yet deeper horizons in old growth soils. Variation with depth in allocation of distinct photosynthetic compounds may be a strategy for forests possessing deep, well-networked root systems. Organic acid exudation may be an important strategy for deriving mineral nutrition in deep horizons where reactive mineral surface area is greater and geogenic nutrients are relatively more available (Brantley et al., 2012). In contrast, sugar exudation fueling microbially-mediated organic matter decay and microbially-derived carbonic acid may offer more nutrition in shallower horizons, where organic matter concentrations and heterotrophic microbe abundances tend to be greater (Fierer et al., 2007). Regardless of whether EOC represents OAs or substrates for microbes to mineralize to CO2, we emphasize that less pronounced pumping of EOC into deep soil horizons by agricultural and regenerating forest systems represents deep soil modifications of soil weathering agents.
The in situ CO2 concentrations are similar to those reported elsewhere (Selker et al., 1999; Hasenmueller et al., 2015), and their depth distributions are similar to those measured and modeled near the current study (Oh and Richter 2005). Those in the current study suggest that land use history is an important determinant of belowground CO2 production (Figure 8). The lack of deep roots in agricultural plots drives a decline in soil [CO2] in non-surficial horizons compared to forested plots, a strong indication that root loss limits photosynthate pumping into deep soils and curtails the development of an acid weathering front. More nuanced but equally important, the consistent pattern of greater growing season [CO2] in old-growth compared to regenerating forests, particularly at 300 cm, also suggests that old-growth hardwood forests generate greater quantities of CO2 during the growing season than regenerating forests. If robust across the landscape, these data suggest that ~70 y of forest regeneration can bring soils’ CO2 depth distributions to an intermediate state between that of agricultural and old-growth forest soils. If we consider the old-growth forest soils as a baseline, the data speak both to the strong degree of recovery in soil CO2 imparted by ~70 y of forest regrowth following agricultural abandonment and, simultaneously, to the meaningful legacy of historic land conversion on in situ [CO2] deep in the subsurface, even after decades of forest regeneration.
Microbially-driven δ13C-CO2, but not microbial CO2 efflux or decay rates, reflect land use
We observed no difference in lab-incubated CO2 efflux rates across land uses (Figure 9a). The contrast between this lack of a land use effect on microbial CO2 production rates and an evident effect of land use (agricultural vs. forest) on in situ [CO2] suggests the land use effect on in situ CO2 may be primarily driven by root respiration, and less so by microbes. However, microbes served as a sensitive indicator of land use history in a more nuanced fashion. δ13C values of ex situ (i.e. microbially) respired CO2 at all depths suggest that substrates mineralized to CO2 in old-growth forest soils were comprised of C3 photosynthate to a greater degree than those in regenerating forest soils, which likely supported corn in the past (Figure 9b). This phenomenon is consistent with EOC pools in old-growth forests being more 13C-deplete at some depths than in regenerating forests, and with inferences derived from [EOC] and EOC/SOC ratios that greater root density in old-growth forests drives a degree of photosynthate pumping deep into soil profiles. We cannot know why, if [EOC] are higher at greater depths in old-growth forest soils compared to regenerating forest soils, we did not observe correspondingly higher CO2 efflux rates. Microbial activities may have been inhibited by OAs (Rusu et al., 2016) if they, and not sugars, comprised a meaningful proportion of EOC.
We also did not observe any significant effect of land use on any of the exo-enzymes assayed, in spite of activity declines with depth consistent with other studies (Stone et al., 2014; Loeppmann et al., 2016) and strong visual cues that APase activities were lower in agricultural soils than in forested (Figure 10). Previous work in these regenerating forests reports tree P demand being satisfied in part by both organically-derived P and Ca-associated P likely derived from historic P fertilization and liming (Richter et al., 2006b). Data in the current study indicate that rates of organically-bound P release are similar between forest types, hinting that microbially-released, organically-bound P may serve an equally important role in forest growth in both old-growth and regenerating stands.
Loss of deep roots drives changes in the deep soil biogeochemical environment. The duration of these changes is dictated by the centennial timescale over which roots regenerate. We specifically demonstrate that persistent, frequent disturbances imposed aboveground can alter deep root densities, microbial community composition, the sources of microbially-generated CO2, EOC abundance and contributions to the SOC pool, and in situ [CO2], all well beneath zones of highest rooting densities. We also demonstrate the value of microbial community composition data as integrative signalers of biogeochemical environments, and specifically that the abundance of some classes of the bacterial phylum AD3 may be useful as an index of root density where root biomass measurements are impractical. Though some differences between root depth distributions and the biogeochemical phenomena linked to them may be related to differences in plant species (i.e. pine vs. hardwoods), we emphasize that plant species variation with ecosystem succession is an integral feature of ecosystem regeneration. We thus highlight the importance of top-down, biotic processes as drivers of the soil biogeochemical environments that promote soil development deep within Earth’s critical zone. Because frequent disturbance of formerly deeply-rooted ecosystems is a hallmark of human activities, we suggest that further exploration of a lithosphere-centric complement to Ruddiman’s atmospherically focused hypothesis – the Ruddiman-CZ hypothesis – seems a worthy task, across diverse land conversion scenarios. Such work can result in more accurate parameterization of Earth system models, and can help quantify the influential role of humans in the evolution of Earth’s critical zone in the Anthropocene.
Data Accessibility Statements
All DNA sequences are publically available through National Center for Biotechnology Information BioProject PRJNA395599 (https://www.ncbi.nlm.nih.gov/bioproject/).
We are grateful to Will Cook and Paul Heine for assistance in the field, to Dr. Bruce Barnett, Dr. Roxane Bowden, and Carl Heroneme for assistance in the laboratory, and to the Kansas Biological Survey for logistical support.
This work was funded by the National Science Foundation grants EAR-1331846 and DBI-1262795.
The authors have declared that no competing interests exist.
Contributed to conception and design: SAB, DdeBR and PLS
Contributed to acquisition of data: SAB, CAL, SB, KM, ZB, EH, RS, RF, DdeBR
Contributed to analysis and interpretation of data: SAB, DH, PLS, CAL, SB, KM, DdeBR
Drafted and/or revised the article: all
Approved the submitted version for publication: all