Relative effects of land conversion and land-use intensity on terrestrial vertebrate diversity

cSAR model

We used the numerical cSAR model¹⁶ to calculate native species loss of four taxonomic groups (mammals, amphibians, reptiles, birds) caused by 45 LU types that were mapped onto a reference 5 × 5 arcmin grid (we also call individual grid cells landscapes in the following) of the global land area excluding Greenland and Antarctica. Calculations were based on (a) gridded LU-intensity and LU-type information (see below), (b) effects of LU-intensity on species richness derived from recently published meta-analyses^5,21, and (c) information on species distributions and habitat affiliations from IUCN and Birdlife International databases^41,42. For presentation of results, we aggregated the calculated effects of the 45 LU-types into those of six broad LU-types (cropland (30 annual crop types); pastures (non-grassland converted to grassland); grazing land (natural/ near-natural areas with livestock grazing); builtup (sealed areas); plantations (11 permanent crop types plus timber plantations), and forests (natural/ near-natural forest under forestry); see Supplementary Data 2 for details).

In the below formulae, we use the following indices: g = taxonomic group, n = grid cell, b = broad LU-type. We calculated the total number of native species losses for each taxonomic group g and grid cell n as

Here, ({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{rm{pot}}}}}}}) is the potential species richness in pristine ecosystems, ({{{{{{rm{A}}}}}}}_{{{{{{rm{n}}}}}}}^{{{{{{rm{cur}}}}}}}) is the pristine ecosystem area where no LU occurs (in m²), ({{{{{{rm{A}}}}}}}_{{{{{{rm{n}}}}}}}^{{{{{{rm{pot}}}}}}}) is the grid cell’s terrestrial area (in m²), ({{{{{{rm{h}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) is the affinity parameter, ({{{{{{rm{A}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) is the area of the LU-type, and ({{{{{{rm{z}}}}}}}_{{{{{{rm{n}}}}}}}) is the grid cell’s SAR exponent taken from ref. ⁴³. The model’s components are described below.

Potential species richness ({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{rm{pot}}}}}}})

We defined potential species richness of a landscape as the number of species for which the area-of-habitat (AOH) under pre-human or pristine conditions overlap the landscape (here referred to as native species and native AOH). Following ref. ⁴⁴, we used range maps of all mammal, reptile and amphibia species provided by the IUCN⁴⁵ and bird species by Birdlife International⁴⁶ databases to calculate gridded species richness via, first, overlapping each species’ range polygons with a 5 × 5 arcmin reference raster, second, constraining the resulting list of species per raster cell to those adapted to the pristine ecosystem(s) of these raster cell as defined in ref. ²⁵, and, third, constraining the resulting list of species by each species’ elevational range, also provided by the IUCN⁴². Here, we are interested in the total historical range of extant species⁴⁶ and hence included all parts of the range where the species were indicated as (i) Extant, Probably Extant, Possibly Extinct, Extinct and Presence Uncertain, (ii) Native and Reintroduced, and (iii) Resident or present during the Breeding Season or the Non-breeding Season, in the cited data sources.

We first rasterized each species’ range polygons using the raster and fasterize packages in R⁴⁷. Second, for each terrestrial grid cell in our reference raster, we created a species list by extracting each species’ gridded range using the velox package in R. Third, we ascertained that each cell’s species list contained only terrestrial species by excluding species which exclusively have aquatic habitat affiliations. The species’ habitat affiliations were directly taken from the IUCN and Birdlife databases^42,46. Fourth, we removed species from this cell’s species list which, according to the IUCN, are not affiliated with that cell’s pristine ecosystem. We therefore manually assigned the habitats distinguished in the ICUN habitat affiliation scheme to one or several of the 14 broad ecosystem types distinguished and mapped in ref. ²⁵ (Supplementary Data 4). The maps in the referenced study “approximate the original extent of natural communities prior to major land-use change”⁴⁸ and, hence, represent pristine ecosystems or potential vegetation types. Fifth, we excluded species whose elevational range did not overlap the elevational range of the grid cell using the GMTED2010 dataset (www.usgs.gov). These refinement steps were taken because species’ range maps usually deliver coarse-scale extent of occurrence rather than AOH information⁴⁴. Finally, we counted the species identities in each grid cell as ({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{rm{pot}}}}}}}). The species lists created in this step were also used for later steps, referred to in the appropriate sections.

Areas of pristine ecosystems (({{{{{{rm{A}}}}}}}_{{{{{{rm{n}}}}}}}^{{{{{{rm{cur}}}}}}}),({{{{{{rm{A}}}}}}}_{{{{{{rm{n}}}}}}}^{{{{{{rm{pot}}}}}}})) and LU-types (({{{{{{rm{A}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}))

The potential pristine ecosystem area ({{{{{{rm{A}}}}}}}_{{{{{{rm{n}}}}}}}^{{{{{{rm{pot}}}}}}}) is defined as the cell’s entire terrestrial area (excluding water bodies as defined by the land mask of the HYDE 3.2.1 database⁴⁹). As the area of pristine ecosystems currently found in each grid cell (({{{{{{rm{A}}}}}}}_{{{{{{rm{n}}}}}}}^{{{{{{rm{cur}}}}}}})), we used the proportion of ({{{{{{rm{A}}}}}}}_{{{{{{rm{n}}}}}}}^{{{{{{rm{pot}}}}}}}) marked as wilderness and non-productive/ snow areas as described below. The area of each of the 45 LU-types within each grid cell (({{{{{{rm{A}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}})) was extracted from respective land cover and LU maps applying the approach outlined in ref. ⁵⁰ – with 2010 as year of reference wherever possible (Supplementary Data 2).

Following ref. ⁵⁰, builtup land, total cropland (including annual and permanent crops/ plantations), permanent pastures (areas used as pastures for more than five years) and rangeland (available in the two sub-categories natural and converted) extents were taken from the LU database HYDE 3.2.1⁴⁹ which was adapted to include rural infrastructure areas by assigning 5% of each grid cell’s cropland area to builtup land. We then split the total cropland cover into areas used for 41 different annual and permanent crops by integrating data from the Spatial Production Allocation Model (SPAM) for 2010^51,52 and adjusting them to cropland extent in the data from ref. ⁵⁰. To comply with the IUCN habitats classification scheme⁴², some of these crops were grouped into the plantation category (permanent crops), while the remainder was grouped into the cropland category (annual crops; see Supplementary Data 2 for details).

Wilderness areas were derived from the combination of human footprint data, i.e., a spatially explicit inventory of human artefact density available for 1993 and 2009^53,54 and intact forest landscape data for 2000 and 2013⁵⁵. Core wilderness areas without human use were defined as having a value of zero human footprint and, in forests, being part of an intact forest landscape⁵⁵. Within forests, the additional category of peripheral wilderness was introduced for areas where either only zero human footprint is recorded, or only an intact forest landscape exists.

The area remaining in each grid cell after allocating the above land cover types represents area covered by used forests and other land with mixed land uses⁵⁶. Hence, in addition to the approach in ref. ⁵⁰, forests were split into deciduous and coniferous forests based on the description of the ESA CCI land cover categories⁵⁷. This distinction was necessary for the differentiated allocation of wood harvest (see below). A further refinement was applied by identifying plantation forests, defined as areas in non-forest biomes converted to forests for forestry and areas in forest biomes converted to non-native forest types⁵⁸, which were linked to the IUCN habitat class plantations (Supplementary Data 2).

As in ref. ⁵⁰, the remaining area not allocated to any of the land cover or LU types above is denoted as “other land, maybe grazed”⁵⁶. These lands, typically treeless or bearing scattered tress, were allocated to converted grasslands on areas that potentially carry forests or to natural grassland on areas where the potential vegetation would not consist of forests²⁵.

To arrive at the six broad LU-type aggregates compatible with the IUCN and Birdlife habitat affiliation schemes^42,46 and PREDICTS categories²¹ (needed for quantifying LU-intensity effects, see below in section “Affinity parameter” for details), we rearranged and aggregated the described LU layers as needed (see Supplementary Data 2 for an overview). (a) Builtup remained as described above. (b) Cropland was defined as annual crops, covering the respective 29 SPAM categories plus fodder. (c) Pastures were defined as areas where pristine ecosystems were converted to grasslands and includes permanent pastures and converted rangelands from HYDE 3.2.1⁴⁹, plus those parts of “other land maybe grazed” located in forest²⁵. (d) Grazing land was defined as natural or near-natural areas where grazing occurs and includes natural rangelands from HYDE 3.2.1⁴⁹, plus 50% of each grid cell’s open forest area and 25% of each grid cell’s peripheral wilderness area, the latter two assumed to be only occasionally grazed and hence given low grazing intensity (see below), plus those parts of “other land maybe grazed” located in non-forest²⁵. (e) Forests were defined as forests where forestry occurs and includes 100% of each grid cell’s closed forest area, 50% of each grid cell’s open forest area, and 25% of each grid cell’s peripheral wilderness area, the latter two assumed to be only occasionally used for forestry and given low intensity (see below). (f) Plantations were defined as areas where pristine ecosystems were converted into plantation-like LU and include the 11 SPAM categories representing permanent / plantation crops, plus used forests identified by ref. ⁵⁸ as plantations (see above). As stated above, these aggregated broad LU-types were needed to align the different LU categorizations used in the different data sources with each other. The effects on biodiversity were then calculated on each of the 45 LU-types and afterwards aggregated to the six broad LU-types to give a better overview.

Continuous LU-intensity indices (({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}))

We constructed continuous LU-intensity indices ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) for each of the 45 LU-types based on gridded management descriptors¹⁵. For this purpose, we used two different sets of intensity indicators (called Set 1 and Set 2) to compare and combine their impact on predicted species loss. We used two indicator sets to account for the multidimensional nature of LU-intensity^9,12 and to include a wide range of available data products. For an overview of which data products and assumptions went into the individual sets, please refer to Supplementary Data 2.

Set 1

Set 1 is taken from the human appropriation of net primary production (HANPP) framework, a socioecological indicator basically describing the LU mediated extraction of biotic resources in the context of global biogeochemical cycles²³. We used the ratio of HANPP_harv to NPP_pot as a systemic metric to assess LU-intensity^12,22, with HANPP_harv being harvested or extracted biomass and NPP_pot being NPP of potential natural vegetation, i.e. the vegetation existing under current climate conditions in the hypothetical absence of LU²³. The ratio HANPP_harv/NPP_pot relates harvest to the productivity potential of the land where the harvest takes place and is, thus, robust against geographic differences in natural productivity. As it is related to energy availability in ecosystem food chains, it may be linked to the species-energy relationship, the strongest correlate of spatial biodiversity patterns at larger scales⁵⁹.

For calculating NPP_pot, LPJ-GUESS⁶⁰ version 4.0.1 was used in its standard configuration but with nitrogen limitation disabled and forced by the CRU-NCEP climate data^61,62 aggregated from 6-hourly to monthly fields.

HANPP_harv of all LU-types except builtup was calculated based on the FAOSTAT database by principally accounting total biomass flows via conversion and expansion factors as outlined in ref. ⁶³. As a special case, HANPP_harv of built-up was assumed to be half of the actual NPP, which was defined as 1/3 of the potential vegetation in ref. ⁶⁴. This results in a constant intensity on built-up land of ~17% of NPP_pot.

HANPP_harv of permanent and non-permanent crops was spatially downscaled following 40 permanent and non-permanent crop-specific production patterns from the Spatially-Disaggregated Crop Production Statistics Database (SPAM⁵²), merging minor SPAM categories such as “robusta coffee” and “arabica coffee” to ensure consistency with FAOSTAT reporting. Additionally, we added the LU-type fodder, which was downscaled following NPP_pot patterns.

Harvest of natural and plantation forest is reported by FAOSTAT in the four categories industrial roundwood, wood fuel, and coniferous and deciduous. We allocated industrial roundwood harvest to closed forests, while we split wood fuel harvest in proportion to productivity between closed and open forests, independently for deciduous and coniferous forests, respectively. For Set 1, we assumed forestry harvest to follow the patterns of forest NPP_pot⁶⁵. These intensity definitions were used for both natural and plantation forest.

Reported harvest on grazing land and pastures was allocated following patterns of aboveground NPP accessible for grazing as reported in ref. ⁶³. Following the assumption that systems with low natural productivity allow for a lower maximum harvest than systems with high productivity, we assigned a maximum harvest intensity of 40% at a level of accessible NPP of 20 gC/m² and increased this linearly to a maximum grazing intensity of 80% at 250 gC/m². Such, harvest was concentrated on grazing land and pastures with high productivity. In cases where the calculated national grazing land and pasture harvest demand surpassed NPP availability on grassland, we used information on fertilization rates on grassland⁶⁶ to either adjust NPP or harvest data: NPP was boosted in countries where more than 5% of overall fertilizer consumption was applied to grasslands, while countries where no relevant fertilization of grasslands occurred, the reported harvest demand was reduced accordingly, assuming it will be met from other sources. This intensity definition was applied to both (natural) gazing land and (converted) pastures.

Set 2

For the LU-intensity indicator Set 2, we used published data from different sources. For cropland we used the input metric nitrogen application rates (in kg N/ha of cropland)^12,22, available for 17 major non-permanent crops^67,68. For crops from the SPAM categories (see above) not covered by these data, we used the within-grid-cell area-weighted average of other crops in the same cell. For areas designated as cropland in our data (see above) but not in the available N application data, we assumed national average values of the respective crop.

For pastures and grazing land, we used gridded livestock information⁶⁹. We used information on the typical weight per animal to calculate livestock units⁷⁰ and aggregated the data for all ruminant species (buffalo, horses, cattle, sheep, goats). This data on livestock numbers per grid cell was then divided by land area per grid cell to arrive at livestock densities, which were applied to the extent of grazing land and pastures. Please note that this dataset contains information on the number of livestock (per species group) per area in a grid cell and thereby differs from the grazing intensity metric applied in Set 1⁷¹, as grazing animals may be fed from other sources than grassland⁷².

For builtup, we aggregated a 1 km built-up area density map for 2014⁷³ to the target resolution of 5 arc min and used it as is as intensity indicator.

For natural and plantation forest, we used the same data as described above for Set 1, but we assumed forestry harvest to follow another pattern. We calculated the difference between potential and actual biomass stocks⁷⁴ and allocated forestry harvest within each country according to these patterns, i.e., the share of national forest harvest allocated to a forest cell corresponds to its share in the national difference between potential and actual biomass stocks.

Scaling of LU-intensity indices

For the purpose of applying linear functions on species richness loss caused by LU-intensity (see below, affinity parameter), we scaled each LU-intensity indicator to values between 0 and 1, with 0 being no intensity (hypothetical) and 1 being the intensity threshold above which an increase of intensity causes no further increase of species loss. This threshold is not necessarily the highest recorded value of an intensity indicator, as effects may be regionally variable. We therefore winsorized some LUI indicators to that intensity threshold before scaling them (dividing by this threshold). These thresholds were defined as follows.

In Set 1, maximum intensity was assumed to be reached at harvesting 100% of NPP_pot on cropland.

In forests (natural and plantation), maximum intensity was derived from ref. ⁷⁵, which limits sustainably harvestable aboveground biomass in forests to 30% of NPPpot. In concordance with the HANPP framework, we included the belowground biomass destroyed by forestry using biome-specific factors⁷⁶.

On grazing land and pastures, maximum intensity was defined as removal of all NPP accessible for grazing. This considers only the aboveground and non-woody parts of NPP_pot. The maximum removable aboveground share was estimated as 50% of NPP_pot, and the proportion of non-woody vegetation was estimated as 30% (in closed-canopy land cover types) or 100% (on open land cover types)⁷¹. HANPP_harv/NPP_pot was assumed to be at its maximum intensity level when the maximum level of grazing intensity, as described above, was reached. The resulting thresholds are in line with literature^77,78, and assume that maximum intensities will be reached faster in systems with low natural productivity.

In Set 2, for all crop types (permanent and non-permanent) except legumes, N application rates were capped at 150 kg N/ha, i.e., we assumed that 150 kg N/ha was the maximum LU-intensity on cropland, beyond which no further species richness loss occurs, i.e., after which an increase of N application rates causes no further increase in species loss based on ref. ⁷⁹. For legumes, under the assumption that they need less N fertilizer due to their N-fixing capabilities, we assumed the following cap values, based on information provided in ref. ⁸⁰: beans and lentils at 110 kg N/ha, chickpeas at 100 kg N/ha, soybean at 70 kg N/ha and cowpeas, pigeon peas and other pulses at kg N/ha 90.

For pastures and grazing land, maximum intensity was defined as the per biome 80^th percentile of livestock-density.

The intensity of builtup area was not winsorized.

Affinity parameter (({{{{{{rm{h}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}))

The affinity parameter ({{{{{{rm{h}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) can be regarded as a LU-intensity dependent weighting factor for the area of each of the 45 LU-types used here. For low affinity, i.e., a small fraction of native species is left due to LU, the area of this LU-type (({{{{{{rm{A}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) in formula 1) is down-weighted, resulting in higher species loss ({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{rm{loss}}}}}}}) (and vice versa). The affinity parameter consists of two terms, (a) ({{{{{{rm{r}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}), the fraction of species affiliated with a given LU-type, and (b) ({{{{{{rm{f}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}), the fraction of ({{{{{{rm{r}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) that remains when LU-intensity (({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}})) rises to a particular level.

The fraction of species affiliated with a certain LU-type under minimal ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) (({{{{{{rm{r}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}})) is based on the habitat affiliation information taken from the IUCN Red List API⁴⁵ and BirdLife data⁴⁶ cross-tabulated with our mapped LU-types (Supplementary Data 2). We calculated ({{{{{{rm{r}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) by dividing the number of species affiliated with a certain LU-type (({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}^{{{{{{{rm{pot}}}}}}; {{{{{rm{LU}}}}}}}})) by the number of native species expected in this cell under pristine ecosystem conditions (({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{rm{pot}}}}}}})) as

Please note that for the two unconverted broad LU-types grazing land and forests, respectively (see above), we assumed no land conversion prior to its use, leading to ({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}^{{{{{{{rm{pot}}}}}}; {{{{{rm{LU}}}}}}}}={{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{rm{pot}}}}}}}). We further assumed that the whole fraction of LU-type affiliated species ({{{{{{rm{r}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) are present in a given LU-type as long as ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) is minimal (here < 0.17, see below). If ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) is higher, ({{{{{{rm{r}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) will be further reduced according to functions which vary across LU-types (({{{{{{rm{f}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}})). These functions were derived from estimates of linear mixed effects models published in a recent analysis of fractional species richness declines in local assemblages. The authors used the PREDICTS database and distinguished four broad LU-types and three categorical intensity levels (Minimal, Light, Intense)^5,21. The four broad LU-types from this study were assigned to the six broad LU-types we used in our aggregation (see Supplementary Data 2). Note that while cSAR calculations were performed on LU maps that distinguished all 45 LU-types, the intensity-driven loss functions were identical to all types belonging to one of the four categories distinguished in PREDICTS, e.g., for all 30 annual crop types assigned to the cropland category or for all 11 permanent crop types plus timber plantations assigned to the plantation category (see Supplementary Data 2 for details).

Here, we assume that species loss associated with a particular LU-type under Minimal intensity is equal to ({{{{{{rm{r}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}), and that ({{{{{{rm{f}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) for the categories Light and Intense is the proportional SR decline relative to Minimal,

with s_n,b the fraction of species loss caused by the respective intensification step (from Minimal to Light or Intense, respectively), ({{{{{{rm{Est}}}}}}}_{{{{{{rm{int}}}}}}}) the published species loss estimates for each of the three categorical intensity levels and ({{{{{{rm{Est}}}}}}}_{{{{{{rm{Minimal}}}}}}}) the respective published estimate for Minimal (Supplementary Fig. 5). To create a continuous relationship between ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) and ({{{{{{rm{f}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}), we assigned these three categorical intensity levels to fixed continuous ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) values (with ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) scaled between 0 and 1) and linearly interpolated between them. The definition of the three intensity levels in the literature does not allow an exact assignment along the continuous ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) gradients and we took the most parsimonious assumption that Minimal is between the 1^st and 33^rd (i.e., at 0.17), Light between the 34^th and 66^th (i.e., at 0.5), and Intense between the 67^th and 99^th (i.e., at 0.83) percentile of the scaled intensity gradient. We then placed these s_n,b values in the respective coordinate system (see Fig. ED5) and linearly interpolated between the values to yield an intercept and slope for each interpolation step. This procedure allowed us to calculate as

For the extrapolation below ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) = 0.17 (i.e., below Minimal), we decided to keep the relative species richness constant because the highest possible ({{{{{{rm{f}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) must be 1, in accordance with the IUCN habitat affiliation scheme. In other words, we assumed that all species adapted to a certain LU-type according to the IUCN habitat affiliations can thrive in the respective ecosystems as long as the LU-intensity is minimal. For relative intensity levels >0.83, we argue that extrapolation outside the measured intensity range is uncertain, and that an increase in LUI above 0.83 (i.e., Intense) might not necessarily result in even stronger effects on SR. See Supplementary Fig. 5, which illustrates the results of these considerations and shows the continuous effect of ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) on SR used in this study.

The affinity parameter ({{{{{{rm{h}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) was then calculated as follows and inserted into formula 1 (cSAR model).

Species loss caused by LU-intensity (({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{rm{loss}}}}}}{{{{{rm{int}}}}}}}))

In order to calculate the relative impact of LU-intensity on species richness, we re-ran the model with ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) = 0 in all grid cells and LU types, thereby effectively setting ({{{{{{rm{f}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) = 1 and ({{{{{{rm{h}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}={{{{{{rm{r}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}). The results of this model can be considered as delivering the land conversion effect without any possible enhancement by intensification. In addition, we designed a hypothetical, back-of-the-envelope intensification scenario where ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) = 1 in all grid cells and LU-types.

The contribution of intensity to the species richness loss was then calculated as

With ({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{{rm{loss}}}}}}; {{{{{rm{conv}}}}}}}}) being the results of the ({{{{{{rm{LUI}}}}}}}_{{{{{{rm{n}}}}}},{{{{{rm{b}}}}}}}) = 0 model and ({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{rm{loss}}}}}}}) from Eq. (1).

Native area-of-habitat loss of individual species

The cSAR model calculates by how many species the native species pool is reduced in response to LU in each 5 arcmin grid cell. However, it does not identify the individual species lost. To estimate each species’ native AOH loss, we randomly drew the predicted number of species lost from the native species pool of each cell.

First, we rounded the number of species lost as calculated by the cSAR model to the next integer for losses from both conversion (({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{{rm{loss}}}}}}; {{{{{rm{conv}}}}}}}})) and intensification (here taken as ({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{rm{loss}}}}}}}-{{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{{rm{loss}}}}}}; {{{{{rm{conv}}}}}}}}), see section above: “Species loss caused by LU-intensity”). To avoid rounding all values below 0.5 to 0, and, hence, to underestimate low levels of species loss, particularly in species-poor regions, we used a two-step rounding routine. First, prior to actual rounding, we randomly decided whether a number is rounded to the next higher or lower integer, with the likelihood of either decision depending on the decimal number’s (positive or negative) distance to 0.5 (i.e., the decimal number gave the likelihood of rounding up). Second, we took the species list used to generate ({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{rm{pot}}}}}}}) (see above under potential species richness) and modified it to either contain only species affiliated or unaffiliated with each LU-type, yielding two species lists for each grid cell and LU-type, respectively. The list of species affiliated with a particular LU-type was then used to select species predicted to get lost due to intensification, while the list of species not affiliated with it was used to select species lost due to conversion.

From each grid cell, we then randomly drew as many species from these lists as determined by the rounding routine above, considering each LU-type and whether the number of lost species was caused by intensification or conversion. However, in each cell, each species could only be drawn once, independently of whether it was affiliated with several LU-types. As a consequence, the order in which LU-types are considered when drawing species is relevant for the outcome of the calculation. For instance, species simultaneously unaffiliated with cropland and affiliated with natural forest may never be drawn in response to intensification of natural forest if losses due to conversion into cropland are always handled first. Therefore, we randomly iterated the sequence by which LU-types were considered, i.e., the order of LU-types, in the random draw routine in each of 100 repeated runs.

We repeated the random-draws 100 times to yield a representative sample and processed the resulting 100 lists of species-per-cell losses in the following way. For each of the 100 runs, we summed the areas of all cells each species was drawn from, i.e., predicted to be lost, across all LU-types and within individual LU-types, yielding 100 area sums per species (one per run). From these 100 areas, we calculated the mean and the 0.025^th and 0.975^th quantiles as 95% confidence intervals (CIs). The means and CIs were then divided by the species’ global AOH (sums of cell areas in native range), thereby yielding the proportional global AOH loss attributable to current LU in general, and to different LU-types or land conversion vs. LU-intensity in particular.

Description/ presentation of results

All cSAR model calculations were based on global land use maps that distinguish 45 LU-types as described above. For the sake of simplicity, we present results aggregated to the six broad LU-types cropland, pastures, natural grazing land, built-up, plantations, and forests (natural/ near-natural forest under forestry; see Supplementary Data 2). All calculated SR decreases are expressed in percentage losses relative to ({{{{{{rm{S}}}}}}}_{{{{{{rm{g}}}}}},{{{{{rm{n}}}}}}}^{{{{{{rm{pot}}}}}}}).

Summary statistics mentioned in the text and Supplementary Data 1, 5 and 6 were calculated as follows. Global, biome-wide and nation-wide average species losses due to conversion, LU-intensity or both were calculated as cell-area weighted means across all cells with native terrestrial vertebrate species either excluding or including wilderness areas (which, for this purpose, are defined as cells where the sum of all LU area equals 0). The percentual land area exceeding a certain threshold of calculated SR decline were calculated by dividing the area sum of all cells exceeding that threshold by the area sum of all cells with native species excluding wilderness.

Differences among average AOH losses (across all taxonomic groups) mentioned in the text and Supplementary Data 3 were modelled using generalized linear models assuming a binomial distribution (proportional AOH loss between 0 and 1), each species’ mean AOH loss (mean of 100 random draw runs) as response, and either (a) IUCN categories, (b) land use types, or (c) taxonomic group as predictor variables. Differences between predictor variable levels were then alculated by multiple comparisons via p-values adjusted with the Tukey method. A p-value of < 0.05 was taken as statistically significantly different. We used the function glm from base R⁸¹.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com