Introduction

Forests are under pressure globally, largely attributed to human-induced climate change with rising temperatures, shifting precipitation patterns, and an increase in extreme heat events and droughts1,2,3,4,5,6. Human activities threaten forests through fragmentation and degradation, thereby endangering crucial and unique habitats7,8,9,10. Yet humans depend on forest ecosystems for carbon storage, temperature regulation, soil fertility, resource use, and recreational space11,12. Protected areas represent a key measure to protect and conserve essential ecosystem services and habitats13,14,15. They play a major role in preserving various dimensions of biodiversity across different scales16, mitigating climate change13, and safeguarding ecosystem functions17. One of the objectives outlined in the recently adopted Kunming-Montreal Global Biodiversity Framework 2022 is to expand the coverage of protected areas (Target 3, 18). With the establishment of more protected areas, it is imperative to evaluate the effectiveness of existing protected areas and their associated outcomes. However, this assessment should not rely solely on measuring protection coverage18,19. The effectiveness of protected areas has been investigated and defined in various ways, for example referring to the achievement of specific conservation targets within these areas20 or equating it to management effectiveness21. Here, we define protected area effectiveness in terms of their contribution to sustaining ecosystem health and functioning.

The International Union of Conservation of Nature (IUCN) established a set of categories as a comprehensive framework to designate protected areas with a primary focus on biodiversity preservation22. UNESCO’s biosphere reserves, established under the Man and the Biosphere program, differ from strict nature reserves in that they serve as integrated areas for both nature and human activities, acting as model areas for sustainable development23,24,25. This globally applied framework spans 748 reserves, typically organized into three zones with varying levels of management24,26. The core, buffer, and transition zones range from stricter protection to sustainable use of resources and the integration of education and research. 174 of these reserves have known forest cover. While these protected, forest-covered areas are widely diverse, they are all united within the world network of biosphere reserves. The effectiveness of biosphere reserves in creating learning places for sustainable development needs to be assessed following a transdisciplinary research agenda. Although these are increasingly utilized27, this analysis is only the beginning of understanding biosphere reserves’ complex structures and objectives. We explicitly address the functionality of forest ecosystems, and not biodiversity nor the services biosphere reserves provide for humans. If biosphere reserves in forest biomes are ecologically effective, their contribution to biodiversity conservation should be measurable. In fact, biodiversity plays into the functioning of ecosystems and is, therefore, indirectly integrated into the outcomes of this study28. Ecosystems inside these model areas for sustainable development would then be in a better condition than those outside. To date, the global analysis of the ecological effectiveness of protected areas has depended on comparable measures based on globally available border data29. Regional studies based on remote sensing data investigated species richness30 or functional diversity31 in a single biosphere reserve. To our knowledge, this is the first time that the ecological effectiveness of forested UNESCO biosphere reserves has been tested globally using remotely sensed data.

Satellite remote sensing imagery stands as the sole source of information presently available for understanding and investigating ecosystem functioning on a global scale, offering relatively high temporal and spatial resolutions32. We build upon the definition of ecosystem functions from a remote sensing perspective as the direct or indirect dynamics of ecosystem processes, such as primary productivity and evapotranspiration32, and the ability of those ecosystems to adapt to internal and external changes33. Employing multiple proxies for primary productivity helps identify different patterns of ecosystem functioning across a wide range of forested biomes (Fig. 1)34.

Fig. 1: Graphical interpretation of the forest ecosystem function proxies used in this study.
figure 1

Variations in vegetation vitality assessed via gross and net primary productivity as well as via the enhanced and normalized difference vegetation index; Variations in temperature regulation capacities via hot day land surface temperatures. Variations in precipitation and water regulation via evapotranspiration and variations in ecosystem extents and state via forest patch connectivity.

Full size image

In this study, we analyzed multiple proxies of ecosystem functions in all forest biomes found in biosphere reserves as pertinent proxies of forest ecosystem function32,33,35. The chosen proxies are derived from a review of important, reliable, and currently available remotely sensed indicators of ecosystem functioning and were further filtered to all proxies relevant to forest ecosystems32. We used the GLAD tree cover product, established as a dataset extensively used for its global coverage, high resolution at 30 m, and high accuracy, based on multiple years for pixel-wise value validation35. We derived the annual sums of the preprocessed MODIS Aqua gross and net primary productivity, annual maxima of MODIS Aqua vegetation indices (the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI), and annual sums of the MODIS Terra net evapotranspiration. We calculated the temperature regulation capacity based on hot days in the MODIS Aqua Land Surface Temperature product. We created forest patch connectivity products using Thiessen polygons on the Global Forest Cover Change dataset. We expect these analyses to provide information on biosphere reserves’ effectiveness, although we are aware that the conclusions derived from such an approach can be limited in the absence of a counterfactual analysis, which was not possible in this instance. To assess the effectiveness of biosphere reserves in supporting forest ecosystem functions worldwide, we compared their ecosystem functioning with their surrounding areas based on the seven proxies at a 1 km resolution (Fig. S1). Furthermore, we filtered the data collection and investigated changes over time in two-time steps from 2010 to 2016 and 2017 to 2022.

Results and discussion

Effectiveness of biosphere reserves

According to the available border data, at least 242 (32%) of the global biosphere reserves are located in forested biomes, and 174 have known forest cover. Among these, 119 biosphere reserves had accurate border information and sufficient tree cover to assess forest ecosystem functions. We defined tree cover based on the GLAD tree cover product as a pixel with >30% at 30 m resolution and >50% at 1 km resolution. Biosphere reserves were deemed to contain sufficient tree cover for analysis when the land cover share of tree-covered areas exceeded 20% inside the biosphere reserve (Fig. S1, Tab. S1). We filtered for biosphere reserves, established before 2010 to find potential changes after designation. We extracted remotely sensed ecosystem functioning proxies within biosphere reserves and in their surroundings from 2010 to 2022, leading to nearly 8 million single observations/pixels at 1 km resolution. We defined the surroundings via a buffer using the natural logarithm of the size of each area multiplied by 500 to account for the varying size of the investigated biosphere reserves (Fig. S6). The buffer was then bound to a rectangular format for faster processing (Fig. S1). We modeled forest ecosystem functions inside biosphere reserves against their surroundings and the effects of place (biome) and time on seven proxies (Fig. 2, Tab. S2). Our approach to defining the effectiveness of biosphere reserves in maintaining multiple forest ecosystem functions (E) can be read as follows:

$${E}_{{br}}=\left({P}_{{in}1} > {P}_{{out}1}\right)\wedge \cdots \wedge \left({P}_{{in}7} > {P}_{{out}7}\right)$$
Fig. 2: Model predictions and standard deviation for selected forest ecosystem function proxies inside and outside the 119 biosphere reserves.
figure 2

Models are depicted against the respective observed median for two times 2010–2016 (dots) and 2017–2022 (triangles). Predictions above 0 indicate higher proxy values inside than outside. Predictions below 0 indicate the opposite. Colors refer to biomes. Modeled inside/outside difference of A gross primary productivity in kg*C/m2, B forest patch connectivity, C evapotranspiration in kg/m2/year, D hot day land surface temperature in °C. Note that y-axes for C evapotranspiration and D hot day land surface temperature are reversed to visualize the effectiveness in water captivity and cooling on top like for A gross primary productivity and B forest patch connectivity.

Full size image

P1–7 are the seven proxies of forest ecosystem functions modeled with their corresponding location per biome for two-time steps and their corresponding pixel-wise forest cover share. Pin and Pout define the inside and surroundings of each biosphere reserve respectively. An effective biosphere reserve (Ebr) is defined as one where all proxies have higher modeled forest ecosystem functions inside as compared to their surroundings. The classification as an effective or ineffective biosphere reserve depends on whether only small differences (e.g., 0.1 “better” inside than outside) or large differences have been modeled. To assess the effect of the selected biosphere reserves on forest ecosystem functions, we modeled each proxy against the “Inside/Outside” forest cover share, biome, and the two-time steps, with the biosphere reserves as random factors using the linear mixed effect model with the highest AIC. The explained variance for these seven proxy models was for fixed effects R2 = 0.11–0.38 and for random effects R2 = 0.91–0.99 (Tab. S3). To specifically assess “Inside/Outside” patterns, we modeled each proxy’s difference between the inside and outside against the median forest cover share, biome, and the two-time steps, with the biosphere reserves as random factors. The proxy models with the inside-outside comparison as a fixed effect rather than as an explanatory variable showed for fixed effects R2 = 0.21–0.82 and for random effects R2 = 0.5–0.91 (Table S4, Fig. 3). Of 119 biosphere reserves, only 18 (15%) show higher ecosystem functions inside than outside when considering all functions (Table S2, Fig. 2).

Fig. 3: Forest ecosystem functions in biosphere reserves.
figure 3

A World map forest biomes and UNESCO biosphere reserves. Forest biomes include boreal forests/taiga, mangroves, temperate broadleaf and mixed forests, temperate conifer forests, tropical & subtropical coniferous forests, tropical and subtropical dry broadleaf forests, tropical and subtropical moist broadleaf forests (Ecoregions 2017, Resolve); biosphere reserves: blue = 297 available border data; of those: yellow = 17 established after 2010; red = 119 forest covered and included in this study. B World map of ecosystem function effects in 119 biosphere reserves. The biosphere reserves served as random effects for the three modeled proxies: hot-day temperature, connectivity, and gross primary productivity with fixed effects of forest cover share, inside or outside the biosphere reserve, biome, and two periods of time 2010–2016 and 2017–2022. Each proxy’s random effects were scaled, multiplied by 255, and transformed into RGB values. The colors follow the triangle legend and show the corresponding effect sizes of the three proxies. Smaller windows show the same effects for Central America and Europe.

Full size image

Biome-specific patterns were observable for all proxies, especially for gross primary productivity and evapotranspiration (Fig. 2A, C). When cutting the output of modeled effectiveness inside and outside for each proxy in 10 groups, the average of biosphere reserves close to zero for all proxies is 79 biosphere reserves (Fig. S9). Consequently, 79 biosphere reserves do not show substantial differences between inside and outside on average (Fig. 2). Given that only 66% of the areas exhibit minimal differences from their surroundings, it remains to be investigated if those reserves have a positive effect on their surroundings or no effect at all. Modeled primary productivity and connectivity were higher inside than outside for 65 and 93 biosphere reserves, respectively (Fig. 2A, B). Hot day temperatures and evapotranspiration were predicted to lower inside than outside for 98 and 76 biosphere reserves, respectively (Fig. 2C, D). Substantial temporal changes (from dots to triangles) in proxy medians were rare except for forest patch connectivity, which decreased substantially for single biosphere reserves (Fig. 2B). Biome-wise temporal changes for connectivity occurred inside the biosphere reserves in the conifer and boreal temperate forests with −0.02 and 0.03 reduced modeled connectivity respectively (Fig. S4). This shows that conditions were maintained over time for single proxies and single biosphere reserves: for example, the Fontainebleau Biosphere Reserve, close to Paris, has a rather low level of connectivity to other forested patches in the landscape, but connectivity was still found to be higher inside the reserve than outside. The same applies to two coastal biosphere reserves, Sian Kaan in Mexico and Georgian Bay in Canada, which have higher connectivity inside than outside. The Xiriualtique Jiquitizco in El Salvador is a coastal biosphere reserve as well: while for this biosphere reserve, gross primary productivity is slightly higher inside, all other proxies (evapotranspiration, hot-day temperature, and forest patch connectivity) show higher values outside the biosphere reserve (Fig. S3). One explanation for the observation that individual biosphere reserves do not perform better inside than outside in terms of forest functionality in all proxy results may be related to regional ecological and climatic characteristics and differences in management. We anticipated and confirmed higher connectivity of forest patches within 85 biosphere reserves than outside them (Table S5). Modeled primary productivity (gross and net primary productivity, the enhanced and normalized difference vegetation index) was higher inside than outside for most biosphere reserves (with n = 65, 66, 58, and 91 of 119, respectively) and increased slightly from 2010–2016 to 2017–2022 for all but the forests in the tropical dry broadleaved and the boreal biome. This general increase might be attributed to the general trend of global greening36. We detected the highest primary productivity and the highest difference to the surroundings in the Selva el Ocote Biosphere Reserve in Mexico, a tropical, moist, broadleaved forest area (Table S5). The modeled effect of lower evapotranspiration within biosphere reserves was more pronounced in tropical and subtropical broadleaved forests than in temperate broadleaved forests (Fig. 2C).

Forest ecosystem function proxies in and around biosphere reserves generally decreased over time but showed minor biome-specific differences: predicted primary productivity increased slightly in tropical broadleaf and temperate coniferous forests and decreased slightly in tropical coniferous forests. In all biomes, hot day temperatures increased, and connectivity decreased. Changes were more pronounced outside than inside biosphere reserves, which can be interpreted as an indication of effectiveness (Table S4). Apart from the assumed higher functionality of temperature and water regulation within biosphere reserves, a general increase in temperatures through global warming was confirmed. The year of establishment is no indicator of a better-functioning biosphere reserve. We plotted the modeled effect size of the 119 biosphere reserves for gross primary productivity, connectivity, and hot-day temperature against the designation year of each reserve (Fig. S8). Other than a positive trend of gross primary productivity for recently established reserves, again potentially attributed to the general trend of global greening36, no pattern was identified.

Global distribution of effective biosphere reserves varies across biomes

As of now, a complete dataset of the world network of biosphere reserves remains unavailable. Among the 474 datasets we were able to collect, 119 contain sufficient forest cover to assess ecosystem functions (Fig. 3A, Methods).

We selected three proxies to display the modeled effect size of each biosphere reserve (Fig. 3B, Tables S4 and S6). In temperate biosphere reserves, forest patch connectivity is a higher influencing factor than gross primary productivity or temperature: high and low connectivity differences with their surroundings are more pronounced than in other forest biomes (Fig. 3B Europe, washed green and brown dots). The fragmentation of temperate forests is caused by intensive forestry, transport routes, and secondary causes from rising temperatures, such as pests, and continues to weaken the forest ecosystem functionality 35. Tropical and subtropical biomes showcase higher random effects for hot day temperatures than for primary productivity or forest patch connectivity (Fig. 3 Central America and Central Africa, violet dots). This effect can be interpreted as a decreasing resilience due to rising temperatures and decreasing water availability in tropical forests6.

Across biomes, we predicted for only 18 of 119 biosphere reserves a higher primary productivity and connectivity, lower temperatures on hot days, and lower evapotranspiration compared to their surroundings in 2017–2022 (Table S5). When reducing the number of proxies from seven to only three (Gross primary productivity, temperature, and connectivity), the number of effective biosphere reserves increases to 42 (35%) (Fig. 3B). Based on these findings, we advocate for multiproxy assessments, always considering the most parsimonious dataset and model for robust interpretation and derivation of conclusions.

Some smaller reserves displayed more pronounced random effects towards temperature (Fig. 4), suggesting that larger reserves generally exhibit greater temperature regulation capacities, likely attributable to their higher forest cover shares37. Similarly, the random effects of forest patch connectivity are higher for most of the smaller biosphere reserves. Connectivity in larger biosphere reserves might be more balanced.

Fig. 4: Modeled random effects of 119 biosphere reserves for three forest ecosystem function proxies.
figure 4

Model with fixed effects of forest cover share, inside or outside the biosphere reserve, biome, and two periods of time 2010–2016 and 2017–2022. Each random effect was scaled, multiplied by 255, and transformed into RGB values for three ecosystem function proxies, namely connectivity, hot-day temperature, and gross primary productivity. Each circle represents one biosphere reserve, and the circle size relates to the area size of each reserve. Ternary plot created with the Ternary-package in R60.

Full size image

Higher tree cover constitutes healthy forests in biosphere reserves worldwide

Higher tree cover shares exert a significant impact on ecosystem functions, as shown by all forest-modeled ecosystem function proxies in our study (Fig. 5).

Fig. 5: Modeled ecosystem functioning proxies for forest cover shares in 119 biosphere reserves.
figure 5

A Gross primary productivity, B net primary productivity, C enhanced vegetation index, D normalized difference vegetation index, E hot day land surface temperature, F evapotranspiration, and G forest patch connectivity. Modeled observations are displayed in blue for forests inside biosphere reserves and in red for surrounding forests. Ribbons represent 0.95 confidence intervals. A model with fixed effects of forest cover share, inside or outside the biosphere reserve, biome, and two periods of time 2010–2016 and 2017–2022.

Full size image

We defined tree cover shares as the percentage of tree cover per pixel based on the GLAD tree cover data (Fig. S1, Table S1). A linear relationship between forest cover shares and modeled proxies emerged: as forest cover shares rose, so did the primary productivity, while temperatures declined, evapotranspiration decreased, and forest patch connectivity increased. This trend holds true across all forest biomes (Fig. S7). Forest cover share is part of the linear mixed effect model and influences all proxies significantly, as does, to a lower extent, the difference between the inside and outside (Table S4). The tropics exhibited a more pronounced increase in primary productivity with higher forest cover shares compared to temperate latitudes. Modeled gross primary productivity and the Enhanced Vegetation Index (EVI) were higher outside at the highest forest cover shares, which is likely related to ecosystem functionality on a landscape level. The EVI’s standard deviations were more pronounced than those of other proxies, perhaps partly explained by its sensitivity to greenness in tropical forests. In tropical and subtropical biomes, biosphere reserves exhibited higher primary productivity compared to their surroundings, whereas the reverse holds for temperate forests. In boreal forests, temperatures were nearly the same inside and outside. Biosphere reserves in the boreal biome indicated the lowest primary productivity values and decreased over time. This could be related to the observed boreal forest biome shift with less productivity assessed38. One anomaly was the three biosphere reserves with tropical and subtropical dry broadleaf forests for which higher evapotranspiration outside than inside was predicted. Modeled forest patch connectivity increased with increasing forest cover shares throughout the biomes, except for tropical and subtropical dry broadleaf forests and mangroves (Fig. S7). Generally, higher connectivity did not lead to higher primary productivity. Larger forest blocks are typically relevant for hosting larger species populations and are also essential for reducing edge effects and providing microclimatic and ecohydrological buffering (e.g., 36). This does not imply that small forests do not have very high conservation values but rather indicates that both small and large forest areas need to be targeted by conservation strategies39.

Measuring the effectiveness of biosphere reserves and other area-based conservation instruments

We highlight the impact of proxy selection on the results, emphasizing the need for multiproxy analyses to obtain a comprehensive understanding of forest ecosystem functioning. Remote sensing offers comparable, freely available global datasets, and the spatial and temporal resolution of these datasets will increase even more in the future32. We appreciate the easy access and large computational capacities of cloud-computing platforms such as the Google Earth Engine. Our research offers insights relevant to optimizing sustainable forest management strategies in biosphere reserves and potentially other forested protected areas worldwide.

While there is general agreement that the effectiveness of conservation and sustainable development requires a proxy-based assessment, the indicators proposed to date for assessing forest conservation are relatively simple. Rather than data on management action and outputs such as forest areas under formal protection or certified management (e.g., CBD’s Global Biodiversity Framework, Target 10, 18), which only indicates the intention to improve forest management, we urgently need outcome indicators. These must allow us to assess, with high temporal and spatial resolution, whether management promises are being kept and translated into ecosystem health. In this context, we need to make better use of available remote sensing data, which can be used to analyze not only tree cover but also the health status and trends of forest ecosystems. The Ad Hoc Technical Expert Group (AHTEG) on Indicators for the Kunming-Montreal Global Biodiversity Framework40 will need to ensure that the indicators to be used do not lag behind rapidly evolving technological capabilities. This includes the integration of scenarios of ecosystem functionality developments in the absence of any form of designated protection. We show that a global approach using observable proxies is effective in providing insights into regional characteristics. However, our study did not take country-specific governance or management regimes into account. Future studies will show if and how social and societal aspects influence the effectiveness of biosphere reserves. Additionally, improved monitoring at finer spatial scales is recommended for local management.

In any case, we suggest that a multi-proxy methodology focusing on ecosystem performance, as implemented in this study for UNESCO biosphere reserves, could be applied to other protected area frameworks. This approach could also be adapted to other ecosystems and land cover types or by adding a level of detail and distinguishing between forest types. A basic requirement is reliable and up-to-date georeferenced data on area-based interventions. In this context, we regret that there is no official source for the boundaries of all UNESCO biosphere reserves, which limits their evaluability.

In some reviews, biosphere reserves are excluded from the assessment of effectiveness41,42,43, concerning doubts as to whether they actually represent effective conservation action or merely bureaucratic labels44. Here we would like to emphasize the role that biosphere reserves could play as modern instruments for safeguarding the functioning of the biosphere - following a convivial conservation approach45, under a type of governance that includes people and their needs. However, UNESCO biosphere reserves cannot just aim at being model areas for participatory management without benefiting ecosystems. Particularly in the context of the rapidly advancing climate crisis, human well-being will increasingly depend on the health of ecosystems. Any conservation action that does not have a measurable positive impact on ecosystem functioning, including our proposed proxies and beyond, will lose its credibility and justification.

Methods

Summary

The main datasets based on satellite imagery used were the GLAD forest cover 2010 product46 and seven different preprocessed products of NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS) with derived forest ecosystem functioning proxies, all accessed and processed through the cloud-computing infrastructure Google Earth Engine and the connected large repository of geospatial datasets47. The primary vector dataset is a collection of border information from the world network of biosphere reserves. All datasets are described in Table S1. For a spatio-temporal analysis at 1 km, we collected a 20 percent random sample for each of the seven forest ecosystem functioning proxies inside and outside for two time spans 2010–2016 and 2017–2022.

Forest definition and area selection

The global GLAD tree cover dataset 2010 at 30 m resolution (https://glad.umd.edu/dataset/global-2010-tree-cover-30-m) served as the basis for our definition of forest. We filtered tree cover at 30 m with a threshold of ≥ 30 % and transformed the data to a binary tree cover/no tree cover format. We used the binary tree cover at 30 m to define tree cover at 1 km with a threshold of ≥50%48. This tree cover dataset 2010 at 1 km with values from 50 to 100 % coverage was used in the analyses to define and select forested biosphere reserves (example: Fig. S1).

The set of UNESCO biosphere reserves border data was acquired from the World Database on Protected Areas49, a European Collection of biosphere reserves’ border information50, and directly from local biosphere reserves administrations.

We excluded biosphere reserves that (1) were available as point information only, (2) were established after 2010 since we wanted to detect changes over time under designated areas, (3) had not only less than 20% forest cover inside but also less than 10% forest cover in their surroundings (example: Fig. S1), (4) were outside of forest biomes, namely tropical and subtropical moist and dry broadleaf and coniferous forests, temperate broadleaf and mixed and conifer forests, boreal forests/taiga, and mangroves. The selection criteria are comparable to standard practice (e.g., 14).

Out of 748 biosphere reserves, we obtained geospatial border datasets for 474. Among these, 308 were established before 2010. Of those 174 had correct border information and lie within forest cover biomes, with 119 having sufficient forest cover. Due to the wide range of area sizes of all biosphere reserves, the outside, which is the surrounding area, was defined with a buffer calculated with the natural logarithm of the size of each area multiplied by 500 (Fig. S6). The buffer was then bound to a rectangular format for faster processing. This way, each biosphere reserve is surrounded by a rectangle proportional to the size of the reserve where pixels are defined as outside. We did not exclude designated areas, such as other protected areas or biosphere reserves, from the outside area since we wanted to depict the real-world characteristics of both the inside and outside biosphere reserves. Nevertheless, we accounted for this fact as a source of uncertainty and extracted the percentage of protected area cover in the surroundings of the biosphere reserves as a moderating variable.

Proxies of forest ecosystem functions

The selected proxies follow the overview of relevant, convincing, and potentially available remotely sensed ecosystem functioning indicators suitable for forest ecosystems32. The proxies are forest connectivity, evapotranspiration, the Normalized Difference and Enhanced Vegetation Index, gross and net primary productivity, and land surface temperature. Forest connectivity for each biosphere reserve and its surroundings for both time steps was produced by (1) creating tree cover 2016 (excluding tree cover loss before 2016) and tree cover 2022 (excluding tree cover loss before 2022) datasets, (2) computing Thiessen polygons around forest fragment centroids, (3) computing the ratio of each Thiessen polygon size and the corresponding forest cover size, which results in values ranging from 0 to 1, (4) sampling 10% inside and 10% outside, and (5) transforming to table data. Except for forest connectivity, all proxies were processed as raster data for each biosphere reserve and its surroundings by (1) if required, being scaled to 1 km, (2) masked to tree cover 2010 at 1 km, (3) temporally filtered to means per pixel for the time steps 2010–2016 and 2017–2022, (4) defined/masked as inside or outside, (5) sampled 10% inside and 10% outside, and (6) transformed to table data. The sampling was conducted using the stratified random sampling function in the Google Earth Engine, with the same randomization seed for all reserves. The number of fetched pixels depends on the size of the biosphere reserve and corresponds to 20% in total for the reserve and its surroundings. The per-pixel formula for both time steps 2010–2016 and 2017–2022 for each proxy was defined as follows: Forest connectivity 2016 and 2022, mean of annual sums of evapotranspiration in kg/m²/year, mean of yearly maxima of the Normalized Difference and Enhanced Vegetation Index, mean of the annual sums of the gross and net primary productivity in kg*C/km2, mean of days ≥30 °C of land surface temperature in °C (as a proxy for the cooling function of forests37), mean of the annual sums of precipitation in mm/year. For a detailed product and processing overview, see Table S1.

Area characteristics

For each biosphere reserve, elevation and slope from the global multi-resolution terrain elevation data, GMTED 201051 were masked from water areas and scaled to 1 km. We did not use the commonly used SRTM data, since we required data above 60 °N and the overall accuracy of GMTED data is good52. The water bodies distribution was based on the MODIS Terra land water mask at 250 m and the water mask 201053 was selected and scaled to 1 km.

The land use land cover change data from the MODIS land cover type yearly data at 500 m54 was processed by selecting land cover 2010, 2016, and 2021, comparing 2010 with 2016 and 2016 with 2021 pixel-wise and upscaling to 1 km (if >50 % land use land cover change in 1 km = 1). The burned area data from the MODIS Burned Area Monthly Global 500 m55 indicating the day of a year a fire happened was processed by selecting the time series for each time step (2010–2016 and 2017–2022), transforming values to 1 if burned and 0 if not burned and scaling to 1 km. The metadata for each biosphere reserve, based on the described datasets above, contain the year of establishment, biome and ecoregion information56, tree cover percentage for 2010, 2016, and 2022, water cover percentage for 2010, burned area percentage for 2010–2016 and 2017–2022 and protected area cover percentage 202349. For a detailed product and processing overview, see the Supplementary Material (Table S1).

Statistics and linear mixed effect models

The dataset generated in the Google Earth Engine encompasses 7,938,752 values. They correspond to 119 biosphere reserves and their forest cover inside and outside for seven ecosystem functioning proxies and two-time steps. We used the R programming language for all statistics with RStudio 2023.06.1 + 52457. To investigate if ecosystem functions in biosphere reserves differ from forests in their surroundings over time and worldwide, we used a linear mixed effect model framework, specifically the lmer—function, which is part of the lme4—package58, integrating the nested structure of our data into those models. All forest ecosystem functioning proxy data is continuous. We test multiple models for each proxy to identify the best fit. The biosphere reserves served as random effects. We fitted as fixed effects forest cover percentage, biome, time steps, and the inside and outside to explore which combination best describes the behavior of each proxy. The Aikake Information Criterion (AIC) was used to select the best model for our dataset. The general structure of all models using the maximum likelihood where:

  1. (1).

    proxy ~ Inside_Outside, random = ~ 1|Biosphere_Reserve

  2. (2).

    proxy ~ Inside_Outside *Forest_Cover_2010, random = ~ 1|Biosphere_Reserve

  3. (3).

    proxy ~ Inside_Outside * Forest_Cover_2010*biome, random = ~ 1| Biosphere_Reserve

  4. (4).

    proxy ~ Inside_Outside * Forest_Cover_2010*biome*time, random = ~ 1| Biosphere_Reserve

  5. (5).

    proxy ~ Inside_Outside * Forest_Cover_2010*biome+time* Inside_Outside, random = ~ 1| Biosphere_Reserve

    To model the difference between the inside of the biosphere reserves and their surroundings, we build on the same model structure and used the median proxy difference of inside and outside:

    1. (1)

      proxyDifference_Inside_Outside ~ Forest_Cover_2010_median*biome*time, random = ~ 1| Biosphere_Reserve

Uncertainty and limitations

In the course of this analysis, we investigated various sources of uncertainty. Since the number of biosphere reserves per biome is not balanced, we modeled the data with and without biomes with few numbers (n = 3). Since no significant differences were observed, the biomes with small numbers of biosphere reserves remained in the dataset. We modeled and mapped the standard deviation of the biosphere reserves as random factors (Fig. S4). We compared the standard deviation of each model with the size of the biosphere reserves (Fig. S5). Other data on the number of observations per pixel per proxy per time step, population density, overlapping protected areas, water coverage, elevation, and slope are available upon request. Other sources of uncertainty that were not investigated and that may influence the results were globally different governance and management regimes and forest types. In this study, we did not include a counterfactual to evaluate the effectiveness of biosphere reserves. As a result, our findings should be interpreted with caution. Without a counterfactual, it is difficult to disentangle the effects of protection from those attributable to location. While we expect our results to capture certain aspects of effectiveness, we cannot determine the relative importance of each factor. To elaborate on this, we suggest using structural equation models. Future studies could create added value by investigating the interactions between the proxies and other social or societal factors. Our analysis is based on satellite imagery at 1 km resolution with the known influencing factors of cloud cover, atmospheric conditions, and the study-related selection of spectral, temporal, and spatial resolution. Therefore, we cannot conclude on a fine scale or address rapidly changing environmental dynamics. Validating remotely sensed observations with ground-truthing is especially relevant to derive conclusions at the local and regional levels. Local assessments covering socio-ecological aspects, especially when investigating effects in biosphere reserves, are recommendable (e.g., 59).

Reporting summary

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