We selected all Mexican MAB-UNESCO Biosphere reserves (hereafter reserves) located in the Mesoamerican hotspot, which were established before the year 2000 (Fig. 1, Table S1, Appendix 1). We exclude from our analysis the biosphere reserves La Encrucijada and Pantanos de Centla because of the low performance of forest/no forest identification (see below), therefore we include 19 reserves. The reserves show a gradient of forest loss, with the remaining forest cover within the reserves ranging from 36% to 98%. We combined the Lacan-Tun/Montes Azules and Los Petenes/Ria Celestún reserves into two reserve complexes called Lacandona and LPRC, respectively, because they shared borders, same biome, and management.

Fig. 1.

Location of the studied biosphere reserves in Mexico. Colors show the dominating vegetation type in each reserve. All the Mexican territory below the Tropic of Cancer not included in any protected area was considered in our study as unprotected zones.

(0.43MB).

To characterize forest changes inside reserves and unprotected areas (areas in Mexican tropic non included in any PA, see Fig. 1), we performed a supervised classification of Landsat images for the years 2000 and 2020. We used three major classes: forest, no-forest, and water. We evaluated the accuracy of this classification using a set of independent class cover validation points and determined that the mean (±SD) value of the overall satellite image classification accuracy was 94.4% ± 0.02 for the 2000 year and 94.4% ± 0.04 in 2020 (Table S4, Appendix 2). We conducted a cover change analysis to identify the areas with forest loss, forest regrowth, and no change for a period of 20 years, from 2000 to 2020 in the whole studied area.

We also calculated forest loss rate, forest regrowth rate, and forest fragmentation rate within reserves. To calculate forest loss rate (FLR) and forest regrowth rate (FRR) we use the next formula:

Protected areas are often located in isolated places with steeper slopes, higher elevations, and poorer soils than unprotected places (Joppa and Pfaff, 2009). Thus, a simple comparison of forest changes between protected and unprotected regions may overestimate the effect of legal protection (dos Santos Ribas et al., 2020). Matching analysis is a quasi-experimental approach that statistically allows comparing sampling units with similar covariates and, therefore, to test the effect of treatment without the effect of covariates (Ho et al., 2011). This technic matches sampling units of treatment and control groups with similar covariates values to test the effect of a treatment. Matching technics have been used previously to assess the effect of PAs on forest loss (Yang et al., 2021), but its use regarding forest regrowth and forest fragmentation in PAs is less frequent (but see Sims, 2014).

We used two different sampling strategies for the matching analysis: the first with one million points uniformly seeded in the entire studied area (both inside and outside reserves) at 1 km of separation, and the second with 1500 circular microlandscapes with 3 km of radius randomly seeded in the same area (see Fig. S1 and S2, Appendix 1). Each sampling point has information regarding whether they are located inside reserves or in unprotected areas, and about the changes in forest cover that occurred during 2000–2020 in its specific location (i.e. forest loss, forest gain, or no change). We used the microlandscapes to document forest fragmentation inside reserves and in unprotected areas using formula 2. Furthermore, each sampling unit had the value of the three covariates that may act as confusing factors on forest cover changes: distance to cities, distance to roads, and agriculture suitability (see Appendix 2). To properly compare forest cover changes in reserves and unprotected areas we only include those samples located within forest cover in the year 2000 for the forest loss data (n = 423,900), and within no forest cover in 2000 for the forest regrowth data (n = 331,623).

We estimated propensity scores for forest loss, forest regrowth, and fragmentation data using binomial error distribution with a logit function where the response variable was the treatment (protected or unprotected). After matching, all standardized mean differences for the covariates were below 0.1 which indicates a good balance between control and treatment samples (see Appendix 3). The matching analysis was performed using the MatchIt package (Ho et al., 2011) of R (R Core Team, 2021).

As underlying drivers of these forest changes, we considered two economic indicators: (i) the distance to major cities (localities with a population larger than 15,000 people) as a proxy of access to markets, and (ii) the non-farm occupation (i.e. the mean proportion of people, of the municipalities surrounding the biosphere reserves, working in the professional, services, and industrial sectors; Table S2, Appendix 1) in the year 2000. Population growth rate (during 2000–2020), population density in the year 2000, and rural settlement density in the year 2000 were also used as indicators of demographic factors. Using a correlation and a collinearity analysis, we proved that these indicators were not correlated (see Fig. S3, Appendix 1). During the first stage of our study, we considered including other indicators, such as unemployment rate, marginalization, human development, and agricultural subsidies. However, they were not included in our analysis because of the high correlation with the selected indicators (r > 0.5, p < 0.05, Variance Inflation Factor (VIF) > 2). Most selected indicators of underlying drivers are common in other studies (Aide et al., 2013; Figueroa et al., 2009). The theoretical relationship with forest loss and details of the calculation of these indicators can be found in Table S3 of Appendix 1 and Appendix 2, respectively.

We tested for the effect of reserves (i.e., protection) on forest loss, forest regrowth, and forest fragmentation. For the first two variables, we used generalized mixed models with binomial error distribution for each one. We included forest change as binary response variables (forest loss/no change and forest regrowth/no change) and the protection condition (inside reserves or unprotected) as a fixed effect factor. To test for the effect of the protection condition on forest fragmentation rate and the number of patches in the year 2020, we used a generalized linear mixed model with a Gaussian error and a Poisson error, respectively. We included the biome type as a random effect. We performed the models with pre-matching and matching data. Matched data models were weighted using the weights automatically calculated by the matching algorithm that accounts for the potential differences in the number of samples in each treatment.

To determine the effect of the indicators of socioeconomic drivers on the forest loss, fragmentation, and forest regrowth rates, we used a multimodel inference approach with linear models (Burnham and Anderson, 2002). Because of the small sample size (n = 17), we limited to 3 the maximum number of terms in the models to avoid overfitting. Therefore, for each response variable, we constructed all models that represent all possible model combinations with 1, 2, and 3 predictors and their pair interaction (i.e., approximately 600 models). To calculate standardized parameter estimates (i.e., slopes), we applied a z-score transformation for all predictor variables. We calculated the Akaike Information Criterion corrected for small samples (AICc), and the model with the lowest AICc value was considered the best model. The percentage of deviance explained by the best model was calculated with square-R (goodness-of-fit). The effect size of each variable was estimated with the model-averaged standardized parameter estimates.

Our findings confirm the importance of reserves as key tools to prevent forest loss and to preserve landscapes with a low fragmentation degree. Nevertheless, we did not find evidence supporting the idea that reserves also reduce the forest fragmentation rate. In other words, biodiversity within reserves cope with a lower fragmentation degree than in unprotected zones, but fragmentation increases at similar rates within reserves and unprotected zones. Since forest loss has higher deleterious effects on biodiversity than fragmentation (Fahrig, 2017), our findings represent less concern news for biodiversity inside reserves, without neglecting the potential negative effects of fragmentation (Ewers and Didham, 2006).

Surprisingly, we did not find support to the hypothesis that reserves promote forest regrowth. This contrasts with other works that suggest a positive effect of protected areas on forest recovery (Borda-Niño et al., 2020). This apparent contradiction may be clarified by our matching analysis. In our study, a simple comparison of protected and unprotected areas, without controlling for environmental covariables (pre-matching condition), supports the expected higher forest regrowth within reserves. However, controlling for environmental covariables obliterated such differences. Therefore, the forest regrowth inside and outside reserves does not seem to depend on the establishment of the reserves per se, but on the distance to cities and agriculture suitability (Fig. S4, Appendix 1). Thus, forest regrowth occurs in areas where agricultural activities are not favored because of the biophysical environment (i.e. steeper slopes, high elevations, and poor soils) or because of economic causes (a large distance to markets promote land abandoning and forest regrowth).

Our findings indicate that the forest loss rate inside the reserves is mainly related to two main underlying drivers, positively to the density of rural settlements and negatively non-farm occupation. These results suggest that the lack of job opportunities in rural areas, beyond those offered by the agricultural sector, accentuates deforestation as the demand for agricultural land increases in the absence of other livelihood options. As reported by previous studies, reduction of forest loss rates is associated with the increase of non-farm labor activities (Hoang et al., 2014; Wunder, 2003), since these allow local people to obtain revenues without the need of causing more deforestation for arable land (Curtis et al., 2018; Kovacic and Viteri Salazar, 2017; Vedeld et al., 2007). In parallel with these studies, we found that reserves embedded in municipalities with greater labor opportunities in non-farm sectors, especially in the industrial activities (e.g., manufactures, electricity, or construction, see Table S4, Appendix 1), endure lower forest loss rates (Goers and Lawson, 2009). Also, a higher human development index was associated with a greater non-farm occupation (Fig. 4). This finding concurs with studies pointing out that non-farm activities allow local people to obtain higher profits and deal with poverty (Haggblade et al., 2007). According to our study, the emerging non-farm occupation may promote a “win-win” scenario with a positive balance between conservation and development, especially considering that an incipient forest transition process in southern Mexico has been identified (Vaca et al., 2012). However, it is important to note that this “win-win” scenario is likely limited to regions dominated by subsistence agriculture, such as the study reserves in Mexico (Bonilla-Moheno and Aide, 2020; Curtis et al., 2018). Additional studies are required in places where large-scale agriculture oriented to the production of commodities occurs, as the diversification of labor activities could have a relatively weak effect on deforestation rates in contexts where agriculture is typically supported by large capitals.

Fig. 4.

Relationship between the non-farm occupation, forest loss rate (a), and human development index (HDI, (b). Each point represents a single reserve, lines represent the generalized linear model and the grey shadow is the 95% confidence interval.

(0.2MB).

As expected, the population growth rate was a strong predictor of forest cover changes, being positively related to forest fragmentation rate, and negatively associated with forest regrowth. Human population size is considered a major driver of forest loss (Aide et al., 2013) and fragmentation (Li et al., 2010). Human pressures such as logging and agriculture are expected to increase with increasing the population size (Ehrlich and Holdren, 1971), which in turn can promote the loss and fragmentation of the remaining forest, and limit the recovery of the forest in degraded lands. Our findings contrast with those of Borda-Niño et al. (2020), which suggest that the forest regrowth rate can increase with increasing (not decreasing) demographic factors such as the density of rural settlements. This counterintuitive finding can be related to the environmental policies after the establishment of the reserves. In Mexico, for example, the creation of reserves often requires the relocation of rural communities that were within the limits established for the reserve, which leads to farmers abandoning their cropland allowing the forest to regrow (Fig. S5, Appendix 1). We infer that a higher forest regrowth occurred especially in reserves that had a higher density of rural settlements in their surroundings. Indeed, previously to the year 2000, when some reserves have not been legally established, the amount of forest loss rate inside reserves was significantly and positively related to the density of rural settlements (R2 = 0.4, P < 0.01, Fig. S6, Appendix 1). Although difficult to support with empiric data, another possibility is that there exist active restoration efforts, as occur in some regions (Juan-Baeza et al., 2015), potentially associated positively with demographic factors (the larger the population, the greater the possibility of restoration).

Finally, we have some caveats about our work to enunciate. Our work assessed the influence of socioeconomic variables measured at the municipality on forest spatial changes, however, there exists the possibility that a finer scale (e.g. locality) outcome more reliable results, nonetheless, this is still a challenge since some variables are not available at fine resolution. In addition, there exist other underlying and proximate drivers (e.g. environmental factors) non-assessed in our work that influence forest spatial changes in protected areas that may need specific hypothesis tests in future studies. On the other hand, since here we assessed forest spatial changes in 20 years, our work may not identify properly early stages of forest regrowth (mostly given by small trees and shrubs) and therefore, underestimates forest regrowth and forest fragmentation.