Long-term change in phosphorus behavior and the degree of P saturation in typical cropland after different fertilization practices

ABSTRACT

Soils with continuous application of fertilizer phosphorus (P) may present high accumulations of P and may be at high risk of P loss. In this study, we investigated the effects of different fertilization strategies on soil P dynamics in three typical cropland soils over broad spatiotemporal scales in China. The soil P speciation (characterized both by Hedley fractionation and ³¹P nuclear magnetic resonance (³¹P NMR)) and soil P sorption characteristics were examined under five fertilization treatments: unfertilized control (CK); chemical nitrogen (N) and potassium (K) (NK); N, P and K (NPK); NPK plus straw (NPKS); and NPK plus manure (NPKM). We found that with an increase of 100 kg ha⁻¹ in the cumulative P budget, the increases in Olsen-P and CaCl₂-P contents ranged in 1.98–7.35 mg kg⁻¹ and 0.46–1.46 mg kg⁻¹, respectively, among the three soils. The recovery rates of P in Hedley fractionation (87–111%) were much higher than those of ³¹P NMR spectroscopy (40–91%). More percentage of organic P was detected in Hedley fractionation than ³¹P-NMR in acidic soil, while more percentage of organic P was detected in ³¹P-NMR than Hedley fractionation in neutral and alkaline soils. Chemical fertilizer combined with manure or straw can decrease the bond strength (k) values compared to chemical fertilizer alone, while change trends in the maximum P sorption capacity (Qm) were dependent on the soil types. Moreover, the degree of P saturation (DPS) was significantly correlated with the cumulative P budget in each site. A change point in the relation between DPS and soil organic carbon (SOC) was found. After the change point, the DPS increased by 0.63–1.06% as the SOC increased by 0.1 g kg⁻¹. Agricultural measures that increased the contents of SOC and available P could increase DPS values effectively, especially manure amendment, which should be carefully considered to reduce P loss risk.

KEY WORDS:

• P speciation
• P sorption characteristics
• degree of P saturation
• long-term fertilization
• cropland soils

1. Introduction

Phosphorus (P) is essential for crop growth; however, it is the most frequently deficient nutrient in many agricultural regions due to its low availability and poor recovery (Zhang et al. 2021). Fertilizer P is easily adsorbed by the soil and converted from highly soluble P to sparingly soluble amorphous or crystalline P, which leads to a low use efficiency (10–25%) of applied P for crops (Balemi and Negisho 2012). To meet increasing food demand, an abundant supply of fertilizer is needed. However, the overuse of P fertilizer leads to excessive P accumulation in soil and causes finite high-quality P rock resource waste or even water body eutrophication (Li et al. 2019; Rowe et al. 2015). Meanwhile, excessive P accumulation in soil can also alter the amounts and forms of P and the characteristics of P sorption–desorption, which in turn has led to changes in the retention and loss risk of soil P (Withers et al. 2016).

The total stock of P is distributed among several fractions that have different bioavailabilities and degrees of incorporation with mineral or organic constituents in the soil. The forms of inorganic P (Pi) and organic P (Po) are largely dependent on the intrinsic properties of the soil type and can also be affected by agricultural practices and microbial activity (Ma et al. 2021). However, many more studies focus on the forms and transformation of Pi than Po due to its direct bioavailability for plants (Weihrauch and Opp 2018). In fact, organic P compounds are also necessary to maintain vital activities. Moreover, Po can be converted to Pi through mineralization by microbial and root-released phosphatases (Shen et al. 2011). Thus, both Pi and Po should be considered to comprehensively understand the cycling of soil P. Sequential fractionation schemes have long been adopted to quantify P forms, which have provided valuable information on soil P dynamics under different conditions, and the Hedley fractionation procedure has been widely used (Condron and Newman 2011). In addition, nuclear magnetic resonance spectroscopy (³¹P-NMR) is a powerful tool to identify organic P species, and it can complement the information collected by chemical sequential P fractionation (Turner, Mahieu, and Condron 2003). In agricultural soil, the forms of P can be quite different after fertilization management (Yan et al. 2017). Therefore, investigation of the dynamic changes in P speciation both by sequential fractionation and ³¹P-NMR following long-term fertilization in contrasting soils is instructive for rational P management.

The P sorption characteristics are useful for soil P retention capacity and environmental P risk assessments. Soil P isotherms are used to understand the mechanisms of soil P retention. The Langmuir equation is suitable for numerous soil types and can determine the maximum P sorption capacity and bonding strength (Wang and Liang 2014; Daly et al. 2015). The degree of P saturation (DPS) is generally estimated as the proportion of available P to the maximum P sorption capacity and has been increasingly viewed as an environmental indicator of the potential loss of P in soil (Kleinman and Peter 2017; Blombäck et al. 2021). A more direct way to estimate the risk for P loss is by measuring the dissolved P extracted by deionized water or CaCl₂ solution. However, these indicators are soil specific and cannot provide information on the additional capacity of soils to retain P. The DPS considers both the available P concentration and the maximum P sorption capacity of soil. Several studies have reported that DPS was strongly positively correlated with soil water-soluble P and labile P (Gatiboni et al. 2021, Fu et al.). In addition, DPS can be changed by land use and fertilization management (Abdala et al. 2012; Fu et al. 2021). For instance, Yan et al. (2017) reported that DPS was significantly higher in manure-treated soils than in soils treated with chemical fertilizers alone or with straw retention. Thus, the control of DPS in agricultural activities is important for ecosystem health, particularly in regions under long-term fertilization.

In the present study, three typical croplands of China are investigated: black soil (Luvic Phaeozem according to the Food and Agriculture Organization of the United Nations (FAO) classification), fluvo-aquic soil (Calcaric Cambisol in FAO) and red soil (Eutric Cambisol in FAO). These soils are found across highly diverse climatic zones, soil types and management practices and represent important agricultural regions and the most nutrient-additive areas of China. The changes in P speciation were determined by chemical fractionation and ³¹P-NMR in these three typical croplands, and the variations in soil P sorption characteristics and physicochemical properties were also considered. We hypothesized that under five long-term (1990–2015) fertilization treatments, i.e., unfertilized control (CK); chemical nitrogen (N) and potassium (K) (NK); N, P and K (NPK); NPK plus straw (NPKS); and NPK plus manure (NPKM), (1) the contents of Pi and Po would decrease in no-P treatments but increase in P-applied treatments owing to different crop P uptake and fertilizer P input; (2) the maximum P sorption capacity and bonding strength would be higher in no-P treatments than P-applied treatments; and (3) the DPS would increase in P-applied treatments but decrease in no-P treatments, which may both affected by P behaviors and soil properties.

1. Materials and methods

1.1. Site description and experimental design

Three long-term field fertilization sites were selected: Gongzhuling, black soil (GZL, 43°30′, 125°48′, located in the middle temperate zone), Jilin, northeastern China; Zhengzhou, fluvo-aquic soil (ZZ, 34°47′, 113°40′, warm temperate region), Henan, central China; and Qiyang, red soil (QY, 26°45′, 111°52′, subtropical region), Hunan, southern China. The cropping systems included monomaize cropping at GZL and wheat-maize double cropping at both ZZ and QY. No irrigation was provided to crops at GZL and QY, while the wheat crop was irrigated two or three times and the maize crop was irrigated once at ZZ. Information on the climatic conditions and initial soil physicochemical properties for the three sites is briefly shown in Table 1.

Table 1. Initial surface soil properties (0–20 cm in 1990) and climate conditions (data are the means from 1990–2015) at three long-term experimental sites (Gongzhuling, Zhengzhou, and Qiyang)

Parameters	GZL	ZZ	QY
Annual mean precipitation (mm)	591	641	1426
Annual mean evaporation (mm)	1409	1808	1435
Cumulative effective temperature (>10)	2800	5169	5600
Annual mean temperature (°C)	6.6	14.7	18.0
Cropping system	Single-cropping, maize	Double-cropping, maize/wheat	Double-cropping, maize/wheat
Soil classification in FAO	Luvic Phaeozems	Calcaric Cambisol	Eutric Cambisol
Soil classification in China	Black soil	Fluvo-aquic soil	Red earth
Soil texture	Clay loam	Light loam	Clay
Bulk density (g cm^−3)	1.2	1.4	1.2
Clay content (<0.002 mm, g kg^−1)	321	134	410
Soil pH (soil:water = 1:2.5)	7.6	8.3	5.7
Organic carbon (g kg^−1)	13.2	6.7	6.6
Total P (g kg^−1)	0.6	0.6	0.5

A randomized complete block design was used in the experiments, and detailed information regarding the long-term experiments is given by He et al. (2015). The following five treatments were assessed in this study: (1) CK, (2) NK, (3) NPK, (4) NPKM, and (5) NPKS. N, P and K were applied as urea, calcium triple superphosphate, and potassium sulfate. The annual fertilization rates are provided in Table 2. For each site, the total N applied (inorganic plus organic) was the same for the NK, NPK, NPKS and NPKM treatments except the NPKS treatment at QY, which received extra N from straw; inorganic P fertilizer was the same for the NPK, NPKM and NPKS treatments; and inorganic K fertilizer was the same for the NK, NPK, NPKM, and NPKS treatments.

Table 2. Annual rates of fertilizer N, P and K application in inorganic and organic forms at the three long-term (1990–2015) fertilization sites

Sites	Treatment^a	Maize		Wheat
		Inorganic N-P-K^b (kg ha^−1)	Organic N-P-K (kg ha^−1)	Inorganic N-P-K (kg ha^−1)	Organic N-P-K (kg ha^−1)
GZL	CK	0-0-0
	NK	165-0-68
	NPK	165–36-68
	NPKS^c	112–36-68	53-6-58
	NPKM^d	50–36-68	115–39-77
ZZ	CK	0-0-0		0-0-0
	NK	188-0-78		165-0-68
	NPK	188–41-78		165–36-68
	NPKS	188–41-78		121–36-68	44-8-86
	NPKM	188–41-78		50–36-68	115–66-92
QY	CK	0-0-0		0-0-0
	NK	210-0-70		90-0-30
	NPK	210–37-70		90–16-30
	NPKS	210–37-70	10-1-28	90–16-30	21-4-37
	NPKM	63–37-70	147–62-82	27–16-30	63–27-36

^aCK: unfertilized control; NK: nitrogen and potassium; NPK: inorganic nitrogen, phosphorous and potassium; NPKS: NPK plus straw return; NPKM: NPK plus manure.

^bInorganic N fertilizer as urea, P as calcium triple superphosphate, K as potassium sulfate.

^cAll of the maize straw was applied at GZL (approximately 7.5 t ha⁻¹) and ZZ (approximately 6.0 t ha⁻¹), while half of the straw (approximately 4.5 Mg ha⁻¹) was incorporated into the soil at QY.

^dThe manure was pig manure from 1990 to 2015 at GZL (23.0 t ha⁻¹) and QY (42.0 t ha⁻¹) and horse manure from 1990 to 1998 and cattle manure from 1999 to 2015 at ZZ (12.9 t ha⁻¹).

Weeds were controlled by hand weeding, and pesticides were applied when needed. Crops were harvested manually with sickles close to (approximately 0.02 m) the ground, and all of the harvested biomass was removed from the plots with little crop residue remaining (except in the NPKS).

1.2. Soil sampling and analysis

1.2.1. Soil characterization

For each treatment, soil samples were collected annually at a depth of 0–0.2 m immediately after the maize harvest. A soil auger (10 cm internal diameter) was used to collect the soil samples from three points in each plot, which were thoroughly mixed and treated as one composite sample. All of the soil samples were air-dried, sieved (2 mm) and stored for analysis. In this study, soils from 1990, 2000 (GZL and QY)/2002 (ZZ), and 2015 were obtained for analysis; these soils were archived samples that had been air-dried prior to long-term storage.

The soil organic C (SOC) was measured by sulfuric acid–K dichromate oxidation (Walkley and Black 1934); the soil pH was measured using a glass electrode (soil: water = 1:2.5). Total P was extracted with H₂SO₄-HClO₄, soil Olsen-P was extracted using 0.5 mol L⁻¹ NaHCO₃ (Olsen et al. 1954), CaCl₂-P was extracted using 0.01 mol L⁻¹ CaCl₂ (soil:solution = 1:5, shaken for 15 min at 25°C), and all P contents were determined using the spectrophotometric method (Murphy and Riley 1962). The contents of CaCO₃, free Fe₂O₃ and Al₂O₃ were determined following the methods of Lu (2000). In addition, the soil specific surface area (SSA) was measured using the Brunauer–Emmet–Teller method using a Gemini VII 2390 surface area analyzer (Micromeritics Instrument Corp., USA).

1.2.2. Sequential P fractionation and solution ³¹P-NMR spectroscopy

To determine the soil P speciation, the sequential extraction procedure of Hedley, Stewart, and Chauhan (1982), as modified by Tiessen and Moir (1993), was used. Briefly, triplicate subsamples of each soil (1 g) were sequentially extracted as follows: shaking in 30 mL of deionized water containing anion exchange resin, followed by extraction with 0.5 mol L⁻¹ NaHCO₃ (pH = 8.5), 0.1 mol L⁻¹ NaOH, 1.0 mol L⁻¹ HCl (dHCl) and hot concentrated HCl (cHCl), after which the soil residue was mineralized with concentrated H₂SO₄-H₂O₂ (300 µL 30 mg⁻¹ soil) at 350°C for 3 h (rate of 4°C min⁻¹). Between each consecutive step, the soil extract was shaken end-over-end in 50 ml centrifuge tubes for 16 h and then centrifuged for 10 min at 25,000 × g at 4°C. The supernatant was decanted for analysis. Both inorganic and organic P (calculated as the difference between the total P content after persulfate digestion and the inorganic P content) levels were determined in 0.5 mol L⁻¹ NaHCO₃, 0.1 mol L⁻¹ NaOH and 12 mol L⁻¹ HCl. All P fractions were detected using the spectrophotometric method (Murphy and Riley, 1962). In consideration of plant availability, P fractions were divided into labile P, moderately labile P, and nonlabile P (Rodrigues et al. 2016). Labile P includes Resin-P, NaHCO₃-Pi, and NaHCO₃-Po. Moderately labile P refers to NaOH-Pi, NaOH-Po and dHCl-P. Nonlabile P includes cHCl-Pi, cHCl-Po and residual-P.

Phosphorus (P-NaOH-EDTA) used in the ³¹P-NMR procedure was extracted by shaking 5 g of air-dried soil with 0.25 mol L⁻¹ NaOH + 0.05 mol L⁻¹ Na₂EDTA solution (soil: solution = 1:20) by shaking on a reciprocal shaker for 16 h at 20°C, and the centrifuged supernatant solution was passed through a 0.45 μm membrane filter. The extract was frozen at −80°C and lyophilized for ³¹P-NMR analysis. Freeze-dried NaOH-EDTA extracts were redissolved in 1 ml of 1 mol L⁻¹ NaOH and 0.1 ml of D₂O and transferred to 5 mm NMR tubes. Solution ³¹P NMR spectra were obtained using a Bruker AVANCE spectrometer (Bruker Ltd., Germany) operating at 243 MHz with a 7 μs pulse (45°), a delay time of 2.0 s, and an acquisition time of 0.25 s. The P compounds were identified by their chemical shifts (ppm) in accordance with Turner, Mahieu, and Condron (2003). The classes of P compounds include inorganic orthophosphate (~5.9 ppm), orthophosphate monoesters (3.0 ~ 5.8 ppm), orthophosphate diesters (−0.5 ~ 2.0 ppm), and pyrophosphate (~-4.5 ppm). The P extracted by NaOH-EDTA (P-NaOH-EDTA) was determined in diluted extracts (1:25) to prevent interference from EDTA. The soil contents of specific P compounds were calculated using the corresponding peak area portions and the P-NaOH-EDTA contents.

1.2.3. P sorption

In the sorption study, solutions with various concentrations of P that were used to define the sorption isotherms were prepared in a 0.01 mol L⁻¹ KCl solution (pH = 7.0) using KH₂PO₄. Soil samples of 1.0 g in triplicate were added to a 50-ml centrifuge tube followed by the addition of 25 mL of 0.01 mol L⁻¹ KCl solution containing P concentrations of 0, 10, 20, 40, 60, 80, and 160 mg L⁻¹. Two drops of chloroform were added to inhibit microbial activity. The suspensions were shaken on an end-over-end mechanical shaker at 25°C for 1 h. The mixture was then equilibrated for 24 h at a constant controlled temperature (25°C). After centrifugation at 5000 × g for 10 min, the supernatant was analyzed for P using the molybdenum blue colourimetric method. The maximum P sorption capacity (Qm, mg kg⁻¹) was calculated using the linearized form of the Langmuir equation:

(1) $\frac{C}{q}=\frac{1}{kQm}+\frac{C}{Qm}$(1)

where C is the P concentration in the equilibrium solution (mg L⁻¹); k is a constant related to bond strength (L mg⁻¹); q = q₀+ q₁, q is the total amount of P absorbed (mg kg⁻¹); q₀ is the original sorbed P (mg kg⁻¹) determined by the P release at 0 initial P concentration; and q₁ is the amount of added P sorbed by the soil and was calculated from the difference between initial and final solution P concentrations (mg kg⁻¹).

The DPS was obtained for all of the samples by the following equation:

(2) $DPS\hspace{0.17em}(\mathrm{\%})=\frac{Olsen\hspace{0.17em}-P}{Qm}\times 100$(2)

where DPS is the degree of P saturation, Olsen-P corresponds to the available P extract by 0.5 mol L⁻¹ NaHCO₃ (mg kg⁻¹), and Qm is the maximum P sorption capacity calculated from the Langmuir equation (mg kg⁻¹).

1.3. Cumulative P budget

Maize and wheat were harvested each year, and the grain yields and straw biomass were determined from the entire plot area. Crop samples were air-dried and then oven-dried at 65°C to a uniform moisture level and weighed. Plant samples were ground (<0.15 mm) and digested with H₂SO₄–H₂O₂, and the P concentration of the extracts was measured using the molybdenum blue colourimetric method (Page, Millar, and Keeney 1982). P uptake at each harvest was calculated by multiplying the grain and straw P concentrations by the grain yield and straw biomass. The annual P budget (from 1990 to 2015) was calculated by the fertilizer P input minus the annual crop P uptake (sum of grain and straw P uptake). The cumulative P budget was the sum of the annual P budgets.

1.4. Statistical analyses

All of the statistical analyses were conducted using SPSS 20.0. The data presented for all the chemical and physical analyses were the mean values of three replicates. A Pearson correlation analysis was used to explore the relationships between the cumulative P budget, P speciation, DPS and soil properties.

2. Results

2.1. The cumulative P budget, soil Olsen-P and CaCl2-P

The cumulative P budget and Olsen-P and CaCl₂-P contents in the three soils after different long-term fertilization treatments are presented in Table 3. The cumulative P budgets were −543〜756, −490〜2100 and −94〜3477 kg P ha⁻¹ among GZL, ZZ and QY, respectively. In the P-applied treatments, the cumulative P budgets were positive (i.e., P input>P output) and ordered as NPKM>NPKS>NPK for all three sites. For the treatments without P application (CK and NK), the cumulative P budgets were negative, with much smaller deficiency values at QY than at the other two sites.

Table 3. The cumulative P budget and the contents of Olsen-P and CaCl₂-P in each fertilization treatment at the three long-term experimental sites

Treatment^a	GZL			ZZ			QY
	Cumulative P budget (kg ha^−1)	Olsen-P (mg kg^−1)	CaCl2-P (mg kg^−1)	Cumulative P budget (kg ha^−1)	Olsen-P (mg kg^−1)	CaCl2-P (mg kg^−1)	Cumulative P budget (kg ha^−1)	Olsen-P (mg kg^−1)	CaCl2-P (mg kg^−1)
Initial	-	10.9 ± 0.8	4.6 ± 0.1	-	7.6 ± 0.4	2.3 ± 0.1	-	13.8 ± 0.8	4.4 ± 0.1
CK10	−135.9 ± 9.5	3.3 ± 0.1	4.3 ± 0.1	−276.6 ± 19.3	3.7 ± 0.1	2.6 ± 0.04	−43.4 ± 3.2	5.7 ± 0.2	4.1 ± 0.1
NK10	−254.3 ± 17.8	5.3 ± 0.2	3.3 ± 0	−239.9 ± 16.7	3.4 ± 0.1	2.1 ± 0.04	−87.7 ± 6.1	5.4 ± 0.2	5.2 ± 0.1
NPK10	83.3 ± 5.8	11.5 ± 0.5	6.1 ± 0.1	347.8 ± 24.3	12.3 ± 0.7	4.1 ± 0.1	310.3 ± 21.7	21.4 ± 1.0	7.1 ± 0.2
NPKS10	140.7 ± 9.8	10.9 ± 0.5	6.0 ± 0.1	473.3 ± 33.1	13.7 ± 0.6	4.3 ± 0.1	339.9 ± 23.7	21.9 ± 1.4	8.3 ± 0.2
NPKM10	409.8 ± 28.6	72.3 ± 3.6	10.3 ± 0.3	1264.7 ± 75.8	40.4 ± 2.1	7.3 ± 0.2	1471.3 ± 73.5	80.4 ± 4.3	13.1 ± 0.5
CK25	−271.2 ± 21.6	2.9 ± 0.1	3.5 ± 0.1	−490.2 ± 39.2	2.3 ± 0.1	3.3 ± 0.04	−77.3 ± 6.1	3.9 ± 0.1	3.7 ± 0.1
NK25	−543.3 ± 43.4	3.4 ± 0.1	2.8 ± 0.1	−426.6 ± 34.1	2.1 ± 0.1	3.0 ± 0.1	−94.1 ± 7.5	4.4 ± 0.2	3.9 ± 0.2
NPK25	127.9 ± 10.2	41.1 ± 2.1	11.7 ± 0.3	559.7 ± 44.7	15.1 ± 0.7	8.6 ± 0.2	871.9 ± 69.7	55.8 ± 2.7	10.3 ± 0.5
NPKS25	240.5 ± 19.2	34.4 ± 1.7	11.0 ± 0.3	893.5 ± 71.4	16.2 ± 0.8	8.9 ± 0.2	915.1 ± 73.2	57.3 ± 2.8	11.3 ± 0.5
NPKM25	756.6 ± 45.3	124.2 ± 5.2	22.8 ± 1.5	2100.1 ± 105	53.0 ± 2.1	15.5 ± 0.3	3476.9 ± 139	141.3 ± 4.5	25.3 ± 0.9

^aInitial: soil samples from 1990 at three sites; subscript 10: soils from 2000 for GZL and QY and 2002 for ZZ; subscript 25: soils from 2015 at three sites.

The Olsen-P and CaCl₂-P concentrations among all fertilization treatments and soil types were in the ranges of 2.1–141.3 mg kg⁻¹ and 2.1–25.3 mg kg⁻¹ and much higher in the P-applied treatments than in the no-P treatments, especially in the NPKM treatment. Both the contents of Olsen-P and CaCl₂-P were significantly positively correlated with the cumulative P budgets among the three sites (P < 0.001) (Figure 1). With an increase of 100 kg P ha⁻¹ of the cumulative P budget, Olsen-P and CaCl₂-P increased by 1.98–7.35 mg kg⁻¹ and 0.46–1.46 mg kg⁻¹ in the order of GZL>QY>ZZ.

Figure 1. Relationship between Olsen-P, CaCl₂-P and the cumulative P budget.

2.2. Soil P speciation determined by sequential fractionation and ³¹P NMR spectroscopy

Soil P speciation can be affected by both the soil type and fertilization treatment, as shown in Table 4 and Figure 2. The Hedley procedure recovered 88–109% (average of 98%), 92–103% (average of 97%) and 87–111% (average of 100%) of the total P (extracted by H₂SO₄-HClO₄ and independent of sequential fractionations) in GZL, ZZ and QY, respectively. The proportions of labile P were smallest in the initial soil samples and were 11%, 5% and 5% at GZL, ZZ and QY, respectively; the proportions of nonlabile P were largest for both GZL (51%) and QY (65%), while the moderately labile P was largest in ZZ (66%). Compared to that in the initial soils, the proportion of nonlabile P increased in the no-P treatment at all three sites, while the proportion of labile P (except for NK at QY) and moderately labile P decreased (Figure 2). In contrast, the proportion of nonlabile P decreased in the P-applied treatment at all three sites, while the proportion of labile P showed an increasing trend. Furthermore, P_it (sum of all Pi fractions and residual P) was dominant in all treatments and soils, ranging from 73 to 93%, whereas P_ot (sum of all Po fractions) was less than 27% in all samples. The average proportion of P_ot among all fertilization treatments and years was slightly higher in GZL (18%) than in ZZ (10%) and QY (11%).

Table 4. The recovery rate and the ratio of Pi to Po obtained from the Hedley procedure and 31P NMR spectroscopy at the three long-term experimental sites

Site	Treatment^a	Total-P (mg kg^−1)^b	Hedley procedure			^31P-NMR procedure
			Sum of P fraction (mg kg^−1)	Recovery (%)	Pit/Pot^c	NaOH-EDTA (mg kg^−1)	Recovery (%)	Pit/Pot^d
GZL	Initial	660	582	88	10	554	84	2
	CK10	530	509	96	6	225	43	1
	NK10	490	502	102	6	207	42	1
	NPK10	540	586	109	6	487	90	7
	NPKS10	600	621	103	4	546	91	2
	NPKM10	780	847	109	5	509	65	10
	CK25	530	495	93	3	234	44	1
	NK25	470	462	98	3	234	50	1
	NPK25	730	697	95	4	415	57	4
	NPKS25	800	781	98	5	560	70	5
	NPKM25	1350	1343	99	6	1066	79	12
ZZ	Initial	640	628	98	8	280	44	3
	CK10	580	569	98	9	232	40	2
	NK10	560	568	101	12	240	43	2
	NPK10	730	672	92	7	362	50	5
	NPKS10	730	726	99	7	425	58	4
	NPKM10	850	822	97	7	636	75	8
	CK25	560	552	99	10	225	40	2
	NK25	550	566	103	10	242	44	3
	NPK25	910	902	99	10	397	44	13
	NPKS25	960	906	94	7	485	51	13
	NPKM25	1120	1048	94	9	641	57	19
QY	Initial	490	502	102	7	427	87	10
	CK10	480	484	101	6	308	64	4
	NK10	480	499	104	6	317	66	5
	NPK10	710	657	93	7	475	67	13
	NPKS10	680	755	111	8	476	70	13
	NPKM10	1040	1112	107	9	936	90	16
	CK25	470	472	101	7	375	80	7
	NK25	460	478	104	6	330	72	9
	NPK25	1140	987	87	9	902	79	24
	NPKS25	1090	1029	94	13	943	86	49
	NPKM25	1760	1756	99	13	1295	74	19

^aInitial: soil samples from 1990 at three sites; subscript 10: soils from 2000 for GZL and QY and 2002 for ZZ; subscript 25: soils from 2015 at three sites.

^bTotal P extracted by H₂SO₄-HClO₄ and independent of Hedley fractionation.

^cP_it represents the total inorganic P contents (sum of all Pi fractions and residual P), and P_ot represents the total organic P contents (sum of all Po fractions).

^dP_it represents the sum of orthophosphate and pyrophosphate contents, and P_ot represents the sum of orthophosphate monoester and orthophosphate diester contents.

Figure 2. The P speciation obtained from Hedley fractionation and ³¹P NMR spectroscopy under different soil types and fertilization treatments. Initial: soil samples from 1990 at three sites; subscript 10: soils from 2000 for GZL and QY and 2002 for ZZ; subscript 25: soils from 2015 at three sites.

The content and proportions of P speciation following procedures of solution ³¹P-NMR analysis are shown in Table 4 and Figure 2. The recovery of P extracted by NaOH-EDTA (P-NaOH-EDTA) ranged from 40% to 91% of the total P (extracted by H₂SO₄-HClO₄) among all of the soil samples. Inorganic orthophosphate concentrations were in the ranges of 81–974, 143–607 and 231–1229 mg kg⁻¹ in GZL, ZZ and QY, respectively, which occupied more than 50% of the P-NaOH-EDTA in almost all samples. In contrast, the inorganic pyrophosphates were less than 4% of P-NaOH-EDTA and were not found in NPKS at ZZ in 2015. Orthophosphate monoesters dominated the Po pool and ranged from 5% to 50% of P-NaOH-EDTA in most of the samples. Orthophosphate diesters were mostly found in the acidic soil of QY but were not detected in the alkaline soil of ZZ. The proportions of orthophosphates in the P-applied treatments were higher than no P treatments, while the proportions of pyrophosphates and monoesters were higher in no P treatments than P-applied treatments in each site. Moreover, P_it (sum of orthophosphate and pyrophosphates) ranged from 38% to 98%, and P_ot (sum of orthophosphate monoesters and diesters) ranged from 2% to 62% in all samples. The average proportion of P_ot among all treatments and years was in the order of GZL (31%)>ZZ (20%)>QY (9%).

2.3. P sorption characteristics

The P sorption isotherm parameters based on the Langmuir model are presented in Table 5. All of the P sorption data fit the model well, and the coefficient of determination values (R²) were greater than 0.949. The values of maximum P sorption capacity (Qm) and binding energy (k) differed among treatments and soils. The average value of Qm was 608.5 mg kg⁻¹ among all treatments and years (ranging in 263.2–833.3 mg kg⁻¹) in QY, which was higher than the average values of GZL (377.2 mg kg⁻¹, ranging in 208.3–625.0 mg kg⁻¹) and ZZ (298.0 mg kg⁻¹, ranging in 212.8–434.8 mg kg⁻¹). The k values were in the order of QY>GZL>ZZ. In general, the Qm and k values were lowest under NPKM at each site; the highest k values were found in the NK treatment in 2015, but the highest Qm values were found in different treatments among the three sites.

Table 5. The P sorption indices and DPS for each fertilization treatment and soil type at the three long-term experimental sites

Treatment^a	GZL				ZZ				QY
	Qm^b (mg kg^−1)	k^c (L mg^−1)	R^2	DPS^d (%)	Qm (mg kg^−1)	k (L mg^−1)	R^2	DPS (%)	Qm (mg kg^−1)	k (L mg^−1)	R^2	DPS (%)
Initial	454.5	0.057	0.989	2.3	434.8	0.027	0.962	1.7	555.6	0.090	0.952	2.5
CK10	625.0	0.050	0.961	0.5	357.1	0.024	0.962	1.0	500.0	0.093	0.978	1.1
NK10	303.0	0.152	0.990	1.8	357.1	0.024	0.973	1.0	714.3	0.070	0.973	0.8
NPK10	312.5	0.097	0.967	3.7	333.3	0.029	0.986	3.7	714.3	0.070	0.984	3.0
NPKS10	263.2	0.062	0.952	4.2	322.6	0.023	0.958	4.3	769.2	0.043	0.978	2.8
NPKM10	294.1	0.030	0.958	24.6	322.6	0.020	0.993	12.5	263.2	0.052	0.986	30.6
CK25	454.5	0.067	0.966	0.6	232.6	0.025	0.952	1.0	526.3	0.143	0.989	0.7
NK25	322.6	0.157	0.963	1.1	243.9	0.030	0.980	0.9	769.2	0.217	0.995	0.6
NPK25	434.8	0.106	0.949	9.5	243.9	0.030	0.951	6.2	714.3	0.156	0.999	7.8
NPKS25	476.2	0.070	0.955	7.2	217.4	0.023	0.985	7.5	833.3	0.135	0.999	6.9
NPKM25	208.3	0.037	0.964	59.6	212.8	0.021	0.966	24.9	333.3	0.095	0.964	42.4

^aInitial: soil samples from 1990 at three sites; subscript 10: soils from 2000 for GZL and QY and 2002 for ZZ; subscript 25: soils from 2015 at three sites.

^bMaximum P sorption capacity.

^cConstant related to bonding strength.

^dDegree of P saturation, DPS = Olsen-P/Qm

The DPS values were in the ranges of 0.5–59.6%, 0.9–24.9% and 0.6–42.4% in the GZL, ZZ and QY soils, respectively (Table 5). Compared with the initial soils, the DPS values gradually decreased over time in the no-P treatments but markedly increased in the P-applied treatments, especially under NPKM. In general, the DPS values were >24.9% under NPKM at each site after a 25-year continual combination of chemical fertilizer and manure, while the DPS values were <9.5% in the remaining treatments.

2.4. Correlations between the cumulative P budget, P speciation, DPS and selected soil properties

Summarizing the P fractions obtained by the Hedley procedure into organic and inorganic forms, it can be seen that there was a linear and significant increase in P_it (sum of all Pi fractions and residual P) with increasing cumulative P budget (P < 0.01), and the change rates (slopes of the linear equation obtained from the P_it response to cumulative P budget) were in the order of GZL (0.57)>QY (0.35)>ZZ (0.17). Additionally, the P_ot (sum of all Po fractions) increased with increasing cumulative P budget, but the changing rates were much lower than the P_it at each site, which was in the range of 0.0086–0.0439 (Figure 3). Furthermore, a positive linear relationship can be observed between the cumulative P budgets and DPS. For an increase of 100 kg ha⁻¹ of the cumulative P budget, the DPS increased by 4.39%, 1.25% and 0.86% in the order of GZL>QY>ZZ.

Figure 3. Relationship between the P_it (sum of all Pi fractions and residual P in Hedley procedure), P_ot (sum of all Po fractions in Hedley procedure) and degree of P saturation (DPS) with the cumulative P budget.

Of the selected soil properties, the highest pH was observed at ZZ, followed by GZL and QY (Table 6). Specifically, the pH decreased slightly over time for each treatment at GZL and ZZ but obviously decreased under NK, NPK, and NPKS over time at QY. The average SOC contents were ordered as GZL>QY>ZZ, and the SOC contents in the P-applied treatments were higher than those in the no-P treatments at all three sites. In addition, the concentration of CaCO₃ was highest at ZZ, while the Fe₂O₃ and Al₂O₃ concentrations were highest at QY. CaCO₃, Fe₂O₃ and Al₂O₃ concentrations differed slightly among treatments and years. The SSA values were higher at QY than at GZL and ZZ. The SSA values were similar among treatments within the same year; however, after 25 years of fertilization, they strongly decreased in GZL and ZZ but not in QY.

Table 6. Selected soil properties in all fertilization treatments and years at the three long-term experimental sites

Items	GZL				ZZ				QY
	Min	Max	Median	Mean^a	Min	Max	Median	Mean	Min	Max	Median	Mean
pH (H2O)	6.0	7.7	7.1	7.0	7.9	8.5	8.2	8.2	4.0	5.8	4.8	5.1
SOC (g kg^−1)	10.5	21.4	13.2	14.1	5.9	11.5	7.4	8.0	6.8	14.3	9.8	9.7
CaCO3 (g kg^−1)	10.5	37.0	17.5	21.3	52.6	77.4	64.0	66.6	4.6	15.2	7.1	8.5
Fe2O3 (g kg^−1)	1.3	2.5	1.7	1.9	0.7	1.2	1.1	1.0	2.3	5.9	4.3	4.0
Al2O3 (g kg^−1)	1.3	1.9	1.5	1.6	0.6	0.7	0.7	0.7	2.1	3.0	2.7	2.6
SSA (m^2 g^−1)^b	5.0	18.4	7.5	11.0	4.3	12.1	8.6	8.5	23.6	29.9	27.5	26.8

^aMean value of measurements in all fertilization treatments and years.

^bSSA: soil specific surface area.

This study focused on the relationships between soil properties and DPS at each site because DPS can reflect both the P concentration and P sorption characteristics. The Pearson correlation analysis showed that DPS was significantly positively correlated with the SOC at all three sites (P < 0.001) (Table 7). Additionally, the DPS was negatively correlated with pH in ZZ (P < 0.05) and negatively correlated with Fe₂O₃ in QY (P < 0.05). Furthermore, the response of the relationship between DPS and SOC was well fitted by the linear-linear model (P < 0.01), with changes of 16.5, 8.7 and 10.2 g kg⁻¹ in SOC content at GZL, ZZ and QY, respectively (Figure 4). After the change point, the DPS increased by 0.63%-1.06% as the SOC content increased every 0.1 g kg⁻¹ in the order of GZL>QY>ZZ. The coefficients of the equation after the change point were 9.8-, 3.4- and 5.7-fold those of the previous equation at GZL, ZZ and QY, respectively.

Table 7. Pearson’s correlation coefficients of the degree of P saturation (DPS) and selected soil properties

		pH	SOC	CaCO3	Fe2O3	Al2O3	SSA	DPS
GZL	pH	1
	SOC	−0.096	1
	CaCO3	0.461	−0.016	1
	Fe2O3	−0.581*	0.196	0.003	1
	Al2O3	−0.631*	−0.092	−0.034	0.828**	1
	SSA	−0.234	−0.169	−0.779**	−0.293	−0.1	1
	DPS	0.193	0.888***	0.199	0.171	−0.091	−0.386	1
ZZ	pH	1
	SOC	−0.758**	1
	CaCO3	−0.908**	0.526*	1
	Fe2O3	−0.691**	0.639*	0.649*	1
	Al2O3	−0.346	0.386	0.093	0.225	1
	SSA	0.747**	−0.308	−0.766**	−0.168	−0.103	1
	DPS	−0.643*	0.893***	0.417	0.483	0.352	−0.229	1
QY	pH	1
	SOC	0.339	1
	CaCO3	−0.196	−0.104	1
	Fe2O3	−0.886***	−0.39	0.519	1
	Al2O3	−0.945***	−0.127	0.321	0.892***	1
	SSA	0.001	−0.325	−0.541*	−0.166	−0.108	1
	DPS	0.452	0.922***	−0.231	−0.563*	−0.332	−0.293	1

Figure 4. Relationship between degree of P saturation (DPS) and soil organic carbon (SOC). The arrow indicates the change point of the response curve.

3. Discussion

3.1. Effects of long-term fertilization on soil P content and speciation

In heavily fertilized soils, the P level usually exceeds the agronomic optimum required for crops, resulting in large P accumulation and changes in P behaviors (Shi et al. 2013; Yan et al. 2017). In this study, the cumulative P budget differed greatly among fertilization treatments and soil types, which may be related to the large variation in the fertilizer P input and crop P uptake. Variation in the P contents can be largely explained by the differences in the P budget in agricultural soil (Rowe et al. 2015). The change in Olsen-P was 1.98–7.35 mg kg⁻¹ with each increase of 100 kg ha⁻¹ of the cumulative P budget across our study sites, which was similar to the results of 1.44–5.74 mg kg⁻¹ around China (Cao et al. 2012). In addition, the increasing rates in both Olsen-P and CaCl₂-P by each 100 kg ha⁻¹ cumulative P budget differed greatly and were ordered as GZL>QY>ZZ. This large difference among the three sites can be explained by the soil physical-chemical properties, texture and temperature (Aulakh, Garg, and Kabba 2010; Cao et al. 2012). The increase values were higher in GZL and QY than in ZZ, which may be due to the high SOC and clay content usually maintaining a high P content (Zhu et al. 2013). In many countries and regions, the dissolved P extracted by 0.01 mol L⁻¹ CaCl₂ (CaCl₂-P) is often used as a direct indicator of the proneness to P leaching (Blombäck et al. 2021). Thus, establishing the relationship between Olsen-P, CaCl₂-P and the cumulative P budget is helpful for soil P management for both agricultural production and environmental protection.

The sequential P fractionation approach developed by Hedley, Stewart, and Chauhan (1982) has been widely used to differentiate labile and stable P fractions in soils, expanding our knowledge of the transformation of inorganic and organic P under different land uses and management practices (Negassa and Leinweber 2010). However, this method can only describe the distribution of soil P over ‘operationally defined’ fractions (Gatiboni et al. 2021). Solution ³¹P NMR spectroscopy enables the identification and speciation of Po classes (Hashimoto and Watanabe 2014). However, this technique is unsuitable for the identification of inorganic P (Pi) species. In summary, there are limitations in providing information about the biogeochemical behavior of P by one method alone. In our study, Hedley fractionation and solution ³¹P NMR spectroscopy were combined to substantially characterize P transformation in typical agricultural soils after long-term fertilization. Generally, the average recovery of all fertilization treatments and years at each site was above 98%, indicating that the Hedley method is suitable for a wide range of soil types (Negassa and Leinweber 2010). In contrast, the average recovery of P extracted by NaOH-EDTA was only 76% in QY, followed by GZL at 65% and ZZ at 50%, which may imply that solution ³¹P NMR spectroscopy is more suitable in acidic soil than alkaline soil (Peng, Fang, and Zhan et al. 2009). Moreover, there were great discrepancies among the contents and percentages of P_it/P_ot determined by these two methods. The average P_it/P_ot was higher in Hedley fractionation than ³¹P-NMR in GZL and ZZ, and the trend was opposite in QY. It can be concluded that more percentage of organic P was detected by Hedley fractionation than ³¹P-NMR at QY, while more percentage of organic P was detected by the ³¹P-NMR procedure than Hedley fractionation at GZL and ZZ. This result could be partly because alkaline extractants (i.e., NaHCO₃ and NaOH at a pH of 8.5) are more efficient at extracting Po in acidic soils of QY than in neutral soil of GZL and alkaline soil of ZZ (Turner et al. 2005). On the other hand, if the NaHCO₃, NaOH or hot cHCl Pi extracts contain significant concentrations of poly and pyrophosphates that are not detected by the molybdenum blue method and thus are included in the Po estimation by subtraction, the overestimation of Po fractions is possible, whereas if the acid reagents induce Po hydrolysis, underestimation of Po fractions can occur (Turner et al. 2005).

In the present study, the black soil in GZL, fluvo-aquic soil in ZZ and red earth in QY were derived from contrasting parent materials and had different initial physicochemical properties. Moreover, different long-term fertilization treatments also strongly changed selected soil properties, as shown in Table 6. These soil properties can affect the contents and fractions of P directly or indirectly (Wang and Liang 2014; Weihrauch and Opp 2018). The organic P was highest at GZL, which may be due to its high SOC content in the black soil; dHCl-Pi was highest at ZZ, which may be related to the high CaCO₃ content in the alkaline fluvo-aquic soil, and NaOH-P was highest at QY, which may be due to the high Fe₂O₃ and Al₂O₃ levels in the acidic red soil. Previous work has reported that organic P is associated with SOC, the dHCl-Pi fraction is largely bound to Ca, and the NaOH-Pi fraction is bound to the exterior of Al and Fe oxides (Tiessen and Moir 1993; De Schrijver et al. 2012).

3.2. Effects of long-term fertilization on soil P sorption characteristics

Fertilizer P application can influence soil P sorption behavior (Pizzeghello et al. 2011; Abdala et al. 2012; Yan et al. 2017). Our results indicated that chemical fertilizer combined with manure or straw can decrease the constant k values compared to chemical fertilizer alone, while changes in Qm were dependent on the soil type. Similar results in paddy soil have been reported by Lin, Wang, and Lin et al. (2011). In contrast, Yan et al. (2017) found that swine manure and rice straw increased the Qm compared to chemical fertilizer, while there were no significant differences in the constant k. These inconsistent results may be attributed to large variations in the type and rate of fertilizer P and the diverse soil properties and cropping systems among these study areas. Both Qm and k can be affected by soil properties and management systems (Daly et al. 2015; Fink et al. 2016; Jalali and Jalali 2016). It has been found that Qm was significantly correlated with the content of iron oxides in subtropical soils (Fink et al. 2016). Daly et al. (2015) reported that binding energies were strongly negatively correlated with soil pH, while positively correlated with extractable aluminum. In addition, organic amendments significantly increased P sorption capacity compared to chemical fertilizer treatment, which may be due to the increased amount of amorphous Fe and Al in organic treatments (Yan al. 2017).

Among the three soils, the mean values of Qm and k among the different fertilization treatments and years were ordered as QY>GZL>ZZ, and the mean pH values were in the order of ZZ>GZL>QY, which may indicate that Qm and k decreased with increasing pH. Quintero, Boschetti, and Benavidez (1999) also found that the P retention and maximum buffering capacities (calculated as Qm×k) in soils decreased as the pH increased from 3 to 7. When the pH value is low, a large number of variable positive charges occur; thus, the adsorption of phosphate radicals is strong. Increasing pH can inhibit P adsorption and decrease the adsorption energies (Abdala et al. 2012; Hou et al. 2018). In addition, the Fe₂O₃ and Al₂O₃ contents were highest in QY and lowest in ZZ, while the CaCO₃ content was the opposite. Wang et al. (2014) found that the P adsorption capacity tends to be relatively strong when there are a large number of iron and aluminium oxides. It has also been found that adsorption energies (k) were small in limed soils (Broggi et al. 2011). Furtherly, soils containing large quantities of clay can fix more P than soils with low clay (Arai and Sparks 2007). Soil organic matter can also increase P fixation by directly increasing the number of adsorption sites or indirectly inhibiting the crystallization of Al, which in turn can increase P sorption (Borggaard et al. 1990). A comparison of the black soil with the fluvo-aquic soil shows that the black soil contains more clay and SOC contents, which may explain why the capacity of P sorption is stronger in GZL than in ZZ.

4.3. Relationships between cumulative P budget, soil properties and DPS and implications for P management

The DPS has been used in many countries to evaluate the risk of P loss from agricultural soils. However, the DPS is determined by indirect methods in most studies, and the P sorption isotherms are not measured because it is laborious and time-consuming (Pizzeghello et al. 2011; Blombäck et al. 2021). In this study, we directly calculated the DPS as the ratio of Olsen-P in the soils to the Qm derived from the P sorption isotherms. Our results showed that the DPS ranged widely from 0.5% to 59.6% among fertilization treatments and soil types, which might be attributed to huge differences in cumulative P budgets and soil properties after different long-term fertilization practices. The DPS was significantly positively correlated with the cumulative P budget (P < 0.01), and the increasing rates were GZL>QY>ZZ. The increases were 17–57 mg kg⁻¹ for P_it (sum of all Pi fractions and residual P) and 0.86–4.39 mg kg⁻¹ for P_ot (sum of all Po fractions) across the study sites with each increase of 100 kg ha⁻¹ of the cumulative P budget. These results indicate that the accumulation of P in inorganic P fractions is more obvious than organic P in soils after a history of fertilization, which is consistent with other previous reports (Gatiboni et al. 2021). Moreover, our results also showed that the DPS values were significantly correlated with several soil properties, such as SOC, pH and Fe₂O₃. SOC is generally considered to be an important factor regulating soil DPS values, but this is largely based on theoretical considerations (Zou, Fu, and Cao 2011). Our linear-linear model demonstrated that there is a change point in the response relationship between DPS and SOC, providing evidence of a remarkable effect of SOC on the DPS in agricultural soils with high organic matter content. In addition, the positive relationship between the DPS and SOC indicates that the potential of P release may increase with increasing organic carbon content. This may be attributed to that SOC can provide organic acid anions competing for the binding sites with PO₄^3-, which may decrease the soil retention of P (Weihrauch and Opp 2018).

Overall, all the CaCl₂-P and Olsen-P contents and the DPS values were significantly positively correlated with the cumulative P budget among the three soils with contrasting physicochemical properties. Generally, the environmental risk of P loss from agricultural land can be estimated considering two aspects: the dissolved P content in runoff and the threshold of available P content or DPS in the soil (Davis et al. 2005). The soil Olsen-P contents in several fertilization treatments in 2015 were higher than 40 mg kg⁻¹ in our study, above which many soils in China are at risk of P leaching (Zhao et al. 2007). The DPS values in all treatments were much lower than 25% among the three soils in 2015, except for the MPKM. A DPS value of 25% is usually established as a critical value, above which soils tend to have a higher desorbable P available for runoff or leaching (Pothig et al. 2010). The soil CaCl₂-P and Olsen-P contents and DPS values were relatively high under NPKM after 25 years of continual fertilization, which indicates that chemical fertilizer combined with manure amendment may cause P loss in these three soil types. In the N-based fertilization plans, only total N applied (inorganic plus organic) was kept at the same level, and the extra P and K brought by manure and straw were ignored. Thus, to prevent the risk of P loss, the amount of manure added should be carefully considered (Pizzeghello et al. 2011). Moreover, more effort should be made to better utilize the large amount of cumulative P in soil for finite P resources and water security.

4. Conclusions

This study found that the contents of Olsen-P and CaCl₂-P and the DPS values were significantly positively correlated with the cumulative P budgets (P < 0.001), indicating that a large accumulation of P in soil can improve P availability but also increase the risk of P loss from agricultural systems. The soil P speciation, P sorption characteristics and DPS values differed greatly among fertilization treatments and soil types. There were great discrepancies in the contents and percentages of soil Pi and Po characterized by Hedley fractionation and ³¹P NMR, which confirmed the need for the combined use of different methods to determine P speciation for the best management of soil P. Chemical fertilizer combined with manure or straw can decrease the constant k values compared to chemical fertilizer alone, while changes in Qm were dependent on the soil type. There is a change point in the response relationship between DPS and SOC. After the change point, the coefficients of the linear equation were 3.4–9.8-fold that before the change point. Agricultural measures that increased the contents of SOC and available P could effectively increase DPS values. In particular, N-based fertilization plans are widely used in China, and manure amendments or straw return should be carefully evaluated because they can add abundant amounts of extra P.

Acknowledgments

We acknowledge all staff for their valuable work associated with these long-term Monitoring Network of Soil Fertility and Fertilizer Effects in China.

Disclosure statement

No potential conflict of interest was reported by the authors.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This work is supported by the GDAS’ Project of Science and Technology Development [2019GDASYL-0103033, 2019GDASYL-0104012], the National Natural Science Foundation of China [41977103] and China Agriculture Research System of MOF and MARA [CARS-170203].