Carbon–Mercury Interactions in Spodosols Assessed through Density Fractionation, Radiocarbon Analysis, and Soil Survey Information

Soils comprise the largest terrestrial pool of C and Hg on Earth, and these elements have critical feedbacks to problems ranging from atmospheric pollution and climate change to public health. Empirical evidence suggests these elements cycle closely in a wide range of soils, but mechanistic studies of their interactions within distinct soil organic matter (SOM) pools and between different soil types are needed. Here, we report findings of a novel approach to investigate C–Hg interactions, primarily in Spodosols, in which we: (i) examined density separated topsoil and illuvial horizons of four contrasting Spodosols, and used radiocarbon to investigate interactions between Hg and C cycling in distinct SOM pools; (ii) assessed broader patterns across Spodosols and other soil orders using USDA soil survey laboratory data. Consistent with other studies, C and Hg concentrations of individual soil horizons were positively related across the four contrasting Spodosols. Carbon and Hg were also positively related in the density fractions comprising individual soil horizons, but radiocarbon analysis revealed fundamental differences in Hg retention in modern, C-rich fractions vs. low-C fractions containing less modern radiocarbon. The lack of significant site-to-site variation in C and Hg across these sites (and Spodosols more broadly), contrasted against significant differences between horizons and fractions, suggests processes controlling C–Hg interactions are consistent across the taxonomic order. Furthermore, significant differences between other soil orders indicate that processes controlling soil formation—as represented by soil taxonomy—can explain differences in C–Hg interactions and their distribution across soils. L.E. Nave* Univ. of Michigan Biological Station Pellston, MI, 49769 and Univ. of Michigan Dep. of Ecology and Evolutionary Biology Ann Arbor MI, 48109

A t spatial levels ranging from individual ecosystems to all of Earth, soils are the largest terrestrial storehouses of C ( Jobbagy and Jackson, 2000) and Hg (Bank, 2012). These two elements play very important, but in many ways disparate roles in soils and ecosystems. On one hand, organic C is a metabolic currency, the skeleton of biomolecules, and an exchange site for elements ranging from nutrients to toxic metals, making it centrally important to many biogeochemical cycles. On the other hand, Hg is a toxic, reactive, bioaccumulating metal with a highly mobile gas phase, having deleterious impacts on organismal, soil, and ecosystem processes at molecular to global levels (Driscoll et al., 1994;Selin, 2009). The closely coupled cycling of these two elements and their links to pervasive problems-from atmospheric pollution to public health-argue for an integrated approach to their study.
At regional down to landscape levels, amounts of C and Hg in soil (whether concentrations or stocks) co-vary closely. Mercury can be naturally more abundant in soils in volcanically active regions and where sedimentary rocks are rich in organic C (Turekian and Wedepohl, 1961). In most regions, however, anthropogenic atmospheric deposition (whether modern or legacy) is a significant, if not always dominant, source of Hg in soils (Fitzgerald and Lamborg, 2014;Schuster et al., 2002). Deposited via dry and wet mechanisms (e.g., Risch et al., 2017), Hg in soils is also "lost" to the atmosphere via volatilization (such as through fire), and to receiving water systems via runoff and recharge (Grigal, 2003). On timescales of centuries, soils are net sinks for Hg, especially when complexed with organic matter (Skyllberg et al., 2006;Yu et al., 2014). Sorption to SOM and the subsequent hydrologic redistribution of soluble organomineral complexes explains some of the distribution of Hg across soillandscape units (Demers et al., 2007;Driscoll et al., 2013). In catenas or toposequences, soils in low-lying landscape positions are typically net accumulators of C and Hg, because in part of in situ accumulation of both elements, as well as inputs from contributing source areas (Nave et al., 2017;Kolka et al., 2011). Because Hg and SOM solubility is strongly influenced by soil pH, and the mobility of both elements depends on soil texture and hydrologic fluxes, these soil chemical and physical factors influence C and Hg distribution (Demers et al., 2013;Mitchell et al., 2009). These factors-soil pH, texture, and hydrology-are important predictors of C and Hg distribution over landscape to regional levels, often through relationships with other factors such as atmospheric deposition, soil clay content, and soil drainage regime or annual precipitation (Obrist et al., 2011;Richardson et al., 2013;Stankwitz et al., 2012). In turn, many of these same factors explain the landscape to regional distribution of SOM, of which C and Hg are constituents (Rasmussen et al., 2018).
The vertical distribution of SOM within soil profiles and its variation across space are often studied to assess models of SOM stabilization. Research on SOM stabilization in mineral soils has advanced considerably in recent years, as conceptual models have become better integrated with laboratory and analytical methods that can directly measure properties of interest (e.g., composition, conformation, and stability of C forms). There is now abundant evidence that SOM persistence is typically not because of inherent recalcitrance or acquired through processes such as "humification"; instead, incorporation into aggregates and adsorption on mineral surfaces are viewed as the dominant processes that result in stable SOM (Kleber et al., 2011;Lehmann and Kleber, 2015). These functionally discrete, physically distinct pools (aggregates and organo-mineral complexes) can be reliably separated by density fractionation (Golchin et al., 1994;Swanston et al., 2005;Crow et al., 2007); when combined with radiocarbon analysis as a measure of C stability (Frank et al., 2012;Trumbore, 2009), these methods reliably isolate SOM pools differing in their proximal sources (e.g., plant detritus, microbial residues), functional composition (e.g., aliphatics, aromatics, alkyl C), and turnover times (Swanston et al., 2005;Schrumpf et al., 2013;Trumbore, 2000).
Conceptual models of SOM stabilization have many similarities to models of soil formation, in part because SOM is critical to so many properties of and processes occurring in soil. The Simonson (1959) model treats soil-forming state factors (sensu Jenny, 1941) as the constraints that influence four soil forming processes: gains, losses, transfers, and transformations. Forming through these four processes, soils are dynamic natural bodies that exhibit continuous variation in their myriad properties. These include properties such as C stocks, Hg depth distributions, and the distribution of SOM among functional pools, and soil taxonomy provides a tractable framework for categorizing variation in these properties and testing how soil formation, SOM stabilization, and C-Hg interactions differ among soils.
Carbon and Hg are closely linked in soils and play different roles in organisms, soils, and ecosystems. Spodosols are particularly well suited for studying C and Hg and their links to soil formation, because their genesis is readily interpreted according to the process model of soil formation (Simonson, 1959), and their mechanisms of SOM stabilization are fairly well known from detailed site-level studies (Kaiser et al., 2002;McFarlane et al., 2013;Schulze et al., 2009). Furthermore, mechanisms controlling Hg speciation and interactions with SOM have been studied in Spodosols across a wide geographic range (do Valle et al., 2005;Gomez-Armesto et al., 2018;Pena-Rodriguez et al., 2014;Roulet et al., 1998). However, few studies meld pedogenesis, geochemical detail, and broad patterns across sites (e.g., Richardson et al., 2013), and as yet no researchers have reported Hg data for soil density fractions. Our goal in the present study was to close this data gap and use radiocarbon as a direct measure of SOM stability to understand how C and Hg interact during their transfers and transformations in Spodosols. Through a detailed cross-site study of four disparate Spodosols, and a broader analysis of the USDA-Natural Resources Conservation Service (NRCS), National Cooperative Soil Survey (NCSS) Soil Characterization Database, we tested five hypotheses: (i) The positive relationship between C and Hg in bulk soils also exists in the density fractions that comprise those soils; (ii) horizons and density fractions in which the C is more modern have higher C and Hg concentrations; (iii) Hg/C ratios narrow in increasingly modern SOM pools (horizons and density fractions); (iv) between-horizon and between-fraction variation in Spodosol C and Hg properties is more significant than between-site variation; (v) Spodosols differ in C and Hg properties and relationships from other soil taxonomic orders based on the pedogenic processes and physical and chemical properties that make them unique.

Soil Sampling
At the NEON sites (BART, DSNY, UNDE), the NEON team used a corer (9100 series Power Probe, AMS Inc.) to collect five pairs (n = 10) of complete soil profiles (i.e., two cores from within each of the five NEON soil array plots). One core of each pair (i.e., five cores per site) from each of these sites are described and reported in this analysis. At UMBS, 15 complete soil profiles were collected (in 30-cm increments) from a 20-ha area using a split-core sampler. Individual horizons from 7 of these 15 profiles were described and density separated for this analysis; density fractions from the two with median bulk soil physical and chemical characteristics were analyzed for radiocarbon. Field sampling of soils at the four sites occurred between February (DSNY) and July (UMBS) 2015, and involved the collection of full soil profiles (into the C horizon) in stainless steel corers with butyrate plastic liners. Upon collection, intact cores (in 1-m increments in their liners) from the NEON sites were capped, placed in a cooler, and shipped on ice to the Oregon State University Core Marine Geological Repository Laboratory, where they were stored (as long as 8 wk) at 4°C until processing. At UMBS, individual genetic horizons were measured (thickness) and described (moist Munsell color, texture class) in the field, then quantitatively removed from their core liners and transferred to 4 mil HDPE bags, which were placed in a cooler. At the end of the day of sample collection, these samples were placed on cafeteria trays and air-dried at ~20 to 25°C for 1 wk, reaching <5% volumetric moisture content during that time.

Laboratory Processing
Throughout this paper, we report data for three pedogenically distinct horizons in each profile: the uppermost non-litter horizon of each soil (Oa, A, or AE depending on site), the uppermost illuvial (spodic) horizon, and the deepest B. First, to characterize the portion of the profile that is most influenced by surface detritus inputs and atmospheric Hg deposition, we report data for the uppermost non-litter horizon of each soil: the Oa, A, or AE horizon (depending on site), to which we will refer as the O/A horizon from this point forward. Because the eluvial horizons of Spodosols are characteristically depleted in SOM, C, and Hg, we do not present data for E horizons. Instead, the uppermost illuvial (spodic, i.e., B1) horizon is the next horizon of interest. Lastly, to contrast the properties of this uppermost illuvial horizon with those of the illivual horizon that is least altered from parent material, our third portion of interest is the deepest B horizon within each profile.
Soil cores from NEON sites (BART, DSNY, UNDE) were carefully cut open to minimize plastic contamination (petroleum-derived plastic can bias radiocarbon measurements) from the core liners. Soils were described according to NRCS procedures (Schoeneberger et al., 2012). Because of the narrow internal diameter (3.5 cm) of the cores some soil characteristics could not be described (i.e., horizon boundaries, very coarse roots, and very coarse soil structures); however, a separate 100+ cm soil pit was excavated at each site in conjunction with the NEON site installation. When in doubt, thick horizons contained in the cores were split into multiple horizons to provide more detailed descriptions of the soils and greater flexibility in later compositing of horizons for density separation and analysis. High resolution photos were taken of the cores, individual horizons were collected, weighed, air dried, passed through a 2-mm sieve to exclude roots and coarse fragments, and subsequently weighed again. After describing the five cores from each site, site-level composite horizons were decided based on the combination of the most common horizons from the cores, and the most representative horizons compared with the soil pit dug on site. Fine earth fractions of individual horizons present across the cores were composited and then split with a rifflebox to ensure homogeneity of sample splits sent to subsequent laboratories. For example, at UNDE, all five of the cores possessed A horizons, while only three cores possessed Bs2 horizons; thus, the final number of discrete samples that were aggregated into a site-level composite for further processing and analyses are n = 5 and n = 3, respectively, for the O/A and deepest B horizons at UNDE (Table 2). At UMBS, air-dried samples were sieved (2-mm mesh) to isolate pebbles and gravel (>2 mm), coarse root (>2-mm diam.), fine root (<2-mm diam.), and fine earth fractions; these fractions were subsequently oven dried (60°C for rock, coarse and fine root; 105°C for fine earth) before making final measurements of mass on each to the nearest 0.01 g.
Oven-dried subsamples of the fine-earth fractions (i.e., bulk soil horizons) from all four sites were density separated according to Heckman et al. (2014). Briefly, 10 g (for O/A horizons) or 30 g (for B horizons) were suspended in sodium polytungstate adjusted to a density of 1.85 g cm −3 to float the free light fraction (FLF), which was aspirated and rinsed in triplicate to remove residual polytungstate. Residual (non-floating) material was sonicated at a rate of 1500 J/g -1 of soil to release light fraction material (<1.85 g cm −3 ) held within aggregates (OLF, the occluded light fraction). After aspirating the OLF and associated supernatant sodium polytungstate off the top, the supernatant polytungstate was poured off, and the remaining heavy fraction (HF) remaining at the bottom was collected and rinsed repeatedly with water and dilute MgCl 2 . Across the four sites, O/A horizons consistently yielded material in all three density fractions; however, the FLF and OLF were often not present in B1 horizons, and were rarely recovered from the deepest B horizons of each site. When trace amounts of material were recoverable from any given density fraction from a specific horizon, analytical preference was for C concentration and radiocarbon ( 14 C) content. To ensure that the density fractionation procedure was not resulting in the loss or gain of mass, C, or Hg, we performed several calculations to compare the properties of bulk soil horizons to the density fractions isolated from them. First, after oven-drying and weighing the density fractions isolated from each sample, we summed the individual fraction masses and computed mass recovery as a proportion of the starting sample. Proportional mass recoveries on an individual sample basis ranged from 0.98 to 1.00. Second, to ensure that neither C nor Hg were lost nor gained during density fractionation, we calculated the weighted average concentration of C and Hg across the density fractions recovered from each bulk soil sample, and compared the weighted averages with the measured concentration of C or Hg (C m ) on that bulk soil sample (see next section for analytical methods) as: Where C m is the measured elemental concentration for the bulk soil sample and C p is the predicted elemental concentration for the bulk soil sample, based on the elemental concentrations (C 1 , C n ) and mass proportions (P 1 , P n ), respectively, of the first through n th density fractions (up to 3 per sample). Comparisons revealed no systematic deviation between predicted and measured elemental concentrations, and regression analyses showed very strong correlations between the C and Hg concentrations for bulk soils and their mathematically recombined density fractions. Specifically, r 2 values between C m and C p were very high for C and Hg (0.998 and 0.994, respectively), with neither model having a slope significantly different from 1:1.

Chemical Analyses
Composited soil samples (i.e., bulk soil horizons) and their density fractions in this multi-institutional collaborative study were analyzed for total C and Hg concentrations, and also graphitized for analysis for 14 C content. Details for each are as follows.
Total C concentrations were analyzed on subsamples of each bulk soil or density fraction at up to three facilities; r 2 values indicated that the methods agreed very closely in their estimates of total C concentration (linear regression r 2 = 0.997, 0.979, and 0.995). Where multiple measurements of C concentration were made on a given sample, we took an average and used the mean for data analyses; where only one was available, that value was used. Two of the three facilities that analyzed total C concentrations utilized elemental analyzers (Virginia Tech and UMBS Analytical Table 2. Key descriptive, physical, and chemical properties of the four soils detailed in the cross-site analysis. Data derive from the number of composited samples of each horizon (indicated in the "n" column). The proportional distribution (as percent) of bulk soil material among the three density fractions (FLF, free light fraction; OLF, occluded light fraction; HF, heavy fraction) is presented following the characteristics of bulk soil C (concentration as per cent and radiocarbon Fraction Modern as a proportion). Laboratories); these were on an Elementar Vario MICRO cube elemental analyzer (Elementar Analysensysteme GmbH) and a Costech Analytical CHN Analyzer Thermo Scientific EA Isolink (Costech Analytical Technologies; Thermo Scientific), respectively. The third facility, the Carbon, Water, and Soils Laboratory at the USDA-Forest Service, Northern Research Station (Houghton, MI), measured the C concentrations of known sample masses directly, using an inline manometer to monitor the partial pressure of CO 2 produced during sample preparation for 14 C analysis.
Preparation for 14 C analysis began by placing each precisely weighed sample (bulk soil or density fraction) into a quartz tube, sealing the tube under vacuum, and combusting the sample at 900°C to convert all C into CO 2 . The CO 2 was then reduced to graphite through heating at 570°C in the presence of hydrogen gas and an iron catalyst (Vogel et al., 1987). Subsequently, graphitized samples were measured for 14 C content by accelerator mass spectrometry at Lawrence Livermore National Laboratory (Davis et al., 1990). Data were normalized to the Oxalic Acid I 14 C standard, followed by background subtraction determined from 14 C-free coal or wood and a δ 13 C correction to account for isotopic fractionation. For this analysis, we carefully considered several alternative units for reporting 14 C values of bulk soil horizons and their density fractions, ultimately selecting to report 14 C values as the fraction of the C in each sample that derives from modern sources (Fraction Modern [Fm]). In other words, Fm values quantify the proportion of C in a soil horizon or density fraction that entered that sample since atmospheric nuclear weapons testing in the 1950s (sometimes termed "bomb carbon"). Readers are referred to Ceballos-Nunez et al. (2018) for discussion of current views on the language of radiocarbon in open systems, and problems associated with metrics such as standard radiocarbon age and mean residence time. Radiocarbon data, including identifier information and analytical errors (standard deviations of n = 5 to7 analytical replicates per sample, all of which were <0.005 Fm) are provided in Supplementary Table T1.
Total Hg concentrations of bulk soils and density fractions were measured at the UMBS Analytical Laboratory on a Milestone DMA-80 automated mercury analyzer (Milestone Inc.) according to USEPA method 7473 (U.S. EPA, 1998). For each sample possessing an Hg concentration, we computed Hg/C as the mass ratio of total Hg to total C (in ng g -1 ). These metrics of Hg were available for all 12 of the bulk soil horizons from the four sites (which were themselves composites of n = 2-7 samples per site); sample masses were insufficient for 5 of the 30 density fractions across the four sites.

Soil Survey data
To compare results from our four intensive study sites to Spodosols more broadly, and compare C and Hg concentrations and relationships in Spodosols to soils of other taxonomic orders, we acquired soil survey information and laboratory characterization data from the NCSS. Specifically, we accessed the August 2014 version of the Soil Characterization Database via the International Soil Carbon Network, and queried it to compile a dataset of pedons (and their constituent horizons) possessing C, Hg, and key physical, chemical, and soil taxonomic information. This dataset consisted of 14,007 horizons from 2181 pedons spanning all soil taxonomic orders; we parsed it into Spodosol (595 horizons from 91 pedons in its final, cleaned condition for analysis) vs. "all-soils" (5409 horizons from 1051 pedons) datasets for separate analyses. Detailed information on data handling and preparation for the NCSS datasets are in the supplementary information.

Data Analysis
We tested our hypotheses using a combination of parametric and nonparametric statistical approaches, preferentially utilizing the former for the more detailed cross-site analyses and utilizing the latter for some of the tests comprising our broader assessment of NCSS data. In the cross-site analysis, most response variables had to be ln-transformed to normalize their distributions for parametric tests. These tests included: (i) simple linear regressions to test relationships between C and Hg concentrations, Hg/C ratios, and Fm values of bulk soil horizons, and density fractions within and across horizons; (ii) two-way ANOVAs to test for significant effects of site and horizon [or density fraction within O/A horizons] on C and Hg concentrations, Hg/C ratios, and Fm values; (iii) one-way ANOVAs to separately test for significant effects of site and density fraction within B horizons. These lattermost tests (the one-way ANOVAs) were necessary because the absence of some density fractions from B horizons of certain sites (i.e., an incomplete factorial design) would not allow for two-way ANOVAs. When ANOVA test results were statistically significant, we used the Fisher's Least Significant Difference multiple comparisons method to identify significantly different groups.
To test for significant differences in C and Hg concentrations and Hg/C ratios across all genetic horizons (and the three key horizons of O/A, B1, and deepest B) in NCSS Spodosol pedons, we used the nonparametric Kruskal-Wallis test of medians (with Dunn's multiple comparisons method) because the response parameters could not be transformed to normality. We tested whether between-site variation was significant in the NCSS Spodosols dataset by running Kruskal-Wallis tests on C and Hg concentrations and Hg/C ratios and using state of origin, taxonomic suborder, great group, and subgroup as categorical predictors (group variables) and using Dunn's test to identify whether groups differed significantly from one another. We also used Kruskal-Wallis and Dunn's tests to analyze differences in the C, Hg, and pH values of soil samples across all taxonomic orders in the "all-soils" dataset. To explore continuous relationships between these variables, we ln-transformed their values to improve normality and ran simple linear regression to test significance of relationships between C, Hg, Hg/C, and pH in Spodosols. To test whether Spodosol C-Hg relationships differed from other soil orders in the "all-soils" dataset, we used best subsets regression to test for the statistical significance of soil %C (ln-transformed), pH (ln-transformed), and taxonomic order (represented by categorical dummy variable coding, with Spodosols as the default) in predicting soil total Hg concentrations (also ln-transformed). We selected the strongest model on the basis of its adjusted R 2 (coefficient of multiple determination), C-p, and partial tand Pstatistics of the individual predictors. In all tests, we set P < 0.05 as the a priori threshold for accepting results as statistically significant.

Patterns and Sources of Variation across the Four Intensive Sites
The four Spodosols investigated in detail contrasted in climate, vegetation, parent material, topography, and time available for soil formation, and in keeping with that variation, they differed in morphology, taxonomy, and physical and chemical characteristics (Tables 1 and 2). Within each horizon, total C and Hg concentrations spanned one to two orders of magnitude across sites, and the percentage distribution of bulk soil material among the three density fractions also varied widely. Only surface horizons (O/A horizons) consistently held material in each density fraction. The proportion of bulk soil material recovered in light fraction pools (FLF and OLF) decreased (as HF increased) with depth at all four sites. Parallel to this tendency, the C held in each horizon became less modern (i.e., Fm values decreased) with depth at all sites.
Across the four sites, C and Hg concentrations were strongly positively related. This relationship was true across bulk soil horizons (Fig. 1A), as well as within the density fractions that comprised O/A horizons (Fig. 1B) and B horizons (Fig. 1C).
Fraction modern values (i.e., the proportion of C within each sample that has been input since the 1950s) ranged from 0.3023 to 1.1472 across the density fractions from all sites and horizons, and exhibited several significant relationships with C and Hg concentrations and ratios (Fig. 2). First, heavy fractions (HF) with larger Fm

Fig. 2. Relationships between fraction modern (Fm) values and concentrations of total C (Panel A), total Hg (Panel B), and Hg/C ratios (Panel C) for density fractions of O/A and B horizons from the four sites. Each point represents a value for a density fraction that was isolated by processing two to seven composite samples of a genetic horizon from one of the four sites; density fractions that did not yield sufficient material for Hg analysis are not plotted in Panels B and C.
values had higher C concentrations ( Fig. 2A; r 2 = 0.55, P = 0.006); in contrast, the consistently high C concentrations of light fractions (FLF and OLF) did not vary significantly across their wide-ranging Fm values. Relationships between Fm and total Hg concentrations showed a different pattern; namely, there was no statistically significant relationship for HF, while light fractions (FLF and OLF) with larger Fm values had lower Hg concentrations ( Fig. 2B; r 2 = 0.64, P < 0.001). Lastly, across density fractions from all sites and horizons, Hg/C ratios narrowed as Fm values increased ( Fig. 2C; r 2 = 0.66, P < 0.0001).
Three of the four properties of interest to this analysis (C and Hg concentrations, Hg/C ratios, and Fm values) differed significantly between horizons (Fig. 3), but none of these properties differed significantly between sites. Specifically, C concentrations were significantly greater for O/A than deepest B horizons ( Investigating variation within the density fractions that comprise O/A horizons revealed that fractions differed significantly in C ( Fig. 4A; P < 0.001) and Hg ( Fig.  4B; P < 0.014) concentrations, while sites did not differ significantly in any properties. Concentrations of both elements were significantly higher in the light fraction materials (FLF and OLF) recovered from O/A horizons than in the HF material. Fractions did not differ significantly in Fm (P = 0.09) nor Hg/C (P = 0.176); site P values ranged from 0.107 for Hg concentration to 0.987 for Fm.
Similar to bulk soils overall and the fractions of O/A horizons, properties of the density fractions comprising B horizons did not differ significantly between sites, but differences between fractions were statistically significant for three of the four properties of interest (Fig. 5). Specifically, mean concentrations of total C (Fig. 5A) and total Hg (Fig. 5B) differed significantly between fractions (P < 0.001 for each element), with concentrations of both being significantly higher in the FLF and OLF than in the HF. Carbon in the FLF had significantly higher Fm proportions than in the OLF and HF ( Fig. 5C; P = 0.014). Mean Hg/C ratios of density fractions did not differ significantly between fractions.

Broader Patterns in the NCSS Spodosols Dataset
Soil horizons from the Spodosols in the NCSS dataset (n = 595) exhibited significant differences in C and total Hg concentrations and Hg/C ratios that were generally similar to those observed for the horizons from the four intensive sites. Significant differences in these properties were apparent between many master horizons (Fig. 6), but especially pronounced for the three key horizons of interest in the cross-site analysis (Table 3; O/A, B1, and deepest B in each of the 91 NCSS Spodosol pedons). None of the three properties (C, Hg, Hg/C) of these three key genetic horizons differed significantly between states, taxonomic suborders, great groups, or subgroups.
Bulk soils (horizons) from the broader population of Spodosols showed a similar, statistically significant positive relationship between C and Hg concentrations (Fig. 7A) as observed for bulk soils from the four sites in the detailed cross-site analysis. Furthermore, pH (which ranged from 3.4-8.4) was a statistically significant predictor of C (Fig. 7B) and total Hg (Fig. 7C) concentrations and Hg/C ratios (Fig. 7D) across Spodosols. These pH-dependent relationships were not detectable for bulk soils for the four detailed study sites, which spanned a narrow, intermediate range of pH values (4.2-5.9).
For soil samples from all taxonomic orders in the NCSS "all-soils" dataset (n = 5409), variation in C concentrations was the single strongest predictor of variation in Hg concentrations (Table 4, t = 27.408). However, the strongest model (r 2 = 0.22; P < 0.001) also included soil pH, and identified significantly different relationships between C and Hg concentrations for most soil orders. That model, which set Spodosols as the default (0) and coded each other order to its own categorical dummy variable (1), identified Aridisols as showing a significant, negative relationship between C concentration and total Hg concentration, Ultisols, Inceptisols, and Entisols as having significantly lower Hg concentrations per unit C, and Alfisols and Andisols as having significantly higher total Hg concentrations per unit C than Spodosols. Mollisols did not differ from other orders.
Examining the median C and total Hg concentrations and pH values of bulk soil horizons from all taxonomic orders ( Fig. 8; P < 0.001 for each) in the NCSS "all-soils" dataset revealed significant differences between orders and provided context for multiple regression results. Specifically, Andisols, Spodosols, Aridisols, and Mollisols had the highest median C concentrations (in descending sequence); these medians were significantly higher than most of the other orders. Aridisols and Mollisols had the highest me-dian pH values (8.2 and 7.5, respectively), while median pH values for Spodosols,Ultisols,and Inceptisols (5.0,5.1,5.4,respectively) were significantly lower than pH values for other orders. In terms of their median total Hg concentrations, Andisols had by far the highest (120 μg kg -1 ), differing from all other orders, while Inceptisols (47 μg kg -1 ), Spodosols and Ultisols (44 μg kg -1 each) followed and were higher than those with the lowest (Aridisols at 18 μg kg -1 and Mollisols at 30 μg kg -1 , respectively).

DISCUSSION
Density fractionation and radiocarbon analysis provide the means to place C-Hg interactions in Spodosols in a contemporary paradigm of SOM stabilization and turnover. Recent analyses of soils in regions where Spodosols are widespread (e.g., Richardson et al., 2013;Yu et al., 2014) have thoroughly described drivers of spatial variation in soil Hg, and through soil property measurements have revealed empirical and mechanistic relationships between factors that control C-Hg interactions within soils. However, even as these and other studies have highlighted the importance of SOM, and suggested that soil-forming processes influence C-Hg interactions, none have utilized the two methods that have arguably done the most in recent years to advance understanding of SOM stabilization and its critical role in soil formation. Furthermore, there have been no explicit tests of soil taxonomy as a set of variables for predicting C-Hg relationships across soils. Superficially, some of our findings using these methods confirm that empirical patterns in bulk soils extend to the discrete SOM pools that comprise them; for example, C and Hg concentrations are positively related in bulk soils as well as within their density fractions. Closer scrutiny of these empirical relationships leads to more nuanced (and mechanistic) interpretations that are rooted in soil-forming processes and the taxonomy that categorizes them.
The most important results of our cross-site intensive analysis using density separation and 14 C analysis are those demonstrating that bulk soils are mixtures of SOM pools that cycle C and Hg very differently. This departs from previous work on C-Hg interactions in Spodosols, which has largely viewed individual genetic horizons (or even depth increments) as the fundamental units of soil. One example is the positive relationship between Fm and C concentration in the organo-mineral HF but not in the LF, which consists of substrates that are C-rich but nonetheless span a very wide range of Fm values. This wide range of Fm values in light fractions (FLF and OLF) indicates that presumably "labile" (C-rich, metabolically favorable) detrital inputs can persist in the matrix through physical incorporation into aggregates where their C is protected from microbial attack (Wagai et al., 2009). It is not stability through inherent recalcitrance, nor stabilization through "humification" during gradual decomposition, but physical protection that causes these C-rich substrates to persist (Swanston et al., 2005;Heckman et al., 2014, Schrumpf et al., 2013. In turn, if concepts such as recalcitrance or humification appear to explain variation in Hg concentration across soil horizons (e.g., Yu et al., 2014), the apparent relationship is likely coincidental. Indeed, one of two mechanisms proposed (but not directly measured) by Yu et al. (2014) to explain wider Hg/C ratios in mineral than organic horizons requires no reference to either of those concepts. Namely, Hg may be directly bound to mineral surfaces. Our findings of wider Hg/C ratios in deeper soil horizons and in the organo-mineral HF support this mechanism, and these are further reinforced by results demonstrating that extractable Fe concentrations explain Hg concentrations in lower (B) but not upper (A) mineral soil horizons (Richardson et al., 2013). Density separated heavy fractions in acidic soils with poorly crystalline mineralogy are characteristically high in Fe and Al hydroxides that complex with C (Kaiser et al., 2002) and also directly with inorganic Hg, especially at the low pH values characteristic of Spodosols (Dmytriw et al., 1995;Schlüter, 1997). Thus, the long-reported empirical relationship between C and Hg in bulk soils may be misleading, because it suggests that these two elements are closely coupled by a generalizable mechanism, even in cases where it does not operate. Rather, the mineral surfaces (e.g., Fe hydroxides) that are so important for SOM stabilization may bind to C and Hg independently of one another. As these sorption processes are not limited to Spodosols, further investigation of the distribution of C and Hg among functional groups and mineral fractions may help resolve the degree to which Hg retention is tied to the fate of C across a wider range of soils.
Another example of bulk soil properties belying differences between the SOM pools that comprise them is evidenced by relationships between Fm values and Hg concentrations. In light fraction pools, which span a narrower and more modern range of Fm values, Hg concentrations decline as Fm values increase. This relationship is absent, or potentially even opposite for HF SOM. To the extent that Fm indicates the age of these C-rich detrital inputs to soil (the FLF and OLF; Kaiser et al., 2002;McFarlane et al., 2013;Schulze et al., 2009), their Fm values reflect recent (~decadal) incorporation, and the negative relationship between their apparent age and Hg concentration may reveal an influ-ence of the Clean Air Act. Specifically, rates of atmospheric Hg deposition have been declining in recent decades (Schuster et al., 2002), so Hg inputs to these C-rich OM pools via litterfall and atmospheric deposition have likely declined over the same timescale. In other words, newer have lower Hg concentrations. Obrist et al. (2011) cited this legacy of high but declining Hg deposition  rates as the likely explanation for wider Hg/C ratios in litter and soils containing organic matter that they assumed to be older, but they also noted an alternative (or supplemental) explanation: Hg may be retained as C is lost during decomposition. The latter may be an important mechanism if significant fractions of C and Hg are not bound directly, but rather cycling independently of one another. Our analyses help clarify these explanations, and suggest that they may be more important in some SOM pools than others. Radiocarbon data do indicate that the mixture of C substrates present in bulk soils becomes less modern with depth, and deeper horizons have wider Hg/C ratios, superficially supporting the legacy effects explanation (Obrist et al., 2011). However, topsoil horizons that hold a substantial share of their C in pools of modern, input-driven, light fraction organic matter are fundamentally different than the deeper illuvial horizons where C is scarce, older overall, and mostly intimately associated with minerals. In these horizons, we find that Hg/C ratios do not differ between density fractions, even though they do differ in their C and Hg concentrations and Fm proportions. These illuvial horizon dense fractions are precisely the SOM pools where independent interactions between metal hydroxides, C, and Hg are most likely to be occurring. This suggests that C and Hg are cycling independently of one another in deeper soils, arguing against overgeneralization of (or overemphasis on) direct interactions between C and Hg in soils. Because the design of our cross-site intensive analysis was opportunistic, spanning four sites that simultaneously varied in multiple soil-forming factors, it is important to question whether the results we report for these sites are representative of Spodosols more broadly. Fortunately, opportunistic though it is, our design is hierarchical: it allowed us to study C-Hg interactions within horizons (i.e., between the SOM pools that comprise them), within profiles (i.e., between their horizons), and between sites (i.e., the 4 intensive sites and the 91 Spodosol pedons from the NCSS database). The overarching inference we gain from this design is that similar processes influence C-Hg interactions within horizons as within profiles, and these are consistent across Spodosols. The empirical evidence for this argument is the consistent relationship between C and Hg concentrations across density fractions and bulk soils in the four-site comparison and the wider NCSS Spodosols dataset, strengthened by the notable lack of significant between-site variation in both datasets. In terms of a mechanistic interpretation, we suggest the processes that transform and transfer C and Hg vertically within Spodosol profiles (between horizons) are analogous to the processes that transform C and Hg between the discrete forms that are isolated by density fractionation. In other words, each bulk soil horizon in a Spodosol profile is functionally similar to one of the SOM pools contained within it. Carbon and Hg are input to, and cycle closely within, the organic-rich surface  horizons, where much of the SOM is held in modern, light fraction pools. These elements may persist or even accumulate in physically protected aggregate SOM pools in surface or upper illuvial horizons, and their stability may be shared through direct bonding. In the organo-mineral heavy fraction, which comprises an increasing proportion of the soil matrix (and SOM) in successively deeper horizons, direct relationships between C, metal hydroxides, and Hg may be expected to be the dominant controls on the transformations (and ultimately the outputs) of C and Hg. In this regard, the factors that control SOM stabilizing processes across depths also control the inputs, outputs, transfers and transformations of C and Hg across depths. Namely, litter, climate, mixing organisms and other top-down factors are most important near the surface of the profile; mineralogy, redox potential, groundwater interactions become more important near the bottom. Models of SOM stabilization and soil formation are conceptually similar, and their utility for explaining C-Hg distributions and interactions in Spodosols extends across other soil orders. Distinct soil orders possess different soil-forming processes, but all soilforming processes fall into the four categories defined by Simonson (1959); in turn, the physical and chemical factors that constrain these processes also control C-Hg interactions. Andisols provide a good example of the transferability of a pedogenic perspective on C-Hg interactions, because some of their properties, soil-forming, and SOM-stabilizing processes are similar to Spodosols. Soils in both orders form in acidic, relatively porous parent materials with strong vertical infiltration, and are abundant in poorly crystalline metal hydroxides that form stable bonds with C and Hg (Parfitt, 2009;Spielvogel et al., 2008). The two orders differ in that Andisols form well-developed hierarchical aggregates (Asano and Wagai, 2014); possessing two highly effective processes of SOM stabilization confers unique physical properties such as low bulk density and strong P sorption that differentiate Andisols not only from Spodosols but from all orders as well (Dahlgren et al., 2004). Lastly, while their moderately lower acidity may limit Hg solubility compared with the strongly acidic Spodols (Schuster, 1991), Andisols' volcanic parent materials may contribute some geogenic Hg in addition to the atmospheric inputs that are dominated by anthropogenic sources (Fitzgerald and Lamborg, 2014). Collectively, these properties of and processes make Andisols unique among soil orders, and as a result they showed the highest concentrations of C, Hg, and Hg per unit C in our analysis. Aridisols, with their high levels of inorganic C, high pH, and typically clay-dominated mineral assemblages occupy an opposite extreme from the acidic Spodosols and Andisols, with their poorly crystalline metal hydroxide mineralogy and high capacity for SOM stabilization. Aridisols are rarely studied in the SOM stabilization literature, likely because they are low in organic C (Heuscher et al., 2005), have limited capacity to incorporate C into physically stable pools (Kelly et al., 2017), and any stability that they do possess likely derives from strong climatic control. In our analysis, Aridisols are unique for being the only order to demonstrate a negative relationship between C and Hg concentrations. This likely reflects that increasing total C concentrations in this soil order are a function of higher inorganic C, and associated with higher pH values that limit Hg mobility and enhance its volatilization (Yin et al., 1996;Schlüter, 2000).
In between the extremes of Andisols and Aridisols lay other soil orders with narrower, yet still significant differences in their pH, C and Hg concentrations. These differences are beyond the scope of this paper, but their alignment along soil taxonomic lines highlights the potential of the pedologic framework for assessing C-Hg interactions across soils. These patterns also suggest potentially fruitful research avenues, for example, exploring C-Hg interactions across developmental sequences of soil orders using density separation to investigate discrete SOM pools. For example, differences in C concentration and pH among Entisols, Inceptisols, Spodosols, and Ultisols (almost all of which were significant in our analysis) do not translate into significant differences in Hg con- centrations in soils across these taxonomic orders. Given that these orders differ in their degree of development and possess very different capacity for and processes of SOM stabilization, differences in their C-Hg interactions are likely. Site-specific factors may emerge as important in future analyses aimed at such questions, but our analysis suggests differences between soil horizons and discrete SOM pools are likely to be detectable nonetheless.

CONCLUSIONS
Based on empirical relationships between C and Hg reported in the literature for bulk soils, we hypothesized that C and Hg would be positively related in soil density fractions. Our analysis confirmed this hypothesis. However, the discrete SOM pools represented by density fractions showed divergent relationships between radiocarbon values, C and Hg concentrations and ratios, indicating that generalizations based on bulk soils do not apply equally to all SOM pools present in Spodosols. We argue for recognition of mineral adsorption as the critical mechanism for retaining C and Hg in soils, potentially independently of one another, especially in the most stable pools. More broadly, our results indicate that soil taxonomy is a useful framework for categorizing differences in C and Hg concentrations and their interactions across soils.