Introduction

In terrestrial ecosystems the availability of phosphorus (P) for plants is often limiting primary productivity (Gradowski and Thomas 2006; Elser et al. 2007; Vitousek et al. 2010). For forests it has been suggested that at sufficient P supply in the soil the acquisition of P provided from chemical weathering of primary minerals may sustain tree growth (P acquiring systems; Lang et al. 2016). However, at P deficient conditions, P supply from the mineral phase may not be sufficient to maintain the physiological needs of trees and there is an increasing need to re-utilize P from organic litter, i.e., to recycle P especially within organic layers and the upper mineral horizons (P recycling systems; Lambers et al. 2008; Lang et al. 2016, 2017). These concepts tend to ignore that there might also be another strategy of trees to cope with limited P supply in the topsoil: instead of or in addition to improving P recycling efficiency, trees might increase P uptake from subsoil P reserves, i.e. the soil below 20–30 cm depth, as well as from saprolite (Uhlig et al. 2020) and weathered bedrock. While nutrient uptake from subsoil has been studied occasionally for elements such as nitrogen, calcium and potassium (da Silva et al. 2011; Maeght et al. 2013), these studies did not consider the role of subsoils for P nutrition of temperate forest ecosystems.

Studies in agricultural cropping systems have reported that plants can cover between 10 and 80% of their P demand from subsoils depending on the extent of their root system (Kuhlmann and Baumgärtel 1991; Kautz et al. 2013). Trees in temperate forests usually have much deeper rooting systems than most agricultural crops (Canadell et al. 1996; Fan et al. 2017), which may even reach into the saprolite and weathered bedrock below, although the amount of fine roots in the subsoil strongly depends on soil chemical properties (Strong and La Roi 1983) and profile depth (Kirfel et al. 2019). P uptake via the fine roots of trees is promoted by symbiotic fungi (Plassard and Dell 2010) and, while specific information on the contribution of subsoil inhabiting fungi to tree P uptake are missing, it has been reported that fungal biomass in the rhizosphere of deep tree roots is several times higher than in the surrounding bulk subsoil (Maeght et al. 2013). It can therefore be assumed that P uptake from the subsoil, and potentially also from saprolite, is of at least equal or of even greater importance in forest ecosystems than in agricultural ecosystems, especially when overall P supply is low.

There is sporadic evidence that trees do not only rely on labile P sources, which are available in the short-term, but also incorporate P from more stable sources such as P associated with Ca-minerals (Niederberger et al. 2017). To account for different P sources, phosphorus in soils can be operationally separated into pools of decreasing extractability and availability to plants according to the Hedley sequential fractionation procedure (Hedley et al. 1982). Within this fractionation scheme, P associated with Ca-P minerals is conventionally assumed to be included in the 1 M HCl extract. These Ca-P minerals are partially derived from the parent material, but might also include secondary minerals precipitated from previously biologically cycled phosphate (Angert et al. 2012; Joshi et al. 2016).

Information on the source of Ca-P minerals included in the 1 M HCl extract cannot be obtained from the sequential fractionation, but become apparent from the oxygen isotopic composition of the phosphate (conventionally expressed in delta notation as δ18OP values in ‰). The δ18OP values of secondary Ca-P minerals differ from those of primary (non-recycled) minerals of the parent material, because biological activity, namely the activity of enzymes, will catalyze oxygen isotope exchange (Liang and Blake 2006; von Sperber et al. 2015). Especially the activity of the ubiquitous intracellular enzyme pyrophosphatase affects the δ18OP values due to the repeated recycling of phosphate ions, by which oxygen atoms will be exchanged between phosphate and the surrounding water (Cohn 1958; Blake et al. 2005; von Sperber et al. 2017). During this exchange, a thermodynamic isotope fractionation occurs, shifting the isotopic composition away from the source signal of primary minerals towards an equilibrium distribution of oxygen isotopes between phosphate and ambient water (Chang and Blake 2015). Thus, an interpretation of δ18ΟP values relative to this expected equilibrium δ18ΟP value allows an assessment whether P in the subsoil has been biologically cycled or not. For the sites studied here with igneous and metamorphic parent material, rather low δ18ΟP values of 0–12‰ are expected (Tamburini et al. 2014), and under the given climatic conditions an equilibration with soil water by enzymatic activity should result in a shift towards higher δ18ΟP values. This stable isotope tracer has been applied successfully in microbial incubation studies (Blake et al. 2005; Stout et al. 2014) and a number of ecosystem studies, e.g., along climatic gradients (Angert et al. 2012; Helfenstein et al. 2018), along a young soil chronosequence located in the Swiss Alps (Tamburini et al. 2012), for arable cropping systems (Amelung et al. 2015; Joshi et al. 2016; Bauke et al. 2018) or for tracing sources of atmospheric P inputs (Gross et al. 2013, 2016), but—to the best of our knowledge—it has not been applied to elucidate P cycling in subsoils of forest ecosystems.

Biological cycling of P in soils will also affect P speciation (Turner et al. 2007; Chen et al. 2008; Wang et al. 2019) and the spatial distribution of P, e.g., in the rhizosphere relative to the bulk soil (Richardson et al. 2009). Compared to the Hedley sequential fractionation more detailed information on P speciation can be obtained by synchrotron-based X-ray absorption near-edge structure (XANES) spectroscopy (Kruse et al. 2015; Prietzel et al. 2016; Koch et al. 2018). XANES spectra provide rather quantitative information on the abundance of specific P-minerals and thus help identify changes in the predominant P species (Kruse et al. 2015), e.g., a shift from primary minerals such as apatite to a greater contribution of organic P forms. In addition to this, elemental patterns in undisturbed soil samples at very high spatial resolution can be determined using the high sensitivity and lateral resolution of nano-scale secondary ion mass spectrometry (NanoSIMS) (Mueller et al. 2013; Steffens et al. 2017; Werner et al. 2017).

The mobilisation of P from the soil depends to a great extent on microbial activity (Gyaneshwar et al. 2002; Hallama et al. 2019). Therefore, potential limitations in microbial growth are a crucial factor determining the efficiency of P cycling in the subsoil as well as the supply of P to the tree roots. Nutrient limitation of microbial growth can be identified from the respiratory response of soil microorganisms to glucose and nutrient amendments (Anderson and Domsch 1980; Scheu 1993). Previous studies in temperate ecosystems have shown that microbial activity is frequently limited by carbon (C) supply, but also by limited availability of nitrogen (N) and P (Kamble and Bååth 2014; Meyer et al. 2018).

To study forest P cycling under conditions of different P stocks, Lang et al. (2017) identified forest sites in Germany, which are characterized by differences in P stocks while stand age, vegetation, climate and topography were similar among the sites. We analysed the biogeochemistry of P cycling in mineral topsoils and subsoils of one site with high total P stocks (Bad Brückenau [BBR]), two sites with intermediate total P stocks (Mitterfels [MIT] and Vessertal [VES]), and one site with low total P stocks (Conventwald [CON]). We expected that δ18ΟP values would be at equilibrium with ambient water at sites with low P contents, which would indicate that the majority of P in this pool was part of the biological cycle, whereas δ18ΟP values at sites with high P contents would be in the range expected from primary minerals, indicating a major contribution of non biologically cycled P. We further expected that the transition zone, where δ18OP values changed from a biological signal to a primary mineral signal (Amelung et al. 2015), would be deeper at low P sites compared to sites rich in P. In addition to δ18ΟP values, we recorded XANES spectra for selected depths, to identify at which depth apatite, as a proxy for non-recycled Ca-P minerals, occurs as one of the dominant P species. We used NanoSIMS to visualize P distribution around roots at different depths, assuming that these microscale patterns might reflect whether P acquisition had occurred or not. Finally, we incubated samples with addition of glucose in combination with N, P or both N and P, to determine whether microbial growth and thereby microbial P recycling was limited by nutrient availability. In combining these methods, we tested the hypothesis that, in addition to a more efficient recycling of organic P in topsoil as suggested by Lang et al. (2017), P-limiting conditions in temperate forest ecosystems will also result in increased utilization of P from deeper soil horizons.

Material and methods

Field site and sampling

The four German temperate forest sites studied here are located in the Rhön (Bad Brückenau, BBR, 11°31′ E; 51°31′ N), in the Bavarian Forest (Mitterfels, MIT, 12°52′ E; 48°57′ N), in the Thuringian Forest (Vessertal, VES, 10°46′ E; 50°36′ N) and in the Black Forest (Conventwald, CON, 7°57′ E; 48°01′ N) and are described in detail in Lang et al. (2017). The sites are situated at an elevation between 809 and 1033 m a.s.l. with annual average temperatures ranging from 4.5 to 7.9 °C and annual average precipitation from 1031 to 1395 mm (Table 1).

Table 1 General site characteristics of the four forest sites

The soils at the four sites under study differ in their total P concentrations (BBR > MIT > VES > CON), while stand age, vegetation, topography and climate are very similar. All soils are classified as Cambisols according to WRB (IUSS Working Group WRB 2015) (Table 1). Topsoil pH varies between 4.1 and 4.6 (Table 1). We collected horizon-specific soil samples of the mineral soil to a depth of approximately 2 m. The organic layers of the forest floor were not analyzed here, as their contribution to overall P cycling has already been elucidated in other studies (e.g., Lang et al. 2017). Samples were collected as composites from three sides of the quantitative soil pits described in Lang et al. (2017). Samples were separated in rhizosphere soil and bulk soil directly in the field, with rhizosphere soil being the soil adhering to the roots after shaking. Samples for microbial analyses were kept frozen until processing. At CON and MIT, additional diesel-powered wireline core-drilling was performed to 20 m depth at CON and 30 m depth at MIT to sample subsoil beyond 2 m depth, saprolite (defined as chemically weathered bedrock with physically intact bedrock fabrics (Becker 1895)), fractured bedrock and non weathered bedrock. The drilling method was described in detail in Uhlig and von Blanckenburg (2019). The interface between saprolite and weathered bedrock was found at 7 m at CON and 17 m at MIT. Non weathered parent bedrock was found at 15 m at CON and 26 m at MIT. Bedrock samples from drill cores were ground, while soil samples were air-dried and sieved to < 2 mm. Analyses in this study were limited to the mineral soil horizons and selected depths in the drill cores.

Hedley sequential P fractionation

Data for the sequential extraction of P fractions for the soil profiles up to 1 m depth were reported in Lang et al. (2017). Here, we additionally analysed samples from the drill cores in MIT and CON. For these samples, P fractions were analyzed in duplicates following the Hedley sequential fractionation protocol according to Tiessen and Moir (1993). In each extract, inorganic P (Pi) content was determined by the molybdenum blue method (Murphy and Riley 1962) using a spectrophotometer (Specord 205; Analytic Jena, Germany). Total P contents (Pt) were measured on respective extracts by inductively coupled plasma with optical emission spectroscopy (ICP-OES; Ultima 2, Horiba Scientific, Jobin Yvon S.A.S., France). Organic P content (Po) was calculated as the difference of Pi and Pt contents. According to Cross and Schlesinger (1995) and Negassa and Leinweber (2009), the soil P-pools obtained by Hedley sequential fractionation were classified by their chemical extractability as labile P (resin-P and NaHCO3-Pi and -Po), moderately labile P (NaOH-Pi and -Po) and stable P (1 M HCl-Pi, concentrated HCl-Pt and residual-Pt).

Oxygen isotopic composition of phosphate extracted by 1 M HCl

For the phosphate O-isotope ratio analyses, soil samples were sequentially extracted with 0.5 M NaHCO3, 0.1 M NaOH and 1 M HCl according to the Hedley sequential extraction procedure. However, to increase the P concentration in the final 1 M HCl-extract, samples were extracted at a soil:solution ratio of 1:10 instead of the ratio of 1:60 indicated in the original protocol of the Hedley sequential extraction. The NaHCO3 and NaOH extraction steps were needed to remove the highly and moderately labile Po forms to reduce the concentration of Po in the 1 M HCl extracts, to avoid a bias in measured isotope data due to hydrolysis of Po during the purification procedure. Inorganic P from the 1 M HCl extracts was then further purified according to the stepwise procedure described by Tamburini et al. (2010), which consists of precipitation as ammonium phosphomolybdate and magnesium ammonium phosphate, removal of cations with a cation exchange resin (Dowex 50X8, 200–400 mesh, Sigma–Aldrich, Darmstadt, Germany) and yields silver phosphate (Ag3PO4) as the final precipitate.

Phosphate O-isotope ratios in Ag3PO4 were analysed by the Group of Plant Nutrition at the Institute of Agricultural Sciences (ETH Zurich, Switzerland). Ag3PO4 was analysed by online high-temperature thermal decomposition using a Vario Pyro Cube (Elementar, Hanau, Germany) connected in continuous-flow to an isotope ratio mass spectrometer (Isoprime 100, Manchester, UK) with precision better than ± 0.4‰, calculated on replicate analysis of standards (Vennemann et al. 2002). For complete thermal decomposition a minimal amount of soot was added to weighed samples (lamp black conditioned, CAS-Nr. 23.00–0097; Elementar, Hanau, Germany). Phosphate O-isotope ratios were calibrated against an internal Ag3PO4 standard (Acros Organics, Geel, Belgium; δ18O = 14.2‰) and two benzoic acid standards distributed by the International Atomic Energy Agency (IAEA) (IAEA 601: δ18O = 23.1‰ and IAEA 602: δ18O = 71.3‰). All O-isotope data are reported in delta notation in permil (‰) relative to the Vienna Standard Mean Ocean Water (VSMOW).

To determine the effect of biological P cycling on δ18OP values at each site, we calculated the theoretical values expected for the complete equilibration of oxygen isotope ratios between water and phosphate based on the equation given by Chang and Blake (2015):

$${\delta }^{18}{\mathrm{O}}_{\mathrm{P}}={e}^{(\frac{14.43}{T }-\frac{26.54}{1000})} \times \left({\delta }^{18}{\mathrm{O}}_{\mathrm{W}}+1000\right)-1000$$

where T is the ambient temperature in [K] and δ18OP and δ18OW are the oxygen ratios of phosphate and water, respectively. Not for all sites δ18OW values could be obtained from GNIP (Global Network of Isotopes in Precipitation) stations at similar elevation. Therefore, reference δ18OW values were calculated based on the models by Bowen and Wilkinson (2002) and Bowen and Revenaugh (2003), resulting in δ18OW values of − 10.0‰ (BBR), − 9.6‰ (MIT), − 9.7‰ (VES) and − 9.1‰ (CON). These values were within 0.5‰ of the values recorded by GNIP stations at similar elevation within a 100 km radius. For temperature we used the mean annual temperatures given in Table 1.

XANES spectroscopy

Phosphorus K-edge XANES spectra were recorded at the Soft X-ray Micro-characterization bending magnet beamline (SXRMB, 06B1-1) at the Canadian Light Source synchrotron, Saskatoon, Canada. Air-dried and homogenized bulk soil samples from selected depths of the four sites BBR, MIT, VES and CON were spread as a thin film of fine powder onto a double-sided carbon tape and mounted onto a copper sample holder before being placed in the vacuum chamber. The spectra were collected in fluorescence yield mode (samples) and total electron yield mode (reference standards), respectively, at photon energies between 2130 and 2200 eV with dwell time of 4 s (samples) and 1 s (reference compounds), respectively. The absolute energy was calibrated using the Ar K-edge around 3205 eV; this means that the main peak position of Al phosphate (AlPO4) was around 2152.9 eV. At least two scans were recorded for each sample. Subsequent data treatment and evaluation such as normalization to the intensity of the incident beam (I0), spectra averaging, background correction and normalization as well as linear combination fitting (LCF) were done using the ATHENA software package (Demeter 0.9.20; Ravel and Newville 2005). Linear combination fitting was done on averaged, normalized spectra in the energy range between − 10 eV and + 35 eV relative to the E0 (zero-crossing of the second derivative) for all possible binary to quaternary combinations of the following reference standards: CaHPO4, CaHPO4·2H2O, Ca(H2PO4)2·H2O, Ca10(PO4)6(OH)2, Mg3(PO4)2·8H2O, Mg2O7P2, MgHPO4·3H2O, FePO4·4H2O, AlPO4, P adsorbed on gibbsite, P adsorbed on goethite, and C6H18O24P6·Na·yH2O (phytic acid sodium salt hydrate) as proxy for organic P (Gypser et al. 2018). The R-factor values were used as a goodness-of-fit criterion and significant differences between fits were evaluated using the Hamilton test (Calvin 2013; p < 0.05) with the number of independent data points calculated by Athena estimated as data range divided by core–hole lifetime broadening. If by adding a further reference compound to a fit the R-factor was not significantly better according to the Hamilton test, the fit with the lower number of reference compounds was chosen. If fits with the same number of reference compounds were not significantly different from each other, fit proportions were averaged.

Nano-scale secondary ion mass spectrometry

Nano-scale secondary ion mass spectrometry (NanoSIMS) was used to acquire images of the spatial element distribution in the rhizosphere region of soil aggregates. Selected rhizosphere samples from 0–6 cm and 60–80 cm depth from MIT and VES were used for this analysis. The intact soil aggregates < 2 mm including roots were air-dried and subsequently impregnated with Araldite 2020 epoxy resin (Mueller et al. 2013). The sectioned and polished embedded samples were coated with gold/palladium by physical vapour deposition under argon atmosphere to avoid charging during the NanoSIMS measurements. Before NanoSIMS analysis, areas of interest for the sequence from the root into the soil matrix samples were chosen using an optical stereo microscope (V12, Zeiss, Germany). The soil aggregate cross sections were analysed along a transect (5–8 measurements) from the root into the mineral soil (see Fig. S1 in supplementary material) in imaging mode at a Cameca NanoSIMS 50L (Cameca, France). The dwell time was 10 ms pixel−1 and 512 × 512 pixels were acquired for the 40 × 40 µm2 ion images. The spatial distribution of 12C, 16O, 12C14N, 31P16O2, 27Al16O and 56Fe16O for each measurement spot was recorded. The distribution of 12C was used as an indicator of the epoxy resins, whereas 12C14N indicated the presence of organic matter and 16O, 27Al16O and 56Fe16O were used as proxies for the soil mineral phase. For data evaluation the counts of secondary ions were normalised to the absolute signal intensity. To visualise the spatial distribution, the normalised ion counts of the 12C14N and 31P16O2 secondary ions were smoothed along the transect from the root into the surrounding soil matrix with a moving-average-algorithm over a distance of 10 µm (window size of 128 pixel) (Kendall et al. 1983). An example of NanoSIMS images and the position of the transect scans is shown in Fig. S1 (Supplementary Material).

Microbial biomass and nutrient limitation of microbial growth

Frozen soil samples were gently thawed at 5 °C for 72 h. Depending on the soil horizon, 2–4 g of the upper (~ 0–30 cm), and 5–8 g of the lower horizons (~ > 30 cm depth) were used for measurements of microbial biomass carbon (Cmic) and nutrient limitation of microbial growth. The oxygen consumption rate of the soil samples was measured at 22 °C every 15 min over 72 h using an automatic O2 microcompensation system (Scheu 1992). Cmic was measured by substrate-induced respiration (SIR, Anderson and Domsch 1978) after adding glucose in aqueous solution (12 mg g–1 dry soil for the upper, and 4 mg g–1 dry soil for the lower horizons). The maximum initial respiratory response (MIRR; µl O2 g–1 dw h–1) was taken as the mean of the three lowest measurements during the initial nine hours. Cmic (µg Cmic g–1 dw) was calculated as \({\mathrm{C}}_{\mathrm{mic}}=38 \times \mathrm{MIRR}\) (Beck et al. 1997). Nutrient limitation of microbial growth was investigated by adding N (as NH4NO3(aq)) and P (as KH2PO4(aq)) along with C (as glucose) in the combination C + N, C + P and C + N + P at C:N:P mass ratio of 10:2:1 (Anderson and Domsch 1980). Nutrients were added to the soil samples in solution with 100 µl nutrient solution per g soil fresh weight. For better infiltration of the nutrient solution into the soils, additional 100 µl H2Odest. were added per g soil fresh weight afterwards. Nutrient limitation of microbial growth was quantified as additional cumulative microbial respiration (oxygen consumption) up to the first maximum respiration rate (AMR, Scheu 1993; Tiunov and Scheu 1999). Data are given as relative increase (in %) of respiration after C + N, C + P, C + N + P addition in relation to the treatment where only C was added. For interpretation of nutrient limitation data, we used the increase of AMR after nutrient addition relative to that of only C addition. If AMR, for example, was greater after C + N addition than after C + P addition, this was considered as an indication that microbial growth was more limited by N than by P.

Statistical analysis

Statistical data analysis was performed using Sigma Stat (version 11.0, Jandel Scientific) or StatSoft (ver. 8). Soil and rhizosphere δ18OP values were compared by paired t-test (p < 0.05), correlations among counts of elements recorded in the NanoSIMS images were tested using Spearman correlation coefficients (p < 0.05). Statistical analyses for evaluation of XANES spectra are explained in detail above.

Results

Oxygen isotopic composition of 1 M HCl-extractable phosphate and P speciation

Inorganic P concentrations in the 1 M HCl-extracts decreased in the order BBR > VES > MIT > CON. At the sites where depths below two meters were sampled (MIT and CON), HCl P concentrations remained mostly constant until 10 m depth at CON but showed strong variation with depth at MIT (Fig. 1).

Fig. 1
figure 1

Concentration and δ18ΟP values of the 1 M HCl-Pi pool from Hedley sequential fractionation for different soil depths at the sites with low (MIT, CON; two left panels) and high (BBR, VES; two right panels) P supply from this pool in the soil profile. The vertical dashed line represents the mean value of the theoretical isotopic equilibrium between phosphate and water calculated according to Chang and Blake (2015). The dotted lines represent the possible range of the isotopic equilibrium due to variations in δ18ΟW and temperature. Error bars indicate the standard deviation of the measured δ18ΟP values

The δ18OP values of HCl-extractable P in the top 60 cm of all soil profiles were above or within the range expected for isotopic equilibrium mediated by intracellular pyrophosphatases (Fig. 1). However, below 60 cm soil depth, we observed strong differences between sites. The δ18OP values continuously decreased until a soil depth of about 2 m at the high P site BBR and at the intermediate site VES, whereas at intermediate site MIT and low P site CON, the δ18OP values remained within the range of isotopic equilibrium. The latter two sites showed δ18OP values below equilibrium only at 3.4 and 4.7 m soil depth, respectively (Fig. 1).

To elucidate whether the δ18OP values were linked to the presence of roots, we also analyzed rhizosphere samples and compared their δ18OP values with those from bulk soil. However, we did not find any systematic deviation between bulk soil and rhizosphere soil δ18OP values (p > 0.05; paired t test). Also, the 1 M HCl extractable P fraction showed similar Pi concentrations for rhizosphere and bulk soil samples, close to the 1:1 line (R2 = 0.98; Fig. S2 in supplementary material).

Linear combination fitting of P K-edge XANES spectra (all sample spectra are presented in Fig. S3 in supplementary material) revealed a dominance of Al/Fe-bearing P phases (Table 2). Iron-associated P tended to decrease with depth, whereas aluminum-associated P did not show such a trend. Organic P was found only in the uppermost soil horizon (6–28%) with low P site CON having the highest proportion. By contrast, no Po was detected in the uppermost Ah horizon of the high P site BBR, but it occurred in the second (BA2) and third (CBw1) horizon. Apatite was identified below 80 cm depth in BBR and VES soil profiles, whereas it was only found in the deep saprolite below 5 m in CON and in the non-weathered bedrock at 29.5 m depth in MIT.

Table 2 Specifications of P (% contribution to total P) in selected mineral soil horizons and parent material as obtained by linear combination fitting on P K-edge XANES spectra. The R factor is given as a goodness-of-fit criterion

Spatial distribution of elements at the soil-root interface

Along the transects from a root into the soil matrix at the two intermediate sites VES and MIT, the distribution of normalised counts of 31P16O2 correlated with organic matter (as indicated by 12C14N) in the Ah horizon (R2 = 0.47 for MIT and R2 = 0.55 for VES), whereas in the subsoil 31P16O2 was correlated neither with organic matter nor with the mineral matrix (as indicated by 27Al16O and 56Fe16O). Normalised 31P16O2 counts in the topsoil of both sites were higher near the root surface up to a distance of approximately 20 µm (Fig. 2) in MIT, while no such microscale pattern was observed in the topsoil in VES and in subsoil of both sites (Fig. 2).

Fig. 2
figure 2

Scanned transects using NanoSIMS images from root in mineral matrix in (a) top (0–6 cm depth) and (b) subsoil rhizosphere samples (60–80 cm depth) on MIT and VES sites. The normalised counts of secondary ions were smoothed with a moving-average-algorithm over a distance of 10 µm

Microbial biomass and nutrient limitation of microbial growth

Microbial biomass was highest in the topsoil and decreased with increasing soil depth at all sites. The low P site CON showed substantially larger microbial biomass at all soil depths compared to the other three sites (Fig. S4 in Supplementary Material).

Microbial activity as indicated by substrate-induced respiration was largest for the topsoils (Fig. 3). However, microbial respiration of the topsoil hardly responded to the additional supply of P: the O2 demand measured after addition of C + P was similar to that measured in response to C addition alone. This was different for the subsoils of MIT and CON, which responded stronger to C + P addition than to C addition alone. The addition of C + N increased microbial respiration in all depths at the high P site BBR, with the greatest increase in the deeper subsoil, indicating limited N availability. The other three sites were hardly affected by the addition of C + N, in fact at the intermediate sites MIT and VES microbial respiration even decreased upon addition of C + N. The addition of C + N + P triggered an increase in microbial respiration at all sites and in all soil depths (except for some depths in MIT); the increase in microbial respiration after C + N + P addition was even greater than the sum of C + N and C + P effects. Again, deeper subsoil horizons showed the strongest response (Fig. 3).

Fig. 3
figure 3

Cumulative microbial respiration (measured as oxygen consumption) up to the first maximum respiration rate over 72 h of incubation. Data are given as relative increase (in %) of respiration after C + N, C + P, C + N + P addition in relation to the treatment where only C was added

Discussion

In a general characterisation of P cycling strategies at the four sites studied here, Lang et al. (2017) observed that total P stocks in the soil profile up to 1 m depth—including the organic layers of the forest floor—decreased among the sites in the order BBR > MIT > VES > CON, although total subsoil P stocks in MIT were slightly lower than in VES. Accordingly, in our study of P cycling in the subsoil, we found that the intermediate site MIT mostly presented similar patterns as the P-poor site CON, whereas the other intermediate site VES was more similar to the P-rich site BBR. This suggests that high P concentration in the forest floor and mineral topsoil does not directly affect processes of P cycling in the subsoil; instead subsoil P cycling is related more to the level of P supply at the respective depth.

Concentrations of P in the 1 M HCl extract, which is considered as a proxy for Ca-bound phosphate (e.g. apatite), increased with increasing soil depth at BBR and VES while they remained constant at relatively low levels at MIT and CON until two meters depth (Fig. 1). This depth distribution matches the results from XANES spectroscopy indicating that apatite was only present in substantial amounts in depths which likely represent saprolite or bedrock. In principle, regolith which developed on felsic bedrock is about 20-times thicker than regolith which developed on its mafic counterpart (Bazilevskaya et al. 2013). As the soils of this study developed on bedrock types with varying amounts of mafic compounds, the regolith thickness is expected to increase in the order BBR (mafic bedrock, basalt) < VES (intermediate bedrock, trachyandesite) < CON (felsic bedrock, paragneiss) < MIT (felsic bedrock, paragneiss). We can thus assume that saprolite and bedrock are located at entirely different depths among the study sites: BBR and VES should be characterized by a rather thin layer of saprolite with bedrock occurring within 2 m from the soil surface, whereas weathering in MIT and CON should have proceeded to much greater depths. Based on these general patterns for the distribution of P stocks within the soil and bedrock profile, we propose that the low subsoil P stocks at MIT and CON would force P cycling to proceed to deeper soil depths as opposed to the sites with higher subsoil P stocks in BBR and VES.

Indications of P cycling can be obtained from the oxygen isotope ratio of inorganic phosphate (δ18OP) in the 1 M HCl extract, even though the HCl-extractable P is usually not considered to be part of the immediate biological P cycle in soils. Nevertheless, as the HCl pool includes secondary minerals precipitated after biological cycling, the δ18OP values of this pool provide an integrated signal of P cycling over longer periods of time (Angert et al. 2012; Joshi et al. 2016). For our sites, δ18OP values in the topsoil of BBR and VES were above or within the range of isotope values imparted by biological cycling. The δ18ΟP values above the calculated equilibrium in the upper profile depths can be explained by the evaporative enrichment of the soil water at surface layers (Allison et al. 1983; Barnes and Allison 1983; Hsieh et al. 1998; Gazis and Feng 2004). When using soil water δ18O values measured at the sites of − 3 to − 5‰ (F. Tamburini, unpublished data), the resulting equilibrium value is at 19–22‰, which is in the range of topsoil δ18OP values observed here. It is therefore likely that deviations from our expected equilibrium in the upper part of the soil profile can be explained by evaporative enrichment of the soil water, which should, however, not affect the soil profile below a depth of ~ 10–20 cm (e.g., Hsieh et al. 1998; Gazis and Feng 2004).

With increasing soil depth and increasing 1 M HCl-Pi concentrations in BBR and VES, the δ18ΟP values continuously decreased until they reached values expected for primary minerals in parent material at depths below 60 cm (Fig. 1), thus indicating decreasing P cycling intensity. These findings are in line with studies on agricultural systems with high supply of P. There, the δ18OP values in the surface soil were also close to equilibrium values and decreased in the deeper subsoil to values expected for parent material, where P apparently had neither been accessed by roots nor by soil microorganisms (Amelung et al. 2015; Bauke et al. 2018). However, we also cannot exclude that biological cycling of P occurred in the subsoil, but that either biologically cycled P remained in the more labile fractions, which were removed by NaHCO3 and NaOH extraction before the HCl extraction step, or that any changes in δ18ΟP values were not great enough to impart measurable changes compared to the large reservoir of unrecycled P.

For sites MIT and CON, by contrast, the concentrations of Pi extracted with 1 M HCl remained low until a depth of two meters and δ18OP values remained close to the range of isotopic equilibrium, which indicates that a major proportion of P in this pool had been biologically cycled. This biological P cycling extended to depths below 4 m. Only, at depths deeper than 13 m δ18OP values reached a value of around 6‰, similar to values found in the subsoil at sites BBR and VES (Fig. 1). These δ18OP values correspond to values reported in the literature for parent material containing apatite (~ 10‰; Tamburini et al. 2012; Amelung et al. 2015). Besides, these values coincide with the detection of apatite by XANES (Table 2). We conclude that the general trend in δ18OP values with depth was in line with the transition from saprolite to weathered bedrock and finally to non weathered parent material as also indicated for the same sites by the concentration of refractory elements such as zirconium, which are inert to chemical weathering (Uhlig and von Blanckenburg 2019). Similarly, recent data from Uhlig et al. (2020) indicate that mineral nutrients are originally sourced from 2 to >10 m depth at CON and MIT. Taken together, these results indicate that at the sites with low subsoil P stocks, biological cycling of P proceeded to several meters more in depth than at sites with higher P stocks.

While the δ18ΟP values indicate the biological cycling of P in the subsoils and even in saprolite of the temperate forests in our study, it is not clear whether biological cycling occurred through the activity of roots or microorganisms or both. We did not analyse the distribution of fine roots in the subsoils of these sites. However, the vertical root morphology is influenced by site conditions and effective rooting depth depends on several factors such as annual precipitation, soil order, texture, and chemical conditions as well as tree species and their age (Strong and La Roi 1983), all of which were relatively similar across all study sites (Table 1). Further, trees on soils with greater profile depth can be expected to produce greater fine root biomass in the subsoil than trees on shallow soil profiles (Kirfel et al. 2019). We therefore assume that tree roots at the deeply weathered sites MIT and CON would explore the subsoil to a greater extent than roots in BBR and VES and that δ18ΟP values in the deeper subsoil thus have been partly influenced by the activity of tree roots.

A comparison of δ18ΟP values in bulk soil and rhizosphere samples did not reveal any differences in P cycling, which is probably due to the low spatial resolution when sampling the rhizosphere soil as soil attached to roots after gentle shaking. Therefore, we analysed element spatial patterns at the µm-scale around roots by NanoSIMS. Here, we found that in the topsoil of both VES and MIT the normalised counts of both 31P16O2 and 12C14N secondary ions along the line scans decreased with increasing distance from the root into soil matrix (Fig. 2). We observed that P was enriched rather than depleted nearby the root, and a majority of P was correlated with the accrual of organic matter. Similar gradients for P concentrations in the rhizosphere have been observed before and could even extend to the mm-scale (Richardson et al. 2009). Such patterns may be explained, on the one hand, when assuming that roots grow into old root channels rather than forming new pathways into soils, as, e.g., observed for arable systems (Athmann et al. 2013). On the other hand, P might be enriched in the rhizosphere due to the mobilisation of P from organic and mineral phases by ectomycorrhizal fungi (EMF) (Richardson et al. 2009; Smits et al. 2012; Spohn et al. 2018). For the subsoil samples, however, this relationship was evident neither at VES (high subsoil P stocks) nor at MIT (low subsoil P stocks). Here, the 31P16O2 counts did not correlate with those of organic matter or with those of oxides and no P accumulation was found along subsoil roots (Fig. 2). However, considering the limited number of roots analysed here, the specific contribution of roots to P cycling in deeper soil horizons should be addressed in more detail in future studies.

Subsoil P cycling could be additionally driven by saprotrophic as well as symbiotic microbes, which may be associated with the roots, as is the case for EMF. The activity of enzymes exuded by soil microorganisms should result in a similar exchange of oxygen isotopes in phosphates during P cycling as by enzymes exuded directly by the roots. Indeed, we found microbial biomass in the deeper subsoil (Fig. S4 in Supplementary Material), but in lower quantity than in topsoils. After the addition of nutrients, microbial activity responded strongly in the subsoils of all sites (Fig. 3). Especially in the lower subsoil, a strong response to nutrient addition might be the result of a spatial separation of microbes and resources in the undisturbed soil profile (Preusser et al. 2017), which was overcome by the addition of nutrients in our incubation experiment. The strong response to C + P addition in the subsoil of CON and MIT indicated P limitation of microbial growth, while microbial activity in the subsoil of VES and BBR showed the strongest response to C + N or C + N + P addition, indicating N limitation for these sites. When comparing nutrient limitations in the subsoils with δ18ΟP values a clear pattern emerges: at sites with high P stocks, especially in BBR, where microbial activity appears to be limited by the availability of N but not P, δ18ΟP values indicate that mineral phosphate has not been biologically cycled in the subsoils. In contrast, microbial activity in subsoils of sites with low P stocks appear to be mainly P limited and the according δ18ΟP values indicated biological cycling of P as is needed for a more efficient utilisation of the scarce P resources at these sites. The findings therefore suggest increasing relevance of microbial processes especially at sites with lower subsoil P stocks.

Conclusions

Our results of δ18OP and XANES analyses indicate that biological cycling of subsoil P resources occurs down to several meters soil depth at sites with low P stocks. In this regard, our data support our main hypothesis that P deficiency in the surface soil not only fosters microbial cycling of P in the organic and upper mineral soil layers but also causes the utilization of P from the deeper subsoil. Our data also indicate that despite limitations of subsoil microbial activity by the availability of P and N, the biological cycling of P in subsoils is primary due to microbial activity and less directly linked to the effect of roots. Our results do not discount the possibility that there is also intensive P cycling in the subsoil at the sites with higher subsoil P stocks, but due to larger reservoir inherited from the parent material these effects are not detectable in our study yet.