Weathering Intensity and Presence of Vegetation Are Key Controls on Soil Phosphorus Concentrations: Implications for Past and Future Terrestrial Ecosystems

Abstract

Phosphorus (P) is an essential limiting nutrient in marine and terrestrial ecosystems. Understanding the natural and anthropogenic influence on P concentration in soils is critical for predicting how its distribution in soils may shift as climate changes. While it is known that P is sourced from bedrock weathering, relationships between weathering, P, and other soil-forming factors have not been quantified at continental scales, limiting our ability to predict large-scale changes in P concentrations. Additionally, while we know that Fe oxide-associated P is an important P phase in terrestrial environments, the range in and controls on soil Fe concentrations and species (e.g., Fe in oxides, labile Fe) are poorly constrained. Here, we explore the relationships between soil P and Fe concentrations, soil order, climate, and vegetation in over 5000 soils, and Fe speciation in ca. 400 soils. Weathering intensity has a nuanced control on P concentrations in soils, with P concentrations peaking at intermediate weathering intensities (Chemical Index of Alteration, CIA~60). The presence of vegetation (but not plant functional types) affected soils’ ability to accumulate P. Contrary to expectations, P was not more strongly associated with Fe in oxides than other Fe phases. These results are useful both for predicting changes in potential P fluxes from soils to rivers under climate change and for reconstructing changes in terrestrial nutrient limitations in Earth’s past. In particular, soils’ tendency to accumulate more P with the presence of vegetation suggests that biogeochemical models invoking the evolution and spread of land plants as a driver for increased P fluxes in the geological record may need to be revisited.

1. Introduction

Phosphorus (P) is an essential, often limiting nutrient in ecosystems across the globe [1,2]. P has been studied in plants, soils, rivers, lakes, and oceans—both modern and ancient—and is at the center of critical questions about Earth’s past and future. Today, concerns around P in terrestrial settings focus on soil fertility and crop yields, as well as on agricultural runoff and its effect on eutrophication [3,4,5]. Climate, vegetation, and weathering and erosion rates are inextricably linked, so in order to improve quantitative constraints on potential P fluxes from soils, each of these factors must be considered. The redox-sensitive metal iron (Fe) is often associated with both P and organic matter/carbon (C) in soils (e.g., [6,7,8,9]). The P-Fe oxide pathway is critical for terrestrial P transport; however, it is currently less well-understood than other aspects of terrestrial P cycling, such as P-C associations. Additionally, while much research linking P, weathering, and soil age has been done, such work has typically done on relatively small scales (e.g., chronosequences in Hawaii) relative to the global scale of the questions that we ask here. Therefore, constraining climate, vegetation, and weathering controls on soil P and Fe is critical for predicting changes to fluxes of P from soils to rivers and beyond, and to improving our understanding of how regional soil fertility may evolve in response to a changing climate.

Prior to the evolution and expansion of land plants, changes in the flux of P from continents to oceans are thought to have been a primary control on marine productivity, and consequently, on the concentration of oxygen in the atmosphere [10,11]. It remains debated whether land plants’ colonization of continents would have increased [12] or decreased [13] that P flux. The better we understand the modern distribution of and controls on P in soils, the better we can model global biogeochemical changes in Earth’s past.

1.1. P, Fe, and Weathering and Erosion

P is released from bedrock via the weathering of P-bearing minerals (primarily apatites, which are Ca minerals relatively resistant to weathering [4]). Weathering depends on all soil-forming factors (climate, biology, topography, time, parent material), but is primarily controlled by climate (e.g., precipitation, temperature, seasonality) and time (surface age, length of exposure [14]). Organic matter also serves as an important pool of P in soils ([4,5]), and vegetation, fungi, and microbial life play an important role in both P and Fe liberation and mobilization through symbiotic (mycorrhizal) relationships [15,16]. The presence of plants, then, changes how both P and Fe are distributed in soils—but the links between landscape-scale soil chemistry, plant functional type, climate, and weathering intensity have not been quantified.

P can be depleted from a soil profile within thousands to tens of thousands of years [4,17,18,19,20]. In terrestrial vegetation (biomass), P’s residence time is ~13 years, and residence in soil is ~600 years without P replenishment by dust [4]. Soil P can also be replenished through dust deposition (e.g., [17,21,22]). The natural progression of a soil, with a slow accumulation of Fe and Al oxides and a loss of mobile cations and nutrients, means that P is typically very low in old and/or intensely weathered soils [19]. Other soil properties, such as clay content, are also linked to weathering and vegetation, and could affect the concentrations of both P and Fe. Clay minerals and clay grain-size fractions can sorb considerable amounts of P and Fe, and clay content typically increases with weathering time and/or intensity [23,24]. Weathering intensity is expected to decrease with latitude due to colder temperatures, drier precipitation regimes, and less vegetation and shallower rooting depths [23,24]. Additionally, because the potential volume of erodible material generally increases as weathering intensity increases (though varies with factors such as slope and uplift rate, e.g., [25,26,27,28,29], weathering is linked to P fluxes not only through direct apatite dissolution, but also through erosion rates, e.g., [30,31]). Understanding how weathering intensity is related to soil P (and Fe) concentrations is important for predicting how their erosional fluxes and transport may change in response to climate change (e.g., [30,31]).

After weathering, dust deposition is an important secondary source of P, particularly for old soils whose P has been depleted through weathering, erosion, and biological use. For example, chronosequence studies in Hawaii have shown that for older soils, dust replenishment of P is critical to maintaining fertility once bedrock sources have been depleted [17,32,33]. Dust is also a critical source of P in the Amazon basin, which receives significant dust fluxes from Africa [18,34,35,36], compensating for P loss in deeply weathered tropical soils. Dust-sourced P is less significant for areas that receive little dust deposition, for younger soils that still have a regular input of P from bedrock weathering, or for shorter timescales than either for significant dust accumulation or timescales of soils losing their bedrock-sourced P (i.e., anthropogenic timescales). Arvin et al. [21] found that dust deposition is an important P source in montane soils, despite their relatively “fresh” bedrock. The balance between weathering rate and dust accumulation is also relevant for determining how important the latter is in soil P [22]. Dust is also commonly a source of Fe oxides, increasing soil fertility [35,37,38].

1.2. P and Fe Kinetics in Soils

P and Fe are chemically associated in soils in part because of P’s affinity to sorb to Fe (oxyhydr)oxides (e.g., [4,9,39,40,41,42]). Orthophosphate, the common ion form of free P in soils, is reactive towards those oxides and becomes immobile (not readily bioavailable) when they associate into neoformed minerals or are sorbed onto existing mineral surfaces. Fe-bound P may be a critical mode of P transport among terrestrial ecosystems and from terrestrial/aquatic to marine environments, as oxyhydroxide-bound P is one of the major bioavailable/exchangeable forms of P in fluvial systems [5].

While we know that Fe oxide-adsorption of P is kinetically favorable in soils and other environments where P is being exchanged (e.g., [43,44,45]), our understanding of first-order controls on Fe concentration and speciation (i.e., Fe³⁺ in oxides versus labile Fe²⁺) in soils could be improved. Examining a full suite of extractable Fe (rather than only dithionite-citrate and/or ammonium oxalate) allows a more nuanced picture of Fe phases in soils, which can ultimately be linked to variable P mobility in different soil redox conditions [43,46,47]. It is possible that regional climate factors (precipitation, temperature) exert control on Fe speciation (i.e., reduced vs. oxidized) in soils via soil moisture. Controls on soil moisture are debated and likely vary from profile to profile [48,49,50,51,52,53], and pH (related to soil pore space saturation; often microbially-mediated) is a primary control on reduction/oxidation of Fe phases in soils (e.g., [54]). Additionally, soil porewaters can serve as an intermediate step between recalcitrance and fluvial suspension for dissolved constituents (e.g., [54,55]). Exploring these relationships is key for understanding the likelihood of P mobilization under different moisture situations.

1.3. P and Fe in the Geologic Record

Geologists are interested in constraining many of the processes described above in ancient soils and ecosystems. Weathered P, transported to the oceans via fluvial networks and continental drainage, is typically invoked as a first-order control on marine primary productivity, which is associated with C sedimentation, atmospheric CO₂ drawdown, and oxygen production (e.g., [10,11,56]). Similar to P, Fe can stimulate marine productivity on short timescales (e.g., [57,58,59,60]), but it is of additional interest because in the past, it served as a P sorption site in marine systems, effectively sequestering P so that it is not bioavailable (the so-called “Fe-P trap”) and limiting marine productivity [10].

With these processes in mind, geologists are interested in continent-to-ocean fluxes and mobility of Fe and P through time. However, the continent-to-ocean transport of P, although central to many biogeochemical models, is usually prescribed and not based on actual changes in fossil soil (paleosol) P values or weathering intensity data. Rather, it has been generalized based on crustal P values and fluxes controlled by large-magnitude changes, such as global glaciations or limited modern observations (e.g., [61]). The value in providing robust constraints on the range and distribution of P in soils, as well as climatic and environmental controls on its mobility and bioavailability, will vastly improve the parameters for these models.

In addition to biogeochemical models, geologists are interested in how the evolution of terrestrial vegetation and mycorrhizae (fungal root symbionts) affected P bioavailability, terrestrial mobility, and fluxes to the oceans. Researchers have argued for both an increase [12,56,62,63] and decrease [13,15] in P fluxes to the oceans following the spread of land plants. Constraining how modern vegetation relates to soil P mobility and distributions will help to answer this question.

Finally, the geochemical composition of fossil soils is used to reconstruct a range of climatic and environmental changes; however, these tools are primarily based on small datasets. This work provides ranges of reasonable values for soil chemistries under known environmental and atmospheric conditions, providing critical background information to improve our paleoclimate and paleoenvironment reconstructions.

1.4. This Work

In the context of these concerns, several big-picture questions remain. Despite much research on P in modern soils, ecosystem (vegetation, anthropogenic impact) and climate (precipitation, weathering) controls on its concentration have not been quantified on landscape to continental scales for generalizable conclusions. Additionally, although P is known to have a strong affinity towards Fe (oxyhydr)oxides in soils, Fe species and soil redox conditions have not been considered as potential first-order controls on landscape-scale P bioavailability. Questions about how climate (precipitation, temperature), soil factors (drainage, soil moisture, clay content), and overall weathering intensity affect soil Fe-P associations needs to be interrogated, to see if Fe species ultimately serves as a mediator for P bioavailability. Moreover, once the distribution of and controls on soil P have been quantified, how do those new constraints affect or constrain our understanding of biogeochemical P models in Earth’s past?

We address these questions, along with hypotheses about P in the geologic record, using a large geochemical database of modern soils (n > 5000). The results can inform models of both past and future terrestrial biogeochemical cycles.

2. Materials and Methods

2.1. Sampling

Soil samples were collected in the field by the authors, and received from the USDA KSSL soil repository, and from contributors (total n = 404). A majority (n = 247) of these samples are B horizons, and therefore, reflect longer-term, stable pools of P (e.g., P sorbed to clays or Fe (oxyhydr)oxides, which we test for here, or recalcitrant organic matter) and Fe. These soils include Entisols, Inceptisols, Alfisols, Ultisols, Histosols, Mollisols, and Aridisols and provide coverage from low (Hawaii) to high (Iceland) northern latitudes (Supplemental Figure S1). Soils collected in the field or from repositories were screened for anthropogenic influence (i.e., on developed or disturbed land), with anthropogenically-altered soils noted as such and binned separately in our treatments of the data here. To complement these soils, we used a United States Geological Survey (USGS) soil geochemistry interactive report [64] and collated geochemical data for the Top 5 cm, A horizons, and C horizons of all the soils included in that work (n = 4857). See Supplemental Table S1 for summary descriptions of these two datasets, and Supplemental Table S2 for complete sample information for the new dataset presented here. Geochemical data for the USGS dataset [64] are available at https://0-doi-org.brum.beds.ac.uk/10.3133/sir20175118.

Combining these soil databases yields excellent spatial, climatological, profile depth, drainage, vegetation (as plant functional type/generalized biome), and soil order coverage (Supplemental Figure S1). While the samples are mostly from the United States, variation in soil forming factors across the continent reflect most of the potential global variation. The tropics are somewhat under-sampled, but data from Hawaii are broadly reflective of those latitudes. While there are species-specific relationships between plants and soil geochemistry, generalized biome is acceptable for the continental scale of questions examined in this work, as well as the comparative nature of our inquiries (e.g., [65,66,67,68]). Differences between plant functional types typically swamp differences within one plant functional type. Additionally, because biomes on sub-continental scales are typical for global biogeochemical models, and because more specific data are not available for the latitudinal scale of this work, generalizing by biome scale/plant functional type is appropriate.

2.2. Geochemical Analyses

All analyses were performed on dry soils within the Sheldon lab, as were samples during the USGS analyses. Soils were dried per USDA APHIS regulations upon import and ground to <70 μm for homogeneity during analyses; samples were stored dry in oxic (ambient room) conditions (during the experimental design, tests of pre- and post-drying, with drying temperatures varying between 25, 50, and 100 °C, showed no significant geochemical differences). Bulk elemental geochemistry was performed at ALS Laboratories in Vancouver, BC, Canada, using a four-acid digestion and Inductively Coupled Plasma Mass Spectrometer (ICP-MS) analysis. External precision was 0.01% for all major elements and 10 ppm for P₂O₅, analyzed at ALS. For samples from the USGS dataset, 1σ sample errors were 488 ppm and 1.39 wt% for P₂O₅ and total Fe, respectively) [64]. The 1σ sample errors for total C and organic C were 4.6% (n.d. for inorganic C [64]).

To analyze Fe speciation in soils, we used a four-step sequential Fe extraction following Poulton and Canfield [69] and Raiswell et al. [70]. This extraction protocol separates operationally-defined pools as ascorbic acid (labile Fe; Fe_asc), dithionite- (crystalline Fe oxyhydroxides; Fe_dith), ammonium oxalate- (highly crystalline Fe oxyhydroxides; Fe_ox), and sodium acetate-extractable Fe (Fe in carbonates, analyzed separately from the sequential extraction; Fe_acet). Initially, we included a chromium-reducible sulfur (CRS [71]) extraction separately to analyze for Fe sulfides (e.g., pyrite), but those concentrations were below detection limit for all non-wetland soils. CRS was limited to wetland soils after the initial exploratory extractions, and raw values were normalized to average soil order densities to account for the large difference in density between Histosols and denser, more mineral-rich soil orders [72]. Sample supernatants were centrifuged and diluted 50×, and Fe was measured via ICP-Optical Emission Spectroscopy and ICP-MS at the University of Michigan.

An internal standard of a local soil, dried at three different temperatures, was included in all runs for inter-run consistency checks and to determine reproducibility. The absolute standard deviation for a sample’s reproducibility depends on the concentration of each Fe species, so the relative percent error ((σ/μ) × 100) is used instead. Repeated analyses show that Fe_asc is reliable within 3.7% of the mean, Fe_dith, within 1.7% of the mean, and Fe_ox, within 7.0% of the mean. Reproducibility on the Fe_acet samples was larger, at 15% of the mean, because the means in many of the samples we analyzed were near machine detection limits (~1000 ppm), where a difference of several hundred ppm yields far larger standard deviations than species with higher concentrations. The Fe_asc pool is considered to be a low constraint because a larger proportion of labile Fe may be in a soil’s porewater (rather than solids). Because here we are primarily interested in relative differences between oxide and labile Fe species, and because species’ concentrations tend to differ by orders of magnitude, we have deemed this approach acceptable. For further reading on nuances in sample preparation and extraction protocols, see two recent reviews by Raiswell et al. [73] and Algeo and Liu [74].

Fe speciation data were filtered for outliers by removing any samples where Fe_spec weight percent > total Fe weight percent (i.e., the concentration of a measured Fe pool was higher than the measured total Fe). These discrepancies (n = 107, 95, 139, 106, for Fe_asc, Fe_acet, Fe_dith, and Fe_ox, respectively) likely arose from the 50×-dilution being insufficient, with high Fe concentrations overwhelming the ICP-OES, leading to erroneous measurements. Future work using Fe speciation on soils should consider using stronger dilution factors for high-Fe soils to circumvent instrument limitations.

2.3. Other Data Collated

Geochemical data for the Top 5 cm, A horizons, and C horizons from 4857 soils around the U.S. are available on an interactive portal [64]. Weathering-relevant cations (Al, Ca, Na, K) and the nutrients of interest (P, Fe) along with latitude and total clay content (clay: A and C horizons only) were compiled. For the USGS dataset, clay content was determined by X-Ray diffraction (for 10 and 14Å clays only). Organic, inorganic, and total carbon data were available for a majority of soils in A and C horizons in the USGS dataset (not all had inorganic carbon). It should be noted that a dry soil sample will have a higher relative weight percent organic carbon due to water loss during drying. Soil orders and drainage capability are not given for soils in this database, but because they span the conterminous U.S., it is a reasonable assumption that they include a variety of soil orders. Clay content as grain-size fraction, parent material, and soil moisture regime were also gathered for the newly-collected dataset from USDA Soil Series reports when available (Supplemental Table S2).

Climate data were gathered from Soil Series pages; when these were not available, the PRISM 30-year averages (1981–2010) were used. If neither Soil Series nor PRISM climate data were available (i.e., for non-US samples), country-specific governmental meteorological data or journal articles were used (see Supplemental Table S2 for details).

The Chemical Index of Alteration (CIA) [75] was calculated to assess the degree of weathering for each soil, which reflects all of the soil-forming factors to some degree, but predominantly soil age and climate [14]. The CIA is thought to reflect the weathering of feldspars to form clay minerals where high values represent more intense weathering. Molar oxide ratios are used (Equation (1)). The CIA of an individual soil may not correspond directly with the CIA values of sediments in rivers (e.g., [76]). However, because our soils’ average CIA values (59 ± 15, 42 ± 32, and 57 ± 21 for A, B, and C horizons, respectively) are within 1σ of the North American rivers’ average from that study (n = 7; mean CIA = 66 ± 8.5, range 53–73), though with larger ranges (0.5 to 99 for all three horizons), we proceed with its use here.

CIA = (Al₂O₃/(Al₂O₃ + CaO + NaO + K₂O)) × 100

(1)

2.4. Statistical Analyses

To examine how the total Fe and P correlate with CIA (weathering intensity), a running average with a 10-CIA-unit window was used. Linear least-squares regressions were used to test for significant linear correlations where appropriate, based on the hypotheses that were being tested. Due to the large size of the USGS dataset, p-values are extremely low (<10⁻¹⁶; essentially 0), meaning it is highly likely that the two variables are not randomly correlated. Strength of correlation should be primarily considered before taking any relationships to be predictive, but for simplicity’s sake, coefficients of determination (R² values) are reported. Principle Components Analysis (PCA) was performed for USGS database soils by horizon (Top 5 cm, A, and C). These analyses are used to test the expected relationships (e.g., that Fe and P will be positively and linearly related), rather than to build predictive models. All of the statistical analyses were performed in Matlab.

3. Results

Complete sample information and all geochemical results are found in Supplemental Tables S2–S5 and are available in the online version of this paper. For the larger USGS dataset [64], we have geochemical, drainage, and vegetation (generalized plant functional type) data, but no soil order or Fe speciation information. Results for soil order, drainage, parent material, and Fe speciation analyses are based only on the newly collected dataset (n = 404).

3.1. Fe and P: Concentrations, Soil Order, and Vegetation

3.1.1. P

The average total P concentration (referred to simply as “P” here onwards) in the continental crust is ca. 870 ppm P, and there is relatively little variation between bedrock lithologies [77]. The maximum soil P concentration in this compilation was >20,000 ppm (>2 wt%; C horizon of an anthropogenically-modified soil), an order of magnitude greater than the crustal average. Mean P concentrations for each horizon studied were depleted relative to the crustal average: Top 5 cm, A, and C horizons means were 660, 632, and 508 ppm P, respectively (USGS data), and B horizons from the new dataset had a mean of 937 ppm P. P concentrations vary spatially throughout continental U.S. (Supplemental Figure S2), with hotspots in a variety of different geologic, climatic, and biologic provinces. P concentrations varied among soil orders (Figure 1), with the highest concentrations in Inceptisols (Figure 1A), a relatively weakly-developed soil type.

The vegetation groups are Barren (lacking vegetation), Altered (i.e., affected by human activity), Forest (no differentiation between deciduous and coniferous), Herbaceous/Grassland, Developed/Cultivated (i.e., occupied or manicured), Shrubland, and Other (mostly including early-succession plants or microbial earths). For USGS samples, protocols dictated that roads, buildings, and industrial sites be avoided [64]. Among vegetation types, there is little variability between P concentrations and vegetation types (Figure 2); while Barren landscapes’ B horizons show higher mean P, this is likely due to small bin size (n = 5; Figure 2C). Overall, variability in concentrations in different vegetation types and soil orders was minimal, with most soils being depleted in P relative to the crustal average. Vegetated soils (both natural and anthropogenically-altered) had higher ranges of P than unvegetated (Barren). There was a weak positive correlation with mean annual precipitation (Supplemental Figure S3A), but the range of P values at a given precipitation amount is typically too large (~1500 ppm) to be of predictive use. This trend could also be explained by related factors, such as vegetation coverage (cover versus bare ground) and weathering intensity, which are associated with precipitation.

3.1.2. Fe

The crustal composition of total Fe (Fe_tot) is more variable than that of P because while P is sourced primarily by apatite-group minerals, Fe_tot can be present in a wide range of minerals and lithologies. The Phanerozoic upper continental crust is estimated to have 3.5 wt% Fe [77]. Fe_tot averages in the soils studied here were 2.1 wt% for Top 5 cm, 1.6 wt% for A horizon, 2.6 wt% for C horizon, and 4 wt% for our dataset (primarily B horizons). Density-normalized Fe_tot varied slightly by soil order (Supplemental Table S3), showing the expected trend of modest loss shifting to modest accumulation during the Alfisol-Ultisol transition (Figure 1B), but otherwise was relatively consistent. Oxisols had the highest values and range in Fe, with some variability among the other soil orders but generally within consistent ranges (Figure 1B).

Fe_tot is generally consistent both within a given horizon and between different vegetation cover types (Figure 3). While most soil orders and most horizons had <3 wt% Fe_tot, values ranged up to >15% (Figure 3; Supplemental Figure S2). The only notable variability in Fe concentrations among vegetation is in Cultivated soil B horizons (Figure 4C), which are depleted, and ‘Other’ soil B horizons that are substantially enriched relative to the other vegetation cover groups. ‘Other’ vegetation contains mostly basalt-parented soils from a limited geographical area (primarily Iceland), so this result is likely an artifact of our sampling. There was no strong correlation between Fe_tot in a soil and mean annual precipitation (R²: 0.2; Supplemental Figure S3B).

3.2. Fe and P: Latitude, Weathering, Soil Order, and Clay Content

Latitude, weathering intensity, and clay content should each be associated with each other as well as with Fe (e.g., [24]). Some of the trends described here could be influenced by a latitudinal sampling bias, with most of the samples coming from the mid-latitudes (Supplemental Figure S1; Figure 4A). In the soils analyzed here, maximum weathering (as measured by CIA) decreases as latitude increases (Figure 4B), while clay content (a weathering product) does not show a strong trend with latitude (Figure 4C), with the exception of a weak mid-latitude (35–40° N) peak in clay content in some C horizons. Weathering trends between soil orders behave as expected, with CIA values increasing from Entisols to Ultisols/Oxisols (Figure 5). Clay content and weathering show a well-behaved and expected pattern of increase that matches the theoretical understanding of that metric [24,75,78,79], with B horizons specifically showing the highest clay content until the highest CIA values are reached (Figure 4D, green points). Clay content and geochemical data were not available for the Top 5 cm subset of the USGS dataset, so that horizon is excluded from those comparisons.

Fe_tot content increases moderately with higher latitudes (Figure 6A), higher weathering (Figure 6C), and clay content (Figure 6E). P in soils behaves less linearly with respect to these variables: P increases moderately with latitude (Figure 6B), but rather than decreasing with weathering, as might be expected based on soil age-P relationships, average P concentrations in all horizons vary little, though there is a general increase in range at moderate weathering intensities (CIA~40–70) (Figure 6D). P concentrations are weakly negatively correlated with clay content (Figure 6F). Fe_tot and P showed positive correlations with one another in all horizons, as expected (Supplemental Figure S4).

3.3. Weathering, Clay Content, and Vegetation

Between vegetation types, there is low variability in clay content, with Forests showing slightly lower values in B horizons than the other groups (Supplemental Figure S5). Forests have the highest weathering in natural (unaltered/not developed) soils, followed by Grasslands and Shrublands (Figure 7).

3.4. Fe and P: Drainage and Soil Moisture

Moderately poorly-drained soils had higher Fe and P concentrations than the other degrees of drainage. Perudic moisture regimes had the highest Fe and P values (Figure 8); however, these samples were dominated by basalt-parented soils in Iceland, which may skew the results. Aside from that, Aquic and Xeric moisture regimes had high Fe, and Ustic and Udic soils had high P (Figure 8). Some trends emerge with relatively high Fe/low P (Aquic) and vice versa (Udic).

3.5. Fe Speciation

3.5.1. Fe Speciation and P

Pyrite/Fe in sulfides, normalized to average soil order density, were above detection by titration in only 22 tested non-wetland soils out of 63 from around the continental U.S. (Supplemental Table S5), so that test was excluded from subsequent analyses. Density-normalized pyrite yields from all wetland sediments analyzed (n = 15) were measurable/above detection limit (>0.1 wt%), as opposed to non-wetland soils (Alfisols, Aridisols, and Mollisols), which typically fell below this limit after being density-normalized. In a non-perennially-waterlogged soil (e.g., not wetlands), it is reasonable to assume that Fe contents in sulfides are negligible.

Contrary to expectations (see Section 4.3 below), Fe_dith and Fe_ox pools did not display stronger or more robust correlations with P (Figure 9) than the other Fe pools. There were no clear correlations between P and any Fe species (for all, R² < 0.01). See Supplemental Table S4 for all Fe species results and Fe³⁺/Fe²⁺ ratios.

3.5.2. Fe Speciation, Precipitation, Soil Moisture, and Drainage

Contrary to expectations (see Section 4.2.2 below), Fe speciation pools showed no strong predictive relationships (all R² < 0.2) with mean annual precipitation (Figure 10), but there were differences between soil moisture regimes. Permafrost soils from the North Slopes of Alaska, high in organic matter and Fe, were dominated by labile Fe (Fe_asc). Aquic and Udic soils had high Fe in carbonates (Fe_acet), and Perudic soils (primarily from Iceland, with basalt parent material) had high Fe_ox (Figure 11). Mean Fe_dith was consistent between soil moisture regimes but has an extremely high range in Aridic soils.

Poorly-drained soils had a high range of Fe_asc and Fe_dith, while well-drained soils had high amounts of Fe_acet (Figure 12A). The other Fe species and drainages were consistent. The role of slope was considered, but local, small-scale topographic relief was not readily determinable for most samples.

3.5.3. Fe Speciation, Vegetation, and Soil Order

Forest soils had higher Fe_acet (Figure 13B) than the other vegetation groups, but Fe species concentrations were consistent otherwise. Oxisols had higher Fe_dith and Fe_ox than other soil orders (Figure 14C), but Fe species had little variability between soil orders otherwise.

3.6. Organic Carbon

Organic carbon (C_org) in A horizons ranged from 0 to 61 wt%, and in C horizons, ranged from 0 to 43 wt% [64]. The average C_org:P ratio for A horizons was 17.9:1, and for C horizons, 18.3:1. Relationships between C_org and CIA, P, Fe, and clay content were weak (Supplemental Figures S6 and S7). There were no strong correlations between inorganic C and Fe_tot, or between C_tot and Fe_tot (Supplemental Figure S7).

3.7. Parent Material

There is slight variability between Fe_tot and P concentrations and parent material, with igneous bedrock and ash/volcanics having the highest median values for Fe, with high P values as well (Supplemental Figure S8). While limestone parent materials show high Fe, there are only two samples in this group and they are not included in this discussion. The parent material does not appear to be a predictive control in either Fe_tot or P (Supplemental Figure S8). However, this could be due to the mixed provenance of many modern soils’ parent material (i.e., alluvium/colluvium, glacial deposits, etc. as opposed to compositionally homogeneous bedrock). A more targeted exploration into bedrock-parented soils, Fe_tot, and P could well show more well-behaved relationships.

3.8. Principle Components Analysis

Results from the principal components analyses (PCA) for the USGS data are shown in Supplemental Figure S10. Principle components 1 and 2 explain ~60% of the variance at a maximum. However, all eigenvalues are very low (<0.6) and data are essentially clustered around the origin, so while there are associations, they are very weak and should be interpreted conservatively. Latitude, P, Fe, and Al were associated with PC2 in the Top 5cm and A horizons, and vegetation and CIA were associated with PC1. Groupings changed in the C horizons, with Fe and Al associated with PC2, but P and latitude associated with C_org along PC1. In all horizons, CIA and vegetation vectors are opposite to each other, and were orthogonal to subparallel to most other vectors.

4. Discussion

In this section, we explore our results with a series of questions focused on constraining controls on Fe and total P in soils, defining their relationships on continental scales, and making inferences about how those relationships may impact terrestrial P fluxes. They are centered around key soil-forming factors, such as climate, vegetation, and weathering intensity, as well as soil redox-specific factors (i.e., precipitation, moisture and drainage, and Fe species). Because a majority of samples come from between 20 and 50°N, the interpretations made here are most robust for those latitudes.

4.1. How do Latitude, Weathering, and Clay Content Associate with P and Fe Concentrations?

Most of the expected relationships between soil order, Fe_tot and P, weathering, and latitude were supported. As expected, weathering generally decreases with latitude, dropping from a maximum CIA of >95 at 30–35° N to a CIA ~ 75 at 50° N (the northern limit of the large USGS dataset). Farther north, some B horizons deviate from that pattern, with CIA values higher than might be expected based on their climatic regime [80,81]. Interpreting the latitudinal trends here should include a caveat for latitudinal sampling bias and limitations of the dataset to, mostly, between 20° and 50° N. While the PCA results (Supplemental Figure S10) do not support a relationship between P and weathering, we interpret this discrepancy as being due to the complexities within the P-weathering relationship rather than negating the observations made between latitude, weathering, and P concentrations because the PCA eigenvalues are so low (most < ±0.2), and essentially, are clustered on the origin, which indicates a not-predictive value. P’s relationship with weathering is complex and due to a variety of factors (climate, time, slope/erosion rate, etc.), and by looking only at the end product—which is what is left in the geologic record for analysis—we must inherently work around those limitations and draw the most robust conclusions possible from limited data on these large scales.

Fe_tot behaves as expected, accumulating in B horizons as soil weathering increases, while P accumulation in all horizons peaks at moderately-weathered soils (CIA ~ 60). The latter relationship is expected because a soil that has experienced moderate weathering and has developed somewhat (e.g., Inceptisol, Alfisol) has had sufficient time to have P mobilized from the bedrock/substrate and biotically cycled, but not so long or with such intense weathering that P is depleted from the substrate and removed via erosion and/or loss of biomass. An older and/or more intensely weathered soil (e.g., Ultisol) will have a substrate more depleted of P and will have lost more P. In terms of terrestrial P transport, a moderately-weathered soil is the most likely to have a large pool of potential P to transport. An interesting next step, building off this work, could be to test these expected correlations by mapping soil P with soil age and collecting fluvial P loads.

P concentrations are, on average, far below the crustal average of 870 ppm P. Even accounting for density differences between a typical soil (ρ_soil ~ 1.62 ± 0.2 g cm⁻³; n = 659) [72] and continental crust (e.g., granite; ρ_crust ~ 2.7 g cm⁻³), all of the soil horizons are depleted in P relative to the crustal abundance. There is some variability between soil orders, but all but a small subset were below the crustal abundance (i.e., recent basalt-parented soils). When considering bulk density, the total Fe concentrations show the opposite, with the soil Fe average of 4.4% being greater than the crustal average of 3.5%, again with some variability between soil orders but overall greater than the crustal average (Supplemental Table S3). This has implications for biogeochemical modeling (see Section 4.7.3). While P concentrations could be expected to generally decrease with increased weathering intensity, a moderately-weathered soil is more likely to be at a mid-developmental stage (e.g., Inceptisol, Alfisol) and supporting vegetation that can cycle P without being P-depleted (e.g., [19]). This mid-CIA accumulation could be reflective of the vegetation’s impact on soil nutrients (e.g., mycorrhizal P mobilization, P recycling through biomass/organic matter decomposition), pointing to a key balancing point in a soil’s lifespan where the bedrock is being weathered enough to maximize fertility without yet being depleted. Forest and Altered/Cultivated soils have the highest CIA values (Figure 7), but P concentrations were relatively invariant between types of vegetation, suggesting that weathering intensity exerts a greater control on soil P than the type of vegetation.

Clay content (both percent of grain-size fraction and abundance of clay minerals) should increase with the degree of weathering a soil has undergone, and that is reflected in the results (Figure 4D). Consequently, there are secondary trends associated with clay content, i.e., clay content peaks at lower latitudes where CIA values are the highest. Fe_tot and clay content show a strong positive correlation, as expected (Figure 6E), but P shows only a weak potential decrease as clay content increases, supporting traditional thinking that P is depleted as a soil is increasingly weathered and/or older (Figure 6F). While it is kinetically favorable for P to sorb to Fe and Al oxides (which accumulate in B horizons as a soil ages), clays (often with surficial Fe-Al oxides) may be equally as important for immobilizing P in soils [82].

Interestingly, our results indicate no positive correlations between Fe oxides (Fe_dith and Fe_ox) and clay content grain-size fraction (Supplemental Figure S9), as would be expected based on soil development patterns (i.e., higher weathering intensity leads to accumulation of clay minerals and Fe oxides), suggesting that while particle size may mediate Fe speciation due to varying degrees of reactivity, it is not a primary control. While the USGS database’s clay content is based on mineralogy (determined by XRD), soils used in the Fe species work had clay content determined by settling, which would include both clay- and non-clay minerals that were clay-sized particles. This could explain at least partially the lack of the expected trend. Further exploration of clay mineralogy-specific correlations (or lack thereof) with Fe speciation pools, in the context of soil redox state, is warranted to explain this deviation from our expectations. Therefore, degree of weathering remains a primary control on soil P concentrations (e.g., [19]).

Although dust deposition is a source of P, it is typically secondary to bedrock weathering, not really applicable on anthropogenic-change timescales, and is not rapid in most of our study area (primarily the conterminous United States). Because it is not separately examined or included in this work, anyone using our results of soil P should take regional dust deposition and subsequent P replenishment into account when possible.

4.2. How Do Precipitation, Soil Moisture, and Soil Drainage Affect P and Fe Concentrations, and Fe Species?

4.2.1. How Does Soil Redox Affect P?

The relationships between soil redox factors (precipitation, soil moisture, and soil drainage) and soil P contents suggest that moisture is not a dominant control on soil P. While there is a weak possible correlation between MAP and P (Supplemental Figure S3A), the range of P concentrations for a given MAP is large (ca. 1500 ppm) and is not particularly predictive. Adding more climate data (e.g., correlating continental climate data with the USGS sample locations) could improve this potential relationship. As with other possible relationships in soils, however, it could be mediated by an additional soil-forming factor (i.e., precipitation is related to weathering intensity, and possibly to the presence of vegetation).

There exist only weak patterns between P and soil moisture or drainage. Udic soils have higher mean and ranges of P than other moisture regimes (Figure 8D), and poorly-drained soils can have higher P than well- or excessively-drained soils (Figure 8B). However, the bin sizes for these soil moisture regime types are relatively small; more data could clarify these trends. A relationship between soil P concentration and moisture regime would not necessarily be expected because soil moisture regimes are defined by the number of saturated days per year [83], but not necessarily if those days are linked by season (i.e., they do not need to be consecutive). While soil redox state would affect the (short-term) phase of P adsorption sites (such as Fe oxides) and possibly the accumulation of organic matter, which can mediate P-oxide interactions [8,9], as well as sequester P in organic compounds, only precipitation and its links to weathering intensity would be expected to exert long-term control on P concentrations.

4.2.2. How Does Soil Redox Affect Fe Concentrations and Speciation?

Our primary hypothesis was that precipitation would exert control on soil moisture and therefore on redox conditions in the soil porewaters, controlling Fe speciation. However, mean annual precipitation did not correlate with total Fe and showed none of the expected correlations with the four Fe species. There was perhaps a maximum precipitation (ca. 2000 mm yr⁻¹) for most Fe in carbonates (Figure 10B), but little else. Turning to a more direct approach of constraining soil moisture—soil moisture regime—revealed one interesting correlation: although variability in total Fe was low across soil moisture regimes (Figure 9C,D), there are some patterns of trade-off between species within a moisture regime (Aquic, Udic). Similarly, variability between soil drainage was low, but a general pattern of moderately-drained soils having higher Fe in oxides and carbonate, whereas labile Fe (reduced) showed no trend. Based on these data, soil moisture, redox state, and associated chemistry are highly heterogeneous even within soil orders, moisture regimes, and drainage ability, and broad generalizations cannot be made.

Iron oxide mineralogy and soil magnetic properties have been linked to regional precipitation [84,85,86,87] and soil moisture [42,43,88,89], but relatively few studies explore the links between soil iron oxide mineralogy and soil order [90,91]. Cervi et al. [90] found parent material and climate to be primary controls on soil magnetic susceptibility. A recent survey of Australian soils explored 471 unaltered topsoils’ magnetic properties [92], which was the only study we could find of a similar scale to the work presented here. Hu et al. [92] found parent material to be the strongest control on soil iron mineralogy, with climate (precipitation, temperature) showing weaker influence and vegetation effects on soil iron being regional.

In addition to using the large soil datasets that are increasingly available, more targeted studies on the influence of soil moisture and drainage on redox conditions—using larger datasets and including more detail of soil texture—could strengthen these trends and reveal additional patterns. Many of the soils included in this study were well- or very-well drained, with a smaller subset being poorly- or very poorly-drained. It is possible that relationships between Fe species and moisture are stronger in soil drainage regimes for which we have effectively too small datasets in this study. Additionally, highly localized slopes (along with the context of grain size composition) likely exert strong controls on soil moisture and redox state, but we were not able to expand this study to include those data. Doing so would help clarify potential relationships between soil moisture and Fe species. Using high-resolution slope data (e.g., fine-scale LiDAR/digital elevation models) could be useful for determining these relationships between local soil redox conditions and geochemistry.

4.3. How Do Fe Species Associate with P?

We expected that Fe oxides, having been shown to efficiently sorb P onto their surfaces in other soils (e.g., [8,9,44,93]), would have a more robust correlation than between either labile Fe and P (in reducing microenvironments) or Fe in carbonates and P. However, no species showed a strong correlation with P (Figure 9). Overall, Fe_tot and P in all soils showed a positive correlation, suggesting that Fe species is not a dominant control on P concentration, and that no single species is more likely to be associated with P, at least on the regional scales examined in this work. Because of the weak relationships between Fe species and assorted soil moisture proxies, as well as a surprising lack of correlation with clay content, there is no simple rule connecting Fe species, soil redox state, and P. However, an increase in any Fe-bearing species should correspond to an increase in P due to the general positive relationship (Supplemental Figure S4). Fe (and Al) oxides in the clay fraction specifically may be most important for terrestrial P transport (both in soils and fluvial systems) (e.g., [39,44,94]). We expect that repeated analyses within an individual profile or on a small local scale (i.e., meters) would show stronger expected relationships, but they cannot be generalized to regional or continental scales.

4.4. How Does Vegetation Affect P and Fe Concentrations?

Biology—here, generalized as dominant vegetation on the soil surface (not considering microbial life and mesoscale organisms)—plays an important role in terrestrial P cycling. We found that Barren soils’ B horizons had more P accumulation than other, vegetated soils’ B horizons, reflecting plants’ effectiveness at moving P from regolith to higher in the soil profile ([15,16]; see Section 4.7.4). The other vegetation types showed surprisingly little variability, but P does tend to increase up-profile, suggesting that once a landscape is vegetated—regardless of specific ecosystem type—P is efficiently mobilized by plants to the upper horizons and recycled and stored there. Additionally, vegetated soils had far higher ranges/maximum P concentrations than unvegetated soils, supporting the importance of plants’ role in P mobilization and accumulation. Barren landscapes lack the rooting systems to mobilize and transport P, hence the increased accumulation in B horizons. Halsted and Lynch [95] found essentially no discrimination in P uptake between C₃ and C₄ plants, which our results further support.

Aside from ‘Other’ vegetation (which was biased by Fe-rich parent material), Cultivated soils’ B horizons were the only ones that showed depletion, suggesting that ecosystem variability is less important than parent material and soil mineral accumulation rates for Fe in soils. Deep tilling in agricultural practice could be an additional factor in lower-soil-profile nutrient depletion, but it is less common than it used to be. Similarly, vegetation type does not appear to be a direct control on Fe species within soils; Forests did have greater Fe_acet, which is likely linked to the defined property of Alfisols of high cation exchange capacity and cation mobility, as well as their tendency to be well-drained.

Overall, the presence of vegetation—rather than differences in the type of vegetation—exerts the most control on P accumulation in soils (Figure 2). Although the mean P values were similar between vegetation types within a soil horizon (with the exception of B horizons in Barren soils), non-Barren (vegetated) soils had higher ranges of P than Barren soils. Forests also have higher CIA values than grasslands and shrublands (Figure 7), suggesting that the evolution and spread of different ecosystems could have changed the distribution of weathering intensity. This has interesting implications for the ongoing debate on land plants’ influence on continental P fluxes (see Section 4.7.4).

Vegetation changes with latitude; the soils and ecosystems represented in this work span tropical and temperate forests, continental grasslands, dry shrublands, and high-latitude, little-vegetated landscapes. This inherently contributes to the observed latitudinal trend in weathering because different plants have different interactions with soil and bedrock. While there are ground-level differences between plants within a plant type (e.g., grasses), the functional difference between plant functional types (e.g., grasses vs. trees) are far larger and more meaningful on continental scales. While at local or regional scales, high-resolution, species-specific studies on weathering and plant-soil relationships can be done, that level of analysis is neither possible nor necessary for continental-scale analysis of vegetation-weathering relationships. Global climate, biogeochemical, and environmental change models rely on continent-scale data; therefore, while locally, species variability within a plant functional type may lead to nuanced plant-soil geochemical relationships, the continental scale and umbrella vegetation type is appropriate for this style of inquiry.

4.5. Is Soil Order Predictive of P and Fe Concentrations?

Soil order is not quantitatively predictive of P and Fe concentrations. Soil order, Fe, and P are linked indirectly through the correlation between Fe/P and CIA; this relationship is stronger for Fe, which shows a clear Fe accumulation increasing with weathering (Figure 6A,C), which is also linked to the Inceptisol to Ultisol soil progression common in chronosequences. Soil order and weathering show a very well-behaved and expected relationship, with CIA values increasing from Entisols to Ultisols/Oxisols (Figure 5). An exception to this trend is high-latitude Gelisols, which have longer formation times and lower weathering rates. These results align with existing soil formation frameworks and support proxies that rely upon their use, such as the CIA-based mean annual precipitation proxy, the paleosol weathering index for mean annual temperature, and Bt thickness ([79] and refs. therein). However, because soil order depends on a variety of factors (i.e., weathering stage, vegetation), quantitative relationships between soil order and either Fe_tot or P concentrations cannot be generalized, and samples within one soil order should not necessarily be used to predict Fe or P concentrations quantitatively for other samples of the same order. This is important to remember when using limited numbers of samples of fossil soils (see Section 4.7.1).

4.6. Implications for P in Modern Soils, Climate Change, and Soil Fertility/Food Security

4.6.1. Soil P, Erosion and Transport, and Human Activity

Aside from the geologically-recent uptick in anthropogenic influence on landscapes, continental scale-vegetation was relatively stable over the past millennium in North America (e.g., [96]). Concern for changing P fluxes in terrestrial ecosystems largely stems from the occurrence of harmful algal blooms in rivers, lakes, and coastal areas, often linked to mass-production agriculture and nutrient-rich runoff. The implications for these fluxes from this work center on land use change and soil degradation, as well as climate change-driven shifts in vegetation (e.g., grassland to barren). P would also be lost in that shift because vegetation loss exacerbates erosion and landscape degradation [97].

Human-induced changes (e.g., industrial agriculture, extensive monoculture, unsustainable farming practices, etc.) can lead to natural P being lost through subsidence and erosion, P being added through fertilizer, and reduced plant P retention if natural, diverse ecosystems are replaced with monoculture (e.g., [98]). Land-use changes typically result in the removal of the upper portion of a soil profile (various, e.g., [99]). Because P tends to accumulate in the upper horizons (O,A), our results suggest that this could lead to a noticeable loss in P and in agricultural soils, increased need for fertilizer (see Section 4.6.3). Land-use changes in vegetated soils could increase the P flux from land to lakes and coastal waters, leading to eutrophication. The spread of Aridisols through desertification could mean more regions experience erosion due to vegetation loss, but perhaps in turn lower weathering due to that loss of vegetation and to drier climatic conditions [100]. As desert regions expand, global dust volumes and distributions will change, which would change which areas see P removal (through increased erosion and P loss through dust) and which see P accumulation through dust deposition. Changes in atmospheric and soil CO₂ concentrations are likely to affect both plant and soil productivity (as well as mycorrhizal efficiency) under climate change, although the extent and duration of changes to productivity are debated [100,101,102,103].

4.6.2. Weathering, Climate Change, and Soil P

In addition to changes in vegetation (ecosystem), climate change will impact weathering intensity and erosion, affecting global soil P reservoirs. As a consequence, soil fertility and food security risk will shift and vary depending on how different regions and biomes respond to climate change [104]. Currently, there is a general trend of intense, rapid weathering at lower latitudes as a consequence of high heat and precipitation. Higher latitudes tend to be cooler and drier, and therefore, weathering intensity and rates are lower [14,24,105]. While low-latitude areas are expected to undergo desertification (decreased precipitation, loss of vegetation), higher latitudes may experience an increase in precipitation and temperature. Therefore, the style and intensity of weathering will change regionally. Based on our results, where the range in P concentrations peaks at mid-intensity weathering (i.e., CIA ~ 60; Figure 5) and decreases in older and/or more strongly-weathered soils (as expected), the climate change expected for higher-latitude regions may result in decreased soil fertility (limiting agricultural production, see Section 4.6.3) and increased risk of eutrophication in surrounding waters. The former effect could potentially be mediated by vegetation, which holds P and slows erosion.

Except for the increase in range at mid-intensity weathering, P concentrations were relatively consistent across weathering intensities, which was unexpected. We had hypothesized that as older soils tend to be depleted in P, a higher weathering value would correspond to lower P. This was not demonstrated by our results. We interpret the lack of correlation as highlighting the significance of soil age—rather than weathering intensity (climate)—in controlling soil P, as some previous literature has suggested (e.g., [19,106]). A nuance that this large dataset could be missing is P speciation, which could still show P phase-specific variability with weathering intensity. The implication for P transport from soils remains unchanged from previous studies, where older soils would likely have smaller pools of P for mobilization. Mapping soil P with soil age on a large scale (such as the USGS dataset used here) would be a valuable new addition to the field and could elucidate the P-soil age-weathering puzzle. Additionally, exploring relationships between soils and weathering with species-level variability (on smaller, local scales) could help refine interpretations based on dominant regional plant functional type.

The difference in weathering intensity/CIA values between Forests and other natural vegetation types (Grasslands, Shrublands), as well as the high CIA values in Altered/Cultivated soils suggests that if more land is converted to farming (cultivated) from natural grasslands (e.g., the Great Plains) or if agricultural intensity increases, weathering rates would increase, depleting the soils’ nutrients more quickly and perhaps increasing P fluxes. To our knowledge, this is the first data-based demonstration of the principle on the continental scale rather than as a model result. Understanding how human activity and agricultural degradation of soils affects weathering and resulting P is crucial for predicting changes in terrestrial P.

4.6.3. Soil P, Plants, and Agriculture

Climate change is also likely to affect plant physiology and biological processes, such as P uptake (mediated by mycorrhizal fungi) and C storage, as a result of changing atmospheric carbon dioxide concentrations (pCO₂). Plants rely heavily on symbiotic relationships with mycorrhizae to mobilize and uptake recalcitrant mineral P and fix nitrogen [102,107,108,109]. Mycorrhizal activity and P uptake can also be affected by soil moisture [110,111], which depends on a number of factors (discussed in Section 2.2) that could change in response to altered land use and climate change. P uptake efficiency also varies between species of symbiont [112], which may respond differently to climate stressors such as pCO₂ and temperature [110,113,114]. Some work suggests that increased atmospheric pCO₂ will be advantageous for plant growth, leading to greater terrestrial biomass and C storage [115,116], though potentially with an upper limit on pCO₂ vs. growth returns [117]. However, observations have shown that increases in pCO₂ during and after industrialization did not always increase plant C richness or biomass growth [118,119,120], and the biomass increase effect can still be mediated by P (and N) limitation in soils [121,122,123] and mycorrhizal species [102,114].

Ecosystem-specific C:N:P (Redfield) ratios vary, but throughout natural marine and terrestrial systems, P is the limiting nutrient [4,5]. For example, while the original marine phytoplankton-based Redfield ratio is 106:16:1, in soils, it is naturally higher at 186:13:1 [124]. The C:P ratios in both A and C horizons (16.7:1 and 14:1, respectively) were far lower than what is observed for O horizons (186:1 [124]) or inland streams/rivers (167:1 [125]), and are closer to ratios for microbial biomass in soils (~50:1 [124]). However, this ratio was derived with total P rather than organic P, and is therefore likely an overestimation. In these samples, the C_org:P_org may be closer to the expected value for soils, but further speciation work is needed. Additionally, it is reasonable to expect that mineral soil horizons (rather than organic-rich horizons, like O) would have lower C_org and nutrient concentrations. For developed soils (mean C_org:P = 15.25), this discrepancy could be linked to deep tilling removing nutrients from subsurface horizons and potentially increasing nutrient fluxes as compared to undisturbed or shallow-tilled soils. For natural soils with the potential for cultivation/deep tilling (e.g., grasslands, mean C_org:P = 20.5), the lower C:P ratio in subsurface mineral horizons could be inferred as diminishing returns for deeper tilling. There were no strong correlations between C_org, C_inorg, or C_tot and variables of interest (CIA, clay content, P, Fe) (Supplemental Figures S6 and S7).

Complicating this response is the fact that increases in atmospheric pCO₂ could increase the rate of bedrock weathering (e.g., [126]), which mediates P availability and biolimitation. Moderate increases in pCO₂ and weathering could lead to higher P availability and greater plant biomass, but if weathering increases too much, P in soils could be depleted while plant growth continues to increase. If high weathering rates lead to faster or premature depletion of P, especially exacerbated by anthropogenic activity (agriculture), crop stresses and regional food shortages could occur. Additionally, demand for P as fertilizer could increase; fertilizer P is ultimately sourced from bedrock, which, on human timescales, is a nonrenewable material [127]. Without developing efficient P recycling methods or adding sources such as bone char (e.g., [128,129]), P could run out in a matter of decades [130,131,132]. Therefore, when considering how climate change will affect agricultural yields, it is important to consider plant physiology and P uptake mechanisms together with mineral/geochemical changes in the soil driven by natural and anthropogenic forces.

Overall, this work aids in predicting changes in terrestrial P fluxes primarily by linking soil P, weathering intensity, and the presence of vegetation. However, in this study, “vegetation” remains something of a black box; the mechanism by which the presence of vegetation affects soil P could be through slowing erosion, increasing biotic P cycling, retaining P via plant/microbial biomass, increasing bedrock apatite weathering, or—most likely—a combination of those factors. Exploring each of these factors at large spatial scales, and to the extent possible pairing those analyses with smaller-scale species-specific soil-plant relationships, is a critical next step in elucidating the potential for P fluxes in different terrestrial regions.

4.7. Implications for P in the Fossil Record of Soils and Its Geologic Use

4.7.1. Heterogeneity and Paleosol Representativeness

Fossilized soils (paleosols), which are present in the rock record as far back as at least 3.0 billion years ago, serve as important windows into Earth’s terrestrial past [133]. They have been used to reconstruct ancient atmospheres (e.g., [134,135,136,137,138]), climates ([79] and refs. therein), and terrestrial biology and biogeochemical cycling (e.g., [139,140]). The geochemical composition of fossil soils is used to reconstruct a range of climatic and environmental changes; however, these tools are primarily based on small modern soil datasets. This work provides ranges of reasonable values for soil chemistries under known environmental and atmospheric conditions, providing critical background information to improve our paleoclimate and paleoenvironment reconstructions. Based on the highly variable geochemical results in this work, we urge researchers using paleosols for these types of reconstructions to be cautious in assuming a single paleosol profile to be representative of an entire landscape or basin. More work on how representative a paleosol profile is, chemically and climatically, should be done to incorporate landscape-scale variability in soils (e.g., [141]).

4.7.2. Paleosol Fe, Atmospheric Oxygen Reconstructions, and Microbial Life

Because paleosols form on the Earth’s surface and in direct contact with the atmosphere, their Fe chemistry has been used to reconstruct atmospheric oxygen levels [134,138,142]. Across a range of soil orders, environments, and climates, non-wetland soils in this work reflect oxidizing conditions (based on Fe³⁺/Fe²⁺ ratios; Supplemental Table S4). While Fe³⁺/Fe²⁺ ratios are no longer the primary tool for reconstructing paleo-oxygen levels, these results indicate that that tool is qualitatively robust in the absence of other observations such Cr or S isotopes.

A concern with using fossil soils for reconstructing atmospheric oxygen is that these soils likely hosted microbial life, which could have affected the redox signal (‘biosignature’) left behind. Indeed, redox-sensitive metals have been proposed as biosignatures because they can be mobilized by biologic processes and redeposited within a soil profile [143]. However, if it is suspected that the soils have variable redox states due either to high moisture input or poor drainage, how reliably might they have recorded the atmospheric oxygen signal as opposed to microscale oxic/reducing conditions due to microbial activity? Modern biological soil crusts (BSC; or, cryptogamic soils) are symbiotic communities of microbes (predominantly cyanobacteria), algae, and fungi; these communities have been proposed as analogues for early life on Earth (and Mars) [143]. Here, we found that the soil redox state being recorded by BSCs was oxic (again, using Fe³⁺/Fe²⁺ ratios and dominant Fe speciation; Supplemental Table S4), despite BSCs’ ability to retain water and temporarily shift to reducing conditions. It is unlikely that the microbial communities would have ‘overprinted’ the atmospheric oxygen signal even if they temporarily held water to cycle nutrients and leave behind biosignatures.

4.7.3. Continental Weathering, Nutrient Fluxes, and the Atmosphere

Geologists are interested in constraining many of the processes described above in ancient soils and ecosystems. Weathered P, transported to the oceans via fluvial networks and continental drainage, is typically invoked as a first-order control on marine productivity, which is associated with organic C and organic-associated and inorganic carbonate sedimentation, atmospheric CO₂ drawdown, and oxygen production [56]. Continental weathering, and therefore P flux, is a central control in many biogeochemical models of Earth’s past. However, the continent-to-ocean transport of P is usually prescribed and not based on actual changes in fossil soil (paleosol) P values or weathering intensity data. Rather, it has been generalized based on crustal P values and fluxes controlled by large-magnitude changes, such as global glaciations [144,145] or limited modern observations [61,62].

Based on our results, continental weathering intensity and resulting elemental fluxes should vary with paleogeography (latitude effect), in addition to changes in climate and landmass area. Clay content through time could be taken into account and paired with CIA to reconstruct actual weathering intensity where those data are available in the rock record [146]. Forests also have higher CIA values than grasslands and shrublands (Figure 7), suggesting that the evolution and spread of different ecosystems could have changed the distribution of weathering intensity. Biogeochemical models should take these differences into account as opposed to assuming uniform behavior spatially.

Continental sulfide weathering in the Precambrian (pre-540 Ma) has been invoked as a major source of sulfate to the oceans, affecting organic matter respiration rates, marine alkalinity and oxidation state, and ultimately the concentration of atmospheric oxygen [147,148,149,150]. At the Great Oxygenation Event (GOE) 2.45 billion years ago, pyrite burial and increased ocean sulfate concentrations are thought to have decreased atmospheric methane through microbial sulfate reduction, lowering the greenhouse effect and leading to global glaciations [151,152]. Our results (Supplemental Table S5) constrain the number of Fe-bound sulfides in modern soils, forming under an oxic atmosphere, and suggest the potential for a ‘false positive’ for a reducing atmosphere if a paleosol is not correctly identified as waterlogged (and therefore reducing) or histic. Although it should be noted that using the mere presence of minerals, such as pyrite or uraninite, in the rock record as reducing indicators [153] is simplistic and no longer widely-used in the field.

4.7.4. Vegetation and P in the Phanerozoic (542 Ma Onwards)

The rise of land plants during the Devonian period (419–349 Ma) has been suggested as a major driver of change to the Phanerozoic P and C cycles, which could have affected atmospheric oxygen [15,56,63,154]. While B horizons in Barren vegetation had increased P retention, the lack of variability between other vegetation domains suggests a somewhat binary response—either soils are vegetated and mobilize/accumulate P, or they are not vegetated and are relatively depleted in P. Only soils with vegetation accumulated high ranges of P, supporting the need for biological mediation of P in terrestrial settings.

5. Conclusions

From this large-scale analysis of modern soils’ physical and chemical properties, we draw conclusions about the relationships between modern soil biogeochemistry and climate, and link those to terrestrial P transport. We also describe implications about past terrestrial biogeochemistry and interpreting the terrestrial rock record.

In modern soils, some of the most well-defined (and expected) relationships were between soil composition and weathering. Latitude, clay content, and weathering correlate as expected, with higher latitudes corresponding with lower weathering intensity and clay content. Contrary to our expectations, average P concentrations are relatively consistent across weathering intensities but show a slight peak at middle weathering intensities, which supports commonly applied models of soil age and P loss. Specifically, it helps to clarify the conflation of time and weathering intensity that can occur, pointing to the importance of soil age as a control on soil P depletion. Maximum terrestrial P transport, then, would still likely occur at a midpoint in the weathering-P spectrum, where a soil is mature enough to have mobilized P from apatite and be linked to high erosion rates, but before soil P is too depleted due to age (which would be a lower P flux overall). Contrary to expectations, P in soils was not strongly associated with Fe (oxyhydr)oxides. The scale of such relationships may be much smaller than continental (i.e., locally or per profile). Additionally, no strong, predictive relationships were present between Fe species and precipitation, soil moisture, or drainage, so predicting P transport based on climate or soil redox properties is not possible. Finally, the presence of vegetation (but not plant functional type or generalized ecosystem type) is important for P accumulation in soils, which has implications for how the rise of terrestrial plants may have changed P cycling through geologic history.

Based on modern relationships between soil P, weathering, and vegetation, the spread of land plants in the early Phanerozoic likely increased P accumulation in soils and on continents, rather than increasing the flux of P to the oceans. This would limit the role of terrestrial-marine P transport in marine productivity. Because of latitudinal variability in weathering intensity, biogeochemical models should take paleolatitude into account when using weathering intensity as a driver for P fluxes. Additionally, models should use density-normalized values for terrestrial P rather than bulk crustal values and should use the range of P in modern soils as a quantitative constraint. Finally, while we did not find the relationships between Fe species and climate that we had expected, the overall Fe signature in soils (including microbially-colonized) was oxic despite analyzing a range of soil moistures and drainages, which supports the use of paleosol Fe/redox geochemistry for paleo-atmosphere reconstruction.