Soil phosphorus fractions vary with harvest intensity and vegetation control at two contrasting Douglas-fir sites in the Pacific northwest

Effects of intensive forest management on soil phosphorus (P) are unclear and may impact long-term site productivity. We assessed changes in P availability over 10 years associated with harvest intensity (bole-only vs. whole-tree harvest) and vegetation control treatments (initial vegetation control (IVC) vs. five years of annual vegetation control (AVC)) using a P fractionation procedure. Fractions were characterized at 0–15, 15–30, and 30–60 cm soil depths in two coast Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco var. menziesii) plantations with strongly contrasting soil properties near Matlock, WA (young soils formed in glacial outwash) and Molalla, OR (relatively old soils formed in igneous residuum and exhibiting andic properties). Al and Fe concentrations associated with short-range order minerals were greater at Molalla than Matlock and generally decreased with depth at both sites. We observed decreases in most total-P and P-fraction concentrations across the three soil depths at the Molalla site. Effects were less pronounced and generally inconsistent at the Matlock site. Decreases in total P and P fraction concentrations were greatest in the AVC treatments at Matlock, but opposite trends were observed at Molalla where decreases were greatest with IVC. There was no difference between harvest treatments on the change in P fractions in most instances, with the exception of the 30–60 cm depth at Matlock where concentrations of some P fractions were maintained or increased with bole-only harvesting. Ten-year responses indicate harvest intensity has limited effects on long-term productivity associated with soil P because of the large size of the soil P pools and the relatively small changes in soil P that occurred with treatment. Decreases in P concentrations with AVC at Matlock and IVC at Molalla were larger than the other treatments and highlight the important role of vegetation in P dynamics following harvesting at these sites.


Introduction
Intensive forest management is characterized by increased fiber removals at shorter rotations coupled with increased cultural inputs such as herbicide application for vegetation control. This type of management is very effective at increasing fiber production (Adams et al., 2005), but also increases the risk of reduced productivity over time because of impacts to soil properties that control resource supply for plant growth (Powers et al., 1990). Intensive management is generally associated with a net increase in nutrient removal from a site relative to less intensive practices, and it influences processes that control nutrient loss and retention immediately following harvest, such as leaching, erosion, and growth of competing vegetation (Vitousek and Matson, 1985;Smethurst and Nambiar, 1995). These direct and indirect effects on nutrient pools have the potential to reduce soil resources such as available phosphorus (P).
The USDA Forest Service Long-term Soil Productivity (LTSP) study was launched in 1989 to elucidate effects of intensive forest management on long-term soil productivity in North America (Powers et al., 2005). Research conducted within the North American LTSP study aimed to answer questions that arise from potential changes in soil productivity and concern over observed or hypothesized management effects on soil quality (Boyle and Ek, 1972;Powers et al., 1990). In the U.S. Pacific Northwest (PNW), LTSP research has focused mainly on soil organic carbon (C) (Powers et al., 1990) and nitrogen (N) limitation (Chappell et al., 1991). Soil P has been less studied despite https://doi.org/10.1016/j.geoderma.2019.04.038 Received 4 September 2018; Received in revised form 18 April 2019; Accepted 23 April 2019 acknowledgement of potential NeP co-limitation and other nutrient interactions (Spears et al., 2001;Vitousek et al., 2010). Additionally, there have been few studies in this region that measured changes in soil P over extended periods of time in response to intensive forest management.
Sequential P fractionation was first developed as a method of observing changes in soil P fractions caused by cultivation of field crops (Hedley et al., 1982). More recently it has been used in research to investigate the role of P biogeochemistry in forested terrestrial ecosystems (Giardina et al., 1995;Compton and Cole, 1998;Spears et al., 2001;Ilg et al., 2009;Celi et al., 2013), which differ from field crops in their longer rotation lengths, management regimes, and perennial plant communities. Fractionation has been used to separate labile, moderately-labile, and recalcitrant forms of P and to differentiate inorganic and organic forms within these classes. There have been many modifications (e.g., addition of different extraction steps, varying concentration of extractants, and use of physical extraction techniques such as sonic baths) of fractionation procedures with new technology and insight, but the concept of separating P fractions based on chemical susceptibility to different extracting solutions remains consistent (Robertson et al., 1999). Although these extractions provide only an approximation for theoretical pools, their usefulness in exploring questions pertaining to development and evolution of P pools over time is well-supported (Walker and Syers, 1976;Crews et al., 1995;Schlesinger et al., 1998).
In the PNW, soil P fractionation has been used to compare availableand other P fractions in association with different tree species and how composition changes available-P fractions in stands with intermixed species including Douglas-fir (Giardina et al., 1995;Compton and Cole, 1998;Spears et al., 2000). Vegetation in general has a large influence on P availability and cycling through conversion of inorganic P into less-available organic and occluded forms (Walker and Syers, 1976), and potentially greater influence when the vegetation is capable of fixing atmospheric N (Houlton et al., 2008). For example, Compton and Cole (1998) demonstrated that plant-available P fractions were generally increased in the presence of red alder (Alnus rubra Bong.), and Giardina et al. (1995) also showed that available inorganic P was increased in the presence of alder. Such an effect may be dependent on soil or other site factors (e.g. rainfall, microbial communities, and invasive plant growth) that influence the forms and magnitude of available soil P or total soil N (Holub et al., 2013;Slesak et al., 2016a).
The relationship between N-fixing plant species and available-P fractions becomes relevant in relation to vegetation control practices, which may alter the composition of plants that stimulate rhizosphere interactions and phosphatase activity (Caldwell, 2006). Plant uptake of available P may also change the equilibrium of labile P, with the nonoccluded organic fractions increasing phosphatase activity and decreasing the reserve fractions of P (Sharpley and Smith, 1985). Sequential P fractionation could provide insight on P biogeochemical processes as affected by harvesting and vegetation control treatments and how those treatments and related effects on vegetation uptake of P affect different P fractions. Inference generated by additional knowledge of P dynamics will yield greater understanding of these treatments and their effect on P pools.
The objectives of this study were to determine if harvest-intensity treatments and vegetation control treatments affected soil P fractions over a ten-year period since harvest and treatment, how the responses may differ between sites with distinctly different soils, and whether or not effects on soil P were associated with foliar P nutrition. We measured concentrations of five P fractions from two contrasting sites supporting Douglas-fir plantations in the PNW, analyzed foliage P from Douglas-fir trees on the same plots, and additionally conducted selective dissolution extractions to further characterize soils and provide context to how mineralogy influences responses to changes in phosphorus fractions, were performed. We expected that more intense management practices, especially when applied in tandem (e.g., greater harvest intensity and more prolonged vegetation control), would have greater reductions in P fractions compared to less intensive practices. The ten-year time span of observation in this study is aimed to yield more applicable assessments of treatment effects on long-term soil productivity in early rotation sites managed for forest production.

Site descriptions
The two sites used in this study, Matlock and Molalla, are affiliates of the North American LTSP network (Powers et al., 2005). These sites were established in 2003 to study effects of biomass removal and vegetation control on nutrient cycling, Douglas-fir growth, competing vegetation, and soil properties. The sites have a climate with cool, wet winters and warm, dry summers. Mean annual precipitation is 240 cm at Matlock and 160 cm at Molalla (Table 1; Slesak et al., 2016b).
The Matlock site is located on the Olympic Peninsula in WA with soils classified as sandy-skeletal, mixed, mesic Dystric Xerorthents on 0-3% slopes forming in glacial outwash (Soil Survey Staff, 2019). The Molalla site is located in the western foothills of the Cascade Range in OR with soils classified as fine-loamy, isotic, mesic Andic Humudepts on 2-40% slopes forming in residuum derived from hard extrusive igneous rock and soft tuff mixed with volcanic ash in the upper part (Soil Survey Staff, 2019). Soil properties of each site contrasted strongly (Table 1). Soils at Matlock are much younger, and more coarsely textured with higher coarse-fragment content, low water-holding capacity, and lower nutrient contents relative to Molalla (Devine et al., 2011). More information on the sites and soil responses at them can be found in Slesak et al. (2011Slesak et al. ( , 2016aSlesak et al. ( , 2016b.

Design and treatment application
Sites were clear-cut in spring 2003, with trees de-limbed at the stump and merchantable portions removed using mechanized ground equipment along paths designed to minimize soil disturbance (Harrington and Schoenholtz, 2010). After harvest, both sites were planted with Douglas-fir seedlings in February of 2004 at 3 × 3 m spacing (Harrington and Schoenholtz, 2010), and a 2 × 2 randomized complete block design with four replicates was installed. Each of the 16 treatment plots was 50 × 60 m (0.3 ha) in size, which included a 10-m perimeter buffer around a 30-× 40-m measurement plot.
Two factors, harvest intensity and vegetation control with herbicides were tested. Harvest-intensity treatments included bole-only (BO) harvesting, where only the merchantable logs (≥12.7 cm in diameter) were removed and whole-tree (WT) harvesting, where merchantable logs and most of the logging debris (tree tops and branches ≥5.1 cm in diameter) were removed. Within these treatments, logging-debris mass was similar between sites. Estimated mass of logging debris retained (mean ± standard error) was 22.5 ± 3.0 and 13.5 ± 3.0 Mg ha −1 in the BO and WT harvests, respectively, at Matlock, and 24.0 ± 2.8 and 13.9 ± 2.8 Mg ha −1 in the BO and WT harvests, respectively, at Molalla (Harrington and Schoenholtz, 2010). Estimated P retained in logging debris was 6.3 (~30% of preharvest Douglas-fir aboveground P) and 4.6 (~20% of preharvest Douglas-fir aboveground P) kg ha −1 in the BO and WT harvests, respectively, at Matlock, and 8.9 (~30% of preharvest Douglas-fir aboveground P) and 4.0 (~15% of preharvest Douglas-fir aboveground P) in the BO and WT harvests, respectively, at Molalla (T. Harrington, unpublished data).
The two levels of vegetation control included initial vegetation control (IVC) for the first year only or annual vegetation control (AVC) during the first five years (Harrington and Schoenholtz, 2010). Plant communities within the study sites were characterized five years after treatment application (Peter and Harrington, 2012). At Matlock, the three most abundant species and their areal percent cover in IVC plots were oxeye daisy (Leucanthemum vulgare Lam.; 28% cover), hairy cat's ear (Hypochaeris radicata L.; 27% cover), and Douglas-fir (8% cover). Scotch broom (Cytisus scoparius (L.) Link.) invaded at Matlock and three additional weed control spray treatments specifically targeting the species were used in years 1, 4, and 5 in all plots and were successful in eliminating much of the Scotch broom (Peter and Harrington, 2012). By the tenth year of the study (2013), the Scotch broom had reinvaded Matlock, and its cover differed significantly between the BO (3%) and WT (10%) treatments, but not between the AVC (5%) and IVC (7%) treatments (D. Peter and T. Harrington, unpublished data). At Molalla, the most abundant species in IVC plots were California blackberry (Rubus ursinus Cham. & Schltdl.; 55% cover), velvet grass (Holcus lanatus L.; 27% cover), and Douglas-fir (12% cover) (Peter and Harrington, 2012). Plots with AVC had similar species but with decreased canopy covers except for the crop tree, Douglas-fir, which had covers of 15.2 and 15.9% for Matlock and Molalla, respectively. By year 10, Douglas-fir cover differed between IVC and AVC as follows: 26.3 vs. 58.8% at Matlock and 67.1 and 80.7% at Molalla, respectively (Slesak et al., 2016b).

Soil sampling
Soil sampling was conducted pre-harvest in winter 2003 and again ten years following harvest in spring 2013. Soils were sampled from depths of 0-15, 15-30, and 30-60 cm at five points per plot with bucket augers. Samples were collected from each of the four corners of each plot within the plot area bounded by Douglas-fir trees and one from the central point of each plot. Samples were then composited in the field by depth. Pre-harvest (2003) and ten-year (2013) soil samples were airdried (35°C) and returned to the laboratory where the soil was sieved to pass a 2-mm mesh and archived.

Foliage sampling and analysis
Current-year shoots were collected from the second-whorl of Douglas-fir trees in fall 2013 for estimation of foliar nutrient concentrations. Shoots were ground in a Wiley mill to pass a 1 mm sieve, and then dried at 65°C to constant mass. Samples were submitted to Central Analytical Laboratory, Oregon State University, Corvallis, OR and analyzed for N with dry combustion (LECO Corporation, CNS 2000 Elemental Analyzer) and for P with ICP-OES (PerkinElmer Inc., Optima 2100 DV) after dry-ashing in a porcelain crucible and subsequent dissolution in 1.0 M HCL solution.

Sequential phosphorus fractionation
Extractions obtained from sequential P fractionation reflect availability of each P fraction based on susceptibility to each extracting solution (Hedley et al., 1982; Table 2). Sequential fractionation was performed on soils sampled in 2003 (pretreatment) and 2013 (10 years post-treatment). Each sample was subjected to four separate extractions and a final digestion, providing five total fractions. Starting with 1.0 g of soil, a 0.5 M buffered NaHCO 3 solution was used to extract weaklyadsorbed inorganic P and organic P that is easily mineralized. A 0.1 M NaOH solution was then utilized to extract inorganic and organic P that is more tightly sorbed to mineral surfaces, sesquioxides, and P associated with amorphous Al-and Fe-phosphates. A sonic bath was used in a second 0.1 M NaOH extraction to aid removal of inorganic and organic P held in internal surfaces of aggregates and mineral surfaces. Both fractions extracted using NaOH are considered slow-turnover P (Trasar-Cepeda et al., 1990). The final extraction used a 1.0 M HCl solution to remove any acid-soluble inorganic P in calcium-bound P, which is considered to be mostly occluded and slow to become available. For total P, a separate 0.25 g subsample was milled to a fine powder and digested using a Cu-catalyzed Kjeldahl digestion with 5 mL concentrated sulfuric acid, 0.15 g CuSO 3 and 1.5 g K 2 SO 4 refluxed at 425°C for approximately 4 h (Taylor, 2000). Extracts and digestates from each P fraction were analyzed for P concentration on an ICP-AES (SPECTRO Analytical Instruments Gmbh, SPECTRO ARCOS ICP-AES, Kleve, Germany). Residual-P was calculated as the difference between total-P and the sum of the other four fractions.

Table 2
Terminology of phosphorus fractions in the sequential extraction (adapted from Cross and Schlesinger, 1995).

P-fraction Description Extractant
Labile P Readily plant-available inorganic P and organic P that is easily hydrolysable and contributes to plant-available P 0.5 M NaHCO3 Moderately labile P P associated with Fe-and Al-phosphates and organic P contributing to plant-available P 0.1 M NaOH Moderately recalcitrant P P held at internal surfaces of soil aggregates and organic P that is more strongly sorbed to soil surfaces 0.1 M NaOH with sonication Ca-bound P Unavailable inorganic P that is Ca-bound in primary minerals or precipitated Ca-phosphate 1.0 M HCl Residual P Non-extractable P, and P that is considered to be slow turn-over Concentrated sulfuric acid and Cu-catalyzed digestion D. G. DeBruler, et al. Geoderma 350 (2019)

Selective dissolution extraction
To further characterize soils at Matlock and Molalla selective dissolution extractions were used to determine soil Al and Fe concentrations associated with organometallic complexes and short-range order soil minerals (Wada, 1989). Soils sampled from the 0-15 cm, 15-30 cm, and 30-60 cm depths in 2003 (pre-harvest) from both study sites were subjected to three different selective dissolution extractions according to the Kellogg soil survey laboratory methods manual (Soil Survey Staff, 2014a). Ditionite-citrate extractions were used to determine Fe and Al complexed with organic matter or associated with poorly-crystallized hydrous oxides and amorphous aluminosilicates (Wada, 1989). Ammonium oxide, although similar to Dithionite-Citrate, additionally extracts allophone. Finally, sodium pyrophosphate was used to primarily extract organically complexed Fe and Al.

Dithionite-citrate
A 0.75 g sample of air-dry soil was extracted with 0.4 g sodium citrate and 25 mL of 0.57 M sodium citrate solution in a 50 mL centrifuge tube. Samples were shaken overnight (12-16 h), then centrifuged at 4000 rpm for 15 min. Extracts were diluted to a 1:20 ratio for Fe and Al using deionized water. Concentrations of Al and Fe were determined on an ICP-OES (Varian Vista MPX; Varian, Inc.; Palo Alto, CA).

Ammonium oxalate
A 0.5 g sample of air-dry soil was equilibrated for 30 min with 15 mL of 0.2 M ammonium oxalate buffered to a pH of 3.0 in an extraction tube before a 30 min vacuum extraction of 10 mL. Finally, an additional 35 mL of the ammonium oxalate solution was added to a reservoir tube secured to the top of the extraction tube and a total of 45 mL were vacuum extracted over the course of 12 h in the absence of light. Extracts were diluted to a 1:10 ratio using deionized water. Concentrations of Al and Fe were determined on an ICP-OES (Varian Vista MPX; Varian, Inc.; Palo Alto, CA).

Sodium pyrophosphate
A 0.5 g sample of air-dry soil was extracted with 30 mL of sodium pyrophosphate solution buffered to a pH of 10 in a 50 mL centrifuge tube. Samples were shaken overnight (12-16 h), centrifuged at 4000 rpm for 15 min. Extracts were diluted to a 1:5 ratio with deionized water. Concentrations of Al and Fe were determined on an ICP-OES (Varian Vista MPX; Varian, Inc.; Palo Alto, CA).

Statistical analysis
Effect of treatments on soil P fraction concentrations for each depth were analyzed using a two-way factorial analysis of covariance (ANCOVA), with pre-harvest values used as a covariate, blocked for random effects, and using restricted maximum likelihood (REML) estimations. Changes from pre-harvest estimates of each P fraction were used as the primary response variable in most analyses, where change was calculated as the absolute arithmetic difference between pretreatment and 10-yr post-treatment estimates. We used the 90% confidence interval from each change estimate to independently determine if changes in P fractions were significantly different from zero for each treatment. When significant treatment interactions were found, multiple comparisons with Tukey's adjustment were conducted to detect differences between treatment means. Correlation analyses were used to explore associations between P fraction concentrations in the combined 0-30 cm depth and foliar nutrient concentrations under the different treatments. Finally, we used a simple t-test to determine significant differences in Al and Fe, as extracted using SDE, between the two sites. We used α = 0.1 for statistical significance because of low statistical power associated with the level of replication and high inherent variability in soil parameters, and conducted all analyses in JMP, Version 10 (SAS Institute Inc., Cary, NC 2014).

Selective dissolution extractions and pre-treatment phosphorus fractions
Al and Fe concentrations extracted using selective dissolution methods were elevated at both sites (Table 3) and soils sampled from the 0-15 cm depth at Molalla meet the ammonium-oxalate extracted Al (%) + ½ Fe(%) threshold (2%) for Andic soil indicators (Soil Survey Staff, 2014b). Significant differences in Al and Fe concentrations at several depth and SDE method combinations were detected. Sodium pyrophosphate extracted a greater concentration of Al at Molalla than Matlock in the 0-15 cm and 30-60 cm depths. Ammonium oxalate extracted a greater concentration of Al but a lesser amount of Fe from Molalla soils compared to Matlock soils in the 0-15 cm and 15-30 cm depths. A greater concentration of Fe was extracted with Dithionite-Citrate at Matlock than at Molalla in the 30-60 cm depth.
Soil P concentrations prior to treatment were generally higher at Matlock than at Molalla, and decreased with depth at both sites with the exception of the residual P fraction, which was relatively constant across depths (Fig. 1). After accounting for differences in fine fraction (< 2 mm) bulk density, total P mass was comparable between sites (3.3 Mg ha −1 at Matlock and 3.0 Mg ha −1 at Molalla; Table 1), which was surprising given the distinct differences in soil properties and productivity between sites (Table 1). Of all the fractions, only the moderately-recalcitrant and residual P fractions were similar between sites (Fig. 1). At both sites, Ca-bound P comprised the smallest proportion of total P concentration, accounting for 4.7% to 5.2% of total P at Matlock and 1.0% to 1.4% at Molalla. The largest P fraction at both sites and time periods was the NaOH extracted moderately-labile P fraction accounting for 60.8% to 63% and 52.4% to 58.9% of total P at Matlock and Molalla, respectively. a Standard error of mean, n = 3. b Probabilities in bold are significant (α = 0.1).

Phosphorus fractionation -ten-year responses to treatment
3.2.1. Total P Ten-year responses at Matlock and Molalla showed primarily main effects of vegetation control on the change in total-P at most depths (Table 4). At Matlock, total soil P in the 0-15 cm soil depth was significantly reduced in the AVC treatment over the ten-year period compared to the IVC treatment where total P did not change (Fig. 2). A similar effect of vegetation control was observed at 15-30 cm depth, but it was only significant in the BO harvest treatment (i.e., significant interaction, Table 4). At the 30-60 cm depth increment, main effects of both treatments were significant, which arose from a combination of positive change in total P in the IVC and BO treatments and negative change in the AVC and WT treatments (Fig. 2).
At Molalla, treatment effects on the change in total-P only occurred with vegetation control and were limited to the 0-15 cm and 15-30 cm depth increments. At both of those depths, there was a significant reduction in total P in the IVC treatment but not the AVC treatment (Fig. 2). Although there was no difference between harvest treatments, both harvest treatments had significant reductions in total P at the 0-15 cm and 15-30 cm depth increments over the ten-year assessment period (Fig. 2).

Labile P
At both sites, changes in labile P were only associated with main effects of vegetation control at depths > 15 cm (Table 4, Fig. 2). Treatment effects were opposite at the sites, where the change in labile P was higher and positive with IVC at Matlock, but lower and negative with IVC at Molalla. A significant interaction between treatments at 0-15 cm followed this general pattern, except that differences between vegetation control treatments were only significant with the BO harvest intensity (Table 4). Estimated differences in labile P between vegetation control treatments ranged from approximately 24 mg kg −1 (15-30 cm depth) to 7 mg kg −1 (30-60 cm depth) at Matlock, and from 3 mg kg −1 (15-30 cm depth) to 2 mg kg −1 (30-60 cm depth) at Molalla.

Moderately-labile P
There were significant main effects of vegetation control on moderately-labile P at both sites at the 0-15 and 15-30 cm depths, and at 30-60 cm at Matlock only (Table 4). As with labile P, the responses were opposite between sites, with significantly positive and greater change with IVC at Matlock, and significantly positive and greater change with AVC at Molalla (Fig. 3). At Matlock, there was also a significant main effect of harvest treatment at 30-60 cm, where responses differed by approximately 70 mg kg −1 between BO and WT harvests. No other harvest treatment differences were found at either site, but there were significant increases in moderately-labile P at 15-30 cm in both harvest treatments at Matlock, and significant decreases at 0-15 cm for both harvest intensities at Molalla.

Moderately-recalcitrant P
For moderately-recalcitrant P, there was a significant treatment interaction at Matlock at the 15-30 and 30-60 cm depths, but multiple comparisons did not reveal any significant differences among treatment combination means and there was no consistent trend in the response (see footnotes, Table 4). There was also no significant change in the moderately-recalcitrant P fraction for any treatment or depth at the Matlock site. In contrast, the change in moderately-recalcitrant P was significantly different between vegetation control treatments at all depths at the Molalla site (Fig. 3). In all cases, the IVC treatment was associated with significant reductions in moderately recalcitrant P with little change in the AVC treatment over the ten-year period. There were no significant effects of harvest intensity on the change in moderatelyrecalcitrant P at Molalla, but the fraction was significantly reduced in the BO treatment (Fig. 3).

Ca-bound P
At Matlock, there was a significant main effect of vegetation control on Ca-bound P at the 0-15 cm and 15-30 depth increments (Table 4). In both instances, change in Ca-bound P was negative and significantly different from zero in the AVC treatment whereas change with IVC was positive (Fig. 4). At Molalla, the change in Ca-bound P was different between vegetation control treatments at 0-15 cm depths only, where IVC caused a significant reduction in the P fraction. With the exception of the AVC treatment at 0-15 cm, all treatments and depths had significant reductions in the Ca-bound P fraction at Molalla, which were approximately 20-30% of pretreatment estimates (Fig. 4).

Residual P
There were no treatment effects on residual P at Matlock, but all treatments were associated with significant reductions at the 15-30 cm depth increment (Table 4, Fig. 4). At Molalla, there were significant main effects of vegetation control at the 15-30 cm depth only. In this instance, the response (IVC > AVC) was opposite of that observed for other fractions, with the AVC treatment having a significant reduction in residual P (Fig. 4). With the exception of the IVC treatment at 15-30 cm depth, residual P was significantly reduced in all treatments at the 15-30 cm and 30-60 cm depths at Molalla.

Foliar nutrition and correlations with P fractions
There were no significant treatment effects on N and P concentrations and N:P ratios in Douglas-fir foliage at either site, and no apparent trends in the response (data not shown). Foliar N concentrations across treatments were 1.28% ( ± 0.03) and 1.44% ( ± 0.06) at Matlock and Molalla, respectively. Foliar P concentrations across treatments were 0.24% ( ± 0.01) and 0.21% ( ± 0.01) at Matlock and Molalla, respectively. Mean foliar N:P ratios were 5.3 ( ± 0.3) at Matlock, and 7.0 ( ± 0.6) at Molalla.
We limited the foliar-P fraction correlation assessment to the vegetation control treatments only (pooled across both harvest treatments). For this assessment, we used a mean P fraction value for 0-30 cm from samples collected at the ten-year post-treatment period. At Matlock, there were no significant correlations between any of the foliar variables (N %, P %, and N:P) and P fractions in the IVC treatment (data not shown), but there was a significant positive correlation in the AVC treatment between the ten-year moderately-recalcitrant P fraction and foliar P concentration (r = 0.82, p = .013). At Molalla, the Cabound P fraction was negatively correlated with foliar N concentration (r = −0.66, p = .074). There were no other significant correlations among any of the P fraction and foliar variables (data not shown).

Discussion
Our results indicate potential for reductions in soil P fractions over time with the more intensive practices of WT harvest and AVC at Matlock. In contrast, most of the treatments at Molalla reduced P fractions except for the AVC treatment, where P fractions were either maintained or increased in some instances. Although decreases were detected, total P concentrations of surface soils after ten years were still  relatively high at both sites (1500 mg kg −1 and 970 mg kg −1 at 0-15 cm at Matlock and Molalla, respectively) in comparison to total P of 120 mg kg −1 for some prairie soils (Caldwell, 2006) and 134 mg kg −1 for loblolly pine (Pinus taeda) stands in the southeastern U.S. (Kiser, 2011). Given the relatively large P reserves at our sites, and the small amount of ecosystem P removed in both harvest treatments, it seems unlikely that harvest intensity will induce P limitations to productivity over the course of multiple rotations. Comparisons of the individual fraction ratios with those of the Walker and Syers (1976) model indicate these sites have nearly exhausted reserves of Ca-bound P, which occurs at nearly 20,000 years in some chronosequences (Crews et al., 1995). Relatively lesser amounts of Ca-bound P at Molalla indicate more weathered soils compared to Matlock. Correspondingly, increased proportions of moderately-labile P and residual P fractions at Molalla further support this.
The large proportion of moderately-labile P at both sites is in agreement with the findings of Yang and Post (2011), who concluded that this fraction was a significant proportion across most soil orders and weathering stages. Additionally, the significant labile and moderately-labile pools may be a result of the presence of high contents of Fe and Al associated with hydrous oxides, amorphous alumnosilicates, and organometallic complexes as determined by our selective dissolution extraction. These minerals and complexes are generally considered to be highly reactive with soil P and likely represent a plant available pool (Saunders, 1964;Hodges and Zelazny, 1980;Borggaard et al., 1990;Guo and Yost, 1998;Harsh et al., 2002). Our SDE analysis showed a significantly greater concentration of Al as extracted by both sodium pyrophosphate and ammonium oxalate at Molalla than at Matlock. If both extractions are associated with organo-metallic complexes and short-range order minerals (Wada, 1989), respectively, we would also expect the Molalla soils to be more retentive of soil P compared with Matlock (Saunders, 1964;Hodges and Zelazny, 1980), especially since reactive Al content would be even higher at Molalla (i.e., higher clay content and lower coarse fragment content at Molalla; Table 1). This, however, does not appear to be directly supported by this study as Molalla soils contained a lower concentration and lower total content of Fig. 2. Absolute change in soil total P and labile P concentrations by treatment, depth, and site ten years after harvest with different harvest intensities (bole-only (BO) vs whole-tree (WT)) and vegetation control (annual vegetation control (AVC) vs initial vegetation control (IVC)). Error bars represent 90% confidence intervals for each least squares mean with significant differences indicated by an asterisk. Note differences in scale between sites for the labile P fraction. soil P fractions. Additionally, there were a greater number of significant decreases in P fractions at Molalla than at Matlock. Only ammonium oxalate extracted Fe concentration was greater at Matlock than Molalla. This suggests some other mechanism, besides the differences in the preharvest concentrations of Al and Fe associated with poorly crystalline minerals, is responsible for the contrasting responses of some of the P fractions to treatments. One possibility is Fe and Al concentrations in post-harvest samples may have responded to the different harvest and herbicide treatments (not measured in this study). Effects on Fe and Al contents have been observed in a long-term (44 years) study where increases in P associated with amorphous forms of Al and Fe, were positively correlated with manure fertilization (Pizeeghello et al., 2014). Additionally, because of the differences in bulk densities at these sites, trends in soil P fraction content with soil Fe and Al content (not explored in this study) may differ from those explored in terms of concentration.
Phosphorus fractions generally decreased with depth at both Matlock and Molalla prior to harvest, which is consistent with observations across a variety of managed systems (Spears et al., 2001;Wang et al., 2007;Maranguit et al., 2017). One exception to this trend is the typical increase in Ca-bound P with depth (Chen, 2003) that is associated with increases in soil pH. Ca-bound P did not increase with depth at our sites even though soil pH did (Table 1), suggesting a different geochemical control in these subsoils compared to other managed systems, such as a greater content of poorly-crystalline minerals. Wood et al., 1984, showed that the HCl could solubilize Fe bound P associated with poorly crystalline minerals. Our SDE analysis showed typical decreases in ammonium oxalate extractible Fe concentrations with increasing depth which is matched with decreases in Ca-bound P with increasing depth.
In addition, the residual-P fraction concentration remained relatively consistent across depths at both sites but the ratio of residual P to total P increased. Compton and Cole (1998) observed a similar trend and suggested if the residual P pool is primarily organic P, it may represent a significant source of labile P available to plants (Johnson et al., 2003;Gahoonia and Nielsen, 1992). Although we did not quantify the proportion of organic P in the residual pool, the general treatment responses indicate that most of the changes in P fractions are Fig. 3. Absolute change in soil moderately-labile P and moderately-recalcitrant P concentrations by treatment, depth, and site ten years after harvest with different harvest intensities (bole-only (BO) vs whole-tree (WT)) and vegetation control (annual vegetation control (AVC) vs initial vegetation control (IVC)). Error bars represent 90% confidence intervals for each least squares mean with significant differences indicated by an asterisk. associated with the relative amount of vegetation present, suggesting that biological controls are highly important to changes in the P fractions and cycling at these sites.
Because most treatment effects were associated with the vegetation control treatment, differences in site-specific responses to treatments are likely associated with differences in vegetative communities among treatments and between sites, and their interactions with physio-chemical characteristics of each soil. For example, Scotch broom cover at Matlock was extensive and significantly greater in WT than in BO throughout the study period (D. Peter and T. Harrington, unpublished data). Scotch broom is an N-fixing invasive plant capable of increasing phosphatase activity and decreasing total P (Caldwell, 2006). At Matlock in the 30-60 cm depth, total P, labile P, and moderately labile P each were lower in WT than in BO, which could have resulted from greater uptake of plant-available P by the rapidly growing Scotch broom. Although isolated patches of Scotch broom were found at Molalla, Scotch broom cover did not vary significantly by harvest intensity or vegetation control treatments. In a small-plot (28 m 2 ) study focused on presence versus absence of Scotch broom, Slesak et al. (2016a) found that the moderately-available P fraction was reduced in the presence of Scotch broom at Matlock, but not at Molalla. Slesak et al. (2016a) reported that plots at Matlock were located in a uniform stand of Scotch broom, whereas plots at Molalla were centered on individual Scotch broom plants. Giardina et al. (1995), found that phosphatase activity and available P were higher under a mixed-species Douglas-fir stand in the presence of N-fixing red alder, which is more similar to the uniform Douglas-fir/Scotch broom mixture at Matlock. Given the above, reductions in P with WT at Matlock are likely driven indirectly by the influence of that treatment on Scotch broom abundance rather than a direct result of increased P removal at time of harvest.
At Molalla, nutrient uptake by vegetation in general may have played a larger role in soil P changes over the ten-year study period. Phosphorus fractions are generally smaller (Fig. 1), and site productivity higher, at Molalla compared to Matlock (Slesak et al., 2016b). This relatively greater imbalance between P supply and demand could lead to the observed P reductions in soil. Increased or maintained P fractions in the AVC treatment supports this, because total vegetation (Douglas-fir and competition) uptake in that treatment would have Fig. 4. Absolute change in soil Ca-bound P and residual P concentrations by treatment, depth, and site ten years after harvest with different harvest intensities (boleonly (BO) vs whole-tree (WT)) and vegetation control (annual vegetation control (AVC) vs initial vegetation control (IVC)). Error bars represent 90% confidence intervals for each least squares mean with significant differences indicated by an asterisk. Note differences in scale between sites for the Ca-bound P fraction. been much lower compared to the IVC treatment. Treatment-induced shifts in community composition could have also contributed to the response at Molalla (Peter and Harrington, 2012).
Maintained or increased labile and moderately labile P at Matlock at the 30-60 cm depths in the BO harvest treatment might be attributed to the direct effect of debris on the amount of P removed during harvest, its indirect influence on belowground decomposition (Devine and Harrington, 2007;Sanchez et al., 2006), or the rate at which plants access P and translocate it to surface soils (Eger et al., 2018). The response may also be associated with the aforementioned Scotch broom effect because the BO treatment had a markedly lower amount of Scotch broom present at the Matlock site (Harrington and Schoenholtz, 2010;D. Peter and T. Harrington, unpublished data). Greater overall P retention with BO harvest or increased P cycling with broom (Caldwell, 2006;McGill and Cole, 1981) would allow for greater stand uptake, but some of the available P would be fixed onto mineral surfaces and lead to the observed increases.
Relative decreases in labile P, coupled with increases in exchangeable Ca measured at the same time as this study (Slesak et al., 2016b) could indicate precipitation of calcium phosphates in these soils (Tan, 1982). Giardina et al. (1995) also found increases in soil-solution Ca 2+ and phosphates, which were correlated with greater concentrations of Ca-bound P. However, we did not measure any increases in the Cabound P fraction. Organic molecules in solution also occlude phosphate compounds and Ca 2+ from dissolution during fractionation (Lindsay, 1979;Giardina et al., 1995). Similar studies to ours (Cross and Schlesinger, 1995;Giardina et al., 1995) suggest that the Ca-bound P fraction, extracted by HCl, may include other HCl-soluble compounds besides primary mineral P and calcium phosphates such as HCl soluble Fe and Al bound P as mentioned previously mentioned in this discussion. Further evidence of this possibility was provided by He et al. (2006), who used an enzymatic assay analysis paired with oxidative autoclaving to show that some soil samples contained significant organic P concentrations in the HCl-extracted fraction. This exemplifies the difficulty in interpreting the degree of availability or occlusion within these fractions (Condron and Newman, 2011).
Decreases in residual P over the ten-year study period were observed across most treatment and depth combinations at both sites relative to the pre-treatment concentrations. These decreases were likely related to shifts in P demand associated with changes in vegetative nutrient demand between pre-harvest and year 10 (Vitousek et al., 2010). The residual P fraction is generally considered unavailable (Cross and Schlesinger, 1995). However, recent studies have shown the residual fraction may be more available than previously reported (Condron and Newman, 2011;Vitousek et al., 2010), especially in slightly weathered soils (Guo et al., 2000). This fraction has been shown to support plant growth in forests and to turn-over significantly during a 65-year period (Hedley et al., 1982;Chen et al., 2000;Chatterjee et al., 2013).
Low correlations between soil P pools and foliar P in Douglas-fir in this study suggest plots with higher soil P do not necessarily supply more P to Douglas-fir trees growing on these sites, which is further supported by lack of treatment effects on foliar N and P. Low N:P in foliage (< 10) at both sites and all treatment plots indicates N limitation is more likely than a P limitation (Güsewell, 2004). In addition, foliar P concentrations at both sites were above levels considered adequate for growth (Ballard and Carter, 1986). Therefore, it is likely that P is not limiting to plant growth in these soils, but risks of P limitation may be higher at the Molalla site relative to the younger Matlock site (Root, 2017). Nonetheless, the Molalla site shows little evidence of developing P deficiencies in the near future compared to some forested ecosystems (Crews et al., 1995;Riley and Vitousek, 1995;Vitousek and Farrington, 1997).

Conclusions
Results of this study identified changes in soil P fraction concentrations at 0-60 cm during the first 10 years of Douglas-fir stand growth, depending on site characteristics, harvest intensity, and vegetation control treatments. We also detected significant differences in Al and Fe concentrations as extracted by sodium pyrophosphate and ammonium oxalate between the two sites. The less-intensive treatments of BO harvesting and vegetation control limited to the first year after tree planting (IVC) at Matlock ameliorated many of the decreases in P fraction concentrations, whereas five years of annual vegetation control mediated decreases in P concentrations at Molalla. The contrasting responses to vegetation control at these sites suggest that vegetation abundance and composition has an overwhelming influence on soil P fraction dynamics following harvest. On-site observations at year-ten included less abundance of the invasive Scotch broom under BO harvest than WT harvest at Matlock. Therefore, on sites where BO is effective in reducing Scotch broom and other invasive species, this treatment may be preferred to VC; as both broom cover and soil P loss were reduced. In contrast, sites such as Molalla, where harvest intensity had little influence on the vegetation community, VC may be preferred. These responses and changes in fractions most associated with organic P suggest that biologic controls of P pools are more important than geochemical ones at these two sites, despite elevated concentrations of Fe and Al extracted with SDE methods. With high foliar P concentrations, low N:P, and the lack of consistent positive correlations between P fraction concentrations and foliar P concentrations, there is little evidence to suggest more intensive harvesting has an effect on foliar nutrition at these sites, and therefore little effect on productivity by way of P nutrition manipulation. Continued analysis of changes in P fractions and cycling in response to contrasting management regimes and across more site conditions is needed to understand if these P pools change over longer periods of time more relevant to ecosystem development and management (i.e., rotation length).