Ecosystem service coproduction across the zones of biosphere reserves in Europe

ABSTRACT Biosphere reserves (BR) balance biodiversity protection and sustainable use through different management restrictions in three zones: core areas, buffer zones, and transition areas. Information about the links between zoning and ecosystem services (ES) is lacking, particularly in terms of the relative roles of natural contributions (ecosystem properties and functions) and anthropogenic contributions (human inputs such as technology and infrastructure) in coproducing ES. This study aimed to: (1) analyse how coproduction of four ES (crop production, grazing, timber production, recreation) differs across the three zones of BRs; and (2) understand which predictors (zoning, natural and anthropogenic contributions, other environmental characteristics) best explain ES provision within BRs. To do this, we collected spatial data on 137 terrestrial BRs in the European Union and on 16 indicators of ES coproduction. We used non-parametric pairwise Wilcoxon rank sum tests to calculate differences in indicators between zones. We used model selection and multiple linear regression to identify predictors of ES provision patterns. Anthropogenic contributions showed most differences between zones, with contributions generally increasing from buffer zones to transition areas. Natural contributions did not, on average, differ between zones, however, for recreation and crop production they decreased from buffer zones to transition areas. ES provision differed between zones only for crop production and grazing, which increased from buffer zones to transition areas. Regression analysis showed that natural contributions are the best predictors of ES provision for all four services. Our results indicate that zoning of BRs has implications for ES coproduction.


Biosphere reserves and ecosystem services
Sustainable management of protected areas and working landscapes is essential for halting biodiversity loss and safeguarding ecosystem services (ES). This goal is explicitly targeted in the United Nations Sustainable Development Goal 15 (Life on Land, United Nations 2015). Similarly, the Biodiversity Strategy of the European Union (EU) aims to protect 30% of land in the EU by 2030 (European Commission 2020), in accordance with the Post-2020 Biodiversity Framework of the Convention on Biological Diversity (UNEP 2020). However, only 10% is slated for strict protection with limited human modification (European Commission 2020), with the remaining 20% under less strict protection. Stringent protection of large landscapes can conflict with the production of provisioning ES for human livelihoods (Oldekop et al. 2016) and of cultural services like recreation. Therefore, integrated conservation strategies are needed to effectively deliver biodiversity protection and ecosystem functioning while allowing for the sustainable use of landscapes (see, e.g. Mace 2014;Palomo et al. 2014b;Kremen and Merenlender 2018).
Biosphere reserves (BRs) appear to be ideal examples of areas that can simultaneously meet conservation targets. BRs have been designated under UNESCO's Man and the Biosphere (MAB) Programme since 1976, with the goal of balancing biodiversity conservation and sustainable use (Bridgewater 2002;Bridgewater and Babin 2017;Reed and Price 2019). An important feature of BRs is the spatial zoning of core areas, buffer zones, and transition areas. Each of these zones has different functions, levels of protection and human modification (Table 1) (UNESCO 2021). Core areas are strictly protected and only permit a very low level of modification. Buffer zones usually surround, or are contiguous to, core areas and allow for limited human management and resource use. Transition areas are the least protected and aim for the sustainable use of ES.
Case studies of protected areas outside BRs have shown that zoning schemes have an effect on ES (Eigenbrod et al. 2010;Geneletti 2013;Kärkkäinen et al. 2020). Previous analyses investigating ES delivery in BRs have focused on a single case study (see, e.g. Plieninger et al. 2013;Bieling 2014;Palomo et al. 2014a;Lowell 2017;Negev et al. 2019). To the best of our knowledge, no research on the role of zoning in BRs on the provision of different ES has been conducted across a large number of BRs. This is partly due to a lack of data on the spatial zoning of BRs. As BRs are not an IUCN-listed category of protected area, there are gaps in the World Database on Protected Areas (WDPA, UNEP-WCMC) regarding BRs.

Ecosystem service coproduction
Given that BRs consist of adjacent zones with differing intensities of management and a focus on sustainable ES provision (UNESCO 2016), they are well suited for research on how different levels of protection (i.e. zoning) influence ES coproduction in socialecological systems (Figure 1). ES coproduction occurs through the combination and interaction of natural contributions, defined as biodiversity and ecosystem properties and functions; and anthropogenic contributions, defined as human inputs, such as technology, infrastructure, and knowledge (Table 1) Palomo et al. 2016). Both natural and anthropogenic contributions directly drive ecosystem service provision. Natural contributions include ecosystem properties and functions underpinned by biodiversity (Remme et al. 2014) and intermediate services (Lamothe and Sutherland 2018) that contribute to ES provision and use, such as natural soil fertility or biological pest control (Figure 1). Anthropogenic contributions can include the management or conversion of ecosystems, such as turning forest into agricultural land (Remme et al. 2014) and thus raising the ES supply capacity of ecosystems generally (Lavorel et al. 2020). Or they can include the mobilization of ecosystem functions by harvesting or accessing them, such as by building hiking paths to enable recreation (Lavorel et al. 2020).
Coproduction has been described conceptually Palomo et al. 2016;Rieb et al. 2017), and was prominently featured in the framework of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES, Díaz et al. 2015). However, few studies empirically disentangle anthropogenic and natural contributions to ES (see, e.g. Lavorel et al. 2020) and, to our knowledge, none investigate the role of coproduction on ES provision in conservation areas. As zoning can be seen as a proxy for different management restrictions across BRs, the magnitude and relative importance of both natural and anthropogenic contributions to the provision of different ES might vary across these zones (Table 1). For example, some anthropogenic contributions, such as the use of inorganic fertilizers and pesticides in agriculture, or heavy machinery in forestry, potentially conflict with biodiversity conservation and the provision of regulating and cultural ES (Storkey et al. 2012;Newbold et al. 2015). Thus, it remains unclear how Table 1. Definitions of key terms used in this study.

Coproduction
The combination and interaction of natural and anthropogenic contributions that lead to ecosystem service provision Palomo et al. 2016).

Anthropogenic contributions
Human inputs that directly drive ecosystem service provision including technology, infrastructure and knowledge. We analyse two types of anthropogenic contributions: ecosystem management, and harvest and access (Lavorel et al. 2020). Natural contributions Ecosystem inputs that directly drive ecosystem service provision, i.e. ecosystem properties and functions underpinned by biodiversity and intermediate services that underlie the provision of ecosystem services, such as soil fertility or pollination.

ES provision
The amount of ecosystem services generated per unit area derived from coproduction. Note that we do not distinguish potential and actual use. Commonly used synonyms for ecosystem service provision are ecosystem service supply or yield. Service providing area (SPA) The extent of land cover or land use classes in which a particular ecosystem service can be provided (e.g. forest land cover classes for timber production). Biosphere reserve An area designated under the Man and Biosphere Programme of the United Nations Educational, Scientific and Cultural Organization that aims to reconcile sustainable use and the protection of biodiversity and ecosystems through the use of a zoning system. Core area Strictly protected zone of a biosphere reserve focussing on the conservation of all aspects of biodiversity, usually located in the interior of a biosphere reserve. Buffer zone Area of a biosphere reserve which allows limited management and resource use, usually surrounding or adjacent to a core area.

Transition area
Area of a biosphere reserve in which sustainable management and resource use is practiced, usually located outside the buffer zone.
different BR zones perform in terms of ES coproduction and provision (Coetzer et al. 2014).

Study aims
In this paper, we analyse data on 16 ES indicators across 137 terrestrial BRs in Europe to achieve two main goals. First, we quantify differences in ES coproduction (ecosystem management, harvest and access, natural contributions) and provision of three material ES (crop production, grazing, timber production) and one non-material ES (recreation), to understand how these vary across the three types of BR zones. Second, using statistical models, we explore how patterns of ES provision across BRs are influenced by different drivers, including zoning, natural and anthropogenic contributions, and spatial and environmental characteristics. Regulating ES are not included in our analysis due to challenges in finding appropriate indicators for anthropogenic contributions to regulating ES across European BRs. With our study, we thus cover two aspects that have barely been studied quantitatively so far: the coproduction of ES by nature and humans and the effect of zoning on the provision of ecosystem services. We hypothesise that, for provisioning ES, the levels of anthropogenic contributions and final ES provision will increase from core areas (lowest levels) to buffer zones (intermediate levels) and transition areas (highest levels) as management restrictions decline ( Figure S1 in Supplementary Material). We expect natural contributions to increase from core to transition areas for grazing and timber production but to decrease for crop production and recreation.
We assume that soil biomass productivity, as an indicator of natural contributions to grazing and timber production, will be lower in core areas compared to transition areas, as strictly protected areas are usually not established on highly productive sites (Joppa et al. 2009). Pest control, as an indicator of natural contributions to crop production is expected to be higher in core areas, as less modified habitat likely provides better habitat for species that provide pest control. For recreation, we predict the highest provision in the buffer zone. This is because we assume that the buffer zone has the best combination of accessibility (i.e. recreation infrastructure is highest in the buffer zone) and natural beauty (highest in the core area and lowest in the transition area) (cf., Braat and ten Brink 2008;García-Llorente et al. 2012).

Collection of biosphere reserve zoning data
We considered all 186 BRs in the EU and the UK that belonged to the UNESCO World Network of BRs (UNESCO 2020) as of May 2020 as potential study areas. We obtained spatial zoning data on each BR from a variety of sources (Supplementary Material  Table S1). These included online databases with data at global, national, and regional levels. We also conducted web searches on regional BR networks and specific BRs using local languages for better results. Wherever online datasets were unavailable or out of date, we obtained missing data by contacting individual administration authorities, including BR administrations, regional authorities and national agencies. When no digital data were available or data were incomplete, such as point data only or missing data on zones, we digitized zone maps available on BR websites. We excluded BRs in Croatia, the Canary Islands, Madeira, Azores and overseas departments of France (e.g. Guadeloupe) due to insufficient coverage of the spatial indicator data for the selected ES.
We focused on terrestrial BRs. We also included the terrestrial portion of mixed terrestrial and marine BRs if all three zones were present in the terrestrial portion. To distinguish terrestrial from marine areas, we removed areas within each BR classified as 'Coastal lagoons' or 'Sea and ocean' in the 2018 CORINE land cover dataset (EEA 2020). This process yielded a total of 141 BRs. We then rasterized zone data into a 1 km grid, matching the reference grid of the European Environment Agency (EEA 2017a). We excluded BRs which had no core area, buffer zone or transition area after rasterization due to small patch sizes from the analysis (n = 4). The final BR dataset thus included 137 study sites. We digitized and edited BR zoning geodata in ArcGIS 10.8 (ESRI).

Spatially explicit indicators
We collected spatial data on indicators of anthropogenic contributions, natural contributions, and ES provision per area (henceforth: ES provision) for four services: crop production, grazing, timber production and recreation (Table 2). Two types of anthropogenic contributions were included: one related to ecosystem management, and another related to ES harvest and access (adapted from Lavorel et al. 2020). Indicator selection was driven by data availability, as we required contiguous, comprehensive, and consistent data covering the entire study area (EU and UK). When necessary, we prioritized indicators that could more accurately represent anthropogenic or natural contributions or ES provision over using the most recent data. Based on these criteria, we were able to assess the coproduction of three provisioning ES (crop production, grazing, and timber production) and one cultural ES (recreation) ( Table 2; Details in Supplementary Material S1).
We defined the extent of service providing area (SPA) for each ES by identifying the specific CORINE land cover classes that potentially provide each ES. This included cropland for crop production, pasture for grazing, and forest classes for timber production. The CORINE classes for potential recreation areas included forests, water bodies, mountains, coastal systems, and grassland, as these land cover classes are known to be of high recreational value (Paracchini et al. 2014;Martín-López et al. 2018). Full details are provided in Supplementary Material Table S2. Indicators for each ES were restricted to those areas that overlapped with the relevant SPA for each ES. We resampled all ES indicators to 1 km resolution, matching the EEA reference grid (EEA 2017a). We conducted indicator resampling and processing in ArcGIS 10.8 (ESRI).
To analyse ES coproduction, we calculated the median value of each indicator of ecosystem management, harvest and access, natural contributions and ES provision for each BR zone (n = 411, as there are 137 BRs and three zones in each BR). We chose medians of spatial indicators instead of means to reduce the influence of outliers. Calculating a single median value for each BR zone also reduced effects of spatial autocorrelation between pixels.

Statistical tests of differences in ecosystem service coproduction between biosphere reserve zones
We conducted statistical analyses in R, version 4.0.2 (R Core Team 2020). We used the 'raster' package (Hijmans 2020) for spatial calculations. We applied zonal statistics to calculate the mean and standard deviation of the proportion of SPA for each ES across the three zones of all 137 BRs.
To test for statistical differences in each indicator between core areas and buffer zones and between buffer zones and transition areas, we used nonparametric pairwise Wilcoxon rank sum tests. To account for differences between BRs, we paired zones of each BR. In total, we performed tests for 32 comparisons: two comparisons (core area vs. buffer zone and buffer zone vs. transition area) for each of the 16 indicators.

Predicting ecosystem service provision in biosphere reserves
We used multiple regression analyses and model averaging to assess the importance of different explanatory variables on ES provision for each of the four ES. Explanatory variables included indicators of anthropogenic and natural contributions (Table 2) and zone as a categorical variable (core area, buffer zone and transition area). We also explored an alternative model with two zones: combined core and buffer zone (1) and transition area (2): see Supplementary Material, Table S5. Further environmental or zoning covariates (see details in Supplementary Material S2) encompass terrain ruggedness (index between 0 and 1) to account for the accessibility of the terrain (Levers et al. 2016;Amatulli et al. 2018), proportion of SPA, and a categorical variable that described whether a BR met recommendations with respect to zone proportions (German National MAB Committee 2007; Austrian National MAB Committee 2016). Along with individual effects of each of the above variables, we also included two-way interactions between each coproduction aspect with zone, recommended spatial proportions of zones and terrain ruggedness.
ES provision was log-scaled for grazing, timber, and recreation to achieve normality of residual distributions. All numeric explanatory variables were scaled by dividing them by their standard deviations, and centred by subtracting their means through use of the 'scale' function in R. To account for multicollinearity of explanatory variables, we calculated the squared generalized variance inflation factor (GVIF^(1/(2*Df))) for all predictors and for each model. The ecosystem management indicator for grazing (and its interactions) was the only predictor with a squared GVIF^(1/(2*Df)) >10 (Marquardt 1970) and was thus removed from the grazing model. For the other models, we did not remove any predictors. We performed model averaging of the best models using Akaike Information Criterion (AIC) with the 'MuMIn' package in R (Barton 2018). Averaging was conducted for all models with an AIC delta <2 compared to the best model.

Spatial characteristics of biosphere reserves
The mean area of the 137 BRs is 154,624 ha (standard deviation 245,972 ha, median: 74,115 ha); the smallest BR has an area of 5,896 ha (Paúl do Boquilobo, Portugal) and the largest has an area of 2,395,591 ha (Bassin de la Dordogne, France) ( Figure 2). Core areas have the smallest mean proportion of the total BR area (12%), followed by buffer zones (36%) and transition areas (53%) ( Table 3). Across BRs, from core areas to transition areas, crop production and grazing SPAs made up an increasing proportion of the area of zones and, in contrast, timber production and recreation SPA proportions decreased from core to transition areas (Table 3).

Effect of biosphere reserve zoning on ecosystem service coproduction
Differences between the three zones were confirmed in 17 of 32 pairwise Wilcoxon rank sum tests (Figure 3 and Figure S2 in Supplementary Material). Differences were significant for the comparison of buffer zones vs. transition areas for 11 indicators and for core areas vs. buffer zones for 6 of the 16 indicators.
Overall, we found some support for expected differences in anthropogenic contributions to the four ES ( Figure 3). Indicators of ecosystem management for all four ES had higher values in transition areas than in buffer zones ( Figure 3). Ecosystem management indicators were also higher in buffer zones Table 2. Spatially explicit indicators for ecosystem service coproduction: anthropogenic contributions (ecosystem management (EM), harvest and access (HA)), natural contribution (NC) and provision (P) of four ecosystem services.
x compared to core areas for crop production and recreation, but not for grazing and timber production. Furthermore, harvest and access indicators for grazing and timber production were higher in buffer zones than in core areas, and higher in transition areas than in buffer zones. Recreation infrastructure density was higher in buffer zones than in core areas, but we did not detect a difference between buffer zones and transition areas ( Figure 3). We found weaker support for expected differences in natural contributions to the four ES ( Figure 3). For recreation, as expected, natural contributions were higher in core areas than in buffer zones, and higher in buffer zones than in transition areas. We also found that pest control, the natural contribution for crop production, was higher in buffer zones than in transition areas. For natural contributions to grazing and timber production, we could not detect differences between zones at a 95% confidence level. However, soil biomass productivity, the natural contribution to grazing, tended to be higher in buffer zones than in core areas (p = 0.06).
Our analyses revealed few significant differences between zones for ES provision indicators. We did not find provision differences in these indicators between core areas and buffer zones for any of the four ES, but we did find some differences between buffer zones and transition areas ( Figure 3). As expected, crop production and grazing provision were higher in transition  Table S1 in Supplementary Material for data sources for each biosphere reserve. Table 3. Mean area of 137 biosphere reserve (BR) zones, mean zone proportions in relation to total biosphere reserve area, and the proportion of service providing area (SPA) per zone. Service providing areas encompass all land cover classes in which the respective ecosystem service can be provided. Standard deviation values are provided in parentheses.

Zone
Mean area (ha) Mean proportions of zones to total BR area (%) areas than in the buffer zones, but there was no difference for timber production or recreation between these zones. In contrast to our expectation, recreation provision tended to be higher in core areas than in buffer zones (p = 0.07, Supplementary Material Table S4).

Explanatory variables for ecosystem service provision in biosphere reserves
Natural contributions were significant explanatory variables of ES provision for all four ES (Table 4) in our statistical models. Natural contributions were positively associated with the provision of crop production, grazing, and timber production, but negatively associated with the provision of recreation. However, for recreation, the negative association of natural contributions was smaller in BRs that met the recommended areal proportions of zones (interaction effect). We also found that, for some ES, anthropogenic contributions were significantly associated with ES provision. Ecosystem management showed a significant positive association with the provision of timber production, and the harvest and access indicator showed a significant negative relationship with the provision of crop production. We found no significant influence of zoning on any of the four ES provision values in the model. We did observe that other variables were significant predictors of ES provision, however. Terrain ruggedness was negatively associated with grazing and timber production, and positively associated with recreation provision. For timber production, the combined core area and buffer zone model indicated a positive relationship with ecosystem management (human-induced land alteration), but this relationship weakened with increasing terrain ruggedness (Supplementary Material Table  S5). Unsurprisingly, we found that a higher proportion of pastures or forests leads to higher provision of grazing or timber, respectively. When a BR met the recommended areal proportions of zones, this had a positive association with the provision of crop production and grazing and a negative association with timber production and recreation. For grazing, the model indicated a significant positive association with harvest and access for BRs with the recommended areal proportions of zones.
The explanatory power of the models was highest for grazing (R 2 range 0.49-0.50) and lowest for recreation (R 2 range 0.07-0.09) ( Table 4).

The role of biosphere reserve zoning in ecosystem service coproduction
Most of the significant differences we observed in ES coproduction and provision indicators occurred when we compared buffer zones with transition areas. While the tested differences matched our expectations in more than half of the comparisons (17 of 32), our results also showed that, in some cases, coproduction indicators are not different between zones. This was mainly for comparisons between core area and buffer zone and for indicators of natural contributions and ES provision. From the pronounced differences in ecosystem management as well as harvest and access between zones, we can conclude that zoning is effective for human interventions. Transition areas consistently showed higher values of management pressure than core areas and buffer zones for indicators such as human-induced alteration of net primary productivity (HANPP luc ) or road density in forests. Similarly, an evaluation of the efficacy of zoning in conservation areas in China showed that human disturbances occurred less often inside core zones compared to the surrounding, less protected zones that correspond to BR transition areas (Xu et al. 2016). The relatively limited impact of 'core' vs 'buffer' designation may be due to these two zones receiving similar management. In theory, different zones have different management restrictions, but in practice, they might not, or the restrictions might not be effectively implemented (Hull et al. 2011). The few significant differences for indicators of natural contributions may be explained with weak effects of management restrictions on natural contributions compared to other environmental variables. For example, the natural contribution indicator of timber production and grazing depends among soil on climate and topography, which are not significantly altered by conservation management. The Table 4. Outputs from model averaging for models with Δ AIC < 2 from the best model for total ES provision values (. p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001). Data show the number of averaged models, the range of R-squared values of the averaged models, the number of valid observations and coefficients from the n best models and the standard error (SE) of averaged coefficients. No output exists for predictors labeled with NA, either because they were removed from the model because of a squared GVIF^(1/(2*Df)) >10 or because of model selection. effectiveness of conservation areas also depends on a variety of socioeconomic factors (Geldmann et al. 2019), so that differences in ES provision may not be explained solely by different management practices. For the three provisioning services, we generally observed lower anthropogenic contributions, in terms of ecosystem management or harvest and access, in core areas than in buffer zones. For ES provision, no differences for any of the four ES between core areas and buffer zones were found. This may indicate either that these inner zones are affected by similar management restrictions or that other conditions influence ES provision. Alternatively, it could imply that less intensive management practices are equally capable of producing ES, especially those ES strongly influenced by other factors than management intensities. Our results confirm that transition areas produce higher levels of the three provisioning ES. For example, we found higher crop provision in transition areas than in buffer zones. Transition areas are usually not protected and thus we expected fewer management restrictions for crop production in these zones compared to buffer zones. A similar pattern was found by Eigenbrod et al. (2010) in England, where agricultural production was much lower in strictly protected areas than in the rest of the country.

Main and interaction effects
For the cultural ES, recreation, a common assumption is that moderate land use intensity is attractive and that a certain degree of accessibility is needed for recreational use, as would be expected in the buffer zone of BRs (Paracchini et al. 2014). While we do find that buffer zones contain more recreation infrastructure than core areas, there was no difference in this indicator between buffer zones and transition areas. One explanation could be that we included roads/ tracks that are mostly for agricultural or forestry uses when creating the indicator of recreational infrastructure. While this infrastructure may not be directly intended for recreation, we assumed it is also used for recreation (Turgut et al. 2021). Because the density of tracks is higher in areas with agricultural or forest use, this might have increased the density of recreation infrastructure in buffer zones and transition areas, compared to core areas, and influenced our results. In addition, our data do not allow us to determine the direction of the causal relationship between recreation and infrastructure. The popularity of a site could trigger the construction of recreational infrastructure, so that the infrastructure itself is an important manifestation of recreation . However, for recreation provision, as recorded by the number of geotagged photos on social media, we did not find differences between buffer zones and transition areas, but did observe somewhat higher numbers in core areas than in buffer zones. This pattern that core areas (strictly protected areas) are attractive for recreation was also shown by Eigenbrod et al. (2010).
Previous studies have also identified a positive influence of natural characteristics (number of biotopes and biodiversity) on visitation rates in national parks (Neuvonen et al. 2010;Siikamäki et al. 2015). This underlines the importance of naturalness: high levels of human alteration could lead to a decrease in recreation provision, as cultural services have their optimum in natural and semi-natural areas (Braat and ten Brink 2008). Higher recreation provision in core areas contradicts the design concept of BRs, where the principle is that core areas should focus on protection and not on active use. This may be as core areas in BRs in Europe remain relatively accessible overall. Recreational use, however, is supposed to have a lower impact on protected area ecosystems than provisioning services, such as crop or timber production or grazing (Maxwell et al. 2016).

Explaining ecosystem service provision in biosphere reserves
Our regression analysis indicates that natural contributions were strong significant predictors of ES provision of all four ES. However, we did not find an effect of zoning per se, as a proxy for management restrictions, on the provision of any of the four ES. This is surprising as for two ES we observed higher provisioning services in the transition area than in the core area or buffer zone, when we tested for statistical differences in our first analysis. However, in the alternative regression model, where we combined the core area and buffer zone into a single zone, and compared this to the transition area (Supplementary Materials Table S5), there was an effect of zoning for crop production.
When BRs met the recommended areal proportions of zones, this had a strong significant effect on provision of all four ES. In particular, crop production and grazing were positively influenced when zone proportions met the recommended criteria. As the provision of these two ES was highest in the transition area, we assume that this is mainly due to sufficiently large transition areas (the criteria is of at least 50% of the total area) with large patches of cultivated area. The higher efficiency of large cultivation or timber production areas for provisioning ES may also explain the positive association between the proportion of SPAs and both grazing and timber production.

Crop production
Our model reveals that natural contributions together with harvest and access are significant coproduction factors for crop provision. The positive effect of natural contributions, such as pest regulation, has been acknowledged in other studies (Bommarco et al. 2013;Rega et al. 2018). Surprisingly, we found a negative influence of the embodied energy inputs, such as fertilization or irrigation, on crop production. While this contradicts general findings from some studies (e.g. Lobell 2007;Levers et al. 2016;Schröter et al. 2021), Levers et al. (2016) also found a negative influence of nitrogen input on specific crop types, which highlights the varied relationships between fertilizer application and crop production. Furthermore, human inputs to crop production vary spatially with different levels of intensification across Europe (Levers et al. 2016;Vallecillo et al. 2019), which is possibly not captured in our analyses.

Grazing
Natural contributions (soil biomass productivity of grasslands and pastures) were also a strong positive predictor for grazing provision in our model. While soil, climatic and topographic conditions are input data to the natural contributions indicator (Tóth et al. 2013), primary productivity is used as input data to the provision indicator of grazing (Plutzar et al. 2016). As the association of soil, climatic and topographic conditions with primary productivity is widely acknowledged (e.g. Svoray et al. 2008;Schut et al. 2015) the strong predictive capacity of natural contributions for grazing provision is not surprising. Soil suitability for grassland production has also been shown to positively influence the distribution of dairy cattle in linear regression models (Neumann et al. 2010). The negative influence of terrain ruggedness is in line with findings on the negative impact of slope on grazing for cattle, although sheep and goats can graze in areas with steep slopes (Mueggler 1965;Neumann et al. 2010).

Timber production
For timber production, the model confirmed the importance of natural contributions (soil biomass productivity of forest areas), which is the strongest variable among all standardised predictors in the model. This finding is in line with research on the importance of soil functions for forest productivity (Burger and Kelting 1999;Lang et al. 2016). The ecosystem management indicator of managementrelated pressures was also a strong positive predictor. This finding is in line with other studies showing that the proportion of non-native or planted conifers and vertical heterogeneity, both expressions of management, are important predictors of timber production (Felipe-Lucia et al. 2018), and that the proportion of plantation species strongly predicts harvest intensity (Levers et al. 2014). The negative association of terrain ruggedness has also been found for forests in Europe (Levers et al. 2014). We also found that timber yield increases as the proportion of the area of a BR zone covered with forest (SPA) increases, which may imply that large forest areas could be more productive than smaller ones.

Recreation
We found a tendency for recreational infrastructure to be positively associated with recreation provision (indicated by geotagged photos). This finding is comparable with studies in national parks that have shown positive effects on visitation rates of trails and roads (Schägner et al. 2016) and of offered recreation activities and trails (Neuvonen et al. 2010). There is also a tendency for core areas to have more recreation provision than buffer zones as the pairwise tests for differences revealed. The model, however, showed a negative association between the absence of disturbance and recreation provision. In other words, more geotagged photos were taken in BR zones with high disturbance levels and therefore low natural contributions. This is interesting, as it shows that areas with higher levels of light at night, such as transition areas, potentially provide more recreational ES provision. This could be because of the higher number of people living in these areas, and thus increased demand for recreation, pointing to the importance of accessibility for recreation in addition to naturalness (Paracchini et al. 2014). We also found a positive effect of terrain ruggedness on recreation provision, presumably because people like scenic hilly areas, which is in line with other findings (Weyland and Laterra 2014;Schirpke et al. 2018). However, as the overall predictive capacity of the recreation model is very limited, these results should be interpreted cautiously.

Limitations of the study and challenges for future work
Our study represents, to the best of our knowledge, the first effort to quantify and analyse ES coproduction with large-scale spatial datasets of ES coproduction and provision indicators. However, we had a limited ability to measure coproduction of the selected ES across broad spatial scales with finegrained data. The lack of spatial data on land use intensity is a recognised problem (Kuemmerle et al. 2013), and also limited the number of ES that could be included in this study. Coproduction of ES is particularly relevant for regulating ES, such as water and flood regulation, freshwater supply and pollination, due to considerable anthropogenic contributions (Palomo et al. 2016). Unfortunately, there is a lack of large-scale, spatially explicit data for these services across Europe. Additionally, some ES may have occurred outside the service providing units (SPA) that we necessarily assigned them to. The land cover classes used in our SPA definitions constitute classes where the respective ES are generally, but not exclusively, provided. Especially for recreation, this delimitation is likely to be imperfect. There was also a temporal mismatch between some of the ES indicator data and spatial zoning data of newer BRs. In other words, older ES indicator datasets might not reflect spatial patterns of BRs that were designated recently: Of the 137 analysed BRs, 18 were designated after 2015, which is the latest reference year of the indicator data, with the exception of the OSM data. We encourage future studies to take into account the potential effect of time on human management actions and ES coproduction. The age of the current zoning scheme of a BR might affect differences in ES coproduction indicators. However, the extent of different zones of BRs can be revised and adapted after a BR is designated (Price et al. 2010), and data on these zoning changes are not consistently available.
While we believe our selected indicators meet the criteria of credibility (valid representation, embedded within a conceptual framework, agreement by experts) and feasibility (data availability) (van Oudenhoven et al. 2018), we recognize that these indicators do not represent all aspects of ES coproduction or provision. For example, the level of light at night has limitations as an indicator of the level of disturbance, as it does not reflect other disturbances such as pollution, noise, or historical land clearing. Another example is that the modelled livestock density data, which we used as an indicator for harvest and access, assumes that strictly protected areas (IUCN categories Ia, Ib, II or III) are unsuitable for grazing, and hence designates a livestock density of 0 for these areas. Fortunately, only a negligible proportion of the area of the studied BRs overlaps with such strictly protected areas. Thus, we have some confidence that we do not substantially underestimate the actual levels of grazing activity in BRs. This is further supported by our finding that the mean of the median values for grazing SPA in core areas was different from 0 (mean: 21.7, interquartile range: 4.7-27.2). A similar caveat exists for the data on managementrelated pressures on forest ecosystems, as the degree of protection (whether or not it is in a Natura 2000 site) is an input into this proxy (Hengeveld et al. 2012). The proportion of the area of BR zones covered by Natura 2000 sites differs between the zones: more than 90% of the core areas but less than 20% of the transition areas (data not shown).
In addition, terrain ruggedness is included indirectly in the ES provision indicator of grazing (Plutzar et al. 2016), as slope is used to derive the suitability of a location for grazing. As terrain ruggedness and slope are highly correlated (Amatulli et al. 2018), this might lead to circular reasoning in the regression models. The negative effect of terrain ruggedness on the provision of grazing, as obtained in our model, should thus be considered with caution.
With regard to the indicator of recreational infrastructure (harvest and access), a spatially heterogeneous mapping effort of the underlying OSM data results in varying data completeness between countries (cf. Barrington-Leigh and Millard-Ball 2017). We assume that this has a negligible impact on our results, as we work with pairwise comparisons within each BR for which the variance in data completeness is small. Moreover, the EU has a relatively complete coverage of roads in the OSM data (Barrington-Leigh and Millard-Ball 2017).
All of the above mentioned data limitations could partially explain why the chosen spatial indicators do not reflect differing management restrictions of zones and hence did not differ significantly between BR zones in nearly half of the comparisons (15 of 32) in our analysis. Thus, to better assess and understand ES coproduction across spatial scales, more efforts to provide spatially explicit and continuous data on ES are needed (Vaz et al. 2021).
Finally, other physical or socioeconomic variables not included in our analyses could be as important to ES coproduction and provision as zoning. Future studies could compare BR data on ES coproduction and provision with systematically selected control areas in order to assess effectiveness of zoning (dos Santos Ribas et al. 2020).

Implications for operationalisation of the biosphere reserve concept in policy and management
To date, comprehensive knowledge of the zoning and spatial configurations of BRs has been limited. From our analysis of the average proportions of zones, we can confirm they are consistent with the general design concept of BRs, but also that only 41% (57/ 137) fulfil recommendations for areal proportions of zones as recommended by German and Austrian authorities. We also find an expected increase of the land use categories cropland and pasture from core areas to transition areas, and a higher share of forests (and other more natural land cover classes) in the two inner zones. These findings are in line with the BR design concept, according to which protection is the focus of core areas and buffer zones, and sustainable use is the focus in transition areas.
While every BR should have a management policy or plan (according to the Statutory Framework of the World Network of BRs, UNESCO 1996), UNESCO (2021) provides only limited guidance on management, recognising the different and varied situations of BRs in different countries. Thus, a challenge in finding general patterns of ES coproduction across BRs is that there are no common guidelines on how much human modification or land use intensity is allowed in the different zones. While such standardisation might be desirable from a scientific point of view, a diversity of natural and cultural landscapes and land use practices is also part of the BR 'philosophy' (Ferreira et al. 2018). Additionally, the overlap of different national and EU conservation areas (e.g. nature reserves, landscape protection areas, Natura 2000 areas) presumably varies across the BRs. These varying overlaps might affect management restrictions in different BR zones. A detailed analysis of spatial overlaps of the three zones with protected areas could be part of a future study.
Zoning seems particularly effective for anthropogenic contributions, i.e. direct measures of human modification, which are generally lower for the provisioning services in the inner zones. As these anthropogenic contributions can have negative aspects on other ES and biodiversity (Kleijn et al. 2009;Storkey et al. 2012;Newbold et al. 2015;Beckmann et al. 2019), zoning can be seen as a useful conservation tool. However, whether the existence of a transition area with minor protection benefits ES and biodiversity in the core area remains an open question (Xu et al. 2016). When evaluating the operationalisation of BR zoning, we also need to be careful about attributing the observed differences in natural and anthropogenic contributions and ES provision to the establishment of the zones, per se. In particular, core areas have to be already protected under national legislation, i.e. they may have been established because of pre-existing conditions (e.g. high level of natural contributions). Hence, our analysis did not allow us to assess effectiveness in terms of isolating the causal impacts of BR zones on ES outcomes (see, e.g. Sims and Alix-Garcia 2017 for an impact assessment of BRs). Teasing out causal impacts requires time-series data (before and after BR establishment and zoning), as well as data on control sites (e.g. similar sites where BRs were not designated) which was beyond the scope of this study. As current reporting and data on effectiveness of BR management are limited (Ferreira et al. 2018), this implies a need for regular monitoring and updating efforts of spatial layers across the EU in order to measure the effectiveness of BRs for ES provision relative to other sites. We see a need for standardised ES assessments across zones and BRs (Vasseur and Siron 2019). For instance, the periodic review process of BRs (Price et al. 2010) could include ES and information on management restrictions across zones in a standardised way, in line with the Lima Action Plan for 2016 to 2025, which states that BRs are to be 'recognized as sources and stewards of ecosystem services' (UNESCO 2016, p. 6). According to this plan, BRs should identify important ES and report on their quality and quantity in the periodic reviews. Such information could enable more detailed analyses on ES (coproduction) in the future.

Conclusion
Across 137 terrestrial biosphere reserves in the EU, we found that the inner zones focusing on conservation (core areas and buffer zones) show lower levels of anthropogenic contributions to the four selected ecosystem services than zones focused on use (transition areas). This could imply that zoning is effective in limiting management intensity in the inner zones, in line with the general design concept of zoning in the BRs. While natural contributions and ES provision are also expected to differ between zones through different forms of management, we did not find clear patterns in our indicators. Overall, the majority of expected differences were confirmed, but we also found mixed results with respect to particular aspects of coproduction for single ES. This could be due to data limitations described above. Or it could point to a lack of standardised management restrictions across BRs. However, different management might be appropriate given the diversity of contexts, landscapes, and management practices within BRs. Model selection and averaging revealed that natural contributions positively related to ES provision for all four services. This is in line with other findings on the importance of natural contributions for ES provision and points to the essential need of conserving ecosystem functioning and biodiversity. Surprisingly, zoning per se was not related to the provision of any of the four ES in the regression analysis.
Our study provides the first insights into the relationship of BR zoning and ES coproduction based on a unique and newly compiled large-scale European dataset of terrestrial BRs. Future studies could look into comparisons of BRs with other areas and perform temporal analyses, preferably based on data collected on a BR level in line with the Lima Action Plan of the Man and Biosphere Programme of UNESCO. This study is the first step towards a greater understanding of the role of natural and anthropogenic inputs into the provision of ES in a large set of BRs across Europe. While there is significant scope for future work, our results indicate that BRs can provide a model for balancing natural and anthropogenic contributions to ES, and a model for sustainable use of natural resources more generally. granted to MS. AB and MFL gratefully acknowledge the support of iDiv funded by the German Research Foundation (DFG-FZT 118, 202548816). We thank Mick Wu for statistical support. We thank Berta Martín-López for participation in the workshop and helpful comments on an earlier draft of this paper. We thank UNESCO-MAB and UNEP-WCMC for their useful advice on zoning data collection. We are grateful to all national agencies, organisations and biosphere reserve administrations, who helped through providing zoning data, in particular, the Spanish Organismo Autónomo Parques Nacionales and the German Federal Agency for Nature Conservation (BfN). We thank Carlo Rega (JRC) for providing data.

Disclosure statement
No potential conflict of interest was reported by the author(s).