Calibration of Near-Infrared Spectra for Phosphorus Fractions in Grassland Soils on the Tibetan Plateau

Abstract

Soil phosphorus (P) is essential for plant growth and influences biological processes. Determining the amounts of available P to plants has been challenging, and many different approaches exist. The traditional Hedley sequential extraction method and its subsequent modification are applied to determine different soil P forms, which is critical for understanding its dynamics and availability. However, quantifying organic and inorganic P (Po & Pi) in different extracts is labor-intensive and rarely used with large sample numbers. As an alternative, near-infrared spectroscopy (NIRS) has been employed to determine different P fractions at reasonable costs in a short time. This study aimed to test whether the analysis of P fractions with NIRS is an appropriate method to disentangle the effects of P limitation on high-altitude grassland ecosystems, particularly with fertilizer amendments. We explored NIRS in soils from the grassland soil samples on the northern Tibetan Plateau. First, we extracted the P fractions of 191 samples from the Haibei Alpine Meadow Ecosystem Research Station at four depth increments (0–10 cm, 10–20 cm, 20–40 cm, and 40–70 cm), including nutrient additions of nitrogen (N) and P. We compared the results of the Hedley extraction with the laboratory-based NIRS model. The fractionation data were correlated with the corresponding NIRS soil spectra; the coefficient of determination (R²) of the NIRS calibrations to predict P in P fractions ranged between 0.12 and 0.90; the ratio of (standard error of) prediction to the standard deviation (RPD) ranged between 1.07 and 3.21; the ratio of performance to inter-quartile distance (RPIQ) ranged from 0.3 to 4.3; and the model prediction quality was higher for Po than Pi fractions, and decreased with fertilizer amendment. However, the external-validation results were not precise enough for the labile P fractions (RPD < 1.4) due to the limited number of samples. The results indicate that using NIRS to predict the more stable P pools, combined with Hedley fractionation focusing on the labile P pool, can be a promising approach for soils in alpine grasslands on the Tibetan Plateau.

1. Introduction

Phosphorus (P) is essential for all life forms, especially plant growth [1,2,3]. It is one of the plant’s essential nutrients that influences biological processes, and a deficiency of P, when combined with nitrogen (N), can impede plant growth in terrestrial environments around the world [4,5]. P in soil occurs in two distinct pools: inorganic P (Pi, minerals, and mineral-bound P) and organic P (Po, organic matter, and organic-bound P) [6]. Forms of soil P range from ions in solution to very stable inorganic and organic compounds, which typically include: primary P minerals (apatite, strengite, variscite); secondary P minerals (Ca−, Mg−, Fe−, or Al-phosphates); P adsorbed on the edges of clay minerals; mineral- and organic-associated P; and dissolved P in the form of H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻ [2]. The last two parts are the primary forms that can be assimilated by plants and microbes [6]. Soil P plays a vital role in determining ecosystems’ structures, functions, and processes [5,7,8]. It is essential to distinguish between Pi and Po, since Po, although not directly available for plant uptake, plays a vital role in plant nutrition in forest ecosystems [9,10]. In alpine environments, soil P and its bioavailability also influence the establishment of timber lines [11]. In addition, there is increasing interest in the role of soil P and other nutrients, such as soil N and carbon (C), in alpine grassland ecosystems that are vulnerable to climate warming and nutrient enrichment [12,13,14,15].

In cold and high-altitude regions such as the Tibetan Plateau—which covers an area of more than 2.4 million km² and has an average altitude exceeding 4000 m above sea level [16], with young soils that are at the beginning of their formation [17,18]—weathering coincides with the onset of primary succession; Pi is slowly released from parent material as the primary source [19]. Concerns have been expressed that the ecosystems and plant productivity on the Tibetan Plateau may suffer from P and N limitations [20,21]. Related research shows that nutrient application in the meadow could increase primary plant productivity. It suppresses topsoil microbial activity [22], especially for P enrichment with N; meanwhile, P addition more efficiently increases the N uptake capacity of grasses, since N is an essential limiting nutrient for plant productivity [23]. On the other hand, P is essential for the organization of alpine meadow soil fungal communities; however, the addition of P may reduce soil fungal diversity and long-term ecosystem stability in fragile alpine environments [13]. Furthermore, factors such as overgrazing and fertilization of grassland; the foundation of settlement; and road and railway construction interfere with the fragile ecosystem of the Tibetan Plateau [24,25,26,27]. When the temperature increases in the Tibetan grassland, total N, total P, and available P in the soil decrease, and the coupling between available N and P is lost [28]. This effect has also been shown in high-altitude tree-line systems south of the Tibetan Plateau in the Himalayas [11]. In contrast, Huang et al. [29] reported that the distribution of alpine grasslands on the north Tibetan Plateau did not show an upward shift despite rapid climate warming having occurred from 2000 to 2014. In summary, it would be valuable to explore the role of P in the response of the meadow ecosystem on the Tibetan Plateau to climate change and nutrient enrichment; moreover, the determination of soil P fractions that reflect its bioavailability is also an essential basis for evaluating the influence of soil P on the alpine grassland ecosystem.

Kruse et al. [30] reviewed advanced methods currently used in soil P research. Nevertheless, for a long time, extraction methods such as the Olsen bicarbonate method have been developed to analyze the different P fractions [31]. These methods focus only on the soil labile-P but do not quantify Pi and Po into fractions of different plant availability. The sequential extraction of P after Hedley et al. [32], which was further modified by Tiessen and Moir [33], has become standard in ecosystem research. This method provides seven Pi and four Po fractions. These P fractions are usually grouped into pools of distinct plant availability as follows [9,34]:

(a)

Labile P, which is considered to supply the short-term P demand of plants, including Resin-P, NaHCO₃-Pi and Po;

(b)

Moderately labile P, which can be transformed into labile P forms, including NaOH-Pi and Po;

(c)

Stable P, including HCl-P and residual-P, which hardly contributes to bioavailable P.

However, compared to other methods, Hedley fractionation has a high workload that requires technical expertise, and high costs in the laboratory, making it time-consuming, error-prone, and unsuitable for routine soil analysis [35,36]. Therefore, it would be helpful to determine Hedley P extraction by utilizing a less expensive and time-saving technique. A promising approach is the measurement of near-infrared spectra (NIRS) as proxies for P fractions in the soil fractions with reasonable accuracy [9,37,38,39,40]. In contrast to traditional soil chemical, biological and physical methods, NIRS allows the measurement of large sample numbers in a relatively short time and does not require complicated sample pretreatment [41]. However, there have been no studies focusing on the use of NIRS to predict organic and inorganic Hedley P fractions in meadow soils on the Tibetan Plateau.

The specific objectives of this study are (i) to adapt the Hedley method modified by Tiessen and Moir [33] to the specifics of the grassland soils of the Tibetan Plateau, and (ii) to examine whether P fractions determined by the Hedley method can be predicted from the spectral properties of soils determined by NIRS as applied by Niederberger et al. [9], which serves as a reference method.

2. Materials and Methods

2.1. Soil Samples

Soil samples were provided by the Northwest Institute of Plateau Biology, Chinese Academy of Sciences, and were taken from the long-term field experiment at the Haibei Alpine Meadow Ecosystem Research Station in the northeastern part of the Tibetan Plateau (37°37′ N, 101°12′ E). The mean annual temperature is 1.1 °C, and the mean annual precipitation is 485 mm at 3.200 m a.s.l [42]. Thirty-six plots (6 m × 6 m) using the randomized block design were set up with 6 repeats and 6 treatments. Fertilizers (urea and triple superphosphate,

Table 1. Design of the nutrient addition experiment on the alpine grassland since 2012 (N: nitrogen, P: phosphorous).

Treatments	Fertilizer	Amount (ha−1 year−1)
Control	No	0
P	Triple Superphosphate (TSP)	50 kg
NP	Carbamide CO(NH2)2 + TSP	50 kg P + 100 kg N
N25	Carbamide CO(NH2)2	25 kg
N50	Carbamide CO(NH2)2	50 kg
N100	Carbamide CO(NH2)2	100 kg

) were evenly distributed by hand on the soil surface in the plots, after sunset (for higher moisture), during the growing season in early June, July, and August, each year since 2011 [43]. The buffer strip between each block was 2 m wide, and that between each plot was 1 m wide. Each fertilization treatment was divided into two parts: with and without grazing by yak. The soil samples were taken from four depth increments (0–10, 10–20, 20–40, and 40–70 cm) in 2016 from 4 repeats. All samples were air-dried and sieved (<2 mm), and 191 samples (one N25 treatment sample was missing) in total were analyzed.

2.2. Hedley Fractionation

Soil samples were air-dried, sieved to a fine fraction (<2 mm), and extracted in different fractions, according to Hedley et al. [32], Alt et al. [44], and Niederberger et al. [9]. An amount of 0.5 g from each sample was weighed into a 50 mL centrifuge tube for extraction. The fractionation steps are given in

Table 2. Experimental procedure for the Hedley fractionation steps and assignment of Hedley P fractions to soil P pools.

P Pools and Fractions		Extraction Procedure [9,43]	Properties and Bonding Forms of Pi and Po in the Fractions [45]
Labile P	Resin-P	Anion-exchange resin in resin bag, 0.5 M HCl	Mainly Pi, marginal Po; biologically most available P form; adsorbed on the surface of crystalline compounds.
	NaHCO3-P	0.5 M NaHCO3, pH 8.5	Highly labile P; Pi likely to be plant-available, associated with Fe and Al oxides; Po easily mineralized.
Moderate P	NaOH-P	0.1 M NaOH	Moderately labile P; Pi associated with Fe and Al oxides; Po involved in slow transformation processes.
Stable P	HClconc.-P	HClconc. 85 °C	Very stable Pi; covers P in primary minerals; Po in very stable pools, eventually also derived from particulate organic matter.
	Residual-P	0.5 M H2SO4	Highly resistant and occluded P forms.

. For the regular Hedley fractions, the amount of organically bound P in the 1 M HCl extracts was considered negligible [33] and, therefore, not analyzed. The concentrated HCl (HCl_conc.) fraction was applied.

All Pi fractions were determined by continuous flow analysis (CFA, SEAL Auto Analyzer AA3, SEAL Analytical GmbH, Norderstedt, Germany). Total NaHCO₃-P and NaOH-P were determined by inductively coupled plasma–optical emission spectrometry (ICP–OES, Optima 5300 DV, PerkinElmer, Waltham, MA, USA). To evaluate differences in P fractions between the treatments and soil depth increments, one-way ANOVA and independent-sample t-tests were used.

2.3. Near-Infrared Spectroscopy

All soil samples were scanned with an integrating sphere measured by diffuse reflectance using a Fourier transform near-infrared reflectance spectrometer (Tensor 37, with a resolution better than 0.6 cm⁻¹, a wavenumber accuracy better than 0.01 cm⁻¹ at 2000 cm⁻¹ and a photometric accuracy better than 0.1%; Bruker Optik GmbH, Ettlingen, Germany). The sample was inserted into a glass sample cup with a flat bottom, placed on the scanner, and measured as it rotated. Each spectrum consisted of 64 independent scans. Five replicate measurements were made for each sample, and a mean spectrum was created over the entire spectrum from 12,000 to 3800 cm⁻¹ wavenumbers, i.e., 800–2500 nm wavelength with a resolution of 16 cm⁻¹. The spectral data and model development were processed with the spectroscopic software OPUS/QUANT (Version 7.5, Bruker Optik GmbH, Ettlingen, Germany, 2014). The data pre-processing of a combination of the first derivative and vector normalization (SNV) in the software was used to find the best fit with the value of the respective component (Bruker Optik User Manual version 6, 2006). Corrected partial least square (PLS) regressions then proceeded to regress the information contained in the spectrum of the calibration sample with the P fraction values, which is a popular chemometrics method in NIR analysis [9].

Calibration was performed with leave-one-out cross-validation for the small data sets to determine the accuracy of the NIRS models [46]: one sample was removed from the data set and validated against the remaining subset. This process was repeated until every sample was used once for validation. For each data set, 80% of the samples were used to build the models with cross-validation, and 20% were set aside for external validation to test the prediction [9]; these were chosen by an automatic function in the software (Bruker Optik User Manual, version 6, 2006). Different treatments of N and P were considered in the calibration, but not in the external-validation due to the limit of the samples. No outliers were removed in the model development process, since new outliers were detected after removing the old ones. The differences in performance of the model were compared using the following parameters: coefficient of determination (R²) (Equation (1)); root mean standard error (RMSE) (Equation (2)); the ratio of (standard error of) prediction to the standard deviation (RPD) (Equation (3)); and the ratio of performance to inter-quartile distance (RPIQ) (Equation (4)):

$${\mathrm{R}}^2=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\hat{\mathrm{y}}}_{\mathrm{i}}-{\overline{\mathrm{y}}}_{\mathrm{i}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-{\overline{\mathrm{y}}}_{\mathrm{i}}\right)}^2}$$

(1)

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\hat{\mathrm{y}}}_{\mathrm{i}}-{\mathrm{y}}_{\mathrm{i}}\right)}^2}$$

(2)

$$\mathrm{RPD}=\frac{\mathrm{standard} \mathrm{deviation} \mathrm{of} \mathrm{analyzed} \mathrm{data}}{\mathrm{RMSE}}$$

(3)

$$\mathrm{RPIQ}=\frac{\left(\mathrm{Q}3-\mathrm{Q}1\right)}{\mathrm{RMSE}}$$

(4)

in which $\hat{\mathrm{y}}$ = predicted value; $\overline{\mathrm{y}}$ = mean measured values; $\mathrm{y}$ = measured values; $\mathrm{n}$ = number of samples with $\mathrm{i}$ = 1, 2, …, n; Q1, Q3 are the upper bounds of the first and third quartiles of the measurements, respectively; and $\mathrm{Q}3-\mathrm{Q}1$ represents the interval responsible for 50% of the population around the median [47]. The range of RPIQ might represent the soil attributes’ values better than the standard deviation, since the distribution of soil attributes is usually skewed [47,48]. However, the RPD value has been adopted extensively when using NIRS, and several limits have been suggested and applied [9,40,49]. Hence, we included RPD as one of the evaluation criteria, and bias—the mean of residuals—did not occur very often [47]. Generally, a robust model achieves high R², RPD and RPIQ values and a low RMSE [9,48,50,51]. The performance of the models can be judged as follows:

RPD: (i) ≥2.0, suitable for predicting P in soils; (ii) 1.4–2.0, limited; and (iii) ≤1.4, not reliable [40];

RPIQ: (i) >2.5, excellent; (ii) 2.0–2.5, good model with predictive ability; (iii) 1.7–2.0, good model; (iv) 1.4–1.7, fair model that needs improvement; and (v) <1.4, very poor predictive ability [51].

All other statistical analyses were conducted using Origin v.2021 and Excel 2019.

3. Results

3.1. Hedley P Fractions

No significant differences in grazing treatment were found in Hedley fractionation, so we kept the data volume support for modeling, and the final data of P fractions were averaged. The average total P concentration by Hedley fractionation of all the samples covered a wide range from 381 to 950 µg g⁻¹, and significant differences in total P between treatments were found after N and P amendment. In comparison, each P fraction ranged from 0 to 460.5 µg g⁻¹ (Table S1). The complete P fractionation results at different soil depths are displayed in

Table 3. Soil P contents (µg g⁻¹) from each soil depth for seven Hedley fractionations in different treatments with N (N25: 25 kg nitrogen, N50: 50 kg nitrogen, N100: 100 kg nitrogen ha⁻¹ year⁻¹), P (50 kg phosphorous), and their combination NP (100 kg nitrogen + 50 kg phosphorous).

Fractions	Labile P (µg g−1)			Moderate P (µg g−1)		Stable P (µg g−1)
Treatments	Resin-P	NaHCO3-Pi	NaHCO3-Po	NaOH-Pi	NaOH-Po	HClconc.-P	Residual-P
0–10 cm
Control	20.9 (1.1)	19.0 (1.9)	21.4 (2.1)	28.1 (2.4)	177.5 (2.5)	219.0 (6.8)	86.8 (3.6)
P	183.8 (22.3) ***	112.4 (8.5) ***	68.5 (6.7) ***	82.1 (7.5) ***	181.2 (4.4)	332.3 (14.7) ***	80.5 (1.9) *
NP	139.5 (12.7) ***	96.4 (5.4) ***	49.4 (7.7) **	84.0 (3.6) ***	196.2 (5.0) **	326.1 (19.3) ***	83.5 (3.8)
N25	19.1 (1.7)	23.4 (2.4)	24.0 (3.4)	28.1 (2.0)	178.6 (3.7)	252.9 (18.9)	85.0 (2.9)
N50	14.8 (3.0)	23.2 (1.8)	25.5 (3.0)	27.1 (2.3)	172.8 (4.5)	230.2 (7.4)	79.3 (3.3)
N100	13.6 (1.9) **	20.5 (2.3)	21.4 (2.6)	27.6 (2.3)	179.4 (8.4)	230.8 (8.0)	83.6 (3.0)
10–20 cm
Control	11.1 (1.1)	6.6 (0.7)	14.4 (1.2)	11.5 (0.5)	142.8 (4.1)	203.2 (5.7)	81.0 (3.9)
P	31.2 (5.7) **	24.0 (5.3) *	24.2 (4.3)	20.9 (3.2) *	150.0 (7.7)	257.3 (16.2) *	75.8 (2.5)
NP	29.3 (2.6) ***	21.4 (2.1) ***	20.0 (2.3) *	24.8 (1.9) ***	166.3 (5.8) **	222.8 (7.1) *	72.8 (6.8)
N25	12.0 (1.4)	9.0 (1.0)	16.5 (1.5)	12.4 (1.6)	158.7 (5.7) *	232.7 (14.9)	74.7 (1.7) *
N50	10.8 (1.2)	9.3 (0.7) **	15.0 (1.4)	12.1 (0.8)	143.2 (4.8)	227.5 (8.3) *	75.5 (1.8) *
N100	10.5 (0.8)	8.3 (0.8)	14.0 (0.9)	12.7 (1.1)	153.4 (6.0)	215.6 (7.8)	74.7 (2.1)
20–40 cm
Control	8.6 (1.1)	2.7 (0.4)	8.8 (0.9)	4.6 (0.3)	96.2 (4.2)	273.8 (10.2)	67.1 (2.3)
P	15.5 (2.3) *	10.0 (2.0) **	14.1 (1.6) *	7.2 (1.2)	102.2 (8.1)	298.5 (23.1)	68.6 (3.5)
NP	13.7 (0.9) ***	8.4 (0.7) ***	12.9 (1.4) *	7.0 (0.6) **	106.3 (7.4)	269.4 (12.6)	68.0 (1.8)
N25	7.6 (1.0)	4.2 (0.5) *	10.0 (1.4)	4.7 (0.7)	92.7 (9.4)	288.6 (20.4)	68.3 (2.4)
N50	9.0 (1.1)	4.7 (0.9)	10.9 (1.4)	5.1 (0.5)	93.9 (7.4)	274.9 (16.4)	69.9 (3.4)
N100	8.2 (1.4)	3.7 (0.3) *	9.2 (0.9)	5.5 (0.9)	93.3 (7.2)	283.5 (12.0)	65.3 (2.6)
40–70 cm
Control	7.9 (0.9)	1.4 (0.3)	3.0 (1.1)	3.1 (0.3)	28.9 (8.1)	408.8 (20.3)	43.7 (3.5)
P	9.7 (0.8)	5.3 (0.7) ***	4.8 (0.8)	3.9 (0.4)	34.4 (5.2)	403.8 (11.0)	48.2 (2.7)
NP	9.4 (0.6)	3.7 (0.7) *	4.7 (1.3)	3.5 (0.6)	57.0 (20.6)	395.4 (15.1)	51.4 (1.9) **
N25	6.0 (0.9)	2.6 (0.5)	3.7 (0.7)	2.7 (0.3)	36.8 (14.0)	398.3 (19.6)	48.4 (3.0)
N50	6.4 (1.0)	1.6 (0.4)	2.0 (0.8)	2.9 (0.3)	19.6 (3.7)	437.7 (6.6) **	39.5 (1.7) *
N100	6.5 (1.1)	1.6 (0.3)	3.4 (0.8)	3.0 (0.3)	28.6 (3.7)	420.0 (11.6)	45.3 (2.6)

Asterisks denote significant differences among the treatments at p = 0.05; * < 0.1; ** < 0.01; or *** < 0.001. Values are the mean ± SE, df = 7; treatment with N25 at 10–20 cm, df = 6.

, separately summarized for P pools (labile P, moderate P, and stable P). The highest P concentrations were found in the stable P pool for each depth increment (Table 3). Different reactions of P fractions between treatments and soil depths were detected. In the treatments with P application (P and NP additions), Pi contents (Resin-P, NaHCO₃-P, NaOH-P, and HCl_conc.-P) at 0–10 cm and 10–20 cm exceeded those without P input (Table 3). However, the difference between means was more distinct in those fractions in which the 0 to 10 cm depth had higher p values than the other layers. Notably, the contents of Po fractions (NaHCO₃-Po and NaOH-Po) only increased significantly with the addition of the NP combination. The residual fraction was the only part in which no highly significant differences (p < 0.1) between treatments were found.

The proportion of P in each fraction varied distinctly from different soil depth increments, highlighting some effects more clearly (Figure 1 and Table S1). Similarly, changes occurred mainly in the topsoil (Figure 1a,b). P addition, after several years, increased the proportion of Pi at 0–20 cm. The proportion of moderate and stable P pools also increased with increasing depth. For the labile P pool, the contribution of Resin-P and NaHCO₃-Pi was significantly increased by P addition (6.3% to 10.2% compared with 1.9% to 2.7%), which led to no change in NaHCO₃-Po. NaOH-Pi in the moderate P pool also increased with P addition (4.8–4.9% to 2.6%), but the contents of NaOH-Po were stable (Figure 1 and Table S1). Intriguingly, the contribution of Po to total Hedley-P was higher in the treatments without P addition. Overall, stable P formed the central part of the total P concentration for all nutrient amendments and soil depth increments; meanwhile, phosphate addition significantly changed labile and moderate Pi’s contributions.

3.2. NIRS Models

The quality measures of the NIRS models varied to a great extent over the cross-validation (calibration) data set with different treatments (

Table 4. Quality parameters for NIRS model calibration in the cross-validation process for all Hedley fractions with different treatments.

P Pools		Labile P			Moderate P		Stable P
P Fractions and Treatments		Resin-P	NaHCO3-Pi	NaHCO3-Po	NaOH-Pi	NaOH-Po	HClconc.-P	Residual-P
Control	R2	0.52	0.82	0.67	0.86	0.95	0.79	0.67
	RPD	1.44	2.33	1.73	1.99	4.43	2.24	1.74
	Bias	0.137	0.023	0.344	0.605	0.517	8.15	−0.417
	RMSECV	4.09	3.21	4.47	3.97	12.9	39.7	10.9
	RPIQ	1.5	3.1	2.1	3.4	5.9	2.7	1.9
P	R2	0.08	0.42	0.43	0.53	0.86	0.06	0.66
	RPD	1.04	1.32	1.33	1.46	2.74	0.99	1.71
	Bias	−2.3	−0.164	0.096	0.0252	−3.94	8.15	−0.268
	RMSECV	74.9	34.6	20.2	23	21.5	70.5	8.36
	RPIQ	0.9	1.3	1.1	1.5	3.5	0.7	2.0
NP	R2	0.59	0.62	0.37	0.74	0.75	0.40	0.26
	RPD	1.57	1.63	1.26	1.97	2.04	1.29	1.17
	Bias	0.33	0.397	0.528	0.414	0.865	−0.989	−0.668
	RMSECV	35.7	23.5	16	16.6	30.4	57.8	13.7
	RPIQ	2.1	2.5	1.6	2.9	3.0	1.0	1.4
N25	R2	0.48	0.81	0.54	0.87	0.91	0.70	0.80
	RPD	1.39	2.28	1.48	2.75	3.35	1.83	2.22
	Bias	−0.045	−0.15	0.3	−0.347	−4.45	1.49	−0.179
	RMSECV	4.43	4	6.3	3.9	18.8	43.9	6.83
	RPIQ	1.2	2.1	1.9	3.3	3.8	1.5	2.4
N50	R2	0.02	0.65	0.57	0.77	0.91	0.82	0.69
	RPD	0.92	1.69	1.53	2.1	3.35	2.36	1.8
	Bias	0.422	0.209	0.294	0.893	2.02	5.24	0.241
	RMSECV	6.1	5.2	6.35	4.87	17.9	38.7	9.58
	RPIQ	0.6	1.7	1.4	2.1	3.7	2.6	2.2
N100	R2	0.10	0.76	0.70	0.80	0.90	0.71	0.72
	RPD	1.06	2.03	1.82	2.26	3.25	1.85	1.88
	Bias	0.424	0.0426	0.15	0.203	−1.23	−0.862	0.723
	RMSECV	4.21	3.93	4.26	4.54	18.7	45.9	8.44
	RPIQ	1.3	2.8	2.3	3.0	4.0	2.1	2.3

R²: coefficient of determination; RPD: ratio of (standard error of) prediction to the standard deviation; RMSECV: root mean squared error of cross-validation; RPIQ: ratio of performance to inter-quartile range.

). N and P amendments significantly reduced the performance of the models, especially the addition of P, which had higher RMSECV (RMSE of cross-validation) but lower RPD and RPIQ compared with the control. Among all the treatments with N and P, better calibration models achieved for NaOH-Po than others such as Resin-P (R² = 0.02–0.59; RPD = 0.92–1.57; RPIQ = 0.6–2.1; RMSECV = 4.21–74.9), NaHCO₃-Pi (R² = 0.42–0.81; RPD = 1.32–2.28; RPIQ = 1.3–2.8; RMSECV = 3.93–34.6) and NaHCO₃-Po (R² = 0.37–0.70; RPD = 1.26–1.82; RPIQ = 1.1–2.3; RMSECV = 4.26–20.2). In particular, when nutrient amendments were not considered, the fraction of NaOH-Po (R² = 0.90, RPD = 3.21, RPIQ = 4.3), HCl_conc.-P (R² = 0.79; RPD = 2.19; RPIQ = 2.4), and Residual-P (R² = 0.78; RPD = 2.16; RPIQ = 2.6) in the regression model performed better in calibriation compared with other fractions (

Table 5. Quality parameters for NIRS model calibration in the cross-validation process for all Hedley fractions.

P Pools		Labile P			Moderate P		Stable P
P Fractions		Resin-P	NaHCO3-Pi	NaHCO3-Po	NaOH-Pi	NaOH-Po	HClconc.-P	Residual-P
Cross-validation	R2	0.12	0.31	0.44	0.48	0.90	0.79	0.78
	RPD	1.07	1.2	1.34	1.39	3.21	2.19	2.16
	Bias	0.0314	−0.00023	0.0233	0.0975	0.0348	−0.611	0.0418
	RMSECV	42.1	22.1	11.3	15.2	18.8	38	7.51
	RPIQ	0.3	0.8	1.2	2.3	4.3	2.4	2.6

R²: coefficient of determination; RPD: ratio of (standard error of) prediction to the standard deviation; RMSECV: root mean squared error of cross-validation; RPIQ: ratio of performance to inter-quartile range.

). Overall, the best calibration models were achieved for the control data without nutrient amendments (Table 4), and especially for the NaOH-Po, HCl_conc.-P and Residual-P fractions. In addition, the RPD values were higher for cross-validation than for external validation for all seven Hedley fractions—confirming that external validation requires a larger dataset (Table S2)—which ranged from 1.16–2.29, and RMSEP (RMSE of prediction) ranged from 9.97–59.8. According to the judgment of performance [40,51], all the P fractions, except Resin-P from the non-NP addition samples (control), showed better predictive ability in cross-validation, as well as NaOH-Po in all treatments (Table 4) and in external-validation (Table S2). Notably, the RMSECV of stable P fractions was generally higher than others.

We compared the results of the treatments of P addition, NP addition and the Control set after grouping the Hedley fractions into different P pools. Better predictions with R² were shown when no P inputs were considered (Figure 2). The moderate and stable P fractions could effectively be predicted (R² > 0.86), while only the result of the labile P pool in the control was acceptable, albeit barely (Figure 2g, R² = 0.72). The Po and Pi fractions from labile and moderate P pools were also compared to assess the model. For the Po fractions, R² was higher than for the Pi fractions (Figure S1). The prediction accuracy decreased when the P and NP additions were applied.

4. Discussion

4.1. Influence of Nutrient Amendments on Soil P Pools

Research on nutrient addition to different natural ecosystems is becoming a key issue worldwide [52,53,54]. In tropical forest ecosystems with highly weathered soils, atmospheric P deposition can lead to P fertilization effects through direct input. It can accelerate P release from organic matter and increase P availability [55]. It has been reported that in subtropical forest ecosystems, N addition significantly increased the content of available soil P, but decreased the content of moderately labile P [56]. However, fertilization and plant species composition on the Tibetan Plateau could play a substantial role in P pools in grassland soils. In alpine grasslands, the relative availability of soil N and soil P has been indicated as an important determinant in regulating plants’ N and P contents, controlling plant metabolic and growth rates [57]. As shown for some legumes, N additions could also significantly increase the aboveground net primary production (ANPP) [43]. The addition of P and N during the experiment at the Haibei station also significantly reduced fungal species richness and changed the fungal community composition [13]. Overall, N addition can affect the P fractions in three ways: First, N addition to ammonium compounds can lower the soil pH and promote the dissolution of recalcitrant Ca-bound Pi into moderately occluded Pi and available P fractions [58]; second, N also promotes the secretion of microbial extracellular phosphatase, which enhances the mineralization of Po, thereby increasing the leaching of available P [59]. In addition, plants and mycorrhizal communities are also the driving force behind soil P cycling [60]. Our results partly verify that Po contents were increased by the second way, with NP addition at a 0–20 cm depth (Table 3), indicating the P-limit on the research area. In addition, the P contents of Pi and Po increased significantly under the P additions, but there were no notable changes in the proportions of Po fractions (Table 3). Plant-species composition must be considered to estimate plant access to soil P, as Pätzhold et al. [45] showed in long-term grassland experiments in Germany. P and N additions also significantly influenced the labile and moderate P fractions (Figure 1 & Table 3). Regarding the proportion of each fraction, our study showed values of 8.7% for plant-available labile P pools, and 31% for the moderate P pool (28.4% for Po in this pool) for the control samples (Table S1). In addition, most of the soil P was present in stable P forms (all fractions made up almost 65% of the total P, Table S1), which can be attributed to the formation of the relatively young soil on the Tibetan Plateau [17]. This pool also differed at least between depths. However, significant differences were observed between the P ratios of all the pools compared to different soil depths.

The distinct differences between all depths of the stable P pool could probably be explained by the large C content and soil development difference. In soils with less organic matter, the Pi fractions can have equal importance in controlling Resin-P [61]. Furthermore, strong correlations between the NaHCO₃ and NaOH fractions have been found, suggesting that these fractions do not represent separate pools [33]. The contribution of individual fractions to total P in soils varies substantially between grassland, pasture, arable land, and forest. Significant differences were observed mainly in moderately labile P and stable P fractions. As shown in Figure 1 and the P content for the control samples in our study, the increase in P pools after nutrient amendment is inconsistent with findings from other studies worldwide (Table S3). However, recent research has pointed out different results to determining bioavailable soil P (labile P) using Hedley fractionation. Wu et al. [62] provided better results for labile P fractions using this method than synchrotron-based X-ray absorption near edge structure (XANES) spectroscopy. Hou et al. [63] highlighted the value of soil P dynamics and availability worldwide to predict the response of terrestrial ecosystems to global changes. On the other hand, Negassa & Leinweber [64] reviewed the limit of Hedley fractionation for different land uses. Klotzbücher et al. [65] argued that this method could not provide information to understand P bioavailability, because the sequential extraction using individual steps could not be related to a specific mineral source or its binding strength. In addition, the labile P measured by Hedley fractionation should not be generally defined as plant-available P because most of the P pool is immobilized by microbes [19]. We suppose this difference is mainly related to whether P fertilizer is applied. In summary, in our experiment, the results of Hedley fractionation can provide an indication of the amount of bioavailable P in grassland soils.

4.2. NIRS Models

Chodak [66] and Cécillon et al. [67] emphasized that P calibration rarely performs well in soil, with an R² value between 0.4–0.5 only; this is because P has no direct vibrational response in NIRS. Therefore, its proper prediction depends on the correlation with the soil components absorbed in the spectral range of wavelengths, such as total C and N, particle size, and water content. Our calibration models of individual P fractions achieved relatively satisfactory results for some, but not all, individual P fractions (Table 4 and Table 5). It is obvious that nutrient amendments increased the data redundancy (higher heterogeneity) significantly, resulting in reduced accuracy (Table 4). Generally, the precision and reproducibility of the results of the Hedley fractionation method are of crucial importance for the quality of predictions of NIRS models [9]. All the fractions of the labile P pool (Resin-P, NaHCO₃-P) performed poorly in the cross-validation (Table 4 and Table 5), and fractions of moderate and stable P pools performed better. Richter et al. [68] hypothesized that organic acids slowly displace Pi because many naturally occurring acids, such as citric acid and oxalic acid, compete with phosphates for sorption sites. The poor external-validation results were due to an insufficient sample size and complex external conditions, i.e., heterogeneous sets. The results of NIRS for P determination allow approximate quantitative predictions or the grouping of samples into different categories based on RPD and RPIQ values, at least under our experimental conditions [9,40,49,51]. We attribute the higher RMSE values (Table 4, Table 5 and Table S2) to the significant differences caused in the fractionation and the high concentrations. For many agricultural and engineering applications (such as soil fertility assessments), it is usually sufficient to classify soils as basic test values without requiring an accurate estimate of soil properties. To use NIRS to determine P fractions in grassland-soil-monitoring studies, the precision of the reference values and NIRS measurements has to be improved; the grinding of soil samples is a possibility for the development of improved NIRS measurements. Overall, soil particle sizes influence the spectral absorbance, introducing an error for P measurements [69,70,71]. Additionally, the homogeneity of soil samples, which have lower data redundancy, provides a better quality of P prediction in the NIRS model [9]. Thus, it is necessary to classify the soil samples appropriately.

In addition, we note that our model showed little variance in performance with inorganic and Po fractions, indicating that both P forms are bound to primary soil properties detected directly by NIRS (Figure S1). However, this reproducibility could not provide information about the distribution of Pi and Po within individual fractions. The regressions of the NaOH-Po and HCl_conc.-P fractions exhibited the best performance of all the regressions, and showed similar statistical values for both validations. According to present knowledge, P or phosphate are not excited through NIR light waves to reflect in a particular part of the NIR spectrum. Hence recognizable excitation of P depends on organic or mineral compounds. It can be hypothesized that the separation of the HCl_conc. fraction in Pi and Po does not give reliable estimates of these P fractions. In general, the information on the P forms extracted using the NaOH and HCl_conc. steps are vague, describing these forms only as recalcitrant Pi and Po [33]. Overall, from the Hedley results, it can be interpreted that the effect of depth is more durable than the impact of nutrition addition in the data set. Chang et al. [72] reported that the ability to predict levels of extractable cations varied with the extraction method without explaining the differences; however, spectral changes associated with soil P content are usually associated with organic matter or crystal water. The quality of the Po fraction model is related to the underlying mechanism that produces near-infrared spectroscopy, because organic compounds are more easily irradiated than inorganic compounds [9]; this is because the organic compounds have characteristic infrared absorptions such as carbon-hydrogen bonds [73], while the P compounds correlated with organic matter in soil are detectable in NIRS [74]. In contrast to the Po fraction, the relationship between the Pi fraction and other soil properties detectable by NIRS is insufficient to predict Pi using NIRS [40]. If the model is improved, NIRS can help to evaluate the data obtained with the Hedley method, which means that the model based on the large sample numbers could be more precise. However, even with the results of this study, some trends are shown; for example, although the parent material and nutrient content influence the P status of the soil, other factors become more important.

As noted in previous studies, the labile P pool (Resin-P, NaHCO₃-Pi, and Po) could be recognized as a plant-available form [75,76]. Our results show that with P fertilization, the Resin-P fraction, in particular, increases strongly (Figure 1); this can only be inadequately predicted in the NIRS model. It can be concluded that a combination of NIRS modeling and chemical extraction for the labile P pool, such as an anion-exchange resin, provides the best results overall. Additionally, Reijneveld et al. [77] showed that a combination of chemical extraction using CaCl₂ and NIRS gave the best results compared to conventional soil chemical techniques.

In summary, we showed that the model prediction quality was better for Po than Pi fractions and the moderate, stable P pool than in the labile P pool (Table 4). Since the P content in the nutrient addition experiment was mainly concentrated in the Resin-P fraction of the topsoil, the modeling of Hedley P fractions with NIRS may be a promising method for determining plant-available P from soils with similar properties. Combining the measurement of Pi using Hedley fractionation and using the measurement of the more stable P pool with NIRS also appears promising in the natural grassland without nutrient additions.

5. Conclusions

Hedley fractionation was adapted to the specifics of the grassland soils of the Tibetan Plateau, and showed reliable results of soil available-P for plants. Additionally, the prediction of the Hedley fractions by the NIRS spectra gave satisfactory results for the moderate and stable P pools, except for the labile fraction. When P was added to the soils by fertilizer application, the labile Resin-P pool increased from below 3% in total to rates of to 18%, which was accompanied by a decrease in the model quality. Therefore, NIRS modeling may suit large-scale experiments requiring substantial spatial and temporal replication of soil samples within similar soil physical and chemical properties in combination with chemical soil extractions—such as the resin-P fractionation or the CaCl₂ extraction for labile P pools in soils—to lower the costs and increase the value of soil analysis results.

This study yielded basic P data and a specific description of Hedley fractionation and NIRS modeling for alpine grassland soils on the Tibetan Plateau. It will support the continuation of possible future studies that will help to investigate the responses of soil to environmental factors, such as climate warming and the stability of the alpine ecosystem. Hence, its implementation in extensive inventories would significantly improve data gathering. Further research should improve the accuracy of NIRS models for predicting Hedley precision in determining P fractions and their relationships to other soil parameters.