Influence of Soil Properties on P Pools and Its Effect on Forest Productivity in Mediterranean Calcareous Soils

: Although soil phosphorus is essential for vegetation growth it is not always found in labile forms, hampering its absorption by plants, and is limited in forest ecosystems. This study explores soil P state and availability in calcareous soils, determining which processes affect its different pools and which soil parameters influence forest productivity of a Mediterranean pine species. We used a sequential fractionation method to determine P fractions in the soil according to their lability and their organic or inorganic nature. Those fractions were related to different soil and climatic parameters and to the site index in 32 Pinus halepensis plots of the National Spanish Forest Inventory. Soil texture, carbonates content, water retention capacity and organic matter substantially affected P fractions in the calcareous soils studied. Membrane extractable P was mainly influenced by organic matter-related parameters while the predominant P fraction in these soils, the primary P, was linked to the carbonates content. The biological mineralization processes played a key role in the soil in the studied plots.


Introduction
Phosphorus is essential for life as it is involved in many vital processes in organisms [1]. Soil available P is usually low in terrestrial ecosystems [2] and is considered one of the most limiting nutrients for plant growth and thus, for forest productivity ( [3,4] among others). Although there is abundant total P in the soil, it is not always found in available forms due to its low solubility, strongly dependent on the pH of the soil, its slow diffusion, and its high fixation in soils, hindering its absorption by plants.
During decomposition of plant residues, microorganisms release P via mineralization, and capture P via immobilization, controlling the transformations of this element between organic and inorganic forms and forming a biological subcycle in the general cycle of soil P [5,6]. On the other hand, precipitation-dissolution and adsorption-desorption processes control the abiotic transfer of P between the solid phase and the soil solution [5,7].
Organic matter is key for soil structure and for the regulation of nutrients quantity and cycles [8], thus playing a major role in the availability of P in forest soils. Organic matter is a direct source of P organic forms that are transformed into available P forms for vegetation through the mineralization and the release of microbial biomass [3,9]. Soil microbial activities produce different enzymes, including phosphatase, which can substantially affect the transformation of P organic forms into inorganic or more labile and soluble P forms.
Rainfall and temperature highly affect the mineralization and immobilization rates in the soil, as well as soil moisture which is also key for phosphatases production and for microorganisms' survival [10][11][12]. The mineralization-immobilization processes and the different P pools are also influenced by soil pH in contrasting ways [11,[13][14][15]. Additionally, the parent material, the pedogenesis and the resulting soil type (order), as well as soil texture, affect the fractionation, transformations and distribution of P [11,16,17].
Numerous studies have indicated that current climate change leads to an increase in N depositions and an increase in summer droughts causing nutritional imbalances and P deficiencies in ecosystems, making these changes more pronounced in the Mediterranean areas [18][19][20]. For this reason, and although general processes involved in the biochemistry and availability of P are well known, it is necessary to clarify which soil parameters mostly affect each process both in specific ecosystems, mainly in Mediterranean ones, and in different types of soil [21].
Pinus halepensis Mill. has a wide natural distribution throughout the Mediterranean region, thus inhabiting a wide variety of bioclimates [22]. Due to its high tolerance to drought stress and its low nutrient demand, this species has been extensively used in afforestation targeting ecological restoration and soil erosion reduction.
Some studies found that forest productivity of P. halepensis in natural stands and plantations under different environmental conditions was closely related to the availability of water [23,24] and the amount of organic matter or elements such as nitrogen, magnesium, or potassium [25]. However, these studies did not consider the availability and state of soil phosphorus as a factor that can affect forest productivity. On the other hand, P bioavailability is regulated by the interaction between different P fractions [26,27] and progresses over time in such a way that the available P at a given time may be low, but in the long term, all soil P might become labile [28]. Thus, in soils poor in P, plants can obtain P in more recalcitrant forms through various strategies [15,29] such as through the exudation of organic anions and phosphatases from the roots that increase desorption, solubilization, or mineralization of P from sources of low availability in the soil, among other mechanisms [30][31][32]. Consequently, determining the dynamics and bioavailability of P requires the separation and identification of the different pools of P in the soil [26,27]. For this, we used the sequential fractionation method of P from [33] modified by [34]. This method is widely used for characterizing forms of P differing in their availability for plants and microorganisms and participating in short and long-term transformations in the soil [1,17,27]. Therefore, the procedure enables distinguishing between soil Porg associated with humic and organic compounds that are easily mineralized by microbes and soil Ping retained by chemisorption [35].
Our objectives were (1) to evaluate soil P status and availability in calcareous soils under Mediterranean pine plantations, (2) to assess which soil processes are affecting the different soil P pools and which are driving P availability, and (3) to determine which soil parameters are influencing forest productivity of the studied stands and whether soil P availability is leading this productivity. We hypothesized that (1) there is a relationship between the membrane extractable P in the studied soils and other labile or non-labile P fractions, (2) soil parameters related to organic matter and carbonates content are the most determining factors for P availability in the calcareous soils studied, and (3) forest productivity in the P. halepensis plantations studied is mainly related to P availability and organic matter-related parameters.

Study Area
This study was carried out with soils sampled in 32 plots belonging to the Spanish National Forest Inventory (SNFI) and located in pure P. halepensis plantations. These stands are in the Castilla y León region (northern half of Spain; Table S1). The main characteristics of these plots are presented in Table 1. Lithology was composed of limestone and marl. Soils showed basic pH (from 8 to 8.9) with clay-loam texture and were classified as Calcixerepts within the Inceptisol order.

Soil Sampling
The first 10 cm of mineral topsoil were sampled in the 32 P. halepensis plots [36] following the method of [37]. Only those first 10 cm of soil were sampled as this layer is the fertile soil and the environmental changes are more strongly reflected there. At each plot, four points located 5 m from the center of the plot were sampled (in N, S, E, and W directions). This sampling method allows obtaining representative samples which reflect the soil variability. At each of these points, one disturbed soil mineral sample was taken. The four disturbed samples from each plot were lumped together to attain a composite mineral soil sample per plot (for more details on soil sampling see [36]. Disturbed soil mineral samples were air-dried and sieved (2 mm) before laboratory analyses (done in duplicate).

Phosphorus Fractionation
Phosphorus sequential fractionation was carried out following the method developed by [33] and modified by [34]. This fractionation was designed to remove progressively less available P with each subsequent soil extraction. It is useful to separate plantavailable or labile forms of P from several refractory P pools in soils ( Figure 1). Furthermore, with this procedure, the organic and inorganic P fractions can be differentiated, which is key to understanding the P cycle in forest soils. The first fraction extracted with anion exchange membranes (PAEM) corresponds to the readily exchangeable inorganic P and that which is easily dissolved from solid phases in the soil [34,38]. Then, sequential extractions were carried out on the residue of this first extraction with progressively stronger extractants. The following extraction was performed with 0.5 M NaHCO3, to determine the highly labile P adsorbed in the soil exchange complex, which can be easily transformed into available forms [39,40]. Using basic extractants such as NaHCO3 or NaOH, used in the next step of the fractionation, implies the extraction of organic and inorganic forms of P. The moderately labile P strongly retained by chemisorption [40], was determined by the extraction with 0.1 M NaOH. This P fraction is thought to be associated with the surface of amorphous and some crystalline Al and Fe minerals and is probably available in the mid-term. Aliquots of NaHCO3 and NaOH extracts were acidified to precipitate extracted organic matter, and inorganic P was determined in each of the supernatants. Total P of these extracts was determined after digestion with sulfuric acid and persulfate in an autoclave at 120 °C [41]. Organic P in each extract was calculated as the difference between total and inorganic P. Moreover, for determining the highly and moderately labile P fractions, active carbon was needed to remove the humic substances before determining the inorganic P of each fraction. The inorganic primary P (PHCl1M) was extracted with 1 M HCl. This fraction has low solubility and might be formed by the inorganic P included in Ca primary minerals [7,40]. Finally, the procedure extracts the most recalcitrant fractions which are the stable and the residual P. These fractions are derived from non-alkali-extractable organic debris [34]. On one hand, the stable P (PHClconc) was obtained by the extraction with concentrated HCl at 90 °C in a water bath and corresponds with P bound in the interior of Fe and Al minerals and apatite [34]. On the other hand, the residual P (PHClO4) was extracted with HClO4 and digested at 230 °C and represents the most stable form of P that is available only in the long term [42]. For P determination in the extracts obtained, the method by [43] modified by [44] was applied, using an ultraviolet/visible spectrophotometer Thermo Genesys 20. Analytical determinations were performed in duplicate on all soil samples. Figure 1. Scheme of the sequential phosphorus extraction procedure used in this study (modified from [43]). AEM: anion exchange membrane; Pi: inorganic P; Po: organic P; HClconc: concentrated HCl. Organic P fractions are indicated with green color.

Soil Parameters
The studied soils have been characterized by the analysis of a high amount of physical, chemical, and biochemical soil parameters (see [12,25,45], for more information on the analyses). We extracted soil parameter data for this study from two databases in [46,47]. We discarded some soil parameters because they were obtained as ratios of other selected parameters, thus providing redundant information, such as the ratio of total organic carbon to total nitrogen or the available soil water (the quotient of field capacity and permanent wilting point). We also discarded the sand and silt content according to international criteria (ISSS) since we used USDA criteria because the particle size considered in this second criterion has greater ecological significance. Considering that the high number of soil parameters selected (32 parameters) and the number of plots studied (n = 32) can greatly hinder the statistical study, the soil parameters were clustered into five groups related (1) to soil physical fertility, (2) chemical fertility, (3) soil organic matter, (4) soil composition and (5) soil enzymatic activity ( Table 2 and Table S2). Table 2. Soil parameters used in this study and clustered into five groups (extracted from [46,47]

Forest Productivity
Forest productivity was estimated using the site index (dominant height at a reference age) which is strongly correlated with wood production [48]. We used the site index estimated by [25] for the same Pinus halepensis plots studied here.
Although our objective was to study the effect of soil P availability and status on forest productivity of the studied stands, we also included all the soil parameters studied in this work and some climatic variables, since they are usually important in explaining site indexes. Climatic variables were selected based on previous studies of forest productivity of this species and other similar species [24,36] indicating as climatic variables that most affect forest productivity those related to water availability. For this reason, the Lang index [49], and the annual water index following Thornthwaite [50] have been included in this analysis (data obtained from [47]). Both indexes classify the area as arid/sub-humid.

Statistical Analysis
To determine the relationships between P fractions, Spearman's correlation coefficients (achieved with the Hmisc package in R software; [51]) were calculated, and principal component analysis (PCA) was performed (in R: prcomp). To evaluate the relationships between P fractions and groups of soil parameters we used Spearman's correlation coefficients and canonical correlation analysis (CCA), performed using CCA and CCP packages in R [52,53]. Therefore, canonical correlation analyses were used to explore the joint association between P fractions and soil parameters related to (a) physical fertility, (b) chemical fertility, (c) soil composition, (d) organic matter and (e) soil enzymatic activity. Kaiser's Measure of Sampling Adequacy (in R: KMO; psych package; [54]) was calculated to test whether the sample size was adequate or not for CCA [50]. This test varies from 0 to 1 and indicates the degree to which each variable in a set is predicted without error by the other variables. Lastly, we used general partial least squares regression (PLS; in R: plsr from the pls package; [55]) to study the relationship between site index, P fractions, soil parameters and climatic variables. PLS is a multivariate statistical technique that combines principal component analysis (PCA) and multiple linear regression which allows us to explain both the variability of the matrix of independent variables and the vector of the dependent variable. PLS is particularly useful when the number of predictor variables is greater than the number of observations and when there are many highly related explanatory variables [56]. Before performing all the analyses, the studied variables were tested for linearity, outliers, and normality (in R: shapiro.test), and those variables not normally distributed were transformed (Table S2). Statistical analyses were performed with R software [57].

Soil Phosphorus Fractionation
The studied soils showed a mean total P content of 289 mg/kg, varying between 116 and 611 mg/kg. Organic P represented around 12% of the Ptotal ( Table 3). As expected, the primary P fraction was the biggest fraction in the soils studied followed by the stable P fraction ( Table 3). The primary P includes P forms in primary Ca minerals, accounting for 47% of total P. Stable P accounted for more than 27% of total P. The smallest P fraction in these soils was the membrane extractable P or readily exchangeable P fraction that accounted for 1.2% of total P. In the principal component analysis of P fractions, three principal components accounted for 74% of total variance. The first principal component was mainly positively correlated to the organic highly labile P fraction, the second principal component was positively correlated to the stable P fraction and the third principal component was positively correlated to the primary forms of P Figure 2 and Table S3). According to the Spearman's correlations among all P fractions, the membrane extractable P was positively correlated to the inorganic highly labile P (Spearman rho = 0.47, p = 0.007) and to the organic moderately labile P (rho = 0.54, p = 0.001). The organic highly labile fraction was significantly correlated with the organic moderately labile P (and with the stable P fraction (rho = 0.39; p = 0.027). The highest positive significant correlations were found between the inorganic moderately labile P and the primary P forms (rho = 0.60, p < 0.000) and between the stable and the residual P forms (rho = 0.69, p < 0.000). Besides, total P was only significantly related to the inorganic moderately labile P, the primary P and the stable P (rho = 0.44, p = 0.012, rho = 0.82, p < 0.000 and rho = 0.49, p = 0.005, respectively).

P Fractions and Soil Composition Parameters
In the canonical correlation analysis performed on P fractions in relation to soil composition parameters, we firstly included six parameters (see Table 2), but the Kaiser's Measure of Sampling Adequacy (MSA) was 0.40 and thus, this analysis was not considered statistically acceptable. Therefore, we decided to eliminate one parameter related to soil texture (i.e., clay, silt, or sand content). Clay content was highly correlated to silt content (r = −0.66, p < 0.000), while it was not significantly correlated to sand content (r = 0.06, p = 0.732). Silt and sand contents were also highly significantly correlated (r = −0.74, p < 0.001). Thus, we decided to keep in the analysis the clay percentage and we checked the individual MSA for silt content and the individual MSA for sand content trying to discard one of these two variables. As the MSA of sand was higher than that of silt (0.39 and 0.38, respectively) we discarded silt content in the analysis. Clay and sand contents give good information about soil texture. The Kaiser's Measure of Sampling Adequacy without silt content (MSA = 0.60) indicated that the analysis was statistically acceptable. So, this CCA was performed using five soil parameters related to soil composition: percentage of clay, percentage of sand according to USDA criteria, carbonates, active carbonates, and gypsum content. Two canonical dimensions were significant (Wilk's lambda; p < 0.001 and p = 0.002, apiece) and presented correlations of 0.87 and 0.81, respectively. The first canonical dimension was mainly influenced by the sand content (r = 0.97) and the primary P (r = −0.91). The second canonical dimension was negatively correlated to carbonates (r = −0.77), active carbonates content (r = -0.69) and the inorganic highly labile P (r = −0.54) and was positively correlated to the residual P (r = 0.54, Figure 3 and Table S4).  Table S2.
Clay content was positively correlated with the residual P but negatively with the membrane extractable P (Table 4). Sand content was highly negatively correlated with the primary P and the inorganic moderately labile P, and carbonates content was negatively correlated with residual P (Table 4). Table 4. Spearman's correlations between P fractions and soil composition parameters ("-" means not significant). Significance: *** p < 0.001, ** p < 0.01, * p< 0.05.

P Fractions and Soil Physical Fertility Parameters
The first canonical dimension in the canonical correlation analysis of P fractions and physical fertility-related parameters was significant according to Wilk's lambda test (p < 0.001) with a canonical correlation of 0.89. Kaiser's Measure of Sampling Adequacy (MSA) was 0.63 and, therefore, this analysis was considered statistically acceptable. The first canonical dimension was negatively correlated to field capacity, permanent wilting point, soil porosity, and to all P fractions, and it was just positively related to bulk density (Figure 4 and Table S5).  Table S2.
Strong positive correlations between primary P forms and field capacity and permanent wilting point were found (Table 5). The inorganic moderately labile P fraction and the membrane extractable P were also positively correlated with the previously referred soil parameters. Porosity only showed significant correlations with the inorganic highly labile P fraction and was the parameter least correlated with the first canonical root in the CCA. Bulk density was significantly and negatively related to the membrane extractable P and to the inorganic highly labile fraction, respectively. Table 5. Spearman's correlations between P fractions and soil physical fertility-related parameters ("-" means not significant). Significance: *** p < 0.001, ** p < 0.01, * p< 0.05.

P Fractions and Soil Chemical Fertility Parameters
The canonical correlation analysis (CCA) between the P fractions and soil chemical fertility-related parameters produced two canonical dimensions that captured significant relationships between the two data sets (r = 0.95, p < 0.001 and r = 0.92, p = 0.011, respectively). Kaiser's Measure of Sampling Adequacy (MSA) was 0.61 and, therefore, this analysis was considered statistically acceptable. For the chemical fertility-related parameters, the first canonical dimension was most strongly influenced by the reciprocal cation exchange capacity (1/CEC; r = 0.81) and the second canonical dimension by the available manganese (r = −0.72), followed by the reciprocal of the available iron (1/Fe; r = 0.62) and the exchangeable magnesium (r = 0.61). For the P fractions, the first dimension was mainly correlated to the membrane extractable P (r = −0.63) and to the organic highly labile P (r = −0.59), and for the second dimension, the primary forms of P were dominant (r = 0.87; Figure 5 and Table S6).  Table S2.
According to the Spearman's correlations between P fractions and soil parameters related to chemical fertility (Table 6), the most labile P forms were positively correlated with the cation exchange capacity and with the amount of available manganese. Moreover, the readily exchangeable P extracted with membranes was also positively correlated with the exchangeable potassium. Organic highly labile P fraction was also positively correlated with the amount of available zinc and iron, as was the residual P. Organic moderately labile P was also positively related to the cation exchange capacity, available zinc and exchangeable potassium. Inorganic moderately labile P was negatively correlated with the available manganese and positively with the exchangeable magnesium, as the primary P fraction that was also positively related to pH and negatively to available iron. Table 6. Spearman's correlations between P fractions and soil chemical fertility-related parameters ("-" means not significant). Significance: *** p < 0.001, ** p < 0.01, * p< 0.05.

P Fractions and Organic Matter Parameters
A biplot of the first two dimensions of the CCA analysis based on P fractions and soil organic matter variables is shown in Figure 6. The explained correlations are 0.92 (p < 0.001) and 0.81 (p = 0.041) for the first and the second canonical correlation, respectively. Kaiser's Measure of Sampling Adequacy (MSA) was 0.78 and thus, this analysis was considered statistically acceptable. The first canonical dimension ( Figure 6 and Table S7) was negatively correlated with all the organic matter-related parameters considered, as well as with all P fractions except to the inorganic moderately labile P and with the primary P, two fractions strongly related to each other. The second canonical dimension was mainly negatively influenced by the easily oxidizable C and the microbial biomass N, and positively by the most recalcitrant P forms (stable and residual P).  Table S2.
We found strong correlations between the membrane exchangeable P and all organic matter-related parameters (Table 7). These correlations were all positive. The stable and residual P forms were just positively correlated to total nitrogen, but the primary P form did not show any significant correlation. Table 7. Spearman's correlations between P fractions and soil parameters related to organic matter ("-" means not significant). Significance: *** p < 0.001, ** p < 0.01, * p< 0.05.

P Fractions and Soil Enzymatic Activity
The first canonical dimension in the canonical correlation analysis based on P fractions and soil enzymatic activity-related variables was significant (Wilk's lambda test; p < 0.001) and its canonical correlation was 0.93. Kaiser's Measure of Sampling Adequacy (MSA) was 0.65 and, therefore, this analysis was considered statistically acceptable. The canonical dimension was negatively influenced by all the soil enzymatic activity parameters ( Figure 7 and Table S8). Besides, most of the phosphorus fractions were negatively correlated with the first canonical dimension except the inorganic moderately labile P and the primary P.  Table S2.
As expected, strong significant correlations were found between soil parameters related to enzyme activity and the most labile and plant-available phosphorus fractions (Table 8). The primary, stable and residual P fractions did not show significant correlations with these soil parameters, as well as the inorganic moderately labile P fraction (highly related to the primary P). Table 8. Spearman's correlations between P fractions and soil enzymatic activity parameters ("-" means not significant). Significance: *** p < 0.001, ** p < 0.01, * p< 0.05.

Relationship between Forest Productivity, Soil P Fractions and Soil Parameters
In order to study the importance of the soil parameters considered in the productivity of P. halepensis stands a PLS analysis was performed. The PLS analysis included the site index as the response variable and all the soil parameters, all the phosphorus fractions and two climatic variables related to water availability as predictive variables. The PLS resulted in the selection of two components that explained 67.5% of the variation of the site index (Table S9). These results showed that the first component was mostly related to the composition of the soil (carbonates, active carbonates and soil texture explained 87.8% of the component variance). Soil pH and the permanent wilting point (PWP) were also important variables in this component as well as the primary P fraction (Table S10). Component 1 only used 11.4% of information from predictors and was able to predict 35.9% of site index variability. The second component of the PLS was associated with variables related to soil water (and soil chemical fertility (25.2% of the component variance; Table  S10). Although climatic variables did not explain a high percentage of the variation in any of the two components, they did substantially improve the site index variation explained by the model (67.5% compared to 47.7% without the climatic variables), remaining also the results regarding to the most influential variables in each component (data not shown).

Soil Phosphorus Status and Availability
Membrane extractable P was the lowest fraction in the studied soils. As expected, we found positive correlations between membrane extractable P and highly labile inorganic P since this fraction includes easily interchangeable P compounds and simulates the action of plant roots in dissolving P minerals [39,40]. Thus, these results indicated that the highly labile inorganic P fraction in the studied soils could act as a source of available P in the short term. However, and unlike other studies (such as [16]), we did not find significant correlations between the most available P and highly labile organic P, even when this fraction is considered a very active fraction of soil organic P [58]. Organic P concentration calculated as the sum of PoNaHCO3 and PoNaOH, represented around 12% of the total P, following previous results in other studies in inceptisols or in Mediterranean forest soils [42,59]. The most refractory forms of P (PHClconc and PHClO4) were positively correlated with the highly labile organic P fraction, corroborating that the most stable and residual fractions may even be sources of available P in the long term [60].
Primary and stable fractions were the predominant forms of P in the soils studied as observed in other calcareous soils [1,61], among others). Primary P was just correlated with the moderately labile inorganic P and with total P, indicating that in these calcareous soils total P is mostly formed by the P included in Ca minerals. The relationship between primary and moderately labile inorganic P deserves an in-depth study to understand its significance in P dynamics in limestone soils with high carbonates content. Few studies have indicated that 1M HCl extraction removes not only Ca-P but also Fe and Al-related P forms [62][63][64], which could explain this correlation. In fact, Fe, Al, and Mn are sometimes more significant in controlling P solubility in calcareous soils than free lime [65]. However, this possible Ca-P miss quantification by sequential chemical extraction methods is soilspecific (depending on soil properties) and several studies have shown that in calcareous soils the Ca-P and organic P fractions extracted with sequential methods are comparable and consistent with those extracted using other novel methods such as synchrotron X-ray absorption near-edge structure (XANES) spectroscopy (59,63], among others). Total P concentration in the studied soils was low, but in the range found on forest soils [1,17,27] and was also alike to that found in limestone soils of northwestern Spain [61] and in inceptisoltype soils [42].
The principal component analysis allowed us to differentiate three pools of P according to their lability (see, for example, [60,66]). A first pool related to the first axis of the PCA corresponding to both organic and inorganic labile fractions, that is, the short-term available P fractions. The second PCA was associated with the P forms extracted in the last steps of fractionation, which are the most stable forms and available in the very long term. Finally, in the third axis, the medium-term-available forms appeared, which represent a greater proportion of the total P in these soils and are mainly inorganic forms of P extracted with NaOH and 1M HCl.

Relationship between P Fractions and Soil Parameters
Specific characteristics of each soil determine the solubility and bioavailability of nutrients, and soil properties affect differently to the diverse soil P forms ( [17,59], among others). Therefore, soil processes involved in the availability and stability of soil P will vary in soils with unlike characteristics and in different ecosystems. Our results showed that the organic fractions and the P available forms might be affected by organic matter and microbiological activity, while the P fractions linked to Ca seem to be influenced by the composition of the soil, mainly by its sand and carbonates content.
The highly significant correlations of the parameters related to the composition of the soil were established with the most stable fractions of P. Although important, the direct effect of texture on the distribution of P pools in forest soils has been little studied [17,67]. Strong negative correlations were found between sand content and primary P, a fraction that represents the total P content in the studied soils. The negative influence of sand content on primary P and total P concentrations has been previously demonstrated in soils from very different climatic areas around the world [17,35]. This influence could be interpreted as the reduction of surfaces to which phosphorus can bind, such as the surfaces of the clay-humic complex and, in the specific case of limestone soils, the lower content of active limestone, a material with the size of silt and clay which shows great chemical activity in the soil and is often a better predictor of P behavior than the soil total lime content [65].
Soil aeration and water retention capacity are important factors affecting P dynamics because they play important roles in the mineralization and immobilization of this element [11]. However, the significance of these factors in the P cycle remains uncertain. In the present study the effect of soil aeration in P fractions has been considered through bulk density and porosity. Our results indicated that labile P fractions were reduced with increased soil bulk density, i.e., with reduced soil aeration. Lower soil aeration might imply lower soil microbial activities and subsequently reduced mineralization processes. However, contradictory results have been reported in this matter, with higher mineralization rates seen both in good aeration conditions and in bad conditions [10]. The reason why bulk density greatly affects the soil P available forms may be, on the one hand, to the content of organic matter, since the higher the organic matter, the lower the bulk density. On the other hand, it may be due to aspects related to the textural class.
As expected, parameters related to soil water retention were positively correlated with the available and labile P, since soil moisture is essential for soil microorganisms' survival and for phosphatases production [11,68], thus increasing the labile P fractions. However, whether the soil has more or less water depends not only on its retention capacity but also, and more importantly, on the amount of water that reaches the soil, that is, on rainfall. In this study we did not consider climatic or physiographic variables, we cannot determine either the availability of water or the actual aeration level of the soils studied. Therefore, that the studied soils show high field capacity does not imply that they have enough water for the microorganisms to carry out their activities and consequently, increase the amount of labile P in the soil. In fact, the studied plots are under a Mediterranean climate where annual precipitation is low (less than 500 mm) and the summer drought is intense and long-lasting. Consequently, in the studied plots as in other forest systems, the scarcity of precipitation and low soil moisture, together with a low water retention capacity, can greatly reduce the enzymatic activities in the soil, and so reducing the mineralization of P [45,68].
Phosphorus mineralization rates as well as its solubility, fixation, and concentration of its different pools are influenced by soil pH ( [14,15], among others). However, unlike previous studies, our results did not show a strong effect of pH on P pools, which is probably due to the low pH variation among the studied soils (range 8.0-8.9).
The different phosphorus fractions, mainly labile P, were positively related to cation exchange capacity and the amount of exchangeable potassium, which could be explained by the strong relationship of these parameters with the quantity and quality of soil organic matter. The short-term availability of P to vegetation is strongly influenced by biochemical processes that affect organic matter, while its long-term status is generally determined by geochemical transformations [26,27]. In addition, primary P was negatively related to assimilable iron and manganese since when carbonates or active carbonates content in the soil increase, the amount of assimilable Fe and Mn decreases, while primary P bounded in Ca-minerals increases [66].
As expected, the most labile fractions of P both inorganic (PAEM and PiNaHCO3) and organic P (PoNaHCO3) were highly correlated with soil parameters related to organic matter and enzymatic activity. This indicates that in the soils studied the biological cycle of P has great importance in the availability of this element. The availability of P from organic forms may be due to two possible causes: (1) The first is the general mineralization of soil organic matter which is regulated by the energy demand of microorganisms [69].
Results from CCA and Spearman's correlations completed with the group of organic matter-related parameters showed the relevance of this biological P mineralization in the studied soils, following previous results [45]. (2) The second possible cause of the release of P from organic forms could be the specific release of P due to the action of particular enzymes involved (acid phosphatase and, especially in the studied soils, alkaline phosphatase; [70]) which is governed by the specific need for available P. In general, it is considered that the lower the availability of P in the soil, the higher the activity of phosphatase (either acid or alkaline, depending on the soil pH), since this enzyme is released when available P levels are limiting. In the calcareous soils studied, alkaline phosphatase did not show any correlation with assimilable forms of phosphorus. Moreover, it neither showed high relevance in the first canonical axis of the CCA performed with the group of enzyme activity-related parameters. The absence of importance in the CCA and of correlation between the available P forms and alkaline phosphatase could explain that in these calcareous soils the concentration of assimilable P, although low, may be sufficient to satisfy the demand of the vegetation and soil microorganisms, unlike what was found in acid soils in northern Spain [46]. The rest of the enzymes considered in this study (without phosphatase activity) showed highly significant correlations with the labile forms of P, both in the Spearman correlations and in the CCA. The fluorescein diacetate hydrolysis reaction (FDA), which showed a higher weight in the corresponding canonical dimension, reflects all the hydrolytic activity in the soil [71] and is commonly used as an indicator of the general microbial activity in the soil [72]. Furthermore, all the enzymes were related to each other and to the microbiological parameters Cmic, Cmin, Nmic. Microbial populations, besides participating in the mineralization and solubilization of soil P, are also a reserve of organic P that is released when microorganisms die [3,9]. Our results showed that the microbial biomass of C, N and P were significantly correlated with the most labile phosphorus and the microbial biomass C also correlated with organic P. The microbial biomass P in our study represented 2.7% of the total P. The positive correlation found between microbial biomass P and membrane extractable P is consonant with previous studies [58,59] and with the fact that microbial biomass regulates soil P availability through the release of its P content. Therefore, all these results might indicate that the biological mineralization of P in these soils plays a more important role in P liberation than the biochemical mineralization.
Although in the studied soils, organic forms of P represented a low percentage of total P compared to other studies [17,73], a strong influence of total soil carbon (TOC) and easily oxidizable carbon (EOC) on the distribution of P could be expected. Following previous studies, we found a positive correlation between organic P, available P and highly labile inorganic P with soil organic C ( [59,67], among others), since these soluble and labile forms of P come from the mineralization of organic matter [21]. The positive relationship between inorganic P and easily oxidizable carbon may be because the organic matter could act as a source of labile inorganic P in the studied soils [74]. In addition, CCA results showed that the first canonical factor was highly correlated with TOC and the second canonical factor with EOC, indicating that these two parameters might be good predictors of labile P.

Relationship between Forest Productivity, Soil P Fractions and Soil Parameters
Using globally the large set of soil parameters considered in this work (including the different forms of soil P) we explained around 67% of the variance of the site index calculated by [25] for the P. halepensis stands studied. Our results indicated that site index in the studied plots was conditioned by parameters related to soil texture and carbonate content as well as to the primary P, which is the P pool that mostly represents the total P in these plots. Soil water retention-related parameters also explained much of the site index variance in the studied plots. These results are in line with those obtained by other authors for P. halepensis since this species growth (associated with its site index and productivity) is mainly driven by soil water availability considering both its supply and its retention in the soil [23,24]. Surprisingly, climatic variables included in the model which were related to water supply did not influence the site index as much as have been expected. This may be since in these plots, soil composition parameters, mainly the high presence of Ca, largely modulates trees growth. However, the inclusion of climatic variables did considerably improve site index variance explanation. In addition, soil chemical fertility related to the supply of macronutrients was also important in explaining site index. This may be due to high availability of macronutrients as K and Mg could favor resistance to drought since trees with a good supply of these nutrients reduce water loss through transpiration [24,75]. The availability of micronutrients did not affect the forest productivity of these plots while in agricultural soils this factor is one of the most important for crops' growth [76]. Generally, organic matter decomposition affects the productivity of ecosystems, particularly in nutrient-poor forests such as Mediterranean forests [77]. However, the parameters related to organic matter did not seem to influence the site index of the studied plots. This could be due to the release of nutrients from plant litter depends not only on the activity of the microorganisms but also and mainly, on the litter quality and on environmental factors. In our case, the quality of the litter is the same in all the plots: same species of Pinus with the same chemical composition in needles, same strategies for the conservation of nutrients in their tissues, and same decomposition rates [78]. Furthermore, the plots are in areas with very similar environmental conditions so environmental factors would not make great differences in the nutrient release process between plots either.

Conclusions
In the calcareous soils studied, biological mineralization processes played a key role in the P cycle, and therefore, the parameters related to soil organic matter and enzymatic activity mainly promoted the availability of that element. However, the main fraction in these soils (primary P) was linked to the carbonates content and could act as a source of available P in the medium term. The P available concentration, although low, may be sufficient to satisfy the demands of the vegetation and microorganisms in these soils. The site index of Pinus halepensis in the studied plots was not affected by the soil available P, indicating that this species growth and thus its productivity, is not so limited by soil nutrients and other soil parameters, but by water availability, normally scarce throughout its circum-Mediterranean distribution. Advancing in the knowledge about the internal dynamics of forest soils and their relationship with forest productivity is key to providing outstanding information to sustainable forest management.
Supplementary Materials: The following are available online at www.mdpi.com/article/10.3390/f12101398/s1, Table S1: Location of the 32 Pinus halepensis plots studied (UTM Projection in meters; Datum ED50) and soil textural class at each plot. Table S2: Description, units and transformations of phosphorus fractions and the soil parameters used in this study and clustered into five groups (soil parameters extracted from [46,47]). Variables' transformations were done when necessary to achieve normality. Double-line separates P fractions from soil parameters and dashedline separates each soil parameters group. Table S3: Results from the Principal Component Analysis applied to phosphorus fractions data from 32 Pinus halepensis plots studied. Variables with loadings > |0.70| in bold case. Table S4: Factor loadings (simple correlations between each canonical factor and the variables) of canonical correlation analysis performed on P fractions in relation to soil composition parameters. Variables with higher importance in each component in bold case. Table S5: Factor loadings (simple correlations between each canonical factor and the variables) of canonical correlation analysis performed on P fractions in relation to soil physical fertility parameters. Variables with higher importance in the first component in bold case. Table S6: Factor loadings (simple correlations between each canonical factor and the variables) of canonical correlation analysis performed on P fractions in relation to soil chemical fertility parameters. Variables with higher importance in each component in bold case. Table S7: Factor loadings (simple correlations between each canonical factor and the variables) of canonical correlation analysis performed on P fractions in relation to soil organic matter parameters. Variables with higher importance in each component in bold case. Table S8: Factor loadings (simple correlations between each canonical factor and the variables) of canonical correlation analysis performed on P fractions in relation to soil enzymatic activity parameters. Variables with higher importance in the first component in bold case. Table S9: Variance explained by each component and accumulated variance of the matrix of independent variables (X) and the vector of dependent variable (Y). Selected components in bold case.  Funding: This research was funded by Junta de Castilla y León, project number VA096G19 "Ecosystem services of soils under mixed forest stands versus pure stands. Effect of stand type on soil fertility, water retention and carbon sequestration". R.C.M.-S. was supported by the Youth Guarantee EU initiative (job contract).

Institutional Review Board Statement: Not Applicable.
Informed Consent Statement: Not Applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.