Land use intensification increasingly drives the spatiotemporal patterns of the global human appropriation of net primary production in the last century

Land use has greatly transformed Earth's surface. While spatial reconstructions of how the extent of land cover and land‐use types have changed during the last century are available, much less information exists about changes in land‐use intensity. In particular, global reconstructions that consistently cover land‐use intensity across land‐use types and ecosystems are missing. We, therefore, lack understanding of how changes in land‐use intensity interfere with the natural processes in land systems. To address this research gap, we map land‐cover and land‐use intensity changes between 1910 and 2010 for 9 points in time. We rely on the indicator framework of human appropriation of net primary production (HANPP) to quantify and map land‐use‐induced alterations of the carbon flows in ecosystems. We find that, while at the global aggregate level HANPP growth slowed down during the century, the spatial dynamics of changes in HANPP were increasing, with the highest change rates observed in the most recent past. Across all biomes, the importance of changes in land‐use areas has declined, with the exception of the tropical biomes. In contrast, increases in land‐use intensity became the most important driver of HANPP across all biomes and settings. We conducted uncertainty analyses by modulating input data and assumptions, which indicate that the spatial patterns of land use and potential net primary production are the most critical factors, while spatial allocation rules and uncertainties in overall harvest values play a smaller role. Highlighting the increasing role of land‐use intensity compared to changes in the areal extent of land uses, our study supports calls for better integration of the intensity dimension into global analyses and models. On top of that, we provide important empirical input for further analyses of the sustainability of the global land system.


| INTRODUC TI ON
During the 20th century, industrialization and globalization processes have greatly increased humanity's physical imprint on the planet (Arneth et al., 2019;Krausmann et al., 2013;Kummu et al., 2010;Steffen et al., 2015). With the industrial intensification of land use, yields could be raised by crop breeding and increasing external inputs to agricultural production, such as fertilizer, pesticides, and mechanical power, a process often described as the green revolution (Pellegrini & Fernández, 2018). These trends have driven land system trajectories in many locations in the past decades, and they have also shaped the pressures land use exerts on ecosystems and ecosystem processes IPCC, 2019), including biodiversity loss (IPBES, 2019;Newbold et al., 2015;Pereira et al., 2012).
Global analyses of the effects of land use across spatial units require spatial information on the evolution of the extent of land use, but also of land-use intensity. At the global level, the time series of maps of the areal extent of different land-use types are available (Ellis et al., 2021;Hurtt et al., 2020;Klein Goldewijk et al., 2017). In contrast, maps on land-use intensity are available only for selected land uses. Three dimensions of land-use intensity have been discerned recently Kuemmerle et al., 2013): (i) input intensity, referring to the amount of (agricultural) inputs per unit of land (e.g., fertilizer application rates), (ii) output intensity, referring to the amount of biomass outputs from the land-use activities (e.g., crop or forestry yields), and (iii) system-level intensity, i.e., the extent to which land-use alters properties of land systems such as water availability or ecosystem stocks and flows. At the global level, spatial data on the evolution of (i) and (ii) are available; however, they are largely limited to croplands, respectively, specific crops. Consistent maps of land-use intensity dynamics across land-use types are lacking.
We contribute to filling this research gap by producing consistent global maps of the spatial evolution of land use and land-use intensity during the last century (1910 to 2010). We employ the framework of the human appropriation of net primary production (HANPP), which measures the scope of human intervention into ecological carbon flows by quantifying the difference between the potential NPP and the NPP prevailing in an ecosystem after harvest (Haberl et al., 2014). HANPP provides a systemic perspective on land use that allows for analyses across all land uses, such as growing different types of crops, livestock grazing, and forestry, in a given spatial unit. HANPP has previously been mapped at the global level for the year 2000 (Haberl et al., 2007), with information on its temporal evolution available at the aggregate level of world regions and for selected countries (Chen et al., 2015;deSouza & Malhi, 2018;Kastner, 2009;Krausmann et al., 2013), and the local level (Niedertscheider et al., 2017;Pritchard et al., 2018;Qin et al., 2021).
Advancing established accounting rules and available data, we present a spatially explicit reconstruction of the temporal evolution of HANPP during the last century , including the following components of HANPP: (1) the NPP associated with society's biomass harvest (HANPP harv ), including NPP lost due to deforestation and (2) the change in NPP induced by land use (HANPP luc ). In addition, we explore the role of changes in the potential NPP (NPP pot ), i.e., the NPP that would prevail in the absence of land use but with contemporary climate, which forms the baseline of measuring human impacts on NPP flows.
Applying the HANPP framework for time-series analyses is well suited to investigate changes in land-use intensity, as it covers two of the three land-use intensity dimensions introduced above, namely output intensity (via HANPP harv per unit land) and a system-level metric, i.e., the effect of land use on overall NPP flows in ecosystems (via HANPP as a percentage of NPP pot ). In addition, the HANPP framework is consistently applicable across different land-use types (such as cropping, livestock grazing, and forestry) and its results are expressed in units that can straightforwardly be linked to ecological and socio-economic processes. Next to investigating changes in land-use intensity, the HANPP framework can also quantify the impact of other major drivers of land system change, namely the effects of e.g., expansion or contraction of areas under human use, and the impact of changes in the efficiency of societal biomass appropriation (e.g., via asking how much biomass harvest is generated per unit of overall HANPP). Using NPP pot as the reference level, the framework is well suited to explore how changes in potential vegetation productivity through environmental change have affected the dynamics of the indicator. In addition, spatially explicit maps of the temporal evolution of HANPP allow going beyond analyzing national-level trends and enable assessing how HANPP has changed across any socio-economically or biogeographically defined region.
We address the following research questions: How have spatial patterns of HANPP and its components evolved during the 20th century? How did the patterns change across the world's major biomes? What were the proximate drivers of these changes, such as the role of land-use expansion, specific elements of land-use intensity, and changes in natural productivity? Finally, how large are the uncertainties of the spatial evolution of HANPP, and which set of input data and assumptions contributes most to these uncertainties? 2 | MATERIAL S AND ME THODS Following Haberl et al. (2014), HANPP is defined as the amount of NPP appropriated by societal activities in a given year through harvest (HANPP harv ) or land-use change (HANPP luc ): (1) HANPP = HANPP harv + HANPP luc K E Y W O R D S global assessment, historical reconstruction, human appropriation of NPP, land-use change, land-use intensity, net primary production HANPP harv is the aggregate of biomass harvested that enters socioeconomic processes (including biomass grazed by livestock) and biomass compartments that are destroyed during harvest but do not enter socioeconomic processes such as unused crop and forestry residues or biomass destroyed through deforestation. Due to high uncertainties in the deforestation term, we present results for this part separately and refer to this compartment as HANPP defo and label the sum of all other HANPP harv flows HANPP harv* , i.e., HANPP luc denotes the difference between the NPP of the potential natural vegetation (NPP pot ) and the NPP of the vegetation prevailing in the studied year (NPP act ). HANPP luc can be positive, if land-use practices lower ecosystem productivity, or negative in cases where they increase ecosystem productivity, for instance, through irrigation in drylands, and is calculated as: We assessed HANPP as the annual flow of carbon and as a percentage of NPP pot . We calculated HANPP components at a resolution of five arcminutes for the years 1910, 1930, 1950, 1960, 1970, 1980, 1990, 2000, and 2010. To tackle the uncertainties in the spatial distribution of HANPP, we modulated six central input datasets and assumptions (Table 1). This way, we calculated a total of 216 global maps of the HANPP indicator for each year studied.
The starting point for the calculations was the compilation of consistent land-use reconstructions for each year, discerning six major land-use types. Each grid cell contains information on the fractional cover, i.e., the share of the respective land use in the total land area of the grid cell, of the discerned land uses. For this, we used state-of-the-art reconstructions of historical trends in land use and land cover (Klein Goldewijk et al., 2017;Ramankutty & Foley, 1999;Ramankutty, personal communication, 2018) and processed them in order to achieve full fractional cover for the sum of all land-use layers. This way, summing up across land uses will consistently equal the total land area of each grid cell. We prioritized the different input maps and reduced fractional cover of the land uses with lower priority, in cases where inconsistencies between data sets exist, i.e., summing up the different land uses would exceed the available land area.
A detailed account of how the land use data were constructed can be found in Supplementary Method Note S1. Using different input data on the distribution of cropland areas, we arrived at three modulations for the land-use reconstructions used in our assessments ( Table 1).
The resulting land-use datasets distinguish six land-use categories at a resolution of five arcminutes: cropland, grazing land, forestry land (i.e., forests that we considered to be used for forestry), infrastructure areas, wilderness areas, and non-productive lands.
For each of these land-use categories, we separately calculated and mapped HANPP and its components and arrived at the overall HANPP by summing across each pixel of these maps. NPP pot served as a reference level from which the human appropriation of NPP was calculated. We considered three modulations of this input dataset ( Table 1). Two of them were derived from the global dynamic vegetation model (DGVM) LPJ-GUESS  and one from the DGVM JSBACH (Lasslop et al., 2020). In all cases, we ran the models without land-use information but with historical climate forcing for the years 1901 to 2015. LPJ-GUESS version 4.0.1 was used in its standard configuration and forced by the CRU-NCEP (Harris et al., 2014;Viovy, 2018) climate data aggregated from 6-hourly to monthly fields. In addition, we performed an LPJ-GUESS simulation with nitrogen limitation disabled to ensure consistency with the previous studies (e.g., Haberl et al., 2007;Krausmann et al., 2013), and comparing these different model runs enabled us to identify the diverging spatial patterns and temporal trends of NPP induced by nitrogen limitation Smith et al., 2014;Wårlind et al., 2014). We used the model output of NPP per year and unit area in terms of grams of carbon per m 2 and year, which was available at a spatial resolution of 30 arcminutes. We built decadal averages of this output around the years included in our assessment, downscaled the resulting maps to five arcminutes to match the resolution of the land-use data, using bilinear interpolation, and removed values below zero. To avoid biases, we only used the two LPJ-GUESS runs in the presentation of the results and included the JSBACH model only in the uncertainty assessment. This is because the NPP values derived by JSBACH are at the very high end of global estimates of NPP, whereas the LPJ-GUESS runs are positioned at the center of the current estimates (Ito, 2011;Yu et al., 2018).
The land-use and NPP pot datasets provided the starting point for downscaling data on HANPP harv and HANPP luc from the national level to the grid. For the first three years in our study (1910, 1930, and 1950), we relied on HANPP harv reconstructions from a previous study (Krausmann et al., 2013) Table 1 summarizes the modulations performed in this study. In total, we calculated 216 versions of HANPP for each year investigated.

| Uncertainty assessment
Note that the main results are based on 144 of this 216 version only (see above). We used all 216 versions to perform a comprehensive assessment of spatial uncertainty and to identify which of the (2) HANPP harv = HANPP harv* + HANPP defo (3) HANPP luc = NPP pot − NPP act TA B L E 1 Overview of the different modulations performed in the study to assess the uncertainty in spatial patterns

Parameter Rationale for modulating input parameters Specific modulations performed
Land use High uncertainties exist in the historical extent of land use. We explore this using different layers for the extent of croplands, which also affects the estimates for the extent of the other land use types (Supplementary Note S1) Cropland reconstruction from HYDE 3.2 (Klein Goldewijk et al., 2017) Cropland reconstruction from Ramankutty and Foley (1999;Ramankutty, personal communication 2018) A combination of the two above, calibrated to match the totals for cropland extent in Krausmann et al. (2013) Productivity of the potential natural vegetation NPP pot NPP pot is the baseline for the HANPP framework. Different global vegetation models provide varying estimates for NPP pot , both in the overall magnitude and in spatial and temporal patterns Global vegetation model LPJ-GUESS with nitrogen limitation of NPP  Global vegetation model LPJ-GUESS without nitrogen limitation of NPP  Global vegetation model JS BACH (Lasslop et al., 2020) (included for the uncertainty assessment only)

Amount of biomass harvest on croplands
Estimates of harvest volumes at the national level are subject to uncertainties. To cover these uncertainties, we rely on published estimates on ranges of HANPP harv* (Krausmann et al., 2018) Low estimate

Medium estimate
High estimate

Allocation of biomass harvest on croplands
Different allocation rules were explored, as historical data on spatial distribution of crop harvest are not available Harvest follows patterns of NPP pot , in line with the previous work (Haberl et al., 2007) Harvest follows patterns of NPP pot and synthetic nitrogen fertilizer inputs (Nishina et al., 2017) Allocation of biomass harvest on grazing lands Different allocation rules were explored, due to large uncertainties on spatial distribution of grazing pressures (Fetzel et al., 2017) Grazing pressure is more concentrated in regions of higher productivity Grazing pressure is more evenly spread out across grazing lands

Extent of land degradation in drylands
The effect of land degradation on NPP is highly uncertain and we used published estimates (Zika & Erb, 2009)  National level livestock feed demand is calculated based on data on livestock populations and feed requirements. Feed availability from non-grazing-land sources is estimated based on feed use statistics and assumptions on non-market feed sources. The difference between feed demand and feed availability from non-grazing land sources is assumed to be HANPP harv* on grazing lands. HANPP luc is assessed based on information on potential vegetation cover of grazing lands and on land degradation in dryland areas. Spatial allocation is based modelling NPP available for grazing per pixel Forestry land (FAO, 2020;Krausmann et al., 2013;Schulze et al., 2012) National level wood harvest statistics are extrapolated to estimate HANPP harv* . Spatial allocation follows NPP pot patterns and patterns of fertilizer application. Spatial allocation is based on NPP pot patterns. HANPP luc is assumed to be zero.
For details please refer to Supplementary Method Note S2. modulated input factors contributed most to the observed variance.
We did this by calculating the variance between modulations of a single parameter and comparing it to the overall variance in each pixel.
We assumed a linear model between the values of HANPP and the different parameter modulations and used an ANOVA to quantify individual and overall sums of squares (SSq). The ratio between the parameter SSq and the overall SSq gave the variance contribution of the parameter. This variance decomposition was performed for every pixel, which allowed us to investigate the spatial differences in the most critical parameters. In addition, we provide additional maps showing the standard deviation of the components of the HANPP framework in the supplementary material ( Figures S3 and S5).

| Decomposing HANPP to investigate the contribution of major land-system change drivers across biomes
To quantify how different drivers contribute to the land-system change throughout time, we employed the additive implementation of the log mean divisia index (LMDI) decomposition approach (Ang, 2005). A mathematical formalization of this approach and an elaboration of its strengths compared to the other index decomposition approaches can be found in Ang (2005). In index decomposition approaches, the variable of interest (in our case HANPP as a share of NPP pot ), is split into meaningful terms, constructing a mathematical identity. The method then quantifies the contribution of each of these terms to the change in the variable of interest. Summing up the individual contributions will give the change in the variable of interest. We employed an approach that decomposes changes in HANPP through the following identity, modified from Gingrich et al. (2015): We performed a decomposition analysis at the level of seven ecological biomes (b, for definitions, see Supplementary Method Note S4) and for the global total, to assess the changes in HANPP% (HANPP as a share of NPP pot ) for each land-use type, lu, individually (i.e., infrastructure land, cropland, grazing land, forestry land, and wilderness). Summing up across land-use types, this method enables us to quantify the following drivers' contribution to overall changes in HANPP%: • Changes in area mix, A, expressed in km 2 . Agricultural expansion typically results in changes toward more HANPP-intensive land uses and will drive HANPP% upwards, contraction of intensely use lands, e.g., cropland abandonment, will decrease HANPP%; • Changes in harvest intensity, HANPP harv /A, expressed in tons of carbon per year and km 2 , i.e., the amount of harvested biomass per unit area. Increases here will increase HANPP%. This driver represents an indicator of the output intensity of land use.
• Changes in HANPP intensity, HANPP/HANPP harv , expressed in tons of carbon per year and km 2 , i.e., HANPP per unit harvested biomass; this value is high if HANPP luc is high, i.e., a large part of HANPP is not harvested biomass. Decreases in HANPP luc will decrease HANPP intensity. Therefore, this measure indicates the ratio between the system-level intensity and output intensity and can be interpreted as (the inverse) of a measure for the efficiency of societal biomass appropriation.
• Changes in natural productivity, here considered as 1/NPP pot , expressed in tons of carbon per year and km 2 , i.e., the inverse of NPP pot ; if NPP pot increases due to changing environmental conditions (e.g., longer vegetation growing seasons in the north or increased CO 2 levels in the atmosphere) HANPP% will decrease. (4) HANPP, its components, and NPPpot for 1910, 1950, 1990, and 2010  Decomposing our results sheds light on how major determinants drove aggregate changes in HANPP (Figure 4). The factor that contributed most to increases in HANPP during the 20th century was an increase in harvest intensity, i.e., HANPP harv per unit area, pointing to the major impact of increasing output-intensity of land use on HANPP dynamics. This effect was the highest on croplands but also substantial on grazing land and forestry areas. The harvest intensity effect was markedly larger in the two more recent periods in all biomes, highlighting how agricultural intensification has become a dominant land-system driver around the world since the green revolution. While the scale of this effect differed across biomes and the three periods, it was the main driver of increased HANPP values across all biomes.

| Presentation of the results
While being the major gross driver, the harvest intensity effect was countered by changes in HANPP intensity, implying a reduction in HANPP luc per unit HANPP harv , partly offsetting the effect of higher harvest per unit area, when looking at overall HANPP levels. These counteracting trends emerge because increasing output intensity diminished the difference between NPP pot and NPP act , especially on croplands. The net effect of changes in harvest intensity and HANPP intensity was universally positive (Figure 4).

| DISCUSS ION
The key finding of our study is that dynamics of land-use intensifi-  Table 2) contributing strongest to the observed variance in the respective pixel (i.e., the contribution of the factor is more than 50%). The crucial relevance of land-use intensity for understanding land-use patterns and their implications is increasingly acknowledged in the literature (e.g., Beckmann et al., 2019;Dullinger et al., 2021;Felipe-Lucia et al., 2020;Friedlingstein et al., 2020;Hong et al., 2021;Kuemmerle et al., 2013). Our results further corroborate these statements and clearly show that without appropriate consideration of the dynamics of land-use intensity, assessments, and models, such as integrated assessment models, will increasingly miss the major changes in land systems and might result in biased and incomplete interpretations Pongratz et al., 2018). For instance, as net rates of land-use change have declined in many parts of the world, ideas like "peak land use," i.e., the outlook that pressures on land in terms of used areas stabilize and eventually decline in the near future (Ausubel et al., 2013;Seppelt et al., 2014), have to be critically reflected upon when considering the increasing role land-use intensity changes play in driving overall land system change.
One consequence of this trend toward an increased dynamic of land-use intensity and a decreasing role of agricultural expansion in many world regions is that HANPP, expressed as the share of NPP pot , has peaked or flattened out in the temperate, sub-tropical, and boreal biomes, where increases in harvest intensity have largely been offset by decreases in HANPP intensity (Figure 4).
At the same time, we find strong continued growth of HANPP in the tropical regions, where agricultural expansion continues and where HANPP luc levels remain substantial ( Figure 1) and have increased strongly in some regions in the more recent past, such as in tropical Africa ( Figure 2). However, for the forest biomes, the level of HANPP in the tropics is still considerably lower than for the sub-tropical and temperate biomes, implying that the overall HANPP rates are becoming increasingly similar if these trends continue. In this context, it is important to note that the increases in pressures on terrestrial ecosystems in tropical regions such as the Amazon, the Cerrado, or Southeast Asia are driven, in part, by the consumption in world regions that have enjoyed recent declines in land-use pressures (Hoang & Kanemoto, 2021;Pendrill et al., 2019). Together with a situation where global conservation interest is also concentrated in tropical regions (Dinerstein et al., 2020;Strassburg et al., 2020), this highlights that future land-use trajectories will have to deal with competing interests and will be increasingly driven by global demand and power relations (e.g., Krausmann & Langthaler, 2019).
The trends of the ever-increasing importance of intensification also raise questions about potential limits to intensification processes. To what extent will it be possible to further increase biomass harvests per unit area and at what costs? On the one hand, demand for land-based products is expected to further increase: food and feed are needed to feed a growing and increasingly affluent global population (Tilman et al., 2011), and increased use of biomass for energy and material use is envisaged in bioeconomy strategies around the world as a central measure to curb dependence on fossil fuels (Dietz et al., 2018). On the other hand, negative consequences of intensification such as detrimental effects of nutrient leaching and pesticide use on biodiversity are becoming more and more evident (Conijn et al., 2018;Glibert, 2017;Sharma et al., 2020). Landsystems that consistently show NPP values higher than NPP pot (i.e., negative HANPP luc ), as is currently the case in, e.g., eastern China ( Figure 1), are typically associated with considerable environmental externalities (West et al., 2014). These trends have led to calls for a more "sustainable" or "agroecological intensification" (Garnett et al., 2013;Thomson et al., 2019;Wezel et al., 2015), but it remains to be seen how such approaches can be operationalized at large scales, without increasing areal expansion of agriculture due to lower outputs per unit of land. Another potential way to reduce pressures on land systems that is increasingly explored relates to lowering the overall demand for land-based products, such as transitions to less resource-intensive diets, including the reduction of food waste (Hallström et al., 2015;Springmann et al., 2018;Tilman & Clark, 2014;Willett et al., 2019).  (Yu et al., 2018), in line with the modeled data. With regard to the robustness of our findings, it is also noteworthy that, for larger spatial scales, Krausmann et al. (2013) have highlighted that the HANPP framework is able to provide robust insights, even in the light of uncertain NPP pot trends.
The modeled and observed increases in natural vegetation productivity have been mainly attributed to increases in atmospheric CO 2 concentrations and longer growing seasons in high latitudes as a result of rising temperature levels, which, in turn, have been attributed largely to anthropogenic greenhouse gas emissions (Lucht et al., 2002;Pachauri et al., 2014;Walker et al., 2021;Zhu et al., 2016). In the context of our assessment, it is interesting to see that increases in NPP pot have occurred simultaneously with strong increases in harvest intensity in many parts of the world. Our estimates for harvest intensity are independent of the used DGVM output for NPP pot data and rely on agricultural and forestry statistics. The observed alignment of these trends calls for a better understanding of how much of the continued increases in output per unit area can be attributed to agronomical practices and how strong the role of such more indirect factors has been, is discussed in the literature on the effects of climate change on crop yields (e.g., Rosenzweig et al., 2014;Wang et al., 2018) and forest productivity (e.g., Pretzsch et al., 2014). Insights into this question are important, as it is far from clear that trends in increasing NPP levels will continue in a similar fashion in the future and dynamics will likely differ substantially in different biomes (Jong et al., 2012;Walker et al., 2021;Wang et al., 2020). In addition, these findings call for a more nuanced representation of such global change processes. For instance, while shifts in natural vegetation are included in current global vegetation models, the present HANPP framework does not consider such shifts. It would be interesting to explore how such changes may affect land-use patterns, as, for instance, shifts from herbaceous to woody vegetation will translate into changes in land management practices and preferences (e.g., Linders et al., 2020), which in turn will imply different HANPP patterns.
Our analysis highlights the added value of a comprehensive perspective on land-use change processes, as offered by the HANPP framework. A detailed look into the dynamics of the individual HANPP components allowed disentangling processes to occurr simultaneously in global land systems. By looking systematically across the different components that determine the overall magnitude of HANPP, we were able to draw up a more complete picture of land-system change, focusing in particular on the role of land-use intensity. Furthermore, through modulating central input data, we are able to provide an assessment of spatial uncertainty patterns of HANPP. While we find relatively low uncertainties across large parts of the globe, we see that the underlying land-use data and data on NPP pot are central for the overall observed uncertainty patterns.
Considering that we only modulated cropland areas (Supplementary Method Note S1) for which alternative historical reconstructions were available, the role land-use data play in overall uncertainties is likely even underestimated in our data. For historical periods, these uncertainties will be hard to reduce. For the most recent past and current patterns, the availability of ever-improving satellite observations has the potential to reduce uncertainties related to land use and land cover data, although the reconciliation of ground-based observation (census statistics) with satellite-derived data remains challenging. Other modulated factors had a lower relevance at the global level but might still be important regionally.
We also find that distinct uncertainty patterns emerge, depending on the unit in which HANPP is expressed ( Figure 5) (Field et al., 1995;Ito, 2011;Roxburgh et al., 2005). Differences in temporal trends are mostly caused by uncertain plant-physiological effects of elevated atmospheric CO 2 and different model assumptions concerning nutrient limitation Wårlind et al., 2014). However, if expressed as the share of NPP pot , the contribution of the NPP pot to overall uncertainties at the global level is much smaller, highlighting that expressed in relative terms the indicator is more robust.
It is also interesting to note that for the year 2000 our estimates are by about 2-4 HANPP%-points lower than those presented in earlier studies (Haberl et al., 2007;Krausmann et al., 2013). As the main reason for this difference, we identified the fact that in our assessment HANPP defo and the effect of ecosystem degradation on HANPP luc was estimated via spatially explicit approaches, which yielded substantially lower estimates for the year 2000. The other components were well in line with earlier assessments, with minor deviations due to refinements in the treatment of permanent cropland and the additional data limitation that FAOSTAT no longer reports the harvest of dedicated fodder crops (see Supplementary Method Note S2).
Our spatially explicit reconstruction of the development of HANPP and its components across the major land uses can provide valuable inputs for future research. For instance, the maps allow for assessing trends across spatial entities that do not align with national borders. We have presented such an assessment across the world's major biomes, but other spatial units will be interesting to analyze as well. With HANPP being closely and mechanistically linked to biodiversity patterns and processes, it will be crucial to investigate how present global biodiversity patterns align with past land-use patterns to address the ongoing challenge of biodiversity loss. Such an analysis could scrutinize the idea that biodiversity patterns only adapt to land-use pressure with a time lag. Additionally, the maps could form a baseline for developing spatial scenarios on future landsystem changes that incorporate the land-use intensity dimension, as well as contributing (along with other information) to identifying areas where further intensification is likely possible with limited negative environmental consequences. Addressing this intensity dimension will be crucial if pressing sustainability challenges in the land system, such as curbing biodiversity loss, are to be tackled in a comprehensive manner.

ACK N OWLED G EM ENTS
We greatfully acknowledge financial support by the Deutsche

DATA AVA I L A B I L I T Y S TAT E M E N T
Complete data for all included years and model runs, along with highresolution versions of the maps presented in Figures 1 and 2