Half of global agricultural soil phosphorus fertility derived from anthropogenic sources

The use of mineral phosphorus (P) fertilizers, often referred to as anthropogenic phosphorus, has dramatically altered the global phosphorus cycle and increased soil phosphorus fertility and crop yields. Quantifying agriculture’s reliance on anthropogenic phosphorus requires estimates of its contribution to agricultural soil fertility. Here we present a model of soil phosphorus dynamics simulating phosphorus availability in agricultural soils for individual countries from 1950 to 2017. Distinguishing between anthropogenic and natural phosphorus pools and accounting for farming practices, agricultural trade and crop–livestock recycling, we estimate that the global anthropogenic contribution to available phosphorus in agricultural soils was 47 ± 8% in 2017. Country-level anthropogenic phosphorus signatures vary according to cumulative fertilizer use and phosphorus availability in soil inherited pre-1950, with negligible influence of the trade of feed and food products. Despite different historical trajectories, we find that Western Europe, North America and Asia are similarly reliant on anthropogenic phosphorus, with nearly 60% of the total available phosphorus of anthropic origin in 2017. Conversely, anthropogenic phosphorus inputs in Africa remained low over the study period, contributing only around 30% of available phosphorus. The unequal reliance of agricultural soil fertility and food production systems on anthropogenic phosphorus resources highlights the need for a fairer management of the world’s remaining phosphate rock resources. About half of the current available phosphorus in agricultural soil globally is derived from anthropogenic sources, according to country-scale simulations of phosphorus dynamics between 1950 and 2017.

The use of mineral phosphorus (P) fertilizers, often referred to as anthropogenic phosphorus, has dramatically altered the global phosphorus cycle and increased soil phosphorus fertility and crop yields. Quantifying agriculture's reliance on anthropogenic phosphorus requires estimates of its contribution to agricultural soil fertility. Here we present a model of soil phosphorus dynamics simulating phosphorus availability in agricultural soils for individual countries from 1950 to 2017. Distinguishing between anthropogenic and natural phosphorus pools and accounting for farming practices, agricultural trade and crop-livestock recycling, we estimate that the global anthropogenic contribution to available phosphorus in agricultural soils was 47 ± 8% in 2017. Country-level anthropogenic phosphorus signatures vary according to cumulative fertilizer use and phosphorus availability in soil inherited pre-1950, with negligible influence of the trade of feed and food products. Despite different historical trajectories, we find that Western Europe, North America and Asia are similarly reliant on anthropogenic phosphorus, with nearly 60% of the total available phosphorus of anthropic origin in 2017. Conversely, anthropogenic phosphorus inputs in Africa remained low over the study period, contributing only around 30% of available phosphorus. The unequal reliance of agricultural soil fertility and food production systems on anthropogenic phosphorus resources highlights the need for a fairer management of the world's remaining phosphate rock resources.
The Anthropocene is characterized by the profound, anthropogenic disturbance of the global biogeochemical cycles 1 . In particular, the global phosphorus (P) cycle has been altered in an unprecedented way, both in time and in space 2 . At the core of these modifications stand phosphate rocks and their use in agriculture, mainly as mineral P fertilizers but also, to a lesser extent, as mineral feed for livestock animals. Both mineral P fertilizers and mineral P feed are derived from industrial treatments-to increase the P solubility of phosphate rocks-and are hereafter referred to as anthropogenic P.
The massive use of mineral P fertilizers combined with a sharp rise in the number of livestock animals-and resulting mineral P feed use and manure production-have globally increased cropland soil available P, although with some variations across world regions 3 . At the local scale, P has been added to agricultural soils through massive application of mineral P fertilizers-sometimes far above crop P uptake-resulting in increasing soil P legacy, together with P deficits in some regions 4 . Indeed, contrary to nitrogen, if over applied, P accumulates in soils, thus building a legacy from past agricultural practices. Large amounts of P have also been transferred from grasslands to croplands through the production, transport and soil application of animal manure 5 . At the global scale, P has been displaced through the extraction and international trade of phosphate rocks 6 . Highly concentrated in a few

Modelling phosphorus fluxes in agricultural soils
To achieve this objective, we developed a model that simulates for each country and with a yearly time step the evolution of agricultural soil P stocks (without distinction between cropland and grassland soils) from 1950 to 2017, accounting for (1) direct P fertilizer inputs to soils, (2) P recycling from the crop-livestock cycling loop and (3) geographic displacement of P through the trade of agricultural products (Fig. 1). We modelled soil available P as the interaction between a labile P pool (LP) and a stable P pool (SP), which is commonly done in the literature [17][18][19] . The labile P pool refers to soil P directly available for plant uptake while the stable P pool acts as a slow-release buffer that could replenish the labile P pool following crop P uptake. Each pool was further broken down into a natural (Nat) versus an anthropogenic (Ant) component. This allowed us to track the behaviour of P inputs according to their anthropic versus natural origin in a similar way as isotopic labelling approaches 20 . For each country, the sizes of the four P pools (LP Nat , LP Ant , SP Nat , SP Ant ) were simulated following the size of the previous years, soil P inputs and outputs, and soil P exchanges between labile and stable P pools (Methods).
In our model, P harvested from cropland and grassland soils (P harvest,Tot ) was computed as a function of the size of the labile P pool (equation (4)). The two parameters (β and γ) involved in the relation between P harvest and the size of the labile P pool were calibrated for each country. Similarly, we calibrated the P exchange function between stable and labile soil P pools (equation (5)). Both calibrations were conducted independently for each country and were performed so that the model outputs could reproduce both P harvest time series places such as Morocco and Western Sahara, phosphate rocks have been redistributed, yet unevenly, to the rest of the world, thus reshaping the spatial distribution of soil P. More recently, the trade of agricultural products has also contributed to displacing large quantities of P worldwide through soil application of animal and human manure derived from traded feed and food products 7 .
This alteration of the global P cycle has led to an overall rise in soil P fertility, thus increasing crop yields and helping to achieve food security objectives 8 . However, agriculture has become dependent on phosphate rocks, an alarming situation due to the progressive exhaustion of the easily accessible remaining reserves of this fossil resource 9,10 . Although it will take centuries for phosphate rock reserves to be depleted, short-term challenges are likely to arise in the next years to decades 11 . The prices of mineral P fertilizers are likely to increase in most world regions, and potential geopolitical conflicts might break out, thus putting at risk the resilience of agricultural systems worldwide 9,12,13 .
Those alterations of the global P cycle and of soil P fertility following anthropogenic P supply have already been highlighted by recent studies 2,4,14-16 . All have quantified P fluxes within the global food system and highlighted a massive use of phosphate rocks to fertilize agricultural soils. They have also reported inefficient management of P resources, due mostly to losses to aquatic systems and over-fertilization in certain areas. However, none has estimated the specific and cumulated contribution of anthropogenic P supply to P fertility of global agricultural soils. Providing such estimate is key (1) to assess the anthropogenic disturbance of the global P cycle, (2) to estimate the past and current reliance of agricultural soils and food production on the finite phosphate rock resources and (3) Fig. 1 | Structure of the model for a given country, with specific focus on soil P compartments. The anthropogenic (Ant) versus natural (Nat) soil P pools are represented in purple and green, respectively. LP and SP refer respectively to the labile and stable P pools. The fluxes in solid lines were explicitly simulated. Conversely, those in dotted lines were not modelled because they do not modify the anthropogenic signature of P pools. As illustrated by the pie charts, each flux was subdivided into anthropogenic (purple) versus natural (green) components. The blue box refers to a given country that interacts with the rest of the world (other countries) through the trade of food and feed products. In addition to soil P pools, we also modelled the anthropogenic signature of all P fluxes incoming to or outgoing from agricultural soils (Fig. 1). We assumed that P harvest had the same anthropogenic signature as that of the soil labile P pool in which plant uptake occurred. Similarly, we assumed that imported food and feed products had the same P signature as that of the soil labile P pool of the countries where they came from. The anthropogenic P signature of animal manure was assumed equal to that of animal intake. By construction and because they are entirely processed from the chemical industry, an anthropogenic signature of 100% was attributed to both mineral P fertilizers applied to soils and mineral P feed supplements given to livestock animals (Fig. 1).
We initialized the size of the P pools by assuming that, although some mineral P fertilizers were sometimes applied to agricultural soils in the first half of the twentieth century 21-23 , anthropogenic soil P pools were negligible before 1950. The sizes of the labile and stable natural P pools were estimated on the basis of the observed crop and forage P harvest in 1950 and by assuming an equilibrium between the two soil P pools (Methods). Hereafter, we refer to the signature of the soil labile P pool both in 2017 and throughout the 1950-2017 period.

Spatial and historical variations in anthropogenic phosphorus signatures
Our results show that the use of anthropogenic P-derived from phosphate rocks-greatly modified agricultural soil P fertility worldwide (Fig. 2). In 2017, the global mean anthropogenic signature of the soil labile P pool was 47 ± 8%, suggesting that almost half of the global agricultural soil P fertility was derived from anthropogenic P supply. This result clearly mirrors the overall intensification growth path of the world agriculture, which has relied on the massive use of mineral fertilizers from the 1950s 24 .
Our findings also highlight large spatial variabilities of anthropogenic soil P signatures among world regions (Fig. 2). North America, Western Europe and Asia displayed the highest anthropogenic signatures, with values of 70 ± 8%, 62 ± 8% and 61 ± 6%, respectively. Central and South America had an intermediate anthropogenic signature of 41 ± 10%, close to that of Eastern Europe (39 ± 10%). In Africa, the soil P anthropogenic signature remained moderate at 32 ± 6%. Finally, Oceania, which includes Australia and New Zealand, had the lowest anthropogenic signature, with a value of 18 ± 11%.
We also found that the anthropogenic signatures of soil P were explained by two main factors: (1) the cumulative application of anthropogenic P-mostly as mineral P fertilizer-over the 1950-2017 period and (2) the initial soil P status in 1950-assimilated to natural P in our model ( Supplementary Information section 8). As a result, the massive mineral P fertilization in North America, Western Europe and Asia since the mid-twentieth century translated into high anthropogenic signatures in those regions. Besides, high initial soil P pools in 1950 in Western Europe mechanistically 'diluted' the massive anthropogenic P supply to this region, which has been three times higher than in North America and Asia (Extended Data Table 1). This explains why, despite very different cumulated mineral P fertilization (Extended Data Fig. 1), the three regions all displayed similar anthropogenic signatures in 2017. Finally, and contrary to the stable P pool, the anthropogenic signature of the labile P pool depended on a third and last driver-the P transfer between the soil P pools (Supplementary Information section 8)which underlines the need to account for soil P dynamics when simulating soil P availability and its changes over time.
We found that the soil P anthropogenic signatures displayed highly variable historical trajectories among countries (Fig. 3), thereby reflecting contrasting and region-specific histories of soil P fertilization (Extended Data Figs. 2 and 3). More precisely, our results showed that the anthropogenic signature of industrialized Western countries such as France, the Netherlands and the United States rose sharply from the 1950s to the 1970s ( Fig. 3) because of the very early and massive application of mineral P fertilizers. However, since the 1980s, the anthropogenic signatures of France and the Netherlands have stabilized around values of 71 ± 7% and 49 ± 10% (in 2017), respectively. Conversely, the signature of the soil labile P pool in the United States was still rising in 2017. The observed stabilization in France and in the Netherlands resulted from smaller mineral P inputs over the past four decades, compensated partly by higher manure inputs with fewer effects on the anthropogenic signatures (Fig. 4). Indeed, when the amount of imported feed products is small compared with that of domestically produced feed-which is the general case-manure simply recycles P from domestic feed and forage back to agricultural soils, without affecting the anthropogenic signature of agricultural soils.
The soil P anthropogenic signatures of Brazil, India and China took off in the 1970s, up to values of 61 ± 9%, 67 ± 8% and 74 ± 6% in 2017, respectively (Fig. 3). Those results reflect the Green Revolution and the sharp increase in mineral P fertilizer use that transformed agriculture and soil fertility of most Asian and South American countries from the 1960s. Interestingly, the signatures of those countries have become similar to and even higher than those of Western European countries in 2017.
Finally, the soil P anthropogenic signatures of most African countries (for example, Morocco and Zimbabwe) exhibited a late and slow take-off, with values remaining low and close to 20-30% in 2017. This is due mainly to small applications of mineral P fertilizers, resulting in strong limitation of crop yields by soil P availability in those countries 25 .

Trade did not modify soil P anthropogenic signatures
Given that the trade of agricultural products displaces annually 3 MtP yr −1 , equivalent to 27% of the P traded through mineral P fertilizers 6 , we would have expected trade to substantially modify the soil P anthropogenic signature of large importers such as China and the Netherlands. Yet we found that the trade of feed and food products had almost no impact on the anthropogenic signatures of the soil P pools. The reasons for that were twofold. First, for 86% of the countries studied, the soil P inputs derived from imported feed and food products (and applied to soils as animal manure and sewage sludge) contributed to less than 5% of the cumulative P inputs over the 1950-2017 period (Fig. 4). Second, for a handful of countries, although the P embedded in imported feed and food products represented an important, indirect P input to their Article https://doi.org/10.1038/s41561-022-01092-0 soils, the anthropogenic signature of these imported products was generally close to that of their domestic agricultural soils (Extended Data Fig. 4), thus resulting in a neutral effect on the soil P anthropogenic signature (Extended Data Fig. 5).
We nevertheless found that, for some countries (such as the Netherlands, Saudi Arabia, Israel, Belgium, Egypt, Lebanon, Portugal and Switzerland), imported feed and food products represented an important input of anthropogenic P to agricultural soils. In these countries, P in animal manure and sewage sludge derived from imported feed and food products accounted for 10-48% of their total soil P inputs over the 2007-2017 period, with around half being from anthropic origin. More generally, we found that, at the global scale, 54% of the P embedded in traded products was of anthropic origin; a value amounting to 1.6 MtP yr −1 , or 8% of the global use of mineral P fertilizers (Methods). Those results illustrate that food and feed trade displaces large amounts of anthropogenic P that need to be accounted for when assessing the reliance of global food systems on phosphate rock resources 12,26 .

Discussion
Thanks to our large -scale and time-dependant modelling approach, our results helped to picture both the geographical variations and temporal evolutions of the countries' reliance on anthropogenic P. Overall, we found that anthropogenic P has contributed greatly to build the current P fertility of most agricultural soils, the average global anthropogenic signature of agricultural soil available P being around 47%. We also found that the discrepancies in anthropogenic signatures among countries were explained by both an uneven use of mineral P fertilizers since the Second World War and inequalities in the inherited agricultural soil P stocks in 1950. In some cases, the use of mineral P fertilizers has even exacerbated the existing inequalities in soil available P in 1950. For example, Western Europe was the region that had fertilized the most with mineral P over the 1950-2017 period while its soil available P stocks in 1950, inherited from biogeochemical background and past fertilization practices, were also among the highest in the world 3 . Our findings also illustrate that, on the other side of the gradient, many countries-especially in Africa-have benefited very little from phosphate rocks for sustaining their soil P fertility, the anthropogenic signature of their soils remaining below 32%. These results question the past distribution of phosphate rock resources, which has been based on regional economic and financial ability to buy fertilizers rather than on the search to optimize global food production 27 .
Complemented by proxies on current P fertilizer use as well as on soil P stocks 12 , our results show that world regions are facing differentiated risks of P deficits. Western Europe has high soil P stocks, but its high anthropogenic P signature shows its strong reliance on mineral P fertilizers, suggesting that future strategy should be based on recycling. By contrast, countries such as India and Brazil are particularly at risk as they use massive amounts of mineral P fertilizers to maintain high yields, but unlike European countries, their soil available P stocks are small 3 . Note that the anthropogenic signature we simulated is rather conservative and will decrease only slightly even if mineral P fertilizer use is drastically reduced. This is because increasing P recycling on agricultural soils affects both anthropogenic and natural P fluxes in similar proportion to that of the soil P anthropogenic signature. This makes our anthropogenic metric suitable to estimate the current reliance on anthropogenic P but not appropriate to analyse how soil P fertility would be affected under future farming scenarios.
Our work also offers a way to quantify the anthropogenic P content of agricultural food and feed products by considering not only contemporary soil P inputs but also the inherited soil P legacy from past anthropogenic P supply. Most studies estimating the P footprint of crop products have focused on the annual P supply to soils under the considered crop 6,7,12 . By labelling the anthropic versus natural origin of P in agricultural soils and in agricultural products, our methodology helps to get estimates of crop P footprints that account for past mineral P fertilization.
To conclude, the quantification of the anthropogenic signature of agricultural soil P sheds a new light on our understanding of the Anthropocene, which we analysed through the spectrum of the global P cycle. World regions have unequally benefited from the use of mineral P fertilizers, thus calling for different desirable trajectories for the future. In the Global North and in some Asian countries, the high anthropogenic signatures reflect their massive use of mineral P fertilizers over the 1950-2017 period but also sketch out the potential opportunity (depending on their soil P stocks) to engage an agro-ecological transition aiming at closing the P loop on a regional or national scale. On the contrary, Africa and to a lesser extent some South American countries have low anthropogenic signatures and thus soil P legacy virtually not reliant upon mineral P fertilizer. In such places, which are often facing high food insecurity issues, the potential for increasing yields through higher soil P inputs is among the highest in the world 27 . In this context, it would seem fair to reorient the distribution of phosphate rock resources towards these regions. Such management would require an unprecedented global coordination. Like many other issues, it raises the question of environmental justice: in a finite world, how and to whom should resources be allocated? Standard deviation values result from the uncertainties in the three calibrated parameters (β, the maximum attainable yield without any P limitation; γ, the ability of crops to extract soil P; and μ SPtoLP , the fraction of P in the stable P pool that transfers annually to the labile P pool); n refers to the number of calibrated triplets per country.  For comparison, P harvested from agricultural soils has also been displayed (blue line). For clarity, we did not show the anthropic versus natural origins of P embedded in sludge (in yellow), although they have been considered in the model. Manure input was split into three categories: (1) manure that came from the consumption of domestically produced feed and forage, (2) manure that came from the consumption of imported feed and forage and (3) manure that came from the consumption of mineral P feed. Note that the y-axis scale differs among countries.

Modelling framework
Our model aims at simulating P fluxes within the agricultural system of each country, as well as P flows between countries through the trade of agricultural products (Fig. 1). Therefore, we simulated not only soil P stocks and dynamics but also the fluxes of P embedded in crop and livestock products as well as P fluxes in animal manure and sewage sludge. The P dynamics in agricultural soils (without distinction between cropland and grassland soils) were simulated by using our own two-pool soil P model, inspired largely by that of Ringeval et al. 23 . Similar models have been used in previous large-scale studies and proved successful to reproduce observed P harvest over long-term periods [17][18][19] . For the specific purpose of our study, we distinguished the natural versus anthropic origin of each P flow and stock, following the methodology of a preliminary study conducted at the French scale 23 . Soil P pools were expressed in kgP ha −1 , while fluxes were all expressed in kgP ha −1 yr −1 . We assumed that for each country, the 0-30 cm soil horizon was relevant to study the soil P dynamics following fertilizer inputs, plant uptake by roots and erosion losses.
For each country, we modelled the size of the soil P pools using equations (1)-(3). We distinguished between labile P pools (LP) and stable P pools (SP) and between anthropic (subscript Ant) and natural (subscript Nat) origin of P within each pool.
Where y refers to the year considered, X stands for the subscripts Nat and Ant, and SL, OF, CF and LO correspond to soil P inputs as sewage sludge, animal manure, mineral P fertilizer and outputs as losses via erosion to water bodies, respectively; Ө refers to the bioavailable fraction of manure and sludge, which was assumed to equal 0.8 on the basis of reported values from the literature 18,23 ; T represents the net transfer of P from the stable to the labile P pool and captures soil P dynamics such as adsorption/desorption, precipitation/solubilization and organization/mineralization processes; P harvest refers to the P in the plant biomass that is harvested from both croplands and grasslands, both as grains for animal feed and human food consumption (either domestically consumed or exported to other countries) and as forage (domestically consumed). The P content of each P pool was also modified following annual net changes in agricultural areas. We assumed that agricultural lands increased over forested areas, the soil P anthropogenic signature of which was assumed to be 0% (see Supplementary Information section 5). P harvest . For each country, P harvest,Tot of crop and forage products was annually simulated on the basis of equation (4). Where β (in kgP ha −1 yr −1 ) refers to the maximum attainable yield without any P limitation, and γ (in ha kg −1 P) depicts the ability of crops to extract soil P; P harvest,Tot and LP Tot refer, respectively, to P exported from agricultural soils either for animal or human consumption (in kg ha −1 yr −1 P) and to the size of the labile P pool (kgP ha −1 ). The subscript 'Tot' refers to the sum of anthropogenic and natural P. A calibration procedure (Calibration procedure and uncertainty analysis) was performed to determine the values of β and γ for each country.
P transfers from the stable to the labile P pool. The net flux of P from the stable to the labile P pool was also explicitly simulated, as described in equation (5).
Where μ SPtoLP (yr −1 ) represents the fraction of P in the stable P pool that transfers annually to the labile P pool. The inverse of μ SPtoLP can also be interpreted as a mean residence time of P in the stable P pool. For the labile P pool, a turnover rate of 0.2 yr −1 -corresponding to a mean residence time of 5 yr-was used as proposed by the literature 18 . Similarly to β and γ, μ SPtoLP was also calibrated for each country, thus counterbalancing the 0.2 fixed rate of transfer from labile to stable P pool. The calibration of μ SPtoLP aimed at considering the soil P dynamic specificities of each country.
Input data. Input data were collected for each individual country. Mineral P fertilizer (CF) data were collected from FAOSTAT 29 for the 1961-2017 period. Data for the 1950-1960 period were estimated from the spatially explicit dataset of Lu and Tian 15 by aggregating values per country. We estimated that, overall, almost 1 billion ton anthropogenic P have reached agricultural soils over the 1950-2017 period. This corresponds to an average global use of 19.5 Mt anthropogenic P annually (calculated over the past 10 yr), well within the range of previous studies 16 . Manure (OF) data were estimated by considering seven animal categories (Supplementary Data Table 2) and by multiplying the number of heads 29-32 by species-dependent P excretion rates 33 . All the manure produced was assumed to reach agricultural soils, without further distinction between cropland and grassland soils. Losses (LO) of P from soils were calculated by multiplying country-specific percentages of soil lost by erosion 34 by the here-simulated P pool sizes. Because P losses from soils may apply to both labile and stable P pools indifferently, we assumed that the signature of P losses was similar to that of the overall soil P pools. Sludge P (SL) production was estimated by focusing on sewage sludge coming from human food consumption only (thereby excluding P release from detergents). We then estimated the fraction of P in sludge that would reach agricultural soils by multiplying each national human P excretion by the fraction of sewage sludge that is treated and by the P removal efficiency of treatment plants 35 . We assumed that Mineral P Feed (MF) represented 5% of the global production of phosphate rock 29 annually 36,37 . The global use of mineral P feed was then allocated to each country on the basis of both their total number of livestock heads and their consumption of mineral P fertilizer (Supplementary Information section 2.5). Harvested P data were used to calibrate our model. P harvest from grasslands was assumed to equal livestock forage P demand at the country scale, assuming no trade of forage between countries. Livestock forage P demand was estimated annually by multiplying livestock number of cattle, sheep and goats by their total energy requirement 38 , by the percentage of forage consumption in their diet 39 and by P concentration of forage. The P harvest from croplands was estimated by considering 31 crop species (Supplementary Data Table 5) and by collecting harvested data from FAOSTAT 29 for the 1961-2017 period. Crop yield data for the 1950-1960 period were estimated by assuming that crop-specific production varied linearly with the country human population as in Bouwman et al. 21 . Finally, crop P harvest was estimated by multiplying crop yields by the species-dependent P concentration, derived from the COMIFER 40 . The trade of agricultural products was considered by including 55 crop products, representing more than 85% of the P traded through food and feed at the global scale (Supplementary Data Table 6). From 1986 to 2017, we extracted the data regarding the imported quantities of feed and food products and the countries of origin from the FOATSAT Detailed Trade Matrix 29 . From 1961 to 1985, the imported quantities were derived from the FAOSTAT Food Balances matrix, and the countries of origin were assumed to be roughly the same as those Article https://doi.org/10.1038/s41561-022-01092-0 from the 1986-1991 period ( Supplementary Information section 4). For feed, imported quantities and their origin for the 1950-1960 period were estimated on the basis of the information at year 1961 and scaled down on the basis of the livestock density. For food, we assumed that the traded quantities were negligible for the 1950-1960 period. Note that the quantities and origin of imported food and feed products were used only to compute the anthropogenic signatures of animal manure and sludge, respectively.

Model initialization.
To initialize the model, we had to determine the size of each soil P pool in 1950. We hypothesized that before 1950 the application of mineral P fertilizer was negligible compared with the use of other fertilizers, mainly animal manure. As a result, we considered that LP Ant (1949) = SP Ant (1949) = 0. In 1950, the LP Ant value was set equal to the inputs of mineral P fertilizer that year. The SP Ant was null in 1950 because, by definition, all mineral P fertilizers are applied to the labile P pool. The size of the LP Tot in 1950 was calculated on the basis of the P harvest that year, using equation (4), possibly leading to significant effects of P harvest in 1950 on the initial soil P pool sizes. The size of LP Nat was then calculated as LP Tot minus LP Ant . Finally, on the basis of equilibrium between the labile and stable P pools, we determined SP Nat (equation (5) with T = 0).
Calibration procedure and uncertainty analysis. For each country, an independent calibration procedure was performed to estimate the values of β, γ and μ SPtoLP parameters. The general idea was to (1) test many triplets of parameters, (2) select the ones that allow the simulations to match P harvest time-series available data ( Supplementary  Information section 2.2) and labile P pool size available data in 2005 ( Supplementary Information section 2.7) and (3) use these sets of selected triplets (also called calibrated triplets) to assess the uncertainty in anthropogenic signatures. The adequacy between the simulation outputs and the available data were quantified through calibration scores (equations (6) and (7)). On the basis of these scores, we selected for each country a set of 'calibrated' triplets (Supplementary Information section 6.2). For each calibrated triplet, we ran our model. The set of obtained anthropogenic signatures was then used to compute a mean and a standard deviation value for each country.
Calibration started first by choosing the range of values to be tested for each parameter (Supplementary Information section 6.1). Note that for β, the ranges were country dependent. Then, from all possible values, we produced triplets ( β i , γ i , μ SPtoL,i ) for which we ran the model and computed score 1 and score 2, as described in equations (6) and (7). From 1,500 to 8,000 triplets were tested for each country.
where N refers to the number of years studied (N = 67), P harvest,obs refers to the mean value of P harvest available data for the period studied, P harvest,obs and P harvest,sim refer to the P harvest time-series available and simulated values and LP obs and LP sim refer to labile P pool size available and simulated values at year 2005; P harvest,obs was calculated from international databases (Supplementary equation (3)), and P harvest,sim was calculated with equation (4). Score 1 calculates the normalized root mean square error between the simulated total P harvest and the P harvest available data over the 1950-2017 period. Score 2 computes the error between the total labile P pool simulated at year 2005 and the total labile P pool available data extracted either from Ringeval et al. 3 or from observations representative of country-scale values for the few countries where observed soil P values were available (Supplementary Information section 6.3). In the main text, only results from the calibration performed with labile P pool available data from Ringeval et al. 3 were presented.
On the basis of these scores, calibrated triplets were selected for each country. When possible, all triplets that led to score 1 + score 2 < 30% were selected, and the corresponding countries were classified in category 1 (62% of countries). This category thus encompassed countries with good fit of model outputs with available data on both historical P harvest and the size of the labile P pool in 2005. For some countries, no triplets were able to match this condition. In this case and when possible, we selected the triplets for which score 1 < 30%. The corresponding countries were classified in category 2 (29% of countries). For such countries, the parameters were thus constrained only by P harvest. Finally, when no triplets led to score 1 < 30%, the corresponding countries were classified in category 3 (7% of countries) and not further studied.
In an effort to calibrate the model with observed data on soil available P pool instead of simulated values from Ringeval et al. 3 , we collected data on soil available P measurements (whatever the chemical extraction used) conducted on agricultural soils (Supplementary Information section 2.7). We found enough data for 24 countries, including European countries 41 , the United States and Canada 42 , Tanzania 43 , Botswana and Yemen 44 . For these countries, we found that using the measured data to calibrate the model instead of the simulated data did not significantly change the anthropogenic signatures of the soil P pools (Supplementary Information section 6.3).

Data analysis
For each model output (size and anthropogenic signatures of the fluxes and soil P pools), we computed a mean and a standard deviation value, reflecting uncertainties in the three calibrated parameters. The global and regional averages of the soil P anthropogenic signatures are mean values weighted by the agricultural area of each country. We also presented detailed results for eight contrasted countries (France, the Netherlands, the United States, Brazil, India, China, Zimbabwe and Morocco) that all belong to category 1. The rationale was to better appreciate the temporal evolutions in soil P inputs and in soil P anthropogenic signatures. We displayed only results for the anthropogenic signatures of the labile P pools since those of the stable P pools were very similar.

Limitations of the work
We quantified only the uncertainties in the soil P anthropogenic signatures related to uncertainties in the three calibrated parameters.
Although not quantified here, other uncertainties and limitations exist, in particular resulting from model structure, related assumptions and input data. We discussed the limitations of the model in the Supplementary Information (section 10).
Regarding input data to the model, we recognize large uncertainties in our estimates of P harvest from grasslands, especially in some countries such as India 45 . We explicitly chose to base our estimates of P harvest from grasslands on forage demand and not on forage production given the large uncertainty in net primary productivity of grassland 46 . Data on the trade of agricultural products were used to determine the anthropogenic signature of animal manure and sewage sludge, when animals and humans consumed imported products ( Supplementary Information sections 3 and 4). For several reasons, we underestimated the global P flows that occurred through the trade of agricultural products. However, the effects on the soil P anthropogenic signatures are likely to be small (Supplementary Information section 10.2). We also hypothesized that all the animal manure produced was available for agricultural soil application, which is highly questionable. Other studies have used global average, but those were too uncertain to be included in our estimate 14 . However, because manure application Article https://doi.org/10.1038/s41561-022-01092-0 on agricultural soils is basically P recycling that conserves anthropogenic signature, the possible overestimation of manure application in our approach is likely to have little effect on the estimated anthropogenic soil P signatures. Finally, in our approach, we did not distinguish between cropland and grassland soils. This avoided speculative assumptions on the split of manure and mineral fertilizer between cropland and grassland soils. Yet it is likely that distinguishing between these two land uses would have given higher anthropogenic soil P signatures for croplands and smaller signatures for grassland soils.
For some countries, we did not find any triplet of parameters that allowed the model to replicate the historical available data on soil P harvest. Therefore, these countries (category 3) were excluded from our analyses. The countries in category 2 were not able to match the constraints on both P harvest and labile P pool size in 2005. For these countries, located mainly in Africa, we recognize strong uncertainties in our estimates of the size of soil P pools and thus in the estimated anthropogenic signatures ( Supplementary Information section 6.2).
Another limitation comes from neglecting mineral P fertilizer use before 1950. Although mineral P fertilizer use was small for most countries in 1950 21 , several European countries already had high fertilization rates at that period. We nevertheless chose to start the simulations in 1950 because of the critical lack of accurate global data before that date. A sensitivity test performed by Ringeval et al. 23 -on which our methodology is based-has shown that our estimates are likely to remain robust. They have shown that running the model for 30 years with soil inputs and outputs related to the first year (~1950) before the transient run showed only an ~10% underestimation of the anthropogenic signature of agricultural soils.
Finally, our definition of anthropogenic P-which referred to the use of mineral P fertilizers and mineral P feed derived from the extraction and chemical treatments of phosphate rocks-is questionable. 'Anthropogenic' should more broadly refer to human practices altering the P cycle. This includes, among other things, application of human manure 47 and guano 48 to agricultural soils and P fertility transfers from grassland to cropland soils. However, in a context of increasing scarcity of phosphate rock reserves, we found it more relevant to quantify the modification of the P cycle related solely to the use of this critical and highly processed resource.

Data availability
All data that support the findings of this study are archived on https://data.inrae.fr/privateurl.xhtml?token=4ddb8501-c41d-4ad6-8a09-5c1ebb8289f9. FAO national statistical data were obtained from https://www.fao.org/faostat/en/. Land-use data from HYDE 3.2 were downloaded at https://easy.dans.knaw.nl/ui/datasets/id/ easy-dataset:74467/tab/2. Observed P-Olsen values for agricultural European soils were requested at https://esdac.jrc.ec.europa.eu/ content/lucas2015-topsoil-data. Observed soil available P values for American and Canadian soils were downloaded from the soil test summary website of the Fertilizer Institute (https://soiltest.tfi.org/ tables). For Tanzania, observed data of agricultural soil available P were obtained at https://registry.opendata.aws/afsis/. For Botswana and Yemen, observed data of agricultural soil available P were obtained from the ISRIC World Soil Information website (https://www.isric.org/ documents/document-type/isric-report-201006-inventory-p-olsen -data-isric-wise-soil-database-use). Information on agricultural soil available P from Ringeval et al. 3 was obtained following a request to the authors. Information on the use of mineral P fertilizer was obtained at https://doi.pangaea.de/10.1594/PANGAEA.863323. Source data are provided with this paper. Except for sludge, the figure differentiates anthropic vs. natural origin of phosphorus (P) inputs. The flux of P in manure (OF) was split in 3 categories: (i) P in manure that originated from the consumption of domestically produced feed (ii) P in manure that originated from the consumption of imported feed and (iii) P in manure that originated from the consumption of mineral feed. For clarity, we did not show the anthropic vs. natural origin of phosphorus embedded in sludge, although they have been considered in the model. Fig. 4 | Anthropogenic signatures of imported feed, food and of the labile phosphorus (P) pool over the 1950-2017 period. Data on the anthropogenic signature of the imported feed were only displayed for the countries in which the contribution of soil P inputs coming from imported feed over the 2007-2017 period represented more than 5% of their total cumulative soil P inputs. Data on the anthropogenic signature of the imported food were only displayed for the countries in which the contribution of soil P inputs coming from imported food over the 2007-2017 period represented more than 4% of their total cumulative soil P inputs. For some years no data are displayed because the country studied did not import any feed nor food. Data on the anthropogenic signature of the labile pool (blue) are presented as mean (line) values ± SD (blurred zone). The n set of anthropogenic signature for each country is derived from running the model with all calibrated triplets.

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01092-0 Extended Data Fig. 5 | Trade effect on the anthropogenic signature of the soil labile phosphorus (P) pool. Each figure displays the anthropogenic signature of the soil labile P pool over the 1950-2017 period, with and without trade of agricultural products. Data are presented as mean (lines) values ± SD (blurred zone). The n set of anthropogenic signature for each country is derived from running the model with all calibrated triplets. In green: results from the main computations, where the trade is considered. In purple: results from computations where the anthropogenic signature of all imported food and feed products was set to that of the importing country. Data are shown only for the countries where P embedded in imported food and feed products contributed to more than 5% of total soil P inputs over the 1950-2017 period.