Introduction

Landscape ecology envisions the landscape as the outcome of the complex relationship between humans and nature (Opdam et al. 2018). Land system science defines land systems, similarly, as the terrestrial component of the Earth system encompassing all processes and activities related to the human use of land, including socioeconomic, technological and organizational investments and arrangements (Verburg et al. 2013a). Land system change is seen as both a cause and consequence of socio-ecological processes, rather than as a sole outcome, or emergent property, of human–environment interactions (Verburg et al. 2015). Land system science and landscape ecology share a lot of concepts of interest: the acknowledgement of human activities as a main actor of land use and landscape change (Lambin et al. 2001; Turner II et al. 2007; Roy Chowdhury and Turner 2019); the importance of addressing multiple spatial and temporal scales (Veldkamp and Lambin 2001; Dearing et al. 2010); the attention for human benefits of nature through the concepts of ecosystem services and Nature’s contributions to People (Wu 2013) and the attention for sustainability solutions (Nielsen et al. 2019). At the same time, landscape ecologists have been putting a strong emphasis on characterizing landscape pattern across spatial and temporal scales (O’Neill et al. 1988; Li and Wu 2004; Wu 2004), and understanding how landscape structure and composition emerge from socio-ecological processes and impact the functioning and performance of ecosystems (Nagendra et al. 2004). Land system science has, in contrast, often focused on the characterization and explanation of land use change at the level of individual pixels (observed from remote sensing), plots or at the level of individual or collective actors (Overmars and Verburg 2006; Manson 2007), with less specific attention for the landscape patterns and structures. Despite common analytical methods, these different foci can lead to different representations of the landscape and its dynamics, having repercussions for interpretation and assessment of drivers and impacts of land use and landscape changes.

While both disciplines acknowledge the role of human use of land, both often represent land use and landscapes in empirical studies mostly by land cover. This is a direct result of the relative ease of observing land cover from remote sensing while land use intensity or land use management is more difficult to characterize and data are often absent (Kuemmerle et al. 2013; Erb et al. 2017). This is especially the case for larger scale assessments (Verburg et al. 2013b). Rather than representing the landscape and its composition, the dominant land cover is used to represent the land system and landscape. Some approaches use a representation that describes the different land cover fractions, but ignoring their spatial configuration. It is well-known that many of the Earth’s landscape are mosaics of different land covers, encompassing novel ecosystems (Perring and Ellis 2013). Although coarse-scale land cover representations are always constrained in thematic resolution, ignoring these important mosaic landscapes is an important omission.

More recently land system representations have been proposed to better capture human use of the land, or some aspects of landscape configuration. For example, Ornetsmüller et al. (2018) distinguished permanent agricultural use from shifting cultivation in a land systems classification for Laos while Debonne et al. (2018) represented smallholders vs large scale land acquisitions explicitly given their different drivers and social and environmental impacts. On a global scale, Ellis and Ramankutty (2008) and Václavík et al (2013) used a wide range of environmental and demographic conditions to identify archetypical combinations of socio-ecological conditions defining land use, while van Asselen and Verburg (2012) used information on land cover, input intensity and livestock to subdivide global land cover into land systems. These land system representations have been created for various purposes and, while making the best use of available data, improved over land cover classifications by including specific aspects relevant to landscape ecology and land system science. However, in spite of this progress, most global assessments for biodiversity such as for IPBES and IPCC are still based on relatively rough and classical land cover representations.

In particular, most large-scale biodiversity modelling (especially those using species distribution modelling, SDM) have so far either focused on the impacts of climate change (Thuiller et al. 2019) or on climate change with some broad land cover classes (Thuiller et al. 2014). This has been done in SDM by simple filtering species distributions by land cover classes (Maiorano et al. 2013; Powers and Jetz 2019) or land use types (Newbold 2018). Another approach in biodiversity assessments is to assign expert based values to land use, based on the observed or expected effect it has on biodiversity, or based on the level of naturalness (Pouzols et al. 2014; Di Minin et al. 2016). Such approaches not only require a well-established relationship between different species and the land use impact, but also neglect some landscape characters and land management that play important roles in biodiversity. While some large-scale biodiversity assessments account for nitrogen application and habitat fragmentation as additional stressors (Schipper et al. 2020), these are independently addressed from land cover information without using an integrated landscape approach.

There were several reasons to ignore land use or land cover in these models that relate species occurrence (or abundance) to environmental variables. The first one was that since land cover and land use are mainly driven by climate, they do not bring much additional information, and rather create some multi-collinearity issues when combined with climatic variables in SDMs (Thuiller et al. 2004). The second reason is that past land cover or land use maps were generally of poor thematic resolution (e.g. forest, urban, cropland, grassland and other) which prevent to estimate the tight association that could arise between a given species and for instance extensive croplands. Still, some recent exercises have shown the interest of using detailed land use classification to predict species extinction at global scale (Powers and Jetz 2019). IPBES (Diaz et al. 2019) also identified land use change as the main driver of biodiversity loss and several papers have argued for inclusion of land use change as part of biodiversity assessments (Titeux et al. 2016; Randin et al. 2020).

A way to possibly advance large-scale biodiversity assessments is to provide a land system classifications map that includes more ecologically relevant variables describing landscape and land management conditions. In this paper we aim to present a continental-scale land system classification for Europe making best use of available data that moves beyond the many assessments of biodiversity that are using rough land cover classes and relationships between land cover and habitat for species (Newbold 2018; Powers and Jetz 2019). We focus our analysis on including aspects of landscape structure and land use intensity. There is sufficient evidence that land use intensity (Beckmann et al. 2019) and landscape structure (Walz and Syrbe 2013; Boesing et al. 2017) are strong determinants of the land use impact on biodiversity. Recent studies showed the impact of forest management types on species richness (Chaudhary et al. 2016) and the differences of accounting for forest management in addition to forest cover change on global biodiversity (Schulze et al. 2020). Other studies indicated the role of agricultural management (Beckmann et al. 2019) and livestock grazing (Zhang et al. 2017; Li et al. 2019), as well as the importance of landscape mosaic for conservation (Harvey et al. 2008) and ecosystem service provisioning (Dainese et al. 2019). While land use studies often characterize urban area by the amount of build-up land (Seto et al. 2012; van Vliet et al. 2019) the structure of urbanization has a strong impact on biodiversity (Malkinson et al. 2018). In peri-urban areas the build-up land cover is often not vast, but the wide-spread disturbance through human activities and infrastructure causes much larger declines in biodiversity than the overall build-up area would suggest (Buczkowski and Richmond 2012; Concepción et al. 2016).

The land system typology presented in this paper aims at incorporating important elements of the landscape and land system relevant to biodiversity. To demonstrate the added value of such representation compared to the often-used land cover classifications, we used bioclimatic variables and the land system map to model the species distribution on nine selected bird species that have wide habitat range across Europe with different vulnerability to land use and landscapes. We hypothesized that the predictive ability of models and variable contribution to the spatial distribution is improved by the new land system typology.

Materials and methods

Classification overview

We used an expert-based land system classification approach based on variables for which there is evidence of the effect on biodiversity. We operationalized land systems as areas that host one or more land use activities at the landscape level having an impact on species occurrence and biodiversity. We focused, consistent with traditional classifications, on the following major land systems groups (Table 1): water and wetland systems, human settlement systems, forest systems, grassland systems, cropland systems, and mosaic systems. These land systems are selected based on common land uses in Europe, are represented by major groups in the continent’s land use and land cover datasets (ESA and UCLouvain 2010; European Environment Agency 2018) and by expert opinions towards their importance to study the impacts of land use on biodiversity. Particularly, the inclusion of mosaic system and the separation of forest mosaics and agricultural mosaics is aimed at better capturing the fragmented forest habitat and agricultural composition that certain forest and agricultural species are sensitive to (García-Navas and Thuiller 2020). We also have permanent crops as a sub-system in the cropland category because they can provide vital habitat for certain species and the management styles for permanent crops can generate distinctive biodiversity consequences (Bruggisser et al. 2010; Winter et al. 2018). The threshold of land cover extent for each land system was determined based on two criteria: it indicates the composition of land cover of the system (e.g., covers the majority of the pixel, or the mosaic combinations are meaningful for certain species), and it captures small habitat area requirements. Residential systems also include those with relatively low fractions of build-up land as the impact of associated activities and infrastructure on species is often proportionally large (Buczkowski and Richmond 2012; Concepción et al. 2016).These land systems were further classified based on different land use intensity indicators. We also aggregated classes with small areas (the full expert-based hierarchical classification procedure is included in Fig. S3). Our expert-based procedure, although tailored to the European conditions and for use in biodiversity assessments, follows common patterns with other global and regional land system classifications (e.g. Ellis and Ramankutty 2008; van Asselen and Verburg 2012; Herrero et al. 2014; Malek and Verburg 2017; Kikas et al. 2018).

Table 1 Overview of land systems
Full size table

We chose the spatial extent of the European Union (EU) with the United Kingdom (EU28 +), Norway, Switzerland, and the Western Balkans (Serbia, Kosovo, North Macedonia, Montenegro, Albania, and Bosnia and Herzegovina). However, we excluded Iceland, Turkey and Europe’s Outermost regions and Overseas Countries (mostly islands not present on the European continent) due to data issues, the analytical interests, and their limited importance for assessing continental species and biodiversity. To the best of our knowledge, this is the most complete coverage of land cover and land use analysis at the European scale while most studies do not cover non-EU countries such as Switzerland, Norway, and the Balkans. Such continuity of coverage is important for modeling species migration and projected changes due to climate change. We refer to this region as the non-EU region and the other parts as EU region.

We operated on a 1-km2 spatial resolution. Although some land cover data may have a more detailed resolution (e.g., Corine at 100 m or GlobCover at 300 m), the pursuit for a finer resolution than 1-km2 is hampered by a lack of the high-resolution data on land use intensity and management. Meanwhile, bio-climatic variables used for the large-scale biodiversity assessment and species distribution modeling are derived from global interpolations across the multitude of existing climatic stations. These interpolations have been made at global scale at 1-km2 and are provided from web portals like WorldClim or CHELSA. Some species may be able to thrive in habitat patches smaller than 1-km2, however this character, at least in some cases, can be captured by the mosaic land systems that we represent in our typology. Therefore, a spatial resolution of 1-km2 is an optimal choice for this study, considering data availability and usage for biodiversity assessment.

We complied and harmonized data of land cover composition and intensity indicators using a Geographic Information System based approach. Data were chosen based on the following principles: (1) they are publicly available, (2) we prioritized the highest spatial resolution and most recent date, (3) data produced at the European scale have priority over global datasets, (4) if the data were a result of downscaling, we ensure that no co-variates are used in the downscaling that would induce multi-collinearity with other variables often used in biodiversity modeling. (5) Global data were only used if previous criteria cannot be achieved and data gaps are apparent in particular regions. We mostly relied on global data for the non-EU region. All maps were resampled to a 1-km2 spatial resolution using the Lambert azimuthal equal-area projection that is frequently used for European maps. The details of data and methods used for each system are discussed below (Overview see Tables 1, S3 and S4, Fig. S3).

Extent of land systems

For land cover extent, we used Copernicus products derived from remote sensing data, such as CORINE land cover 2018 (CLC) (European Environment Agency 2018) and Pan-European High Resolution Layer (HRL) thematic maps (i.e., imperviousness, forests, grassland, water & wetness) (European Environment Agency 2015a, b, c, d). A global glacier database was used for the extent of glaciers (Global Land Ice Measurements from Space (GLIMS) and National Snow and Ice Data Center 2012) (Table S3). We aggregated the HRL thematic maps from 20 m to 1 km to derive the extent of (1) human settlement, (2) forests, (3) water and wetlands, and (4) grasslands. For the remaining land covers, including croplands (annual and permanent) and other land cover classes (e.g., shrubs, bareland and rocks) we used CLC maps to resample.

Classification thresholds based on the share of particular land cover within the 1-km2 pixel were chosen based on expert opinions for each system (Table 1). First, this enabled us to study land cover types that may affect species largely but appear on a smaller scale, such as settlement systems. In addition, we used indicators other than land cover composition to characterize these classes, among which one particular example is the peri-urban and villages that may have less built-up land than cropland and forest cover (van Vliet et al. 2019). Despite the small percentage, these classes may cause profound effects on species activities due to the continuity and irreversibility of human disturbances. Second, compared to most land cover classification that determines the class by the largest land cover, we used a threshold of 70% to classify areas dominated by forest, cropland, and grassland. This is to better describe these classes with homogeneity while defining landscapes with several different land cover classes as mosaics. A pixel that contains multiple land cover components smaller than 70% of the cell is classified as a mosaic system, which is further subdivided by its main land cover type. As we acknowledge that all these thresholds are arbitrary and set to balance the importance of the land system components for biodiversity and the number of complex compound classes, we conducted a sensitivity analysis of the thresholds on the extent of land systems (Supplementary Material). We found only small changes in the overall patterns of land systems upon changing the thresholds within reasonable ranges.

Intensity of land systems

We used a range of metrics to measure land use intensity and landscape characteristics (detailed intensity classification of each system can be found in the supplementary material).

Human settlement systems

Urbanization in Europe undergoes multiple dimensional changes in land cover, land use, and socio-ecological activities (Shaw et al. 2020), and most recent changes are small incremental increases (van Vliet et al. 2019). In comparison to current practices that use “built-up” land or a binary classification of urban and rural area (Seto et al. 2012), we classified three intensity levels for settlement systems: urban core (high intensity), peri-urban (medium intensity), and villages (low intensity) based on the degree of imperviousness and spatial connectivity to urban core. We used the intensity of imperviousness (i.e., each cell has a value from 0 to 100%), namely artificial impermeable cover of soil, as a proxy of urbanization degree, following similar mapping exercises (Linard et al. 2012; van Asselen and Verburg 2012; Demuzere et al. 2019). In addition, we considered the adjacency to urban core areas when distinguishing peri-urban and rural settlements (Linard et al. 2012; Stürck et al. 2018) (Fig. S4).

Forest systems

Forest management and forest types are important to accurately assess biodiversity (Chaudhary et al. 2016; Schulze et al. 2020). We used wood production of Europe (Verkerk et al. 2015) and the probability of finding primary forest (Sabatini et al. 2018) to characterize the intensity of forest use in the EU region (Fig. S5). For the non-EU region, we extracted the forest classes from a global dataset (Schulze et al. 2019). The three classes (i.e., primary, naturally regrown, and planted) were reclassified as low-intensity, medium-intensity, and high-intensity to fit our forest classification (while acknowledging these classes do not fully correspond).

Grassland systems

Three indicators were used to map the intensity of grazing grassland: livestock unit density, nitrogen input, and mowing frequency. We selected three ruminant livestock groups (i.e., cattle, goat, and sheep) from the global dataset on livestock density (Gilbert et al. 2018) and converted them to livestock unit (LSU) density (Eurostat 2013), which is the number of animals equivalent to one adult dairy cow producing 3000 kg of milk annually (LSU/km2). The nitrogen application rate was calculated using a novel downscaling of statistical data by the Joint Research Center, which included mineral fertilizer input, the manure application and deposition by grazing animals. The indicator of mowing frequency represents the average number of annual mowing events in the period of 2000–2012 (Estel et al. 2018). The LSU density and nitrogen application rates were reclassified into three levels as low, medium, and high. This is based on the biodiversity dynamics relevant to nitrogen application rate (Kleijn et al. 2009) and supported by other studies (van der Zanden et al. 2016; Estel et al. 2018). When combining the three indicators, LSU was used as the primary indicator and then complemented by nitrogen and finally the mowing frequency (Fig. S6). For the non-EU region, we could only use LSU as the indicator due to data unavailability.

Cropland systems

Erb et al. (2013) and Kuemmerle et al. (2013) inventoried the conceptual basis of cropland intensity and provided a range of potential input and output indicators. We focused on agricultural input alone, and chose nitrogen input and field size as indicators with high relevance for many species, to quantify the agricultural intensity.

The nitrogen application rate was retrieved from the same source and reclassified into the same levels as used in grassland system. For the field size, we used the field size map created by Tieskens et al. (2017) and reclassified into small, medium, and large fields. We used the nitrogen as primary intensity factor and the field size as secondary indicator (Fig. S7). If a cropland pixel had high nitrogen input it was classified as high-intensity cropland unless the field size was in the small class, and vice versa.

For the non-EU region, we used two global datasets. One is the global nitrogen application rate dataset from EarthStat (Mueller et al. 2012; West et al. 2014). The other one is the global field size map (Lesiv et al. 2019). For Switzerland, we assumed cropland above 800 m as medium-intensity and below as high-intensity.

In addition to the intensity of annual crops, we singled out permanent crops because they have different management styles and effects on biodiversity. The understory between stripes of grapes and olive trees are an important indicator both for habitats and management styles (Bruggisser et al. 2010; Winter et al. 2018). Therefore, the second land cover type identified in LUCAS (Eurostat 2018) was selected, representing the understory vegetation, as a biodiversity relevant indicator of management intensity (Fig. S8).

Mosaic systems

Pixels without dominant land cover types were classified as mosaic systems. Two sub-mosaic systems were classified: agricultural mosaics and forest/shrubs mosaics. Croplands and grasslands without tree cover were aggregated to agricultural mosaic systems. This sub-system was further classified into three agricultural intensities: low, medium, and high. Two indicators, LSU and nitrogen application rate from the same source as above, were used. LSU was used as the primary indicator and nitrogen was the secondary indicator (Fig. S9).

The forest/shrub mosaics include mixed land cover types with some portion of forest/shrubs. This is because forest fragments, remnant trees, and small parts of shrubs and bare can serve as valuable habitats for many species. These mixed compositions increase landscape connectivity and preserve potential for restoration (Harvey et al. 2008; Horák et al. 2019). Forest and shrubs were combined as one in the mosaic systems, for their similarity in serving as habitat in a heterogeneous land system. Classes with small areas were assembled into the next class that are most resembling their land cover composition. In total, four forest/shrubs sub-mosaic systems were defined: forest/shrub and cropland, forest/shrub and grassland, forest/shrub and bare land, and forest/shrub and mixed agricultural land.

Using land system map for species distribution modeling

To demonstrate the effect of adding land system and landscape characters for biodiversity assessment, we calibrated a set of species distribution models using an ensemble approach (Guisan et al. 2017). We considered nine bird species with different sensitivity to land use management: three species sensitive to forest management (Columba oenas, Dendrocopos major, Lophophanes cristatus), three species sensitive to agricultural land-use intensity (Alauda arvensis, Motacilla flava, Tyto alba), and three generalist species (Corvus corax, Hirundo rustica, and Phoenicurus phoenicurus) (see S.I. for more details). For each species we retrieved all available occurrence points from the iNaturalist repository (iNaturalist.org 2020) and kept only those of the species that were validated by experts. For species resident in the study area we retained all occurrences collected during the breeding period; for migratory species we retained only locations collected during the breeding period and falling within the breeding extent of occurrence (data on breeding range from Maiorano et al. 2013). Information on the ecology and distribution of the species was obtained from the Birds of the World (Cornell Laboratory of Ornithology 2020).

We compared the ensemble models calibrated with the land system map to models with land cover maps. We hypothesized that the land system map would improve the model’s predictive ability and indicate higher variable contribution (i.e., the relative contribution of the environmental variables to the predicting model) as compared to the use of solely land cover information, especially with species that have specific habitat requirements. We used two land cover maps as comparison: one is an integrated Copernicus Global Land Cover layers (Global Land Cover) (Buchhorn et al. 2019) that contains 10 land cover types (i.e., closed forest, open forest, shrubland, herbaceous vegetation, herbaceous wetland, moss & lichen, bare/sparse vegetation, cropland, built-up, snow & ice, permanent water bodies), and the other one is a map created through reclassifying our land system map to the seven dominant land cover classes (Land Cover EU) (i.e., water/ice, settlement, forest, cropland, grassland, shrubs, bare).

To calibrate all models, we also considered a set of bioclimatic variables. Among the 19 bioclimatic variables available in the Chelsa climate repository (Karger et al. 2017b, a), we selected six variables (i.e., bio2, 4, 8, 10, 15, 19, Table S6) that have low collinearity (Variance Inflation Factor test < 4). We ran the ensemble for every species using three sets of environmental variables respectively: (1) climate variables only, (2) land cover/system maps only (i.e., Global Land Cover, Land Cover EU, and Land System EU), and (3) combined climate with each of the land cover/system maps. The three land cover/system maps were used as categorical data.

The ensemble modelling procedure was developed as follows. We combined five state-of-the-art statistical models: generalized linear models (GLM), generalized additive models (GAM), boosted regression trees (BRT), Random Forest (RF), and maximum entropy modeling (Maxent). Since these models required pseudo-absence (PA) data, we randomly drew 10,000 PA across the whole Europe. This procedure was repeated three times to account for the stochasticity of the PA generation. For each of these three datasets (presence-PA), we then generated three calibration-evaluation sub-datasets for cross-validation. The calibration data (a 70% random part) were used to calibrate the models, and the remaining 30% to evaluate them. Again, this procedure was repeated three times to account for the stochastic procedure. To summarize, for each of the five statistical models, nine models were run and evaluated (45 models in total for each single species).

Models’ predictive ability were evaluated using the true skill statistic TSS, which measures model performance on presence-absence data and, unlike Kappa, is independent of prevalence (Allouche et al. 2006). The ensemble model was then made of all calibrated models with a TSS > 0.3. We further measured the TSS of the ensemble model on the evaluation data, and then used the ensemble to predict species’ probability of occurrence across Europe (Marmion et al. 2009).

We measured the importance of the variables using a standard permutation procedure where each variable is randomly permuted before predictions (while the others are retained as they are). The difference between the original prediction and the one with the variable permuted gives a measure as the importance of the variable (i.e., Pearson correlation). The more different are the predictions, the more important is the variable. This was done for each species, for each variable and repeated three times.

Models and ensemble procedure were performed with the BIOMOD2 package in R (Thuiller 2003; Thuiller et al. 2009). See R scripts in Supplementary Materials.

Results

Current land systems

Figure 1 presents the spatial distribution of land systems across Europe, while Fig. 2 presents the proportion of the different land systems across the four main geographic regions in Europe. Additionally, we show particular land systems and regions in more detail (Fig. 3).

Fig. 1
figure 1

Distribution of land systems in Europe

Full size image
Fig. 2
figure 2

Area division of land systems in different regions (north, west, east, and south) of Europe. The center of each circles indicates the region where the land systems are summarized. The first circle (inner) indicates the area share of the main class, and the second (outer) circle indicates the area share of sub-system classes

Full size image
Fig. 3
figure 3

Mosaic and non-mosaic land systems in Europe and sample regions. Left map shows mosaic land systems resampled to a 2-km resolution for visualization. Right panels: A and B show major land systems without mosaics, while C and D present mosaic systems without major land systems. A Linear clusters of peri-urban and urban core classes (medium and intensity settlement systems) in Po Valley, surrounded by villages, high-intensity cropland and different intensity forest on the mountain; B Clusters of medium (peri-urban) and high (urban core) intensity settlement systems, surrounded by high-intensity grassland; (C) Forest/shrubs with mixed agriculture mosaics in Portugal and Spain. (D) Low-intensity agricultural mosaics in Romania with forest/shrub and cropland and grassland mosaics

Full size image

High-intensity settlements, conceptualized as urban centers, have an average imperviousness degree around 42%. The medium and low intensity settlement classes, which can be understood as peri-urban and villages, have an average imperviousness degree of 12% and 10% respectively, with the rest covered by cropland (37 and 36% respectively) and grassland (23 and 24% respectively). All three settlement classes cover 7.1% of the continent, which is higher than the estimation in other studies (e.g., 1.8% urban built-up by Levers et al. 2018). The peri-urban class is the dominant among three settlement classes in the Western and Southern Europe (Fig. 3A), while village landscapes cover the largest area among the three residential land system classes in Eastern Europe.

Forest systems are characterized by the high tree cover of over 70% often mixed with small areas of grassland and cropland. The high-intensity forest has an average annual wood production of 36.5 m3/km2, and 15.2 m3/km2 and 4.6 m3/km2 for medium and low-intensity forest respectively. Forest systems are the largest land system in Europe with a total of 1.6 million km2 area covering 32.3% land surface. Nonetheless, forests of low-intensity level are the smallest within the three intensity levels, accounting for 30.3% of the total forest classes and 9.8% of the total European land surface. The low-intensity forests are mostly located in the mountain areas of Europe, including the Alps, Pyrenees, the Carpathian Mountains, and in the Scandinavian Mountains. However, large areas of forests in Western Europe are used with high-intensity, particularly in Germany, France, and southern Sweden.

Grassland systems cover about 6.2% of the total European land surface, most of which are in Western Europe. Most high-intensity grassland areas are concentrated in the west and featured by different intensity indicators. For instance, in Ireland high-intensity grassland is characterized by high LSU density and high mowing frequency, whereas the Netherlands has high LSU density and nitrogen application rates (Figs. 3B). In the large areas in central France this class is featured by different levels of nitrogen input but in general frequently mowed and with moderate and high density of livestock. Other grassland classes are scattered in Southern Europe and Eastern Europe, and little grassland in the Northern region.

The annual cropland classes are characterized by a high average of 84% cropland cover and some forest and grassland. The low and medium intensity levels of annual cropland classes have an average nitrogen input of 87 kg/ha and 157 kg/ha respectively. The average nitrogen input of the high-intensity class is as high as 338 kg/ha, almost four times the nitrogen value of low-intensity annual cropland. Cropland systems are the largest systems covering about a third of Europe, except the northern region. In Western Europe, only 1.1% area coverage is classified as low-intensity cropland, while medium and high-intensity cropland both cover more than 10%. Most croplands in Eastern Europe are low to medium-intensity featured by relatively lower nitrogen applications than the rest of Europe. In addition, the majority of extensive permanent crops are found in Southern Europe with an area of 75,968 km2 while the majority of intensive permanent crops are also found in southern Spain and along the coastlines of Italy.

Mosaic systems are diverse systems characterized by a mixed low to medium coverage of cropland (6%-50%), grassland (10%-48%), and forest (6%-44%), without substantial built-up areas. Mosaic systems cover 14.5% of Europe. We distinguish two subsystems: mosaics with forest/shrub (11.0%) and agriculture mosaics (3.4%). The agriculture mosaics are mostly clustered in Western Europe and predominantly as high-intensity agriculture mosaics in Normandy and Brittany in west France (Fig. 3—left panel). In the east of Europe there is also some agriculture mosaics but mostly as low-intensity (Fig. 3D). The forest/shrub and grassland mosaics show a clear pattern along the main mountain ranges, such as the Alps (gaps in Fig. 3A), the Massif Central in central France, and the Balkan Mountains. In Southern Europe, there is only 1.3% of low-intensity agricultural mosaics and the rest of the mosaic land systems are forest/shrub mosaics with cropland (6.6%), grassland (3.3%), mixed agricultural land (1.3%), and bare land (0.4%). This is because of the common agroforestry landscape in Spain and Portugal, where forestry, grazing animals, and crop cultivation are found simultaneously in the same landscape (Fig. 3C).

Land system map contribution to species distribution modeling

Performance of species distribution models with climate and land input

Overall, the models with both climate variables and land cover/system maps have the highest predictive ability (Fig. 4), except for two species (i.e., Lophophanes cristatus and Corvus corax) of which the highest TSS appear in models with only climate variables. Models with the second highest TSS are calibrated with climate variables only and no land cover/system maps. Models with the lowest TSS appear to be using land cover/system maps only and no climate variables. Interestingly, climate variables are visibly dominant in these models, while the inclusion of land cover/system maps with climate variables only marginally improves TSS. However, for models without climate variables, the TSS when using the land system classification remains the highest. For five species, the land system map can improve the model’s predictive ability up to a useful TSS (> 0.3) while models with land cover maps only generate poor TSS (< 0.3) hence not reported in the figure. For the remaining species, the gap in predictive ability between models calibrated with land system map and land cover maps can be as large as 50%, with TSS increasing from 0.34 to 0.53 (e.g., Columba oenas), and 38% from 0.30 to 0.42 (e.g., Tyto alba).

Fig. 4
figure 4

True Skill Statistic (TSS) of selected bird species using ensemble of species distribution models calibrated with different environmental variables. Models with TSS < 0.3 are not reported in the figure. Species are ordered as: sensitive to forest management (Columba oenas, Dendrocopos major, Lophophanes cristatus); sensitive to agricultural land-use intensity (Alauda arvensis, Motacilla flava, Tyto alba); and three generalist species (Corvus corax, Hirundo rustica, and Phoenicurus phoenicurus)

Full size image

Species’ sensitivity to climate and land maps

The variable contribution of different land cover/system maps (Fig. 5) shows a clear pattern. For almost all nine species, bio 4 (i.e., Temperature Seasonality) is the most important variable. Its importance, however, is reduced when land cover/system maps are used in the model. The reduction of its importance has an average value of 17.9%, ranging from 6.1% in Tyto alba to 25.5% in Motacilla flava depending on species.

Fig. 5
figure 5

Variable contribution to species distribution in models calibrated with climate variables and land cover/system maps. Models with only land cover/system maps are not reported in this figure

Full size image

The land system map always has the largest contribution to the spatial distribution of modeled species among the three land cover/system maps. The difference of variable importance between land system map and land cover maps range from 4.9% to 21.6%. Despite the dominant explanatory role of climate variables, land system map sometimes can account for more than 30% among all variables contributing to the spatial distribution, such as Dendrocopos major, Phoenicurus phoenicurus, and Lophophanes cristatus, and even 50% for Motacilla flava (Fig. 5). Particularly for species Dendrocopos major and Motacilla flava, land system map becomes the most important variables, even more so than any climate variables.

While for species sensitive to forest management (three species as group 1 in the first row in Fig. 5) and for generalist species (the last row in Fig. 5), the land system map always demonstrates noticeably higher contribution than the other two land cover maps. This is, however, less profound for the species that are sensitive to agricultural land-use intensity (indicated as the middle row in Fig. 5), particularly for Tyto alba that only shows 5% more importance with land system map comparing to the global land cover map.

Discussion

Including landscape and land management in land use maps

In this study, we have shown a pragmatic approach to synthesize available data into a land systems classification that captures a number of aspects central to landscape ecology that are poorly represented in traditional land cover classifications: the land management intensity for grasslands, arable lands and forests; the importance of land cover mosaics that can be regarded as special land systems in which the functioning is dominated by the mosaics rather than by the individual land covers; and, the role of human disturbance through distinguishing settlement systems in which built-up land cover is far from dominant but human disturbance is playing a key role. We have shown that in an ensemble SDM exercise the land system maps consistently have added value as compared to using land cover classifications. At the same time, data availability at the large scale, the urge to keep the classification simple and close to well-known classification systems, limits the depth to which landscape structure and functional characteristics are included in our land system classification.

Beyond the landscape characteristics that are included in our map, there are many other aspects such as spatial extent of patches, heterogeneity, and connectivity that might be important to ecological function, not captured in our land system classification. The fragmentation of forest patches is one example that could be important for certain species. In the supplementary material we reported a spatial-context analysis of including fragmentation of forests in the classification. The results show that, to a limited degree, the inclusion of such a feature can further improve the predicting and explaining power of species distribution modeling (Table S2, Figure S1 and S2). Many landscape structure and pattern characteristics are largely linked to specific ecological processes. Therefore, instead of including multiple characteristics in a generic land system classification and ending up with too many classes, land system classifications should be adapted to the specific purpose and include those landscape characteristics relevant to the ecological processes studied.

SDM results interpretation

The species distribution modeling results show a clear trend with the land system map consistently having added value as indicated by the TSS and useful information when compared to other land cover maps. This is particularly true when the model considers only land cover/use data and no climate information. The gain in model predictive ability by adding a land cover/system map to climate variables is limited, and land system map does not always contribute to the highest TSS. However, the land system map holds important sources of explanatory power compared to the other two land cover maps as indicated by the higher variable contribution to the models. The limited improvement of TSS by adding land information to climate variables may be explained by that the European-wide bird species distribution is largely determined by climate (Thuiller et al. 2004), and among which the most prominent one is the temperature seasonality (bio4).

Particularly problematic are the results for the species sensitive to agricultural intensities. Our results of both model performance and variable importance indicate the improvement is most limited in this group of species. This may be explained by a more nuanced sensitivity to land use intensity and management than is represented by our classes using common intensity level thresholds that, however, are not optimal for all species. On the other hand, the species tested here are bird species which actively move in the landscape. Species may still be present in those pixels that have low suitability caused by land use intensity but with suitable pixels nearby. Some species may have very specific habitat requirements so that the choice of species also matters. Since we only modeled nine species as an example of the application of the new land systems classification, no large inferences can and should be made based on these examples. The use of species presence data obtained from iNaturalist may also bring some uncertainty in the analysis due to the way in which points were reported and stored.

Comparison with other land cover or land system maps in Europe

Quantifying the intensity of human land use activities and integrating with land cover and land use data remain a challenge for land system scientists. Several efforts have improved the representation of the spatial pattern of human-environmental interactions, including the global representation from Ellis & Ramankutty (2008) and van Asselen & Verburg (2012), the classification of agro-silvo-pastoral mosaic systems in Mediterranean region (Malek and Verburg 2017), grassland management (Estel et al. 2018) and the archetype maps for Europe (Levers et al. 2018). Each of these efforts selected sensible indicators and methods that operate at the corresponding scales and serve to the project purposes. For example, the classification of European archetypical patterns is based on 12 land cover and land use intensity indicators, including a previous version of nitrogen application rate data used in this study. However, rather than an expert-based classification an automated clustering technique (i.e., self-organising maps, SOMs) that calculates the similarities across indicators was performed. van der Zanden et al. (2016) have compared automated and expert-based classification methods and concluded that although large agreement between the two methods exists, sometimes the automated method may overlook small differences in certain categories that are vital for landscape function based on the experts’ opinion. Furthermore, expert-based classification allows a consistent classification for present and future projections of the land system while upon change an automatic classification would also require an update of the classification itself as it is representing the statistical structure of the data. Although our map shares some characteristics with previous maps, it is the first time that a classification system was developed based on those landscape and land use management properties that are important to biodiversity assessment and species distribution modeling.

Robustness and uncertainties

It is challenging to assess the performance of any land system classification. First, the underlying data used for classification include inherent errors and uncertainties. For instance, large areas in the coastal areas of northern England and Ireland were classified as peat bogs in CLC map, hence in our map they were aggregated into wetlands. Some of these areas may be used as pastures in reality and should be classified as grassland systems. Another example is the natural grassland class in CLC, some of which are used as grazing land in mountain areas during part of the year. Furthermore, uncertainty within the intensity indicators may be more pronounced because they are mostly downscaled from statistics or based on different proxies and assumptions. For instance, the nitrogen application rate we used in this study was based on CAPRI (Common Agricultural Policy Regionalised Impact) model result. CAPRI allocated nitrogen input in NUTS2 region to 1-km2 cells based on environment conditions and crop types. The uncertainty from this downscaling work was therefore inherited in our land system map.

Second, the definition of land systems and land use intensity should be evaluated relative to its potential use. Therefore, a traditional accuracy-oriented validation practice will not make sense as a classification is anyhow simply a way to partition variation in a dataset rather than a novel dataset in itself.

In addition to the inclusion of indicators, there is little evidence revealing the (non-)linear relationships between different indicators and species richness. Therefore, some of our thresholds are arbitrary for the continuous variables. In some cases, we used top 25% quantile (e.g., wood production) that is statistically meaningful. In other cases, we used an arbitrary value (e.g., 50kgN/ha) as the threshold to keep consistent with previous studies. We conducted a series of sensitivity analysis on the measure of intensity of land systems. Doing so indicates to what extent the choice of threshold influences the land system classification outcomes, which in our case shows only marginal differences (results in Figs. S10, S11).

Conclusion

The new land system classification, with the most complete coverage of Europe at 1-km2 resolution, resulted in 26 classes of seven major and mosaic land systems. In particular, the forest/shrub mosaic systems describe fragmented forest habitats that are crucial for species while the different land management intensity classes can show profoundly different impacts on biodiversity. The spatial distribution of these classes shows distinctive patterns across Europe: vast low and medium intensity land systems in Eastern Europe, more developed and intensively used land in Western Europe, a more heterogeneous land systems in Southern Europe, and large areas of forest but mostly under medium-intensity in the north. Although designed for Europe, this practice can also be used as an example for other continental and global scale land system classifications that aim for biodiversity assessment or other purposes.

The representation of landscape characteristics in this new land system classification provides crucial information for biodiversity studies. Demonstrated by the species distribution modeling results, the predictive ability of models and variable contribution across different species using the land system map shows added value as compared to traditionally used land cover maps. Although we only analyzed a small set of species, the results presented in this study are promising, given evidence of the importance of counting landscape and land system characters when assessing biodiversity at large-scale. This way both land system science and landscape ecology can better take stock of their human-environmental systems approach in improving the way they represent landscape characteristics in further assessments.