Nature ’ s contributions to people and biodiversity mapping in switzerland: Spatial patterns and environmental drivers

Changes in climate and land use represent significant risks of biodiversity loss globally, affect ecological stability, impact nature ’ s contributions to people (NCP, i


Introduction
Ecosystems are facing severe threats worldwide, compromising the supply of many life-supporting goods and services to people.These threats, driven by human activities, are currently posing a greater risk of global extinction to more species than at any previous time (IPBES, 2019;Canadell & Jackson, 2021).The erosion of habitats and biodiversity leads to the decline, alteration, and shift in the provision of nature's contributions to people (NCP) (Janssen et al., 2021;Walsh et al., 2016), thereby endangering the stability of human livelihood.The challenges of preserving biodiversity and NCP thus require multifaceted approaches that address social, environmental, and economic aspects (Biber-Freudenberger et al., 2020;de Queiroz-Stein & Siegel, 2023).Numerous attempts at protecting biodiversity have failed to slow down the rate of biodiversity loss (IPBES, 2019;Mace et al., 2018), emphasizing that conservation efforts should now be planned holistically, understanding and considering both the interdependence of species as well as the necessity to accommodate human needs (Mace, 2014;Schlaepfer & Lawler, 2022).Indeed, biodiversity carries a major role in the supply of multiple NCP (Bastian, 2013;Harrison et al., 2014;Rey et al., 2023), thus the importance to consider both for conservation initiatives.In this context, the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES) introduced multiple values of nature in its assessments, encompassing nature's values (i.e.biodiversity), nature's contributions to people (NCP; i.e. ecosystem services), and good quality of life (IPBES, 2018).
Evidence reporting the state of these values of nature is essential for designing policies, implementing strategies and practices, and driving behaviours to help achieve the Sustainable Development Goals (UN DESA, 2016), the Aichi Biodiversity Targets (Convention on Biological Diversity, 2016) and the Paris Agreement on Climate change (United Nations Environment Programme, 2015).
In particular, spatially explicit information across large scales reveals patterns that allow moving beyond site-specific understanding to design policies and institutions that work across administrative borders, and focus on ecosystem functionalities and the long-term supply of their services.Mapping NCP and biodiversity plays an important role in supporting effective management processes (e.g.: Beier et al., 2008;Daily et al., 2009) and landscape planning (e.g.Frank et al., 2012;Honeck et al., 2020;Koschke et al., 2012).It also facilitates the assessment of trade-offs and synergies among different NCP and biodiversity components (Ramel et al., 2020;Sylla et al., 2020;Yan et al., 2021;Zhao & Li, 2020), which is crucial for optimizing resource allocation in conservation initiatives (Aryal et al., 2021;Bennett et al., 2009) and to better understand the different relationships between biodiversity and NCP.In addition, mapping of NCP and biodiversity allows, through clustering in bundles, the identification of their joint supply providing a spatial representation of their co-occurrence (Karimi et al., 2021;Qiu & Turner, 2013;Zhang et al., 2022).Bundles were defined by Raudsepp-Hearne et al. (2010) as "sets of ecosystem services that repeatedly appear together across space or time", and allow to pinpoint the key drivers of their distribution.Furthermore, understanding which drivers impact the distribution of NCP and biodiversity bundles is necessary to predict their future evolution (Meacham et al., 2022;Spake et al., 2017).
Biodiversity underpins ecosystems, which in turn provide NCP.Overall, reported relationships between NCP and biodiversity tend to be positive, but are shown to depend on location, specific indicators, and scale of analysis.Global studies (worldwide or continental level) have generally identified positive relationships (i.e.synergies) between biodiversity and NCP, particularly in regulating NCP (Cimatti et al., 2023;Ricketts et al., 2016;Soto-Navarro et al., 2020).For example, it was shown that higher levels of biodiversity can increase NCP provision of carbon sequestration, water quality regulation, and protection from organisms detrimental to humans (Cardinale et al., 2012;Hooper et al., 2005;Tilman et al., 2014).However, more diverse results are shown in regional assessments.For instance Crouzat et al. (2015) showed a positive spatial correlation between plant diversity and crop/wood production in the French Alps, while Kong et al. (2018) showed a negative relationship (i.e.trade-off) between biodiversity and material NCP in China.In addition, although research on functional relationships between biodiversity and NCP is scarce, positive relationships are generally found (Finney & Kaye, 2017;Isbell et al., 2018;Waldén et al., 2023), but not necessarily (e.g.Kleijn et al. (2015) and Winfree et al. (2015), for pollination).This emphasizes the need to integrate NCP and biodiversity considerations in conservation strategies (as shown in Peru by Móstiga et al., 2023) and across relevant scales (Chaplin-Kramer et al., 2022).This is particularly applicable in Switzerland, where the landscape has been extensively shaped by human activities, necessitating the promotion of coexistence between humans and nature due to their spatial proximity.Switzerland exhibits commitment and interest in integrating NCP in nature conservation.In particular, the nationwide biodiversity strategy explicitly states the conservation of ecosystem services in its overall objective (FOEN, 2012).Investigations of NCP in Switzerland have been undertaken through several research studies, with diverse approaches, scales, and resolutions.At the national level, Braun et al. (2018) mapped four NCP indicators, revealing synergies among the majority of them.Jaligot et al. (2019b) identified distinct and complementary bundles of NCP supplied by different Swiss cantons.
Additionally, researchers have investigated regional scales, providing insights into the spatial and temporal dynamics of NCP.For instance, these studies show that the distribution of NCP supply has not been constant over time in the Vaud Canton (Jaligot et al., 2019a), and that potential shifts are expected in the future due to land use and climate change in mountainous regions (Briner et al., 2012).Moreover, they show the possibility of integrating NCP with biodiversity in regional spatial conservation planning, underscoring that excessively prioritizing specific NCP could increase threats on biodiversity (Ramel et al., 2020).Stritih et al. (2021) also highlighted possible changes in NCP provision by mapping risks and uncertainties associated with NCP supply in mountain forests.
Yet, a comprehensive assessment of NCP and their relationship with biodiversity was still lacking in Switzerland.Our study bridged this gap by mapping a set of 15 NCP indicators, as well as a biodiversity indicator based on the modeled distribution of 1482 threatened species.Our specific objectives were to: 1) Create high-resolution models and maps for NCP supply indicators and one biodiversity indicator in Switzerland.2) Assess the relationships − trade-off or synergies − among NCP and between NCP and biodiversity.3) Identify bundles of NCP and biodiversity and their main drivers.
These analyses enhance our understanding of the relationship between NCP and biodiversity at a national scale, offering insights of their spatial patterns across Switzerland.Moreover, these spatial outputs can prove valuable within Switzerland in the context of multi-level planning (Albert et al., 2016).With this work, we aim to contribute to more efficient landscape management and effective conservation strategies.

Study area
Switzerland is a landlocked country located in Europe whose landscape is largely influenced by human activities.Switzerland is showing a tendency towards more settlements, urban areas, and forests and fewer agricultural areas (FSO, 2021a).One prominent feature of Switzerland is its large elevation gradient, spanning from a minimum altitude of 193 m to the Alps' highest summit at 4634 m.The country has been divided into five distinct ecoregions referred to as production region, as identified by the national forestry inventory (Fig. 1, FOEN, 2020a): the Plateau, Jura, Pre-Alps, Alps, and Southern Alps, each with specific ecological and topographic characteristics.

NCP mapping
We mapped the supply of 15 nature's contributions to people (NCP) indicators and one biodiversity indicator guided by the IPBES' classification (IPBES, 2018).Eight regulating NCP, four material NCP, three non-material NCP, and one biodiversity indicator were investigated (Table 2).We produced spatially explicit outputs in the form of raster maps with a 25 m resolution for all of Switzerland.The methodology used for mapping NCP indicators was based on different approaches, reflecting the evolving nature of NCP research (Burkhard & Maes, 2017;Martínez-Harms & Balvanera, 2012).Part of the mapping was based on causal relationships between biophysical information tables derived from primary data (e.g.field surveys) and land use, climatic, and topographic data.We also used regression models for several NCP, for example using niche-based models (Lavorel et al., 2017), which rely on the landscape suitability for a selected set of species associated with certain NCP (Rey et al., 2023).See 2.2.1 to 2.2.15 and Appendix A. Input data used for the mapping is listed in Table 1.

Habitat quality (HAB)
HAB informs on the repartition of natural habitats and their state of degradation.This NCP indicator was computed using InVEST habitat quality module (Natural Capital Project, 2022), which calculates an index value for each raster cell based on a relative habitat score, and its sensitivity and proximity to threats.Habitat scores (ranging from 0 for artificial to 1 for natural; see Table A.2) were set for each land use class using a classification done in Switzerland (GE-21, 2020), based on the urbanity index from O' Neill et al. (1988).Five binary threat layers were generated using land use, population density, and road networks."Urban" threat layer was created selecting settlements and urban land use classes in municipalities exceeding 10′000 inhabitants or densities surpassing 100 inhabitants/km2, while the "Rural residential" layer was based on the same land use classes but for municipalities falling below these thresholds."Primary roads" and "secondary roads" threat layers were generated by rasterizing road network vectors.Habitats sensitivity to threats was defined based on values from the literature (see A.1.1.).

Pollinator abundance (POL)
POL describes the potential for presence of wild pollinators species in a pixel, given the habitat nesting suitability based on land use, and the floral resources available within flying range.Pollinators species were selected based on Kleijn et al., (2015;see Table A.4 for the list).Output map was produced using InVEST pollination module (Natural Capital Project, 2022).

Removal of PM10 by vegetation (AIR)
AIR describes the potential regulation of air quality through removal by vegetation of PM 10 (particulate matter with a diameter ≤ 10 μm).We used a formula of PM 10 removal adapted from Nowak et al. (2006) based on Braun et al. (2018).It estimated the quantity of PM 10 filtrated by vegetation based on the pollutants deposition rates per leaf area surface (Table A.5, values from Remme et al., 2014), the type of land use, and the resuspension rate of PM 10 (see appendix A.1.3.).We coupled land use classes with the dominant leaf-type layer of forests, as well as a digital elevation model, to distinguish broadleaved from coniferous forests.

Carbon stored in biomass (CAR)
CAR estimates the amount of carbon stored in vegetation and soil, thus contributing to the regulation of climate.We used InVEST carbon module (Natural Capital Project, 2022), which attributes a value of elemental carbon stored in each raster cell based on a land use map and a correspondence table containing the mass of carbon stored by land use categories (Table A.6).The values of carbon stored per land use categories were based on the Swiss greenhouse gas inventory of the period 1990-2018(Table 6-4, FOEN, 2020b).The Swiss territory was divided by production region (FOEN, 2020a) and elevation regions (<601 m, 601-1200 m, >1200 m), and each combination (15 total) was computed individually to fit with the values provided by the greenhouse gas inventory.

Nutrient retention by landscape (NR)
NR describes the environment's filtration capability of nitrogen annually, as an indicator of the regulation of water quality.We used InVEST nutrient delivery ratio module (Natural Capital Project, 2022) to estimate the quantity of nitrogen retained by each pixel and thus not reaching the water streams.This was computed using a land use map, corresponding nutrient loads and retention properties, as well as annual water runoff and elevation data.We used information from Jaligot et al. (2019a) for land use biophysical table as well as model parametrization (Table A.7 and A.8).

Sediment retention by landscape (SR)
SR estimates the yearly amount of sediment that is retained by each pixel in the Swiss landscape, thus regulating the erosion of soils.We used InVEST sediment delivery ratio module (Natural Capital Project, 2022), which computes sediment retention using an elevation map, soil properties (erosivity and erodibility) as well as a land use map with corresponding biophysical properties of land use classes.Values of model parameters and biophysical table (Table A.9 and A.10) were based on the study done by Jaligot et al. (2019a).

Protective forests and floodplains (HAZ)
HAZ shows the location of forests and floodplains providing potential protection against natural hazards.The SilvaProtect-CH project provided modeled natural hazard data on rockfall, avalanches, landslides, flood, and debris flow (Losey & Wehrli, 2013).This data was overlaid with forested areas from the Swiss land use map (Table A.11) to identify protective forests potential.In addition to that, alluvial zones representing floodplains (FOEN, 2017) were added to the indicator map.

Pest control species (PC)
PC describes the combined habitat suitability of species identified as predator species to main agricultural pests.To map this indicator, we used a list of 50 predators (Table A.12), based on the study done by Civantos et al. (2012).It includes 2 amphibian, 2 reptile, 34 bird, and 12 mammal species.Invertebrates were not included in this analysis.We modelled the habitat suitability of each of these species and averaged the prediction maps to get the mean suitability value for predator species.The detailed species distribution modelling process is described in section 2.3, and specific model performances are reported in Appendix A.2.

Annual water yield (WY)
WY describes the relative contribution of each pixel to the water yield of the watershed the pixel is on.This indicator was computed with InVEST water yield module (Natural Capital Project, 2022) using annual precipitations, evapotranspiration, soil properties (depth and water availability), and land use.The model was calibrated using data at the watershed level from official hydrological surveys (Schädler & Weingartner, 2002), for 287 available watersheds (Fig. A.1), by modifying the Kc parameter (crop coefficient) in the biophysical table (Table A.13).

Wood provision potential (MAT)
MAT shows an estimation of annual forest growth, as an indicator of the potential supply of wood-based energy and material.We used annual forest growth data from the Swiss greenhouse gas inventory (Table 6-15, FOEN, 2020b) to map this indicator.Values of average wood increment in m 3 (Table A.15) were attributed to each pixel of land use corresponding to productive forest (Table A.14).Switzerland was divided in elevation and production region to use the corresponding value of forest growth (similarly to 2.2.4 -CAR).Pixels located on a slope superior to 110 % were not considered in the analysis as they are not suitable for wood harvesting (Dupire et al., 2015).

Landscape suitability for agriculture (FF)
FF describes the landscape climatic and edaphic suitability to cultivate crops.We selected a list of the most cultivated crops in Switzerland (Table A.16) based on data from the Swiss Farmers Union (USP, 2021).We used the ECOCROP database (FAO, 2007) to extract species-specific optimal growing conditions (monthly precipitations, temperature, and soil pH), similarly to Briner et al. (2012).We modeled optimal growing conditions maps for each species using the "Recocrop" package on R (Hijmans, 2021).The obtained maps represented climatic and edaphic suitability for the selected crops and were aggregated by averaging the value of crop maps.They were then masked to be applied only to food and feed production land use classes (Table A.17).

Medicinal plants (MED)
MED describes the combined habitat suitability of 380 wild plant species identified as having a medicinal potential.A list of medicinal plants was created based on studies from Dal Cero et al. ( 2014) and Rey et al. (2023) (Table A.18).We modelled the habitat suitability of each of these species and averaged the prediction maps.The detailed species distribution modelling process is described in section 2.3, and specific model performances are reported in Appendix A.2.

Landscape suitability for picture-taking (LI)
LI represents the modeled landscape potential for picture-taking linked with nature.We used publicly accessible data from two photo sharing apps for photography (Flickr) and naturalist observations (iNaturalist).We collected geolocation of pictures taken between 2006 and 2021 (Flickr) and 2010 and 2021 (iNaturalist), for a total of 6855 and 3719 observations, respectively.We used automatic image annotation (Schwemmer, 2021) on the geo-referenced pictures to remove pictures not depicting natural elements (e.g.drawings, vehicle pictures), enhancing the general quality of the observations (Fox et al., 2021).We then explained the distribution of these geolocated pictures using a set of environmental predictors (Table A.20), by building a regression model using the "randomForest" R package (Liaw & Wiener, 2002).We used this model to predict the landscape potential suitability for nature picture-taking to the entire Swiss territory.We conducted the analysis with 5-fold cross-validation.For each fold, we assessed the model's performance through area under the curve (AUC), mean square error (MSE) and R-squared.Model performances are shown in Table A.19.

Recreation potential (REC)
REC describes the landscape potential for outdoor recreation.This indicator was mapped similarly to ESTIMAP's "recreation potential" geomatic model done for the European union (European Commission, 2013).Mapping was based on three landscape characteristics: degree of naturalness (DN), natural protected areas (NP), and water components (W).DN was generated by attributing habitat scores to land use categories, as in section 2.2.1 -HAB.NP was created using areas of protected areas of Switzerland, from the TLM_PA item of swissTLMRegio database (swisstopo, 2020b).W was assessed by computing an inverse relative distance to lake coasts, getting the highest value at lake coast and a decreasing value for 2 km.The REC map is a normalized aggregate (sum) of the three landscape characteristics maps.

Emblematic species (ID)
ID describes the combined habitat suitability of 15 species identified as emblematic/iconic, and thus contributing to support cultural identity.The list of species was based on the study from Schirpke et al. (2018) conducted in the Alpine region and shown in Table A.21.We modelled the habitat suitability of each species and averaged the prediction maps to get the mean suitability value map.The detailed species distribution modelling process is described in section 2.3, and specific model performances are reported in Appendix A.2.

Species distribution modelling and biodiversity indicator mapping
We constructed species distribution models (SDMs) for 1482 terrestrial species that have been identified as threatened and classified in the red list by the Swiss government (red list species).The selection of these species for conservation efforts in Switzerland is based on an assessment of the threat level to their populations within the country (FOEN, 2021b).Selected species consist principally of vascular plants (n = 596), arthropods (n = 369), fungi (n = 267), birds (n = 83), and ferns and mosses (n = 62) (see Fig . A.2), and the complete list of selected species along with a detailed description of the method are available in appendix section A.2.
SDMs were built using the N-SDM software (Adde et al., 2023a), which allows the integration of a "global" model, quantifying the species response to bioclimatic conditions across its entire distribution range, with a "regional" model incorporating finer-scale habitat predictors in Switzerland.Ensemble SDMs were computed using the five modelling algorithms available in N-SDM (Generalized Linear Model (GLM), Generalized Additive Model (GAM), Maxnet (MAX), Random Forest (RF), and light Gradient Boosted Machine (GBM)).Candidate environmental predictors used for modelling the distribution of each species were extracted from the "SWECO25" database (v.1.0)(Külling, Adde, et al., 2024) and automatically selected using the "covsel" procedure (Adde et al., 2023b).Model selection and evaluation was achieved using a consensus 'Score' metric averaging the AUC′ (or Somers' D, such as AUC′ = AUC * 2 − 1), the maxTSS, and the CBI (Adde et al., 2023a).Results from the modelling algorithms were mapped over Switzerland on a 25 m resolution grid and ensembled together by averaging the five maps for each species.All information on the parameters used to fit and evaluate the models was stored using the ODMAP protocol (Zurell et al., 2020;appendix B).
Individual species maps were aggregated using the Zonation 5 (v.1.0)prioritization software (Moilanen et al., 2022).Zonation is originally designed to support ecologically based landscape planning, in combining different landscape elements by an algorithm that systematically removes grid cells with the smallest aggregate loss of conservation values in each iteration.By prioritizing areas based on their relative importance for the conservation objective (in our case, red list species), Zonation produces a hierarchical spatial priority map ranging from 0 for minimum priority to 1 for maximal priority (Lehtomäki & Moilanen, 2013).Here, it was used to aggregate the individual red-list species maps produced by the SDMs into a biodiversity indicator (BD).This aggregation was done using the Core-Area zonation 2 (CAZ2) algorithm.The choice of this algorithm was made to strike a balance between achieving high average coverage across all species and capturing the high-occurrence areas of each individual species (Moilanen et al., 2022).

Relationship and bundle analysis
To investigate the relationship between NCP and between NCP and BD, we conducted a correlation analysis to identify potential trade-offs or synergies between the 120 possible pairs of the NCP and BD indicators.We applied min-max normalization to the NCP maps and the BD map to ensure comparability (e.g.Yu et al., 2022) and performed random sampling to select 15,000 points from the study area.To assess the correlations, we conducted Spearman rank correlation (ρ) between each pair of variables (e.g.Zhang et al., 2022).The alpha level for statistical significance was set at 0.05.
We employed a bundle analysis to discern distinct patterns of NCP supply and BD across the landscape (Raudsepp-Hearne et al., 2010).
Here, bundle analysis allows to assess how these elements co-occur and interact over a given area.To do this, we used the k-means clustering algorithm, which has been employed widely in similar analysis (e.g.Cusens et al., 2023;Schirpke et al., 2019;Zhang et al., 2022).To determine the optimal number of bundles, we used the silhouette coefficient method from the "factoextra" R package (Kassambara & Mundt, 2020), by selecting the number of bundles showing the highest average silhouette width.Clustering was done using the "RStoolbox" R package (Leutner et al., 2022) with a sampling of 15,000 points.Built-up areas were masked from the maps used in both the relationship analysis and the bundle analysis.The mask used was based on built areas outlined by the Swiss building zone statistics (ARE, 2022).
To better understand the driving factors behind the identified bundle patterns, we performed a classification analysis using environmental predictors to predict bundles membership, following a methodology similar to that employed by Schirpke et al. (2019).To do that, we selected a set of predictors from the "SWECO25" database (v.1.0)(Table A.25, Külling, Adde, et al., 2024) from a range of topographic, edaphic, habitat type and climatic datasets, similarly to what was done previously in Switzerland in similar large scale environmental clustering study (Lehmann et al., 2010).We extracted values from the bundle map and the predictors using random-stratified sampling, with 10,000 points per bundle category.We excluded predictors displaying a high correlation (|r| > 0.7) (Dormann et al., 2013) and used a random forest classification algorithm using the "ranger" implementation in the "caret" package in R (Wright & Ziegler, 2017;Kuhn, 2008, respectively).We assessed model performance on 30 % of the dataset using several evaluation metrics including Cohen's Kappa, multiclass AUC of ROC, and confusion matrix, where high values of these metrics (>0.8) represent strong agreement between predicted and observed classifications (Fielding & Bell, 1997).
From the model, we identified the most influential predictors based on the mean decrease in accuracy, both for individual bundle classes and for the overall classification model.A higher mean decrease in accuracy for a predictor indicates that its inclusion improves the model's ability to make accurate classification (Archer & Kimes, 2008).We then calculated the marginal effects of the most important variables using the "pdp" R package (Greenwell, 2017) and displayed it with partial dependence plots, which illustrate the relationship between each predictor and the predicted outcome of the classification model, to visualize their individual effect on the prediction.To display the distribution of climatic and land use predictors by bundles, we used density plots for continuous predictors and cumulative barplots for discrete land use/ cover classes.

NCP and biodiversity maps
We produced maps of 15 NCP and one BD indicator, among which distinct spatial patterns emerged (Fig. 2).Pollinator abundance, protective forests and floodplains, sediment retention by landscape, and emblematic species exhibited analogous distributions, primarily associated with mountainous regions (Alps, Southern Alps, Prealps, and Jura), yet limited to moderate elevations.In contrast, other services, such as habitat quality and recreation potential, revealed a more uniform spread across the mountainous areas.Similarly tied to these regions, annual water yield and sediment retention by landscape displayed higher supply at greater elevations.The Plateau region displayed a higher supply of NCP linked to agricultural production: landscape suitability for agriculture, pest control species, as well as medicinal plants.
In contrast, removal of PM 10 by vegetation, carbon stored in biomass, and wood provision potential collectively exhibited similar distributions aligned with forest presence.Lastly, landscape suitability for picturetaking illustrated a relatively consistent and unique distribution throughout Switzerland, closely linked to human infrastructure such as cities and roads.
The BD indicator map displayed scattered patterns covering all regions of Switzerland, with high values of priority score given in Alpine valleys (especially in the Rhône valley), and mean values of 0.62 ± 0.3 for the Jura, 0.41 ± 0.32 for the Plateau, 0.44 ± 0.28 for the Pre-Alps, 0.44 ± 0.3 for the Alps and 0.41 ± 0.3 for the Southern Alps.The mean score metric for all modeled SDMs was 0.88 ± 0.06 (Fig. A.3).
Based on the silhouette coefficient, the optimal number of bundles was found to be four, with an average silhouette width of 0.38 (Table A.26).Fig. 4.A provides a visual representation of the analysis outcome, in which all NCP and BD pixel values are projected to their corresponding bundle centers.Spatial distribution of these 4 bundles across the Swiss landscape revealed uneven patterns.Bundle 1 predominated in steep and mountainous regions, while bundle 2 was prevalent in the Plateau region and Alpine valleys.Bundle 3 displayed a distribution across the Plateau and Jura regions, and bundle 4 was notably dominant in higher-elevation alpine areas, and some Jura regions.
Fig. 4B provides a heatmap displaying the mean values of NCP and BD per bundle and the distribution of land use categories in each bundle.Bundles 1 and 3 primarily included forested areas, encompassing diverse NCP.Bundle 2 was primarily composed of NCP associated with agricultural land, including FF and PC, but also included MED and LI.Bundle 4 was made of a broader range of land use categories, including nonvegetated natural areas such as scree, rocks, and glaciers, as well as agricultural lands.This bundle comprised primarily regulating NCP (HAZ, HAB, SR, WY), along with non-material NCP such as REC and ID.The BD indicator was represented in all bundles but the 4th, with the highest mean value in the 3rd bundle.
Using environmental predictors to explain and predict bundle classification, our random forest analysis accurately classified 86.03 % of the original bundle labels.Evaluation of the model's predictions showed a value of Cohen's kappa of 0.81 and a multiclass AUC of ROC of 0.86.The confusion matrix and per-bundle metrics are shown in Table A .23 and A.24. Predictors displaying the highest correlation (mean decrease in accuracy values) with the spatial distribution of the four bundles were slope, vegetation height, average annual temperature, and light availability (Fig. 5A).The computation of marginal effects of the six main variables is shown in Fig. 5B through partial dependence plots.Bundles 1 and 3 displayed contrasting responses to topographical and climatic conditions (slope, temperature), but shared similar ecological features in term of vegetation composition, especially in the presence of high vegetation, through "beech forests" (Fagus sylvatica) habitat (especially for bundle 3) and "highland coniferous forest" habitat (for bundle 1), as opposed to bundles 2 and 4 which showed a similar response to low vegetation, and a habitat type of "grassland and meadows".
The per-bundle distributions of main identified climatic and land use predictors are displayed in Fig. 6.Continuous predictors were precipitation, temperature, and vegetation height (Fig. 6A-C).The highest mean value for temperature was in bundle 3 (9.16• C), while the lowest was in bundle 4 (2.65 • C).For precipitation, the highest mean value was in bundle 4 (1323 mm), while the lowest was in bundle 3 (984 mm).bundle 3 had the highest mean vegetation height (28 m), whereas bundle 4 had the lowest (4 m) (see Table A.27 for all descriptive statistics per bundle).For discrete land use predictors, "beech forest" predominantly occurred in bundles 1 and 3, "grassland and meadows" in bundles 2 and 4, and "highland coniferous forest" in bundle 1 (Fig. 6D).

Implications for conservation planning
We provided a first comprehensive set of 16 high-resolution maps of indicators of NCP supply and BD, based on the IPBES regional assessment of values of nature (IPBES, 2018).This development of spatially explicit indicators is essential for achieving a more effective balance between ongoing land use changes and ecological integrity, as highlighted in previous research (e.g.Bateman et al., 2013;Grêt-Regamey et al., 2017).
We identified three target levels at which our results can be practically used by stakeholders: cantonal, national, and European.
At the cantonal level, the Swiss government has mandated all cantons to identify, plan and develop their ecological infrastructure (EI) in order to ensure the protection of biodiversity (FOEN, 2021a).This task will benefit from the availability of our NCP in complement to a biodiversity-only approach (as in e.g.Vincent et al., 2019).For example, the canton of Geneva is using NCP indicator maps for the planification of its EI (DETA-DGAN et al., 2018).
At the national level, we provided here the first national-scale assessment for Switzerland, based on a comprehensive list of NCP at high resolution and including an index of biodiversity based on a large number of endangered species, paving the road for a nation-wide spatial prioritization.In practice our maps could be used as an input for a nation-wide prioritization of the EI, as approaches which encompass both NCP and biodiversity indicators for the identification of the EI through weighted spatial prioritization have been tested regionally in Switzerland (e.g.Honeck et al., 2020, Ramel et al., 2020), but not at the national level like done in India (Srivathsa et al., 2023) or at European scale (Kukkala & Moilanen, 2017).NCP maps and knowledge on their synergies and trade-offs have been shown to be valuable to guide the planning and management of the EI at large scales (Chen et al., 2024;Liquete et al., 2015;Ramyar et al., 2020).
Finally, our assessment aligns with the "EU Biodiversity Strategy to 2020" goals, which emphasizes the mission of member states to "map and assess the state of ecosystems and their services" (European Commission.Joint Research Centre., 2020).Switzerland, not being an EU member, was not part of this specific initiative.Nevertheless, the data and insights derived from our assessment can play an important role by inspiring similar iniatives across Europe and supporting the objectives set forward in the European "EU Biodiversity Strategy for 2030" (European Commission.Directorate General for Environment., 2021).

BD and NCP relationship
We assessed the relationships between NCP and BD maps following the "spatial linkage" approach, comparing values of NCP and BD across space (Rey et al., 2022;Ricketts et al., 2016).Our analysis revealed a positive correlation between the spatial distribution of the BD indicator and most of the NCP indicators.This was expected as several studies have shown the positive relationship between biodiversity and NCP supply in different ways (Isbell et al., 2018;Ricketts et al., 2016;Waldén et al., 2023).Notably, our results show the strongest correlation between BD and two NCPmedicinal plants and pest control species − which are niche-based models (Lavorel et al., 2017), derived from SDMs and sharing very similar input data.While other correlations are noteworthily positive, they are relatively modest.Thus, our findings highlight the significant yet partial link between BD and NCP supply in Switzerland.This result is congruent with previous studies, and encourages targeted approaches for conservation of both biodiversity and NCP (Cimon-Morin et al., 2013;Larsen et al., 2011;Naidoo et al., 2008;Xu et al., 2017).

Relationship between NCP, bundles, and drivers
Our findings indicated significant relationships among the spatial distribution of most NCP in Switzerland, which is coherent with what was found in similar studies (Crouzat et al., 2015;Kong et al., 2018;Raudsepp-Hearne et al., 2010;Schirpke et al., 2019).Our results regarding water yield (WY) diverge somewhat from existing studies conducted in the Alpine region.For example, Crouzat et al. (2015) reported a mild yet significant synergy between water quantity regulation and the pest control NCP (PC), and Schirpke et al. (2019) display differing results regarding the relationship between WY and wood provision potential (MAT), carbon storage (CAR) and soil erosion prevention (SR).These contrasts likely stem from variations in scale, geographical region (especially latitude) and resolution used in our analysis.Specifically, highest WY values in Switzerland were most concentrated in places characterized by limited vegetation and species diversity (high altitudes, glaciers, see Fig. 2), and thus appear in a trade-off with NCP linked to forests (CAR, MAT), vegetation in general (SR) and cropland-specific species (PC).
We identified four distinct bundles of NCP and BD for Switzerland in our analyses.This number aligns with Schirpke et al. study (2019), which employed a similar methodology.Dittrich et al. (2017), in the context of a national-scale assessment for Germany, identified eight bundles, by using the same methodology as Crouzat et al. (2015) who found five bundles.This study-specific difference in the number of bundles may be influenced by the spatial extent of the study area.In our case, Switzerland (four bundles) is of a much smaller size compared to Germany (eight bundles).We explored the repartition of bundles based on the approach of Raudsepp-Hearne et al. ( 2010), including NCP and biodiversity components, as recommended by Crouzat et al., (2015).This approach facilitates the overall understanding of NCP and BD distribution.This distribution in each bundle was found to be very diverse, highlighting the multiplicity of roles of supply played by ecosystems.However, each bundle appeared to be supported by specific types of land use, suggesting the vulnerability of their provision to changes in land use, which has been shown in Austria by Schirpke et al. (2023).
Our bundle classification model performed well, allowing us to efficiently highlight the main drivers of bundle's distribution.These drivers are predominantly associated with climate, topography, and  (Sun et al., 2023).Our BD indicator being based on a spatial prioritization, high BD values were not expected at higher elevation, as species diversity decreases with altitude (Sanders & Rahbek, 2012), However, Bundle 4 (linked to high elevations) having the highest habitat quality (HAB) value, remains crucial for biodiversity (e.g. for alpine species, Vittoz et al., 2013), although having low values for our red list species prioritization.
Our bundle approach allows to consider simultaneously NCP and BD in a holistic view (Saidi & Spray, 2018), with one single geographic information layer for the Swiss landscape.Our work can help promote N. Külling et al. social and ecological considerations in nature conservation (Dias et al., 2022), for example highlighting the diverse NCP and BD features provided by one area (as shown in forests by Deal et al., 2012).Here, our results underscore the importance of "beech forest" and "grassland and meadows" as habitats in Switzerland, which are main drivers to the distribution of all four NCP and BD bundles identified.These informations can be leveraged for landscape management decisions that protect both NCP and BD.Indeed, this approach of landscape characterization in bundles has been shown as an efficient way to communicate the numerous functions of the landscape and can be used for the development of informed landscape management plans (Bai et al., 2021;Liao et al., 2023;Malmborg et al., 2021).

Limitations of the study
Our NCP supply assessment does not capture the actual benefits of NCP to people, hindering our ability to identify gaps (e.g. in land use policies) (Mandle et al., 2021).The development of spatially explicit demand indicators (e.g.Schirpke et al., 2019) for Switzerland is necessary to further the connection between human well-being and natural environment.
Large-scale NCP assessments are necessary (Albert et al., 2016;Schröter et al., 2016), but can lack crucial information necessary at smaller scales.For instance, analyses of socio-ecological-technological systems in Geneva canton showed varying drivers and archetype (comparable to bundles) numbers, highlighting a need for more archetypes than at Swiss level (Wicki et al., 2023).A complementary alternative lies in tiered assessments at different spatial scales, allowing to better understand the variety of NCP provided and to design conservation strategies accordingly, with adapted methods for regional and national needs (Grêt-Regamey et al., 2015).
Additionally, though we use a comprehensive list of NCP defined by the IPBES, the choice of NCP and BD indicators are subject to data availability and technical limitations, and model outputs often lack uncertainty measures.This complicates comparability between NCP assessments (Schirpke, Ghermandi, et al., 2023), and we thus advocate for an open access to scripts and data in such assessments, as done here.Finally, although an assessment of NCP and BD was necessary to identify the current drivers of their spatial distribution in the landscape, a look into the future evolution of the identified drivers and bundles is needed to ensure the preservation of NCP and BD supply, and to forecast potential changes (Schirpke et al. (2023)).

Conclusions
In this study, we show significant correlations among the spatial distribution of NCP and BD indicators studied in Switzerland.Although closely linked, NCP and BD are not entirely collinear, emphasizing the importance of informed and targeted approaches in the conservation of both NCP and BD supply.We identified four bundles that represent NCP and BD supply in the landscape, and we showed that they are reliably correlated to a set of environmental predictors.Our findings reveal that climatological drivers, such as temperature and precipitation, along with habitat types (forests and meadows), play a major role in the distribution of these four bundles.This indicates that future global changes will likely have a significant impact on the spatial distribution of BD and NCP supply.Finally, we provide a set of spatially explicit outputs and methods, which could serve as valuable inputs to refine conservation priorities in Switzerland by promoting the integration of NCP, particularly through weighted spatial prioritization.

Fig. 2 .
Fig. 2. Normalized distribution of indicators maps.Color schemes defined by Jenks natural breaks.BD − Red list species, HAB − Habitat quality, POL − Pollinator abundance, AIR − Removal of PM 10 by vegetation, CAR − Carbon stored in biomass, NR − Nutrient retention by landscape, SR − Sediment retention by landscape, HAZ − Protective forests and floodplains, PC − Pest control species, WY − Annual water yield, MAT − Wood provision potential, FF − Landscape suitability for agriculture, MED − Medicinal plants, LI − Landscape suitability for picture-taking, REC − Recreation potential, and ID − Emblematic species.

Fig. 3 .
Fig. 3. Correlogram depicting the relationships between nature's contributions to people (NCP) and biodiversity indicators.Each cell represents the Spearman correlation coefficient (ρ), with non-significant relationships indicated by a cross.BD − Red list species, HAB − Habitat quality, POL − Pollinator abundance, AIR − Removal of PM10 by vegetation, CAR − Carbon stored in biomass, NR − Nutrient retention by landscape, SR − Sediment retention by landscape, HAZ − Protective forests and floodplains, PC − Pest control species, WY − Annual water yield, MAT − Wood provision potential, FF − Landscape suitability for agriculture, MED − Medicinal plants, LI − Landscape suitability for picture-taking, REC − Recreation potential, and ID − Emblematic species.

Fig. 4 .
Fig. 4. A Spatial distribution of the four bundles across Switzerland.4.B (left) Alluvial diagram with flow thickness indicating the distribution of land use categories repartition in each bundle, with 4.B (right) a heatmap indicating bundles mean values for the nature's contributions to people (NCP) and biodiversity indicators.BD − Red list species, HAB − Habitat quality, POL − Pollinator abundance, AIR − Removal of PM10 by vegetation, CAR − Carbon stored in biomass, NR − Nutrient retention by landscape, SR − Sediment retention by landscape, HAZ − Protective forests and floodplains, PC − Pest control species, WY − Annual water yield, MAT − Wood provision potential, FF − Landscape suitability for agriculture, MED − Medicinal plants, LI − Landscape suitability for picture-taking, REC − Recreation potential, and ID − Emblematic species.

Fig. 5 .
Fig. 5.A Predictors importance for random forest bundle classification.The ten most important variables are shown, with the mean decrease in accuracy for the classification model indicated by the dashed black line, and the colored dots representing the mean decrease in accuracy for each bundle.A higher mean decrease in accuracy indicates that its inclusion improves the model's ability to make accurate classification.5.B Partial dependence plots for the six most important variables.Colored lines display bundle-specific responses to the selected variables over their value range.

Fig. 6 .
Fig. 6.A-C Density plots showing the distribution of temperature, precipitation, and vegetation height per bundle.Vertical black line represents the median value.6. D Cumulative barplot showing the relative distribution of land use categories per bundle.

Table 1
Input data used for NCP mapping.

Table 2
List of nature's contributions to people (NCP) indicators, biodiversity (BD) indicator, and codes used in this study.See 2.2.1 to 2.3 and appendix A.1 to A.2 for mapping methods.