Soil parameters affecting longleaf pine (Pinus palustris) site quality in east Texas

Many ecosystems have been degraded through exploitation of their natural resources, or land-use conversion to agricultural and urban use^[1], and restoration is often challenging due to modifications of soils, introduction of exotic invasive species, and lack of adequate resources to adequately conduct the restoration. Site selection is an important step in ecosystem restoration because the original ecosystems may have been greatly altered due to human activities^[2]. The longleaf pine (Pinus palustris Mill.) ecosystems of the southeastern United States are no exception to this degradation. Prior to European settlement, longleaf pine ecosystems occupied vast areas of the southern Atlantic and Gulf Coastal Plain regions of the United States, with approximately 30 million hectares extending between east Texas to Virginia, and stretching as far south as Florida, covering several climatic, physiographic, and many soil types^[3−6].

Longleaf pine was found in a wide range of ecosystems and sites from excessively drained sandhills to poorly drained flatwoods^[7−11]. Longleaf pine was most competitive on the sandier, well-drained sites across the region; however, a relatively frequent low intensity fire return interval, every two to eight years, set by native peoples or from lightning, is regarded as a key factor in historically reducing hardwood and shrub encroachment on most sites where it was found^[12−15]. During the logging and naval stores industry boom of the 1920s, old growth longleaf pine was quickly reduced^[16]; by the mid-1930s, only 10% of the old growth longleaf pine forest remained in east Texas and west Louisiana, but was mostly secondary growth^[17]. This dramatic decline led conservation groups and government agencies to begin conserving the remaining longleaf stands, and also to initiate longleaf pine ecosystem restoration on sites where the ecosystem once existed. However, many challenges exist that hinder this process.

One cause of longleaf pine ecosystem restoration failure is the inadequate consideration of soil suitability for longleaf pine. Soil type can affect the vegetation present, while vegetation can affect the condition of the soils^[18]. Due to these challenges, restoration efforts hypothetically should focus in areas that fit site specific soil/site parameters that support longleaf pine ecosystem restoration with the least management inputs. The objective of this study was to evaluate select soil properties on three soil mapping units (series) currently supporting longleaf pine stands in east Texas and relate these properties to longleaf pine site index.

This study was conducted in the Western Gulf Region of the native longleaf pine range in eastern Texas within portions of the Angelina National Forest (31°2'52.3" N, 94°21'48.96" W) and Sabine National Forest (31°10'56.21" N, 93°43'34.68" W), United States Forest Service Forests and Grasslands of Texas (USA). Sites contained longleaf pine ecosystems before and after the logging boom, and are considered optimum reference sites for possible longleaf pine restoration. All are located on the Catahoula geologic formation, stretching from east Texas to the Mississippi River, that consisted of sandstone, ranging from a few meters to approximately 18 m thick^[19]. As with most of these national forests, recurring prescribed fires on a 3−5 year interval have been used to maintain the site conditions and reduce fuel loading. All sites contained an overstory of longleaf pine, with minimal mid-story or longleaf pine advanced regeneration, and with a variety of herbaceous species and woody plants dominating the understory. The climate for east Texas is described as humid and subtropical, with mild winters with mean low temperatures in January between 2.8 °C to 3.9 °C, with summer temperatures reaching 33.3 °C to 34.4 °C for mean highs in August. Annual rainfall ranges from 1,240 to 1,510 mm with a relatively long growing season^[20].

Sampling locations were located on three different soil series mapping units, exhibiting different soil characteristics on well-drained or excessively-drained soils and ranged in depth and texture to the argillic B horizon: Letney Series (loamy, siliceous, semi-active, thennic Arenic Paleudults), Tehran Series (loamy, siliceous, semiactive, thermic Grossarenic Paleudults), and the Stringtown Series (fine-loamy, siliceous, semiactive, thennic Typic Hapludults).

Plot establishment

Ten, 50 m radius plots were established within each soil series across the two national forests, for a total of 30 plots. Prior to selection, each potential plot was randomly located on relatively pure soil map units determined by soil profile assessments at five points, one point in the center and four in each cardinal direction, 50 m from the center point; verification and identification of the soil series was accomplished using a bucket auger. Any of the points that failed to be consistent with the range of characteristics for the given soil map unit for the site were rejected.

Field sampling

A 10 Basal Area Factor prism was used to determine basal area at each plot center. Site index trees were chosen from the trees recorded with the prism by selecting the six closest trees to plot center that were either dominant or codominant and free of wounds. If six were not recorded by the prism, the nearest suitable trees still within the plot were measured. Annual growth rings from the six trees were quantified from a tree core extracted at DBH (Diameter at Breast Height) to determine age. A laser range finder was used to estimate total height to the nearest 3.05 cm. Total age and height were used in the site index curves developed for longleaf pine^[21]. Soil samples were taken using a bucket auger at plot center to correlate with the collected longleaf pine data. Soil samples were taken from the first three horizons (A, E, and the first argillic B) while individual horizon depths were measured to a depth of 150 cm.

Laboratory methods

Soil textural (sand, silt, and clay) analyses were conducted using the Bouyoucos method^[22] from the A, E, and B horizon samples. For coarse textured soils, 100 g of oven-dried soil was used, while 50 g was used for medium and fine textured soils; each sample was mixed with 100 ml of sodium hexametaphosphate, left for 12 h in deionized water, and then agitated for 15 min. Hydrometer readings were then made at 40 s and at 2 h to obtain total suspended solids. Samples were then poured into a series of sieves dividing the sample into the five sand particle sizes and clay plus silt^[23]: very coarse sand, coarse sand, medium sand, fine sand, and very fine sand with the range in sizes being 1−2, 0.5−1, 0.25−0.5, 0.10−0.25, 0.05−0.10, and < 0.05 mm, respectively. The samples were placed in a forced-draft drying oven at 105 °C until a constant weight was reached, then dry-sieved using a Ro-Tap® Shaker utilizing the same size classifications. Soil samples were dried and weighed prior to sieving and each sand fraction was weighed post sieving.

Bulk density was measured following standard procedures^[22] adjacent to plot center where the other soil samples were collected using a core sampler with 48.25 mm diameter rings, and samples oven-dried at 105 °C until constant weight was achieved and weighed prior and after drying. Field capacity and wilting coefficient were measured using a soil pressure plate apparatus and chambers. Field moist samples were soaked in water for 24 h prior to being placed under the pressure plates at both −31 and −1,500 kPa. Subsamples were weighed moist and then oven-dried at 105 °C to constant weight and then reweighed.

Standard lab methods using an ICP Thermoscientific Analyzer were performed to obtain phosphorus, potassium, calcium, magnesium, nitrogen, organic carbon, and ammonium at the Stephen F. Austin State University Plant, Soil, and Water Laboratory. To obtain pH, a one to two ratio of soil to water using 12.5 g of soil and 25 ml of deionized water method was determined using a pH probe. Electrical conductivity was taken following the completion of the pH using the same prepared sample using an E.C. meter.

Statistical analysis

Analysis of variance (ANOVA) using Proc GLM (General Linear Model) procedure in SAS was used to determine significant differences (p = 0.05). If differences were found among variables, Tukey's mean separation test was then used. Because large set of variables inherently have some correlations, principal component analysis (PCA) was used to summarize all of the variables into unrelated variables (PC1, PC2 ...), and important or significant PCAs were selected to perform regression. The number of principal components evaluated was determined by using randomization in PC-ORD. The top 10 composite variables from each significant PCA were selected and used in step wise regression to determine which variables most influenced longleaf pine site index.

The official descriptions for all three series were: are they are deep, well drained to excessively drained, with some variations in texture, color, and depth of each horizon^[24]. Depth to the first argillic horizon ranged from 23 to 49 cm (mean = 42.5 cm) in the Stringtown series, from 55 to 88 cm (mean = 67.1 cm) in the Letney series, and from 101 to 155 cm (mean = 111 cm) in the Tehran series. The greatest difference is depth to the first argillic (Btl) horizon: Stringtown < 50 cm, Letney 50 to 100 cm, and Tehran Bt1 > 100 cm.

Longleaf pine site indices

ANOVA indicated significant differences for longleaf pine site index (Table 1). Mean site indices for Letney and Stringtown soils were within the USDA Natural Resources Conservation Service (NRCS) range of site indices, but was below for Tehran soils (Table 2).

Soil physical parameters

Within the unweighted soil physical parameters, 12 were significantly different (Table 3). Both depth of A and depth to E on Tehran soils were significantly deeper than Stringtown. As expected, depth to B was significantly different, with Tehran being the deepest and Stringtown being the shallowest. Depth of E was also found to be significantly different, with Tehran being greater than both Stringtown and Letney; depth of B was also significantly different, with Stringtown approximately 73 cm thicker than Tehran, and 31 cm thicker than Letney. Wilting coefficient of the A horizon showed significant differences between Stringtown and Letney soils, with 50% more water held in the Stringtown series (Table 4). B horizon wilting coefficient was significantly greater in Stringtown than the Tehran soils.

Medium sand in the A and B horizons had the highest percent by weight in Tehran soils over the other soils. Medium sand in the B horizon and wilting coefficient of the B horizon were inversely correlated; as medium sand increased, wilting coefficient decreased. As the depth to the first argillic B horizon increased, both total silt + clay and total clay in the B horizon decreased.

Six physical variables weighted by horizon thickness were determined to be significantly different by soil series (Table 4). Field capacity in the E horizon was higher in Tehran soils than the others. Stringtown soils were significantly different from Letney and Tehran soils for field capacity and wilting coefficients weighted by thickness of the B horizon, and Stringtown soils held more moisture at field capacity and at wilting coefficient in the B horizon than Letney and Tehran. A and E horizon organic matter content was highest in Tehran. Organic matter content in the B horizon had the opposite trend, where Stringtown soils were significantly greater than Tehran soils.

Soil chemical parameters

Of the 36 soil chemical parameters not weighted by horizon thickness, exchangeable Ca in the A horizon was the only parameters found to be significantly different; Ca concentration in the A horizon in the Letney soils was significantly higher than in the other two soils.

Weighted by horizon thickness, 17 variables were significantly different (Table 5). Ca weighted by E horizon thickness was not significantly different, but were in the A and B horizons. Organic C in the A horizon was greater in Tehran than in Stringtown soils; and in the E horizon was greater than in Stringtown and Letney soils. The B horizon had the opposite effect, as Stringtown soils contained more organic C than Tehran. Overall, Stringtown contained more total N than Tehran soils, while in the E horizon Tehran soils had more total N; Stringtown had more total N in the B horizon than Letney soils, which had more than Tehran soils. Tehran had more NH₄ in the E horizon, but Stringtown had more in the B horizon than Tehran soils. Tehran had more P in the A and E horizons than Stringtown soils, and more K in the E horizon than Stringtown; Stringtown and Letney soils contained more K in the B horizon than Tehran soils. Stringtown soils contained more Mg in the B horizon than Tehran, and Stringtown soils contained more S in the B horizon than Tehran.

Generally, Stringtown had higher concentrations of nutrients in the B horizon than Tehran soils, although Tehran had higher concentrations in the A and E horizons. Within the A horizon, clay content was highest in the Letney soils which would provide a higher cation exchange capacity. K and Ca within the A horizon which were higher in Tehran and Letney soils; Stringtown averaged lower silt and clay in the A horizon resulting in lower quantities of those nutrients within the A horizon. Total N was highest in the A horizon in the Tehran which also contained the most organic C.

Soil profile nutrients were weighted by horizon depth and then summed for the entire 150 cm soil profile; Ca, Mg, and S were significantly different (Table 6). Stringtown and Letney soils contained more Ca than Tehran, and Stringtown soils contained more total Mg and S than Tehran. Soils with argillic B horizons closer to the surface (Stringtown and Letney) tended to have higher total available nutrient contents than Tehran. Total amounts of Ca, Mg, and S were found to be greatest in the Stringtown soils; Stringtown had the thickest B horizon relative to the 150 cm profile depth, and also had the highest amounts of silt and clay.

Principal Component Analysis ordination

Five variable combinations accounted for approximately 62% of the variation (Table 7) using principal component analysis. PC1 (21% of the variance) was strongly driven by depth to the B horizon, thickness of the B and E horizons, percent silt and clay in the B horizon, total wilting coefficient of the B horizon and entire profile, and total organic matter in the E horizon. PC2 (15% of the variance) was driven by percent medium sand, total sand and silt in the A horizon as well as percent medium sand, total sand, and silt in the E horizon. PC3 (10% of the variance) was driven by field capacity and available water capacity of the A and B horizons, total potential field capacity and available water capacity of the A horizon, total potential available water capacity of the B horizon, and total potential available water capacity for the profile. PC4 (7% of the variance) was driven by field capacity, wilting coefficient, and available water capacity of the E horizon and total field capacity of the entire profile, while PC5 (7% of the variance) was driven by percentage of very coarse sand, coarse and medium sand in the A horizon, percentage of very coarse sand and medium sand in the E horizon, and percentage of very coarse sand, coarse sand, and total clay in the B horizon.

Four significant PCA’s accounted for approximately 63% of the variation among the soil chemical variables (Table 7). PC1 (24% of variance) were concentrations of K, Ca, Mg, and boron in the B horizon, as well as total Mg and boron weighted by depth of the B horizon, and total K, Ca, Mg, and S weighted by depth of the 150 cm soil profiles. PC2 (17% of the variance) was driven by concentration of K, Ca, Mg, S, and Boron in the E horizon, total K, Ca, Mg, S, and B weighted by depth in the E horizon. PC3 (14% of the variance) was driven by total C, P, K, Ca, and Mg weighted by depth in the A horizon, as well as total grams of P weighted by depth of E horizon and total NH₄⁺ and total N weighted by depth of the B horizon. PC4 (8% of the variance) was driven by total C and P within the entire profile, and total N, P, and C in the B horizon.

Seven variables accounted for 67% of the variation (Table 7) when the physical and chemical variables were combined for analysis. PC1 (19% of the variance) was driven by depth to B, thickness of E and B, wilting coefficient of the B horizon, percentage of silt and clay in the B, total potential wilting point of the B, total wilting point in the profile, organic matter in the E horizon, Mg in the B, total N, K, Ca, Mg, S, and B in the B horizon, and total K, Ca, Mg, and S. PC2 (12% of the variance) was driven by field capacity and available water capacity in the A horizon, field capacity and total field capacity of the B horizon, total field capacity in the profile, P, K, Ca, and Mg in the A horizon, NH₄ in the E horizon, NH₄ in the B horizon, total P, K, Ca, Mg and S in the A horizon, and total NH₄ in the profile. PC3 (10% of the variance) was driven by Ca, Mg, and B in the E horizon, Mg, S, and B in the B horizon, total grams of P in the A horizon, and total Mg in the E, and total B in the B horizon. PC4 (9% of the variance) was driven by percent coarse sand, very fine sand, silt and clay, and total sand in the A horizon, percent coarse sand, fine sand, very fine sand, total sand, and total silt in the E horizon, total organic matter in the profile, P in the B horizon, and total C, P, and B in the profile. PC5 (6% of the variance) was driven by concentrations of K, Ca, Mg, S, and B of the E horizon and concentration of C in the B horizon. PC6 (6% of the variance) was driven by bulk density and organic matter in the E and B horizons, concentration of B in the A horizon, concentration of C in the E horizon, concentration of P and C in the B horizon, and total grams of NH₄, and B within the profile, while PC7 (5% of the variance) was driven by clay, wilting point, and total potential wilting point of the E horizon, concentration of Ca in the A horizon, and total Ca in the A horizon.

Regression for site index

Seven variables were the most significant soil physical factors affecting longleaf pine: depth to B, thickness of the E and B horizons, percent silt and clay in the B horizon, wilting coefficient of for the B horizon, wilting coefficient of the profile, and percent organic matter in the E horizon.

The best two-variable model (1) included depth to the B horizon and total wilting coefficient of the B horizon (R² of 0.3984):

Ten soil chemical variables correlated most with site index of longleaf pine: total K, Ca, Mg, and S in the profile, total Mg and B in the B horizon, and concentrations of K, Ca, Mg, S, and B in the B horizon. Only a one variable model (2) best fit the site index (R² of 0.2026):

Combining all variables, the variables most correlated to site index were depth to the B horizon, wilting coefficient of the B horizon, percent silt and clay and depth weight wilting coefficient of the B horizon, the profile weight wilting point of the whole profile, organic matter of the E horizon, concentration of Mg in the B horizon, total N, K, Ca, Mg, S, and boron in the B horizon, and total K, Ca, Mg, and S in the B horizon.

Using step-wise regression, the top variables that affect longleaf pine site index were total N and S in the B horizon, concentration of Mg in the B horizon, total Mg and S in the profile, and wilting coefficient weighted by horizon thickness in the B horizon. These six-variables proved to be the best model (R² = 0.6668). Regression Eqn (3) for site index using these six variables was:

The wilting coefficient was affected by the amount of silt and clay within the profile; as silt and clay decreased, so did the wilting coefficient. Within the unweighted soil physical parameters, field capacity and available water capacity should be affected by this texture correlation; however, it was not found in this study. In fact, available water capacity was highest in the deeper sand soils, suggesting that the pressure plate method used in this study may not have produced reasonable results.

The A horizon, as expected, contained more organic matter than either of the other two horizons (Table 3). The E horizon is characterized as where leaching of humus, silt and clay, and various ions occur, while the B horizon is where accumulation of humus, silt and clay, and various ions occur. In this study the B horizon did contain higher percentages of silt and clay than did the A and E horizons. Within all soils, as depth to the first argillic B horizon increased, the percentage of silt and clay decreased in the B horizon. Conversely, sand increased as depth to the first argillic horizon increased.

Similar to a previous study^[25], the wilting coefficient was influenced by the proportion of B horizon in the entire profile, which had higher wilting coefficients in all soils. Texture and B horizon thickness played a big role due to the inherent ability of fine textured soils to hold more water at the wilting coefficient. However, neither field capacity or available water capacity were statistically different between soils, likely due to the pressure plate system retaining more water than it should have. The B horizon total potential field capacity was significantly highest in the Stringtown soils, likely a result of the increase in silt and clay content in the thicker B horizon in those soils.

The higher percentage of clay may have reduced Ca leaching to the lower profile. The presence of finer texture soils increases cation exchange capacity (CEC), retaining cations. It is unclear why Stringtown had lower concentrations of Ca than Letney and Tehran, or why concentrations of other nutrients were not significantly different. The E horizon total potential field capacity was highest in the Tehran soils, again indicating that the pressure plate method did not produce reasonable results. Finer texture soils have higher CEC, which can result in the presence or ability to hold more cations^[26].

As depth to the B horizon increased, clay content decreased, as did water holding capacity, available water capacity, and wilting coefficient. Texture is inherently related to the amount of water a soil can hold. All three of the soil series in this study had A, E, and B horizons. In this study, texture did not prove to be as important as expected.

Many studies have looked at the relationship between the depth of horizons and site index, with varying results. While no correlation between site index and depth to the first argillic horizon B for longleaf pine was found in east Texas^[27], a negative correlation between site index and depth to finer textured layers in sandy soils for white oak (Quercus alba) was found in Michigan^[28], but the influence began to wane at depths greater than 1.5 m. Site index of radiata pine (Pinus radiata) increased with increasing depth of the topsoil^[29], but evaluated total soil depths in Douglas-fir did not show any significance with soil depths ranging from 50 cm to 100 cm^[30]. No single physical variable had a well-defined correlation with site index for shortleaf pine^[31]; however, depth to the B horizon and texture of that horizon was found to be a good indicator of shortleaf pine on the soils studied, is similar to what our study discovered. Higher soil organic matter resulted in increased site quality for most soils^{[25, 32]} and soil texture was not the only parameter affecting water holding capacity: other factors included OM and soil bulk density along with gravel content and salinity, which also affected water availability.

Ca, N, P, K, and Mg did not affect site index when added to radiata pine^[29]. This partially conflicts with our results, as Ca was the only nutrient that was significant. Growth is often limited in many forests in the southern United States by nutrient availability, promoting fertilization in silvicultural practices. Nitrogen and phosphorous are often considered the most common nutrients that limit growth in southern pine forests^[33]. Similarly, nutrients had a positive correlation with site index in radiata pine^[29], but they did not specify a given depth at which the nutrients were most effective. Our study found that nutrient levels in the B horizon had a strong correlation with site index increase in longleaf pine as did total amounts of nutrients in the profile, which were usually correlated to soils with higher silt and clay concentrations and shallower B horizons. Boron has been shown to positively correlate with growth in pine^[34]; our study also showed boron being positively correlated to longleaf pine site index, but did not differ between soils. Principal component analysis found depth to each horizon as well as certain nutrients within the B and E horizons as important. Higher nutrient levels correlated positively with site index^[29], but they did not specify if this relationship was by horizon or total in the profile. Our study showed low total nutrients within all three soils, and as depth to the first argillic B horizon increased, nutrients available decreased. Stringtown appeared to have greater amounts of nutrients within the B than the Tehran in most situations with Letney soils intermediate in nutrient availability.

Total N, Mg, and S in the B horizon had a positive effect on longleaf pine site index, while concentrations of Mg and S had a negative impact. Indirectly, depth to the B horizon had a correlation to these values as the thickness of the B impacts total available nutrients in that horizon. In addition, as B horizon depth increased, the percent clay within the horizon decreased. The presence of sand in the A and E horizon resulted in a lower cation exchange capacity, allowing for nutrients to leach through these horizons while the B horizon had an increase in clay.

While all of the three soils found on sites in this study historically supported longleaf pine, these results highlight how site index for this species might be driven by variables other than a sandy, well-drained A horizon. The importance of the B horizon depth and the lower ability of the B horizon to allow water to drain was an important variable found in the PCA analysis. Two of the soils in this study had site indices within the NRCS range for longleaf pine, but both were lower than the mean for those soils. The Tehran soil was slightly below the minimum site index for that soil. It could be on that soils, the NRCS underestimated the importance of the depth of the B horizon. In addition, longleaf pine historically on those soils may have benefited more from the short-interval fire frequency than on the other two soils.

Soil physical parameters in the A and E horizons did not appear to greatly influence site index for longleaf pine on these soils in east Texas. However, the depth to B and wilting coefficient of the B influenced site index of longleaf pine on these three soils, which suggest that water availability may play the largest role in affecting site index on these deep, coarse textured soils. Soil chemical parameters in the A and E horizons did not appear significant; however, soil chemical parameters in the first argillic B horizon, as well as nutrient availability in the whole profile did. Soil variables in the B horizon affect site index for longleaf pine the most, while some variables within the whole 150 cm profile also had an effect. This is likely due to the effect of the weighted by horizon thickness of the B horizon had on the total profile because of clay content of the horizon providing for higher available water content, and nutrient storage. Some A horizon parameters showed some slight effect on longleaf pine site index, but this could possibly be due to the amount of organic matter within the A horizon. Each model highlighting different variables reflects the complexity of the interaction of soil variables with site index. Productive forests tend to have soils with favorable physical properties that enhance biological functions. Separating and choosing the most significant of these soil variables can be challenging due to the inherent complexity and interactions among many of them.

Future studies should look at rooting depth within each of these three soils as well as the effect of soils with drainage classes that are known to hold more water. Studies should also consider planting on these three-soil series using the same treatments to how these soil variables affect longleaf pine regeneration.

For those making management decisions on locations with the greatest potential success for longleaf pine establishment, soils with similar A and B horizon characteristics may have the greatest success, along with the use of recurring prescribed fire.

The authors confirm contribution to the paper as follows: study design and manuscript preparation: Oswald BP, Farrish KW, Svehla R; data collection and data analysis: Svehla R. All authors reviewed the results and approved the final version of the manuscript.

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.