Articles | Volume 10, issue 1
Original research article
06 Feb 2024
Original research article |  | 06 Feb 2024

Soil carbon, nitrogen, and phosphorus storage in juniper–oak savanna: role of vegetation and geology

Che-Jen Hsiao, Pedro A. M. Leite, Ayumi Hyodo, and Thomas W. Boutton

Woody-plant encroachment into grasslands and savannas has been globally widespread during the past century, likely driven by interactions between grazing, fire suppression, rising atmospheric CO2, and climate change. In the southernmost US Great Plains, Ashe juniper and live oak have increased in abundance. To evaluate potential interactions between this vegetation change and the underlying soil parent material on ecosystem biogeochemistry, we quantified soil organic carbon (SOC), total nitrogen (TN), total phosphorus (TP), and δ13C of SOC in soils obtained from trenches passing through grassland, juniper, and oak patches on soils lying atop the respective Edwards and Buda limestone formations in central Texas. Soils on the Edwards formation are more shallow and have more rock outcropping than those on Buda. The δ13C values of SOC under grasslands was 19 ‰, whereas those under woody patches were 21 ‰ to 24 ‰, indicating that wooded areas were relatively recent components of the landscape. Compared with grasslands, areas now dominated by juniper or oak had elevated SOC, TN, and TP storage in soils lying atop Edwards limestone. In Buda soils, only oak patches had increased SOC, TN, and TP storage compared with grasslands. Woody encroachment effects on soil nutrients were higher in soils on the Edwards formation, perhaps because root and litter inputs were more concentrated in the relatively shallow layer of soil atop the Edwards bedrock. Our findings suggest that geological factors should be considered when predicting nutrient store responses in savannas following vegetation change. Given that woody encroachment is occurring globally, our results have important implications for the management and conservation of these ecosystems. The potential interactive effects between vegetation change and soil parent material on C, N, and P storage warrant attention in future studies aimed at understanding and modeling the global consequences of woody encroachment.

1 Introduction

Woody-plant encroachment into grasslands, savannas, deserts, and other semiarid and arid ecosystems is a globally extensive land cover change that has been occurring during the past 150 years (Archer et al., 2017; Sala and Maestre, 2014; Stevens et al., 2017). This important vegetation change is driven by several potentially interacting local and global phenomena, including the elimination of naturally occurring fires, chronic livestock grazing, rising atmospheric CO2 concentrations, and climate change (Archer et al., 2017; Rosan et al., 2019; Stevens et al., 2017; Venter et al., 2018; Zhou et al., 2019). Rates of increase in woody cover range from 0.1 % to 2.3 % yr−1 in grasslands and savannas throughout the world (Archer et al., 2017; Deng et al., 2021). These changes in the relative abundances of woody plants and grasses often have significant ecological impacts on arid and semiarid ecosystems. For example, increases in above- and belowground inputs of organic matter derived from woody plants often alter soil carbon (C), nitrogen (N), and phosphorus (P) storage, cycling rates, and stoichiometric relationships (Barger et al., 2011; Eldridge et al., 2017; Finzi et al., 2011; Zhou et al., 2018a).

Little is known regarding how geology might interact with changes in woody-plant cover to influence ecosystem biogeochemistry. Geological factors play key roles in determining fundamental soil properties, including mineralogy, chemistry, and texture, which in turn can be important determinants of soil C, N, and P biogeochemistry (Augusto et al., 2017; Farella et al., 2020). The reactivity of the soil mineral phase can change the stabilization potential of soil organic carbon (SOC) under similar vegetation (Doetterl et al., 2015; Torn et al., 1997), and clayey soils have greater potential to retain mineral-derived nutrients and organic matter across a broad range of soil substrates (Soong et al., 2020). Furthermore, the presence of living organisms (plants and microbes) in the weathering zone of geological materials can influence soil formation and properties through (i) mechanical processes, (ii) changes in soil pH, and (iii) root and microbial exudates that influence hydrolysis reactions and chelation (Anderson, 1988). As plant species can differ greatly with respect to root distribution patterns, root strength, and root exudate production (Canadell et al., 1996; Dietz et al., 2020; O'Keefe et al., 2021), their interactions with parent material may be different and, thus, have important influences on soil development and biogeochemistry.

Throughout the Great Plains of North America, the prevalence of Juniperus spp. has increased during the past century in areas that were previously grass-dominated (Leis et al., 2017; Smeins and Merrill, 1988; Van Auken and Smeins, 2008). Between 1990 and 2010, rates of Juniperus encroachment reported for this region ranged from approximately 30–40 km2 yr−1 in Kansas and Oklahoma to 85 km2 yr−1 in Nebraska (Meneguzzo and Liknes, 2015; Wang et al., 2018a). In Texas, Juniperus cover increased at a rate of 441 km2 yr−1 between 1948 and 1982 (Ansley and Wiedemann, 2008). These dramatic changes in Juniperus cover have been shown to alter key ecosystem properties such as primary productivity, evapotranspiration rates, soil water infiltration and percolation depth, and pool sizes and flux rates of essential elements (Jessup et al., 2003; Leite et al., 2020; McKinley and Blair, 2008; Wang et al., 2018b; Zhang et al., 2023).

In the southernmost portion of the North American Great Plains, the Edwards Plateau spans approximately 100 000 km2 of central Texas and is covered with savannas and woodlands dominated primarily by the woody species Juniperus ashei Buchholz (Ashe juniper) and Quercus virginiana Mill. (live oak) as well as grasses and forbs typical of the southern mixed-grass prairie (Smeins et al., 1976; Smeins and Merrill, 1988; Taylor et al., 2012). Although juniper and oak have been long-term woody components of this landscape, juniper populations in particular have increased markedly during the past century, largely due to reduced fire frequencies caused by chronic livestock grazing and fire suppression policies (Ansley and Wiedemann, 2008; Bendevis et al., 2010; Leite et al., 2020; Wilcox et al., 2007; Marshall, 1995). Woody-plant cover in this region increased from 12 % to 40 % between 1949 and 1990, with over half of this increase attributable to Ashe juniper (Fuhlendorf et al., 1996). Despite these dramatic vegetation changes, little is known about how woody encroachment in the Edwards Plateau region might influence soil storage of C, N, P, and other important elements. In addition, soils of the Edwards Plateau are derived from multiple geological formations, and it remains unknown how woody encroachment might alter nutrient stores in soils derived from different parent materials. At the western edge of the Edwards Plateau, most soils lie atop the Buda and Edwards geological formations. Soils on the Edwards formation are shallow with extensive areas of rock outcropping, whereas those on the Buda limestone are generally deeper and contain carbonate-rich layers on top of the hard limestone (Gabriel et al., 2009; Wilcox et al., 2007).

The purpose of this study was to evaluate the long-term impacts of woody-plant encroachment on the biogeochemistry of soil C, N, and P throughout soil profiles lying on top of the respective Edwards and Buda limestone formations in juniper–oak savannas of the Edwards Plateau in the southern Great Plains. To accomplish this, soil samples were collected to the depth of bedrock from multiple trenches that passed through both grassland and woody-plant communities. We used elemental analyses to quantify pool sizes of SOC, total N (TN), and total P (TP) in soils, and we measured the δ13C values of those same soils to verify woody encroachment and to estimate the relative proportions of the respective grass and woody sources of soil organic matter. Grass species in the Edwards Plateau generally use C4 photosynthesis, resulting in different δ13C values compared with the encroaching woody species, which use C3 photosynthesis (Boutton et al., 1998). Changes in the δ13C signatures at different soil depths can reveal historical shifts in predominant vegetation, possibly due to climatic changes, disturbances, or human activities (Jessup et al., 2003; Zhou et al., 2019). We tested the following hypotheses: (1) woody plants are relatively recent components of grasslands on the western Edwards Plateau; (2) soils under woody plants will have larger pool sizes of SOC, TN, and TP than grass-dominated areas; and (3) vegetation cover will interact with geological substrate to influence SOC, TN, TP, and soil inorganic carbon (SIC).

2 Materials and methods

2.1 Site description

Research was conducted at the Texas A&M AgriLife Sonora Research Station on the western edge of the Edwards Plateau, Texas (3015 N, 10033 W; altitude 670 m). The mean annual temperature was 19.3 C and the mean annual precipitation was 476 mm from 2012 to 2020 (Cooperative Climatological Data Summaries, 2022). Summers are hot and dry with an average July temperature of 30 C, while winters are mild with an average January temperature of 10 C. Most precipitation occurs during May–June and September–October, and precipitation is highly variable between and within years. The Köppen–Geiger climate classification for this region is hot arid steppe (Beck et al., 2018).

The Edwards Plateau is an uplifted and dissected limestone plateau (karst topography) with gentle slopes. Soils in this region are clayey Mollisols with shallow soils on plateaus and hills and deeper soils on plains and valley floors (Gabriel et al., 2009; Wiedenfeld and McAndrew, 1968). The predominant soil map units on plateaus and hills are the Eckrant–Rock outcrop complex and the Prade–Eckrant complex. These include the commonly occurring Harper (clayey, smectitic, thermic Lithic Haplustolls), Prade (clayey-skeletal, smectitic, thermic, shallow Petrocalcic Calciustolls), and Tarrant soil series (clayey-skeletal, smectitic, thermic Lithic Calciustolls), all of which lie atop the Edwards formation (Wilcox et al., 2007). These soils contain large amounts of limestone fragments and limestone outcrops. Depth to bedrock for these soils was generally <0.4 m (Fig. 1 and Figs. S1 and S2 in the Supplement). The clay content in the top 5 cm of Edwards soil is 30 %–40 % and increases with depth to almost 50 % at 20 cm depth (Marshall, 1995). The soil map unit commonly occurring on plains and valley floors is Valera clay, including the Rio Diablo (fine, mixed, superactive, thermic Aridic Haplustolls), Ozona (loamy, mixed, superactive, thermic, shallow Petrocalcic Calciustolls), and Mereta (clayey, mixed, superactive, thermic, shallow Petrocalcic Calciustolls) soil series, which all lie atop the Buda formation (Wilcox et al., 2007; Gabriel et al., 2009). These soils are generally deeper than those lying atop the Edwards formation and contain hard limestone and a caliche layer on top of the limestone bedrock (Figs. 1, S1, S2). The caliche layer, typically light in color, can manifest as either the petrocalcic (Bkkm horizon) or the paralithic (Cr) layer, both potentially cemented or unconsolidated. The petrocalcic layer is comprised of weathered carbonate and cements the B horizon with soil particles (Soil Survey Staff, 2014), while the paralithic layer is a residuum from limestone weathering (Rabenhorst and Wilding, 1986a) (Fig. S3 in the Supplement). The clay content in Buda soil is generally 2 %–5 % lower than that in Edwards soil throughout 0–40 cm profile based on USDA/NRCS data (Official Soil Series Descriptions (Rio Diablo series), 2022; Official Soil Series Descriptions (Eckrant series), 2022; Official Soil Series Descriptions (Valera series), 2022; Official Soil Series Descriptions (Prade series), 2022). Depth to consolidated bedrock for Buda soils was approximately 1.5–2 m, while the depth to the caliche layer was approximately 0.5 m. Although the hard limestone geological formations underlying both the Edwards and Buda soils have contributed somewhat to the formation of these soils, there is considerable chemical and physical evidence indicating that these soils are derived largely from an overlying limestone residuum with distinctly different attributes than the underlying limestone (Rabenhorst and Wilding, 1986a, b; Cooke et al., 2007). The Del Rio Clay, an Upper Cretaceous marly limestone that locally overlies the Edwards limestone, has been proposed as the dominant source of these soils (based on texture, mineralogy, and Nd isotope composition), at least on the eastern portion of the Edwards Plateau (Cooke et al., 2007).

Figure 1Shown are Trench 1 in Edwards soil (a) and Trench 5A in Buda soil (b) at Texas A&M AgriLife Sonora Research Station on the Edwards Plateau, Texas. The beige color in panel (a) is limestone. The yellow-whitish layer in panel (b) is the paralithic layer in which the individual soil particles are aggregated with calcium carbonate.


The study area lies at the southernmost extent of the Great Plains, and the contemporary vegetation of the site is typical of this region. Dominant woody species are live oak and Ashe juniper. Other less abundant woody species include Quercus pungens Liebmann (scrub oak), Juniperus pinchotii Sudw. (redberry juniper), and Prosopis glandulosa Torr. (honey mesquite). Woody species occur commonly in small clusters or as isolated individuals within the surrounding grassland matrix. Dominant grass species include Hilaria belangeri (Steud.) Nash (curly mesquite), Bouteloua curtipendula (Michx.) Torr. (sideoats grama), Aristida spp. (three awns), and Nassella leucotricha (Trin. & Rupr.) Pohl (Texas wintergrass). Grazing in this area was heavy (2–5 ha per animal unit per year) to moderate (6–8 ha per animal unit per year) from approximately 1880 to 2010, but the area chosen for this study has been grazed lightly (9–15 ha per animal unit per year) and intermittently for the past 10 years (Leite et al., 2020). Since 1948, the livestock composition in these grazed areas has maintained an approximate ratio of 60:20:20 with respect to cattle, sheep, and goats, respectively (Marshall, 1995).

2.2 Field sampling

In February 2019, using a backhoe, soils were sampled from three trenches lying atop the Edwards formation and three lying atop the Buda formation. Trench locations were selected to represent (i) the major geological formations underlying the soils in the study area (Edwards and Buda limestone), based on the USDA/NRCS Web Soil Survey (Web Soil Survey, 2022), and (ii) the major vegetation types, including live oak, Ashe juniper, and herbaceous species occurring on the area. Trench lengths and depths were variable, depending on the depth to bedrock and the distribution of bare rock at the surface (Table 1). At 2 m intervals along each trench face, soil samples were collected from multiple depth intervals (0–10, 10–20, 20–40, 40–70, 70–100, 100–120, and 120–150 cm), although depth to bedrock sometimes prevented acquisition of the deeper depths in some trenches, particularly those atop the Edwards formation. Soil depth intervals were selected to represent soil horizon patterns described in USDA/NRCS Web Soil Survey and for comparability with previous studies of the study area. We aimed for the midpoint within each specified depth range for sampling. For instance, samples were taken at 5 cm for the 0–10 cm range, 15 cm for the 10–20 cm range, 30 cm for the 20–40 cm range, and so forth. Each sampling point involved the horizontal insertion of two soil cores (7.6 cm width × 10 cm length) into the trench face. One core was used for the determination of soil bulk density and root biomass, while the other was used for pH, elemental concentrations, and stable isotope analyses. Aboveground vegetation cover within 2 m of each sampling location was mapped and categorized as grassland, Ashe juniper, live oak, or mesquite based on dominant vegetation types. In order to characterize the isotopic composition of plant inputs to the soil organic matter pool, approximate equal amounts of live leaves from three individuals of Ashe juniper, live oak, and the dominant grass species were collected near each trench. Litter from Ashe juniper, live oak, and dominant grass species as well as standing dead grass materials were collected from 0.25 m×0.25 m plots with three replicates close to each leaf sampling location.

Table 1Geological substrates and physical dimensions of soil trenches.

Download Print Version | Download XLSX

2.3 Laboratory analyses

One soil core from each sampling point was oven-dried at 105 C overnight. The core was then passed through a 2 mm sieve to remove rock fragments and retrieve fine roots for quantification of root densities. Soil bulk density was measured and reported on a gravel-free basis (Culley, 1993). The other soil core was air-dried, passed through a 2 mm sieve to remove rocks and large organic fragments, and subsampled for soil pH using a 1:2 soil/0.01 M CaCl2 slurry. Another subsample of the <2 mm fraction was pulverized in a centrifugal mill (Angstrom Inc., Belleville, MI, USA) and used to quantify concentrations of SOC, SIC, TN, and TP as well as the δ13C of SOC. Leaf and litter samples were oven-dried at 80 C for 48 h and pulverized in a centrifugal mill.

Pulverized leaf tissues, litter, and soils were analyzed for C and N concentrations and δ13C values by dry combustion using a Costech ECS 4010 elemental analyzer (Costech Analytical Technologies Inc., Valencia, CA, USA) interfaced via a ConFlo IV with a Delta V Advantage isotope ratio mass spectrometer (Thermo Scientific, Bremen, Germany) at the Stable Isotopes for Biosphere Science laboratory, Texas A&M University. Pulverized subsamples were weighed into tin capsules to measure the percentage of total C and total N. For soils, a second subsample was weighed into a silver capsule and treated with 3 N HCl to remove CaCO3 until no reaction occurred (Nieuwenhuize et al., 1994). This subsample was used to measure the concentration and δ13C value of SOC. Soil inorganic carbon concentrations were computed as the difference between the measurement of total C derived from the nonacid-treated soil and organic C derived from the acid-treated soil. Stable C isotopic ratios were presented in δ notation according to Eq. (1):

(1) δ = R sample - R STD R STD × 10 3 ,

where Rsample and RSTD are the 13C/12C ratio of the sample and standard, respectively. The 13C/12C ratios were expressed relative to the Vienna PeeDee Belemnite (V-PDB) standard (Coplen, 1994). Precision of duplicate measurements was 0.1 ‰ for δ13C. Soil TP was measured using sulfuric acid digestion and colorimetry (Dick and Tabatabai, 1977) at the Soil Testing Lab at Kansas State University, Manhattan, KS.

The fraction of soil C derived from C4 grasses was calculated by the mass balance equation Eq. (2):

(2) δ 13 C soil = δ 13 C grass ( f ) + δ 13 C woody tissue ( 1 - f ) ,

where δ13Csoil, δ13Cgrass, and δ13Cwoody tissue are the δ13C values of soil, C4 grass, and leaf tissues from C3 woody plants (live oak or Ashe juniper), respectively. The fraction, (f), is the proportion of SOC derived from C4 plants and (1−f) is the proportion of SOC derived from C3 plants.

2.4 Calculation of soil C, N, and P densities

As woody encroachment often changes soil bulk density (Hudak et al., 2003; Throop et al., 2012; Throop and Archer, 2008), the quantification of SOC, TN, TP, and SIC densities in different soil layers was performed on the basis of equivalent soil mass (Fowler et al., 2023). In addition, we limited our statistical analyses of soil variables to the 0–10, 10–20, and 20–40 cm depth increments to facilitate comparisons due to differences in maximum soil depth between trenches. The estimation of SOC stocks (Mg C ha−1) was calculated considering SOC concentrations and soil masses in the 0–10, 10–20, and 20–40 cm soil layers along with the trench faces using a cubic spline function (Fuentes et al., 2010; von Haden et al., 2020; Wendt and Hauser, 2013). Although all SOC stocks were calculated and compared as equivalent soil masses, the results were expressed as the 0–10, 10–20, and 20–40 cm soil layers for better clarity. To facilitate comparisons between the unequal sampling depth intervals, we conducted our analyses on SOC density (kg C m−3) using the following expression:

(3) SOC density ( kg C m - 3 ) = SOC stock ( Mg C ha - 1 ) depth interval ( cm ) × 10 .

Total N, TP, and SIC densities were calculated and presented using the same process.

2.5 Data analyses

The effects of geology (Edwards or Buda formations), vegetation type (grassland, Ashe juniper, live oak, or mesquite), and soil depth on soil chemical and isotopic measurements were analyzed using a three-way analysis of variance (ANOVA) with repeated measures in a linear mixed-effects model with restricted maximum likelihood estimations. Trench ID was treated as a random factor. Sampling locations along the trench face and depth were treated as a random factor with repeat measurements. Post hoc multiple comparisons were made using the Tukey adjustment. Dependent variables were tested for the assumption of normality using the Shapiro–Wilk test and then underwent logarithmic or square-root transformations when necessary to meet the assumptions of normality before analysis. The significance of ANOVA outputs was evaluated using an α value of 0.05. All errors were reported as standard errors. All statistical analyses were performed in R (R Core Team, 2013), including the nlme package (Pinheiro et al., 2012) for ANOVA. To visualize soil properties in the horizontal and vertical dimensions along the trench faces, contour maps were developed utilizing all data from all soil depths. Contour maps were created using the local weighted regression surface (loess) function in R.

3 Results

3.1 Soil characteristics

Soil bulk densities were significantly affected by the depth × geology interaction (Table 2), increasing more rapidly with depth in Edwards vs. Buda soils (Fig. 2a). Soil bulk densities were significantly higher beneath grasslands than those under live oak or Ashe juniper. The extremely high abundance of gravel and rock fragments in Edwards soils under live oak and Ashe juniper canopies prevented the measurements of bulk density at depths >20 cm. Soil pH was significantly affected by both soil depth and vegetation type (Table 2), increasing with depth, and it was lower beneath oak than beneath juniper and grass throughout the 0–40 cm profile (Fig. 2b). Fine-root densities were affected significantly by the depth × vegetation interaction (Table 2), with higher densities in soils beneath oak and juniper canopies than in grasslands soils below 10 cm depth (Fig. 2c). Soil inorganic C density expressed, based on equivalent soil mass, was affected by a complex depth × vegetation × geology interaction (Table 2). This interaction was attributable to greater differences between vegetation types on Edwards soils than on Buda soils as well as to the fact that SIC at the 15 and 30 cm depths was substantially higher than at the 5 cm depths beneath grassland than beneath oak or juniper vegetation (Fig. 2d). Lower SIC values beneath oak and juniper on the Edwards soil were related to lower soil pH values beneath their canopies; similarly, lower soil SIC values beneath oak on the Buda soils were also related to lower soil pH.

Table 2Results of the ANOVA testing for the effects of soil depth, vegetation, geology, and their interactions on soil attributes. The abbreviations used in the table are as follows: SOC – soil organic carbon; TN – total nitrogen; TP – total phosphorous; C:NC:N ratio; C:PC:P ratio; N:PN:P ratio; BD – bulk density; SIC – soil inorganic carbon; and “Fine roots” – fine-root (<2 mm) density. Asterisks denote significant changes: p<0.05; p<0.01; and p<0.001. ns represents nonsignificant results.

Download Print Version | Download XLSX

Figure 2Vertical changes in (a) soil bulk density, (b) pH, (c) fine-root biomass density, and (d) soil inorganic carbon (SIC) density beneath grass, juniper, and oak in soils derived from the respective Buda and Edwards formations. The extremely rocky nature of the soil below 20 cm depth beneath oak and juniper canopies precluded measurements of bulk density and root density in the Edwards trenches. Results are given as the mean ± standard error of the mean. Data are plotted at the midpoints of the depth increments.


3.2δ13C values of soils and litter

Litter δ13C values were 26.7 ‰, 24.7 ‰, and 17.1 ‰ under live oak, Ashe juniper, and grass, respectively (Fig. 3a). The δ13C of SOC was significantly affected by vegetation and soil depth (Table 2). The δ13C values of SOC under grasslands were 19 ‰ and 20 ‰ in the 0–10 cm layers in Edwards and Buda soils, respectively (Fig. 3a). Soil δ13C values under woody patches were lower than grasslands on both Edwards and Buda soils, ranging from approximately 24 ‰ to 20 ‰. The δ13C of SOC increased with depth beneath all vegetation types and on both Buda and Edwards soils. The δ13C values gradually increased from the centers of woody patches to their edges, and they reached highest values in the grasslands (Figs. S1, S2).

Figure 3Change in (a) δ13C of litter (above dashed line) and soil as well as change in (b) the fraction of SOC derived from C4 grass beneath juniper and oak in soils derived from the respective Buda and Edwards formations calculated using mass balance. Results are given as the mean ± standard error of the mean. Data are plotted at the midpoints of the depth increments.


The average δ13C value for C4 grasses was 16.3 ‰ ± 0.8 ‰, while the averages for C3 woody plants were 28.5 ‰ ± 0.5 ‰ for live oak and 26.8 ‰ ± 0.4 ‰ for Ashe juniper. These values were applied to a mass balance equation (Eq. 2) to estimate the relative proportion of SOC originating from C4 grasses. The fraction of SOC from C4 grass was approximately 0.5 in the 0–10 cm depth increment and increased up to 0.7–0.8 in the 20–40 cm depth increment (Fig. 3b). Soils under oak had significantly less C derived from C4 grass than juniper on the Edwards soils.

3.3 Soil organic C, TN, TP, and their stoichiometric ratios

The SOC, TN, and TP densities (kg m−3 soil) on an equivalent soil mass basis) in the upper 40 cm of the profile were all significantly affected by the vegetation × geology interaction (Table 2). Soils beneath live oak and Ashe juniper had higher SOC and TN than grasslands throughout the profile on Edwards soils, while only oak had higher SOC and TN on Buda soils (Fig. 4a, c). The cumulative SOC in the 0–40 cm profile on Edwards soils was 23.1 and 21.0 kg C m−2 in soils under juniper and oak, respectively, compared with 15.5 kg C m−2 for the grassland (Fig. 4b). In contrast, cumulative SOC on Buda soils was significantly lower, with 14.2, 10.4, and 10.9 kg C m−2 in soils under oak, juniper, and grassland, respectively. The vertical profiles of TN were similar to SOC (Fig. 4c). Cumulative TN (0–40 cm) was greater in Edwards soils (1.35–1.85 kg N m−2) than Buda soils (0.85–1.12 kg N m−2) and greater under oak vegetation than under juniper or grassland on Edwards soils (Fig. 4d). On Buda soils, oak was significantly higher with respect to cumulative TN (1.1 kg N m−2) than grassland (0.8 kg N m−2). Soils under oak had higher TP densities (kg P per m3 soil) than either grassland or juniper throughout the 0–40 cm profile on Edwards soils (Fig. 4e). The cumulative TP in the 0–40 cm profile on Edwards soils was highest beneath oak (0.254 kg P m−2), followed by grassland (0.210 kg P m−2) and juniper (0.165 kg P m−2) (Fig. 4f). Buda soils were generally lower in TP than Edwards soils, and there was no significant difference in cumulative TP between the three vegetation types in Buda soils.

Figure 4Vertical changes in (a) soil organic carbon (SOC), (c) total nitrogen (TN), and (e) total phosphorus (TP) density (kg m−3) beneath grass, juniper, and oak in soils derived from the respective Buda and Edwards formations. Data are plotted at the midpoints of the depth increments. Cumulative SOC (b), TN (d), and TP (f) values in the full 0–40 cm profile are shown, with significant differences indicated by asterisks: p<0.05, p<0.01, and p<0.001. Results are given as the mean ± standard error. Different letters represent significant differences between vegetation types (p<0.05) within specific geological formations.


Soil organic C concentrations (g C per kg soil) were affected by the depth × vegetation × geology interaction, while TN concentrations were affected by geology and the depth × vegetation interactions (Table S1 in the Supplement). For both SOC and TN concentrations, soils under grasslands were the lowest among the three vegetation types on both Buda and Edwards soils (Fig. S4 in the Supplement). In addition, SOC and TN concentrations were significantly higher on Edwards than on Buda soils. The soil TP concentration was affected by the vegetation × geology interaction (Table S1), with oak on Edwards soils having the highest TP compared with other vegetation types on either Edwards or Buda soils (Fig. S4).

Soil C:N ratios were significantly affected by the vegetation × geology interaction (Table 2), with higher ratios beneath oak and juniper than under grassland on Edwards soils throughout the 0–40 cm soil profile (Fig. 5a). However, on Buda soils, the C:N ratio was greater under oak and grassland compared with beneath juniper throughout the profile. Soil C:P ratios were significantly affected by the vegetation × geology × depth interaction due to (i) larger differences between grasses vs. woody plants in the Edwards than in the Buda soils, (ii) lower values on the Buda soils than the Edwards soils, and (iii) differences in the depth responses between the plant species (Fig. 5b). Soil N:P ratios were significantly affected by the direct effects of vegetation, geology, and depth but were not influenced by any interactions (Table 2). The N:P ratios were generally higher beneath woody plants compared with under grasslands, higher in the Edwards soils compared with the Buda soils, and decreased with soil depth.

Figure 5Changes in the (a) soil C:N ratio, (b) soil C:P ratio, and (c) soil N:P ratio with soil depth beneath grass, juniper, and oak canopies on the respective Buda and Edwards formations. Results are given as the mean ± standard error of the mean. Data are plotted at the midpoints of the depth increments.


4 Discussion

4.1 Evidence of vegetation change

The δ13C values of surface soils (0–10 cm) beneath juniper (21 ‰) and oak (23 ‰) canopies are higher than those obtained for litter (juniper = 25 ‰; oak = 27 ‰) and fresh leaf tissue (juniper = 27 ‰; oak = 29 ‰) obtained for these C3 woody species. This isotopic disequilibrium between soils and current organic matter inputs is even more marked in the 10–20 and 20–40 cm depth increments. This suggests that a considerable proportion of SOC beneath juniper and oak stands was derived from C4 grasses and that this landscape was once a more open C4-dominated grassland. Mass balance calculations based on plant and soil δ13C values (Eq. 2) indicate that approximately 45 %–55 % of SOC at 0–10 cm and 70 %–90 % of SOC at 20–40 cm is derived from C4 grass in soils atop both Buda and Edwards limestone (Fig. 3b). It should be noted that the absolute accuracy of these mass balance calculations could be affected by historical changes in the δ13C values of atmospheric CO2 (Graven et al., 2017), isotopic fractionation during soil organic matter decay (Mainka et al., 2022), and differential decomposition of C3 vs. C4 plant tissues (Wynn et al., 2020). However, these uncertainties are not large enough to alter the conclusion that C3 juniper and oak trees are more recent components of this landscape compared with C4 grasses. Our findings that woody plants are relatively recent components of this landscape are consistent with other studies on the Edwards Plateau and other portions of the southern Great Plains region based on direct observations of plant community changes through time (Briggs et al., 2002; Smeins and Fuhlendorf, 1997; Van Auken, 2008), δ13C values of SOC (Jessup et al., 2003), and sequential remote-sensing imagery (Wang et al., 2018a). A notable increase in woody cover has been observed in Edwards Plateau. From 1986 to 2020, woody cover in these regions increased from 7.6 % to 19.3 %, corresponding to an average annual increase of 0.8 % (Leite et al., 2023). This rate is comparable to observations in South Texas and Oklahoma, where annual increases in juniper cover averaged 0.5 %–1.5 % yr−1 (Barger et al., 2011; Archer et al., 2001; Fowler and Simmons, 2009; Wang et al., 2018a).

4.2 Soil C, N, and P concentrations and densities

Vegetation cover, geology, and soil depth have interacted to significantly alter the horizontal and vertical distribution of soil C, N, and P concentrations and densities within the soil profile (Tables 2, S1; Figs. 4, S4). The most consistently significant factor influencing soil C, N, and P concentrations and densities was the interaction between vegetation type and geology. This interaction was due to larger differences in soil C, N, and P between vegetation types on the Edwards formation compared with the Buda formation. Soil organic C, TN, and TP concentrations and densities were generally higher under oak and juniper canopies than under grasslands within all depth increments (Figs. 4, S4). In addition, the cumulative amounts of SOC, TN, and TP present throughout the entire 0–40 cm depth were always greater in soils atop the Edwards formation compared with those atop the Buda limestone for all three vegetation types (Fig. 4), with the exception that TP under juniper was lower in soil atop the Edwards formation. We speculate that the higher concentrations of C, N, and P in soils atop the Edwards limestone could be attributed to two factors: (1) the higher clay concentration which offers a higher C storage capacity (Six et al., 2002; Basile-Doelsch et al., 2020) and (2) the shallow depth to bedrock (approximately 40 cm) in Edwards soil that constrains root and litter inputs to a limited soil volume. In contrast, depth to bedrock on the Buda formation can be >1 m, providing a significantly greater soil volume in which root and litter inputs can be distributed, thereby reducing C, N, and P concentrations and densities. However, it should be noted that Ashe juniper roots have been shown to achieve rooting depths of 10 m and live oak roots can reach 20 m by exploiting cracks and fissures in the bedrock on the Edwards Plateau (Jackson et al., 1999).

Similar studies in other portions of the southern Great Plains and the western US have also shown that juniper and/or oak encroachment into grasslands results in higher soil nutrient storage (Fernandez et al., 2013; Jessup et al., 2003; McKinley and Blair, 2008; Shawver et al., 2018). Furthermore, our results are broadly consistent with prior studies around the world showing that tree/shrub encroachment into grasslands, savannas, deserts, and other arid/semiarid ecosystems generally results in increased concentrations and pool sizes of soil C, N, P, and other essential elements (Archer et al., 2017; Barger et al., 2011; Blaser et al., 2014; Eldridge et al., 2011; Zhou et al., 2017, 2018b, 2021).

The increases in soil C, N, and P stores under woody canopies following their invasion into previously grass-dominated areas is likely a consequence of multiple factors. First, rates of aboveground net primary productivity in areas encroached upon by woody plants can be as much as 400 % higher than those in adjacent open grasslands at sites across the Great Plains and western North America (Knapp et al., 2008). Although there are currently no comparative data on rates of primary production in open grasslands vs. regions experiencing juniper–oak encroachment in the Edwards Plateau, rates of gross primary productivity have been shown to be 55 % higher following juniper invasion compared with open grasslands across the state of Oklahoma (Wang et al., 2018b). Second, root biomass densities in the top 40 cm of the soil profile in our study area were approximately 1.7 times larger under juniper and oak patches than under grasslands on both the Edwards and Buda limestones (Fig. 2c). Finally, above- and belowground woody-plant tissues have been shown to be biochemically more recalcitrant than grass tissues due to their higher concentrations of suberin, cutin, and syringyl- and vanillyl-lignin subunits (Filley et al., 2008; Boutton et al., 2009). In addition, juniper tissues contain high concentrations of terpenoids (Adams et al., 2013) and oak tissues are enriched in tannins (Mole, 1993); these secondary compounds have the potential to limit the activity of decomposer organisms in the soil environment (Hättenschwiler and Vitousek, 2000). Collectively, the higher rates of above- and belowground organic matter inputs coupled with the biochemical recalcitrance of organic matter derived from woody plants likely favors soil organic matter accumulation and larger stores of soil C, N, and P beneath oak and juniper canopies.

Soil inorganic carbon densities ranged from approximately 10 to 50 kg C m−3 across all samples (Fig. 2d) and comprised, on average, 40 % of the soil total carbon (i.e., SIC + SOC) stores (Figs. 2d, 4a). These values are consistent with prior assessments of SIC in this region (Smith et al., 2014). The SIC densities were significantly influenced by the vegetation × geology × depth interaction, with differences between the vegetation types being more accentuated atop Edwards limestone and with values under grassland declining with depth at a more rapid rate than those under woody plants. The SIC densities were generally lower beneath oak and juniper compared with grassland areas (Fig. 2d). This is consistent with previous studies showing that woody-plant encroachment into grasslands leads to significant reductions in SIC via soil acidification (Liu et al., 2020). We suggest that the lower SIC densities beneath woody-plant canopies may be a consequence of larger pools of SOC and higher root densities under oak and juniper, which generate higher concentrations of soil CO2 via organic matter decay and root respiration. The increased infiltration rates due to a higher root density may export bicarbonate and further enhance carbonate weathering (Leite et al., 2023; Wen et al., 2021). Those higher soil CO2 concentrations would favor carbonic acid (H2CO3) production, which would react with calcium carbonate (CaCO3) to form bicarbonate (HCO3-) that can then dissociate in soil solution and release CO2 from the soil (Ramnarine et al., 2012; Wilsey et al., 2020; Hong and Chen, 2022). We acknowledge that factors such as topographic position and the presence of carbonate-rich horizons may also contribute to complex interactions between soil depth, vegetation, and geology in shaping SIC profiles on Edwards and Buda soils. Our interpretations of SIC profile were based on the assumption that the extent and timing of woody encroachment are consistent across soils atop the Edwards and Buda limestones in the Edwards Plateau. Given that SIC is among the largest pools in the global carbon cycle (940 Pg C; Eswaran et al., 2000) and that SIC is particularly abundant in arid and semiarid regions where woody-plant encroachment is prevalent, SIC loss to the atmosphere following encroachment may be an important impact on the carbon cycle (Liu et al., 2020) and should be considered when evaluating the biogeochemical consequences of land cover changes in these regions.

4.3 Soil stoichiometry

Ecosystem C, N, and P cycles are strongly coupled through the processes of primary productivity, respiration, and decomposition such that a change in the abundance of one element generally results in proportional changes in the other two elements. However, because the P cycle also has a significant inorganic geochemical component, the behavior of P in the soil environment can become less tightly coupled to that of C and N, resulting in shifts in stoichiometric ratios that can have important consequences for ecosystem structure and function (Finzi et al., 2011). In this study, soil C:N ratios were affected by the vegetation × geology interaction (Fig. 5a, Table 2). On the Buda limestone, soils beneath grass and oak had higher C:N ratios (12–13) than those under juniper (10–11) throughout the profile. On the Edwards formation, soil C:N ratios were higher under juniper and oak (12–13) than under grasslands (11). These differences are likely a response to the C:N ratios of the dominant plant species. Prior studies at this site and elsewhere around the western US have shown that the C:N ratios of leaves and roots of woody plants (oak species = 29–34; juniper species = 43–52) are higher than those of grasses (22–23) (Berner and Law, 2016; Marshall, 1995; Pregitzer et al., 2002). Thus, increased soil C and N storage beneath woody plants may be a consequence not only of their biochemical composition (as discussed earlier) but also of their more recalcitrant elemental composition. Soil C:P and N:P ratios were influenced by vegetation cover, geology, and soil depth (Fig. 5b, c; Table 2). These ratios were generally lowest in grassland soils, indicating that woody encroachment increased SOC and TN relatively more than TP. Similar changes to soil C:N:P stoichiometry have been documented following woody encroachment in a subtropical savanna in southern Texas (Zhou et al., 2018a). Soil C:P and N:P ratios were also higher in soils atop the Edwards formation compared with those on the Buda limestone, likely due to higher root biomass densities in the Edwards soils (Fig. 2c) which elevate soil C and N inputs. It is unclear why soil TP does not increase proportionally along with soil SOC and TN where woody plants have encroached. However, ecosystems are generally P-limited and have relatively low P inputs from atmospheric deposition (Mahowald et al., 2008) and the weathering of parent material (Prietzel et al., 2022), and these constraints may limit the ability of plants to acquire P in the same relative proportion as they acquire C and N.

Recognizing that the altered soil nutrient dynamics have far-reaching ecological consequences, the impact of woody encroachment might extend beyond biogeochemical changes. For instance, the modification of microclimatic soil conditions and the suppression of herbaceous diversity under woody canopies could influence broader biodiversity (Archer et al., 2017). Furthermore, the disparity in soil nutrient availability between soils under canopies and nutrient-deprived interspaces (Figs. S1, S2) may escalate the risk of soil erosion in drylands (Puttock et al., 2014; Ravi and D'Odorico, 2009; Wilcox et al., 2022). These considerations, while beyond the primary scope of our current investigation, highlight the multifaceted impacts of woody-plant encroachment on arid and semiarid ecosystems.

5 Conclusions

This study investigated the impact of Juniperus ashei and Quercus virginiana encroachment on soil C, N, and P stoichiometry in mixed-grass prairies on the Edwards Plateau of central Texas, considering the influence of underlying geological variations between soils lying atop two different limestone parent materials – the Buda and Edwards formations. Stable C isotope ratios (δ13C) of soil organic matter revealed that 45 %–90 % of soil C in the 0–40 cm depth interval beneath juniper and oak stands was derived from C4 plants, confirming that these woody plants were recent components of the landscape. The vegetation and geology interaction significantly influenced soil C, N, and P levels, with higher values under juniper and oak canopies compared with grasslands as well as on soils derived from the Edwards formation, possibly due to higher clay content and limited soil volume due to the shallow depth to bedrock (approximately 40 cm). Conversely, the deeper Buda formation (>1 m) allowed more extensive root and litter distribution, resulting in lower element concentrations. Soil C:N, C:P, and N:P ratios were generally higher under woody-plant canopies compared with grasslands, indicating that woody encroachment increased SOC and TN relatively more than TP. While C and N consistently increased – likely because of their close linkage during primary production, respiration, and decomposition – P trends deviated, reflecting influences from geochemical processes. Our results are broadly consistent with prior studies around the world showing that tree/shrub encroachment into grasslands, savannas, deserts, and other arid/semiarid ecosystems generally results in increased concentrations and pool sizes of soil C, N, and P as well as changes in their stoichiometric relationships. Our study also suggests that the magnitude of these changes may be influenced by attributes of the geological formations that underlie the soils. Given the geographic extent of woody encroachment at the global scale, our results suggest that soil conservation practices might need to be tailored according to the underlying geology. Interactions between vegetation change and geology warrant consideration in future studies, and they could play a role in efforts aimed at improving the prediction and modeling of soil C, N, and P storage in grasslands, savannas, and other dryland ecosystems.

Data availability

The data and R scripts presented in this work are available upon reasonable request from Che-Jen Hsiao (


The supplement related to this article is available online at:

Author contributions

TWB developed the research concepts. PAML surveyed and excavated the soil trenches. CJH conducted the fieldwork. AH and CJH carried out soil C, N, and δ13C analyses in the Stable Isotopes for Biosphere Science lab at Texas A&M University. CJH carried the other lab work, statistical analysis, data interpretation, and writing with contributions from TWB, PAML, and AH.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


We thank Butch Taylor and Doug Tolleson for providing background information on land use history and logistical support at the Texas A&M AgriLife Sonora Research Station and Lifei Sun for assistance with fieldwork.

Financial support

This research has been supported by the Texas A&M University Sid Kyle Global Savanna Research Initiative (project no. 513375-41060) and by the USDA/NIFA Hatch project (project no. 1020427).

Review statement

This paper was edited by Jocelyn Lavallee and reviewed by two anonymous referees.


Adams, R. P., Muir, J. P., Taylor, C. A., and Whitney, T. R.: Differences in chemical composition between browsed and non-browsed Juniperus ashei Buch. Trees, Biochem. Syst. Ecol., 46, 73–78,, 2013. 

Anderson, D. W.: The effect of parent material and soil development on nutrient cycling in temperate ecosystems, Biogeochemistry, 5, 71–97,, 1988. 

Ansley, R. J. and Wiedemann, H.: Reversing the woodland steady state: Vegetation responses during restoration of Juniperus-dominated grasslands with chaining and fire, in: Western North American Juniperus Communities, edited by: Van Auken, O. W., Springer, New York, 272–292, ISBN 978-0-387-34002-9, 2008. 

Archer, S. R., Boutton, T. W., and Hibbard, K. A.: Trees in grasslands, in: Global Biogeochemical Cycles in the Climate System, edited by: Schulze, E.-D., Harrison, S., Heimann, M., Holland, E., Lloyd, J., Prentice, I., and Schimel, D., Academic Press, San Diego, CA, 115–137,, 2001. 

Archer, S. R., Andersen, E. M., Predick, K. I., Schwinning, S., Steidl, R. J., and Woods, S. R.: Woody plant encroachment: Causes and consequences, in: Rangeland Systems: Processes, Management and Challenges, edited by: Briske, D. D., Springer, New York, NY, 25–84,, 2017. 

Augusto, L., Achat, D. L., Jonard, M., Vidal, D., and Ringeval, B.: Soil parent material – A major driver of plant nutrient limitations in terrestrial ecosystems, Global Change Biol., 23, 3808–3824,, 2017. 

Barger, N. N., Archer, S. R., Campbell, J. L., Huang, C. Y., Morton, J. A., and Knapp, A. K.: Woody plant proliferation in North American drylands: A synthesis of impacts on ecosystem carbon balance, J. Geophys. Res., 116, 1–17,, 2011. 

Basile-Doelsch, I., Balesdent, J., and Pellerin, S.: Reviews and syntheses: The mechanisms underlying carbon storage in soil, Biogeosciences, 17, 5223–5242,, 2020. 

Beck, H. E., Zimmermann, N. E., McVicar, T. R., Vergopolan, N., Berg, A., and Wood, E. F.: Present and future Köppen-Geiger climate classification maps at 1 km resolution, Sci. Data, 5, 180214,, 2018. 

Bendevis, M. A., Owens, M. K., Heilman, J. L., and McInnes, K. J.: Carbon exchange and water loss from two evergreen trees in a semiarid woodland, Ecohydrology, 3, 107–115,, 2010. 

Berner, L. T. and Law, B. E.: Plant traits, productivity, biomass and soil properties from forest sites in the Pacific Northwest, 1999–2014, Sci. Data, 3, 160002,, 2016. 

Blaser, W. J., Shanungu, G. K., Edwards, P. J., and Olde Venterink, H.: Woody encroachment reduces nutrient limitation and promotes soil carbon sequestration, Ecol. Evol., 4, 1423–1438,, 2014. 

Boutton, T. W., Archer, S. R., Midwood, A. J., Zitzer, S. F., and Bol, R.: δ13C values of soil organic carbon and their use in documenting vegetation change in a subtropical savanna ecosystem, Geoderma, 82, 5–41,, 1998. 

Boutton, T. W., Liao, J. D., and Archer, S. R.: Belowground carbon storage and dynamics accompanying woody plant encroachment in a subtropical savanna, in: Soil Carbon Sequestration and the Greenhouse Effect, edited by: Lal, R. and Follett, R., Soil Science Society of America, Madison, WI, 181–205,, 2009. 

Briggs, J. M., Hoch, G. A., and Johnson, L. C.: Assessing the rate, mechanisms, and consequences of the conversion of tallgrass prairie to juniperus virginiana forest, Ecosystems, 5, 578–586,, 2002. 

Canadell, J., Jackson, R. B., Ehleringer, J. R., Mooney, H. A., Sala, O. E., and Schulze, E. D.: Maximum rooting depth of vegetation types at the global scale, Oecologia, 108, 583–595,, 1996. 

Cooke, M. J., Stern, L. A., Banner, J. L., and Mack, L. E.: Evidence for the silicate source of relict soils on the Edwards Plateau, central Texas, Quaternary Res., 67, 275–285,, 2007. 

Cooperative Climatological Data Summaries:, last access: 14 December 2022. 

Coplen, T. B.: Reporting of stable hydrogen, carbon, and oxygen isotopic abundances, Pure Appl. Chem., 66, 273–276,, 1994. 

Culley, J. L. B.: Density and compressibility, in: Soil sampling and methods of analysis, edited by: Carter, M. R., Lewis Publishers, Boca Raton, FL, 529–539, 1993. 

Deng, Y., Li, X., Shi, F., Hu, X., and Gillespie, T.: Woody plant encroachment enhanced global vegetation greening and ecosystem water-use efficiency, Global Ecol. Biogeogr., 30, 2337–2353,, 2021. 

Dick, W. A. and Tabatabai, M. A.: An alkaline oxidation method for determination of total phosphorus in soils, Soil Sci. Soc. Am. J., 41, 511–514,, 1977. 

Dietz, S., Herz, K., Gorzolka, K., Jandt, U., Bruelheide, H., and Scheel, D.: Root exudate composition of grass and forb species in natural grasslands, Sci. Rep.-UK, 10, 10691,, 2020. 

Doetterl, S., Stevens, A., Six, J., Merckx, R., Van Oost, K., Casanova Pinto, M., Casanova-Katny, A., Muñoz, C., Boudin, M., Zagal Venegas, E., and Boeckx, P.: Soil carbon storage controlled by interactions between geochemistry and climate, Nat. Geosci., 8, 780–783,, 2015. 

Eldridge, D. J., Bowker, M. A., Maestre, F. T., Roger, E., Reynolds, J. F., and Whitford, W. G.: Impacts of shrub encroachment on ecosystem structure and functioning: Towards a global synthesis, Ecol. Lett., 14, 709–722,, 2011. 

Eldridge, D. J., Delgado-Baquerizo, M., Travers, S. K., Val, J., and Oliver, I.: Do grazing intensity and herbivore type affect soil health? Insights from a semi-arid productivity gradient, J. Appl. Ecol., 54, 976–985,, 2017. 

Eswaran, H., Reich, P. F., Kimble, J. M., Beinroth, F. H., Padmanabhan, E., and Moncharoen, P.: Global carbon stocks, in: Global Climate Change and Pedogenic Carbonates, edited by: Lal, R., Kimble, J. M., Eswaran, H., and Steward, B. A., Lewis Publishers, Boca Raton, FL, 15–25, ISBN-10 1566704588, 2000. 

Farella, M. M., Breshears, D. D., and Gallery, R. E.: Predicting drivers of collective soil function with woody plant encroachment in complex landscapes, J. Geophys. Res.-Biogeo., 125, e2020JG005838,, 2020. 

Fernandez, D. P., Neff, J. C., Huang, C. Y., Asner, G. P., and Barger, N. N.: Twentieth century carbon stock changes related to Piñon-Juniper expansion into a black sagebrush community, Carbon Balance Manag., 8, 1–13,, 2013. 

Filley, T. R., Boutton, T. W., Liao, J. D., Jastrow, J. D., and Gamblin, D. E.: Chemical changes to nonaggregated particulate soil organic matter following grassland-to-woodland transition in a subtropical savanna, J. Geophys. Res., 113, G03009,, 2008. 

Finzi, A. C., Austin, A. T., Cleland, E. E., Frey, S. D., Houlton, B. Z., and Wallenstein, M. D.: Responses and feedbacks of coupled biogeochemical cycles to climate change: examples from terrestrial ecosystems, Front. Ecol. Environ., 9, 61–67,, 2011. 

Fowler, A. F., Basso, B., Millar, N., and Brinton, W. F.: A simple soil mass correction for a more accurate determination of soil carbon stock changes, Sci. Rep.-UK, 13, 2242,, 2023. 

Fowler, N. L. and Simmons, M. T.: Savanna dynamics in central Texas: Just succession?, Appl. Veg. Sci., 12, 23–31,, 2009. 

Fuentes, M., Govaerts, B., Hidalgo, C., Etchevers, J., González-Martín, I., Hernández-Hierro, J. M., Sayre, K. D., and Dendooven, L.: Organic carbon and stable 13C isotope in conservation agriculture and conventional systems, Soil Biol. Biochem., 42, 551–557,, 2010. 

Fuhlendorf, S. D., Smeins, F. E., and Grant, W. E.: Simulation of a fire-sensitive ecological threshold: A case study of Ashe juniper on the Edwards Plateau of Texas, United States, Ecol. Model., 90, 245–255,, 1996. 

Gabriel, W. J., Loomis, L. E., and Douglass II, James A. I.: Soil survey of Edwards and Real counties, Texas, Washington D. C., USDA NRCS, 2009. 

Graven, H., Allison, C. E., Etheridge, D. M., Hammer, S., Keeling, R. F., Levin, I., Meijer, H. A. J., Rubino, M., Tans, P. P., Trudinger, C. M., Vaughn, B. H., and White, J. W. C.: Compiled records of carbon isotopes in atmospheric CO2 for historical simulations in CMIP6, Geosci. Model Dev., 10, 4405–4417,, 2017. 

Hättenschwiler, S. and Vitousek, P. M.: The role of polyphenols in terrestrial ecosystem nutrient cycling, Trends Ecol. Evol., 15, 238–243,, 2000. 

Hong, S. and Chen, A.: Contrasting responses of soil inorganic carbon to afforestation in acidic versus alkaline soils, Global Biogeochem. Cy., 36, e2021GB007038,, 2022. 

Hudak, A. T., Wessman, C. A., and Seastedt, T. R.: Woody overstorey effects on soil carbon and nitrogen pools in South African savanna, Austral. Ecol., 28, 173–181,, 2003. 

Jackson, R. B., Moore, L. A., Hoffmann, W. A., Pockman, W. T., and Linder, C. R.: Ecosystem rooting depth determined with caves and DNA, P. Natl. Acad. Sci. USA, 96, 11387–11392,, 1999. 

Jessup, K. E., Barnes, P. W., and Boutton, T. W.: Vegetation dynamics in a Quercus-Juniperus savanna: An isotopic assessment, J. Veg. Sci., 14, 841–852,, 2003. 

Knapp, A. K., Briggs, J. M., Collins, S. L., Archer, S. R., Bret-Harte, M. S., Ewers, B. E., Peters, D. P., Young, D. R., Shaver, G. R., Pendall, E., and Cleary, M. B.: Shrub encroachment in North American grasslands: shifts in growth form dominance rapidly alters control of ecosystem carbon inputs, Global Change Biol., 14, 615–623,, 2008. 

Leis, S. A., Blocksome, C. E., Twidwell, D., Fuhlendorf, S. D., Briggs, J. M., and Sanders, L. D.: Juniper invasions in grasslands: Research needs and intervention strategies, Rangelands, 39, 64–72,, 2017. 

Leite, P. A. M., Wilcox, B. P., and McInnes, K. J.: Woody plant encroachment enhances soil infiltrability of a semiarid karst savanna, Environ. Res. Commun., 2, 115005,, 2020. 

Leite, P. A. M., Schmidt, L. M., Rempe, D. M., Olariu, H. G., Walker, J. W., McInnes, K. J., and Wilcox, B. P.: Woody plant encroachment modifies carbonate bedrock: Field evidence for enhanced weathering and permeability, Sci. Rep.-UK, 13, 15431,, 2023. 

Liu, S., Zhou, L., Li, H., Zhao, X., Yang, Y., Zhu, Y., Hu, H., Chen, L., Zhang, P., Shen, H., and Fang, J.: Shrub encroachment decreases soil inorganic carbon stocks in Mongolian grasslands, J. Ecol., 108, 678–686,, 2020. 

Mahowald, N., Jickells, T. D., Baker, A. R., Artaxo, P., Benitez-Nelson, C. R., Bergametti, G., Bond, T. C., Chen, Y., Cohen, D. D., Herut, B., Kubilay, N., Losno, R., Luo, C., Maenhaut, W., McGee, K. A., Okin, G. S., Siefert, R. L., and Tsukuda, S.: Global distribution of atmospheric phosphorus sources, concentrations and deposition rates, and anthropogenic impacts, Global Biogeochem. Cy., 22, GB4026,, 2008. 

Mainka, M., Summerauer, L., Wasner, D., Garland, G., Griepentrog, M., Berhe, A. A., and Doetterl, S.: Soil geochemistry as a driver of soil organic matter composition: insights from a soil chronosequence, Biogeosciences, 19, 1675–1689,, 2022. 

Marshall, S. B.: Biogeochemical consequences of livestock grazing in a juniper-oak savanna, Texas A&M University, (last access: 29 January 2024), 1995. 

McKinley, D. C. and Blair, J. M.: Woody plant encroachment by Juniperus virginiana in a mesic native grassland promotes rapid carbon and nitrogen accrual, Ecosystems, 11, 454–468,, 2008. 

Meneguzzo, D. M. and Liknes, G. C.: Status and trends of Eastern Redcedar (Juniperus virginiana) in the central United States: Analyses and observations based on forest inventory and analysis data, J. Forest., 113, 325–334,, 2015. 

Mole, S.: The systematic distribution of tannins in the leaves of angiosperms: A tool for ecological studies, Biochem. Syst. Ecol., 21, 833–846,, 1993. 

Nieuwenhuize, J., Maas, Y. E. M., and Middelburg, J. J.: Rapid analysis of organic carbon and nitrogen in particulate materials, Mar. Chem., 45, 217–224,, 1994. 

Official Soil Series Descriptions (Eckrant series):, last access: 16 September 2022. 

Official Soil Series Descriptions (Prade series):, last access: 16 September 2022. 

Official Soil Series Descriptions (Rio Diablo series):, last access: 16 September 2022. 

Official Soil Series Descriptions (Valera series):, last access: 16 September 2022. 

O'Keefe, K., Bachle, S., Keen, R., Tooley, E. G., and Nippert, J. B.: Root traits reveal safety and efficiency differences in grasses and shrubs exposed to different fire regimes, Funct. Ecol., 36, 1–12,, 2021. 

Pinheiro, J., Bates, D., DebRoy, S., and Sarkar, D.: Nlme: Linear and nonlinear mixed effects models, (last access: 29 January 2024), 2012. 

Pregitzer, K. S., DeForest, J. L., Burton, A. J., Allen, M. F., Ruess, R. W., and Hendrick, R. L.: Fine Root Architecture of Nine North American Trees, Ecol. Monogr., 72, 293,, 2002. 

Prietzel, J., Krüger, J., Kaiser, K., Amelung, W., Bauke, S. L., Dippold, M. A., Kandeler, E., Klysubun, W., Lewandowski, H., Löppmann, S., Luster, J., Marhan, S., Puhlmann, H., Schmitt, M., Siegenthaler, M. B., Siemens, J., Spielvogel, S., Willbold, S., Wolff, J., and Lang, F.: Soil phosphorus status and P nutrition strategies of European beech forests on carbonate compared to silicate parent material, Biogeochemistry, 158, 39–72,, 2022. 

Puttock, A., Dungait, J. A. J., Macleod, C. J. A., Bol, R., and Brazier, R. E.: Woody plant encroachment into grasslands leads to accelerated erosion of previously stable organic carbon from dryland soils, J. Geophys. Res.-Biogeo., 119, 2345–2357,, 2014. 

R Core Team: R: A language and environment for statistical computing, (last access: 29 January 2024), 2013. 

Rabenhorst, M. C. and Wilding, L. P.: Pedogenesis on the Edwards Plateau, Texas: I. Nature and continuity of parent material, Soil Sci. Soc. Am. J., 50, 678–687,, 1986a. 

Rabenhorst, M. C. and Wilding, L. P.: Pedogenesis on the Edwards Plateau, Texas: III. New model for the formation of petrocalcic horizons, Soil Sci. Soc. Am. J., 50, 693–699,, 1986b. 

Ramnarine, R., Wagner-Riddle, C., Dunfield, K. E., and Voroney, R. P.: Contributions of carbonates to soil CO2 emissions, Can. J. Soil Sci., 92, 599–607,, 2012. 

Ravi, S. and D'Odorico, P.: Post-fire resource redistribution and fertility island dynamics in shrub encroached desert grasslands: A modeling approach, Landscape Ecol., 24, 325–335,, 2009. 

Rosan, T. M., Aragão, L. E. O. C., Oliveras, I., Phillips, O. L., Malhi, Y., Gloor, E., and Wagner, F. H.: Extensive twenty-first century woody encroachment in South America's savanna, Geophys. Res. Lett., 46, 2019GL082327,, 2019. 

Sala, O. E. and Maestre, F. T.: Grass-woodland transitions: Determinants and consequences for ecosystem functioning and provisioning of services, J. Ecol., 102, 1357–1362,, 2014. 

Shawver, S. E., Brady, J. A., and McGahan, D. G.: Variation of prokaryotic and fungal soil communities across a vegetative transect, Agric. Res. Technol. Open Access J., 15, 1–8,, 2018. 

Six, J., Conant, R. T., Paul, E. A., and Paustian, K.: Stabilization mechanisms of soil organic matter: Implications for C-saturation of soils, Plant Soil, 241, 155–176,, 2002. 

Smeins, F. E. and Fuhlendorf, S. D.: Biology and ecology of ashe juniper, in: Juniper Symposium Proceedings, edited by: Taylor Jr., C. A., Report 94-2, Texas A&M Research and Extension Center, San Angelo, TX, 9–24, 1997. 

Smeins, F. E. and Merrill, L. B.: Long-term change in semi-arid grassland, in: Edwards Plateau Vegetation: Plant Ecological Studies in Central Texas, edited by: Amos, B. and Gehlbach, F., Baylor University Press, Waco, TX, 101–114, 1988. 

Smeins, F. E., Taylor, T. W., and Merrill, L. B.: Vegetation of a 25-Year Exclosure on the Edwards Plateau, Texas, J. Range Manage., 29, 24,, 1976. 

Smith, D. B., Cannon, W. F., Woodruff, L. G., Solano, F., and Ellefsen, K. J.: Geochemical and mineralogical maps for soils of the conterminous United States, U.S. Geological Survey, 386 pp.,, 2014. 

Soil Survey Staff: Keys to soil taxonomy, 12th edn., USDA-Natural Resources Conservation Service, Washington, DC,, 2014. 

Soong, J. L., Janssens, I. A., Grau, O., Margalef, O., Stahl, C., Van Langenhove, L., Urbina, I., Chave, J., Dourdain, A., Ferry, B., Freycon, V., Herault, B., Sardans, J., Peñuelas, J., and Verbruggen, E.: Soil properties explain tree growth and mortality, but not biomass, across phosphorus-depleted tropical forests, Sci. Rep.-UK, 10, 1–13,, 2020. 

Stevens, N., Lehmann, C. E. R., Murphy, B. P., and Durigan, G.: Savanna woody encroachment is widespread across three continents, Global Change Biol., 23, 235–244,, 2017. 

Taylor, C. A., Twidwell, D., Garza, N. E., Rosser, C., Hoffman, J. K., and Brooks, T. D.: Long-term effects of fire, livestock herbivory removal, and weather variability in Texas semiarid savanna, Rangeland Ecol. Manag., 65, 21–30,, 2012. 

Throop, H. L. and Archer, S. R.: Shrub (Prosopis velutina) encroachment in a semidesert grassland: Spatial-temporal changes in soil organic carbon and nitrogen pools, Global Change Biol., 14, 2420–2431,, 2008. 

Throop, H. L., Archer, S. R., Monger, H. C., and Waltman, S.: When bulk density methods matter: Implications for estimating soil organic carbon pools in rocky soils, J. Arid Environ., 77, 66–71,, 2012. 

Torn, M. S., Trumbore, S. E., Chadwick, O. A., Vitousek, P. M., and Hendricks, D. M.: Mineral control of soil organic carbon storage and turnover content were measured by horizon down to the depth at which, Nature, 389, 3601–3603, 1997. 

Van Auken, O. W.: Western North American Juniperus communities: A dynamic vegetation type, edited by: Caldwell, M. M., Heldmaier, G., Jackson, R. B., Lange, O. L., Mooney, H. A., Schulze, E.-D., and Sommer, U., Springer, New York, NY,, 2008. 

Van Auken, O. W. and Smeins, F. E.: Western North American Juniperus communities: Patterns and causes of distribution and abundance, in: Western North American Juniperus Communities: A Dynamic Vegetation Type, edited by: Van Auken, O. W., Springer, New York, 3–18,, 2008. 

Venter, Z. S., Cramer, M. D., and Hawkins, H. J.: Drivers of woody plant encroachment over Africa, Nat. Commun., 9, 1–7,, 2018. 

von Haden, A. C., Yang, W. H., and DeLucia, E. H.: Soils' dirty little secret: Depth-based comparisons can be inadequate for quantifying changes in soil organic carbon and other mineral soil properties, Global Change Biol., 26, 3759–3770,, 2020. 

Wang, J., Xiao, X., Qin, Y., Doughty, R. B., Dong, J., and Zou, Z.: Characterizing the encroachment of juniper forests into sub-humid and semi-arid prairies from 1984 to 2010 using PALSAR and Landsat data, Remote Sens. Environ., 205, 166–179,, 2018a. 

Wang, J., Xiao, X., Zhang, Y., Qin, Y., Doughty, R. B., Wu, X., Bajgain, R., and Du, L.: Enhanced gross primary production and evapotranspiration in juniper-encroached grasslands, Global Change Biol., 24, 5655–5667,, 2018b. 

Web Soil Survey:, last access: 16 September 2022. 

Wen, H., Sullivan, P. L., Macpherson, G. L., Billings, S. A., and Li, L.: Deepening roots can enhance carbonate weathering by amplifying CO2-rich recharge, Biogeosciences, 18, 55–75,, 2021. 

Wendt, J. W. and Hauser, S.: An equivalent soil mass procedure for monitoring soil organic carbon in multiple soil layers, Eur. J. Soil Sci., 64, 58–65,, 2013. 

Wiedenfeld, C. C. and McAndrew, J. D.: Soil survey of Sutton County, Texas, Washington D. C., USDA-NRCS (Natural Resources Conservation Service, United States Department of Agriculture), 1968. 

Wilcox, B. P., Wilding, L. P., and Woodruff, J. M.: Soil and topographic controls on runoff generation from stepped landforms in the Edwards Plateau of Central Texas, Geophys. Res. Lett., 34, 1–6,, 2007. 

Wilcox, B. P., Basant, S., Olariu, H., and Leite, P. A. M.: Ecohydrological connectivity: A unifying framework for understanding how woody plant encroachment alters the water cycle in drylands, Front. Environ. Sci., 10, 934535,, 2022. 

Wilsey, B., Xu, X., Polley, H. W., Hofmockel, K., and Hall, S. J.: Lower soil carbon stocks in exotic vs. native grasslands are driven by carbonate losses, Ecology, 101, e3039,, 2020. 

Wynn, J. G., Duvert, C., Bird, M. I., Munksgaard, N. C., Setterfield, S. A., and Hutley, L. B.: Land transformation in tropical savannas preferentially decomposes newly added biomass, whether C3 or C4 derived, Ecol. Appl., 30, e02192,, 2020.  

Zhang, Q., Boutton, T. W., Hsiao, C.-J., Mushinski, R. M., Wang, L., Bol, R., and Klumpp, E.: Soil colloidal particles in a subtropical savanna: Biogeochemical significance and influence of anthropogenic disturbances, Geoderma, 430, 116282,, 2023. 

Zhou, Y., Boutton, T. W., and Wu, X. Ben: Soil carbon response to woody plant encroachment: importance of spatial heterogeneity and deep soil storage, J. Ecol., 105, 1738–1749,, 2017. 

Zhou, Y., Boutton, T. W., and Wu, X. Ben: Soil C:N:P stoichiometry responds to vegetation change from grassland to woodland, Biogeochemistry, 140, 341–357,, 2018a. 

Zhou, Y., Boutton, T. W., and Wu, X. Ben: Soil phosphorus does not keep pace with soil carbon and nitrogen accumulation following woody encroachment, Global Change Biol., 24, 1992–2007,, 2018b. 

Zhou, Y., Boutton, T. W., and Wu, X. Ben: A three-dimensional assessment of soil δ13C in a subtropical savanna: Implications for vegetation change and soil carbon dynamics, Soil Syst., 3, 73,, 2019. 

Zhou, Y., Taylor, R. J., and Boutton, T. W.: Divergent Patterns and Spatial Heterogeneity of Soil Nutrients in a Complex and Dynamic Savanna Landscape, J. Geophys. Res.-Biogeo., 126, e2021JG006575,, 2021. 

Short summary
Tree cover has increased in grasslands worldwide, with juniper and oak trees expanding in the southern Great Plains, USA. Here, we examine how these changes interact with geology to affect soil C, N, and P storage. Soil concentrations of these elements were significantly higher under trees than grasslands but increased more under trees growing on Edwards soils. Our results suggest that geology and vegetation change should be considered when predicting soil storage in dryland ecosystems globally.