Tropical forest landscapes encompass various combinations of agricultural modification, changing forest cover with the rapid integration with the global commodity market (Deakin et al., 2016). A landscape may exhibit changing land-use practices through a gradient/continuum of agricultural modification and decreasing tree cover. The purposively selected landscapes in Bangladesh, Burkina Faso, Cameroon, Ethiopia, Indonesia, Nicaragua, and Zambia exhibit varying rates and scales of agricultural and forest cover changes (see details for each country in Figure 1; Table 1; Deakin et al., 2016). We applied a standardized set of criteria—such as the level of agricultural pressure (from subsistence to commercial farming), population density, access to markets and roads, and proximity to forest—to identify suitable landscapes in each country (Sunderland et al., 2017). Initially, scoping studies were conducted to ensure that the criteria effectively defined three distinct intensification zones (Zone 1, Zone 2, and Zone 3) representing a gradient or continuum of land use modification in each landscape.

Each landscape demonstrated a transition through the gradient, ranging from areas with the best available forest cover and subsistence agriculture (low intensification/Zone 1) to agroforestry systems (moderate intensification/Zone 2) and finally to monoculture or intensive cash crop systems (high intensification/Zone 3). Zone 1 represents landscapes where subsistence farming is predominant, and people depend heavily on forests. Zone 2 represents areas with intermediate or mixed farming systems, with lower forest access but relatively better access to roads and markets. Zone 3 represents highly intensive, modified landscapes dominated by monoculture agricultural systems, with limited forest access but higher access to markets and roads.

Forest and land cover trends were summarized for each zone using repeat imagery from the Landsat series. To do so, images of surface reflectance from three time periods (1985–1990, 1998–2004, and 2010–2015) were acquired for each site. Images were selected to coincide with the dry season to minimize cloud cover and aid in distinguishing agriculture from other land covers (Table 2). Clouds, shadows, and water bodies were removed using both manual digitization and the provided quality assurance bands. Vegetation indices, including the Normalized Burn Ratio, Tasseled Cap, and Disturbance Index, were calculated for each site, and a combination of vegetation index image thresholding and Maximum Likelihood-based classifications were used to classify land cover for each site. Sufficiently cloud-free imagery was not available for the Indonesian study site for the most recent period, thus 2000 and 2012 normalized top-of-atmosphere Landsat cloud-free composites were used. Land cover classes varied among sites initially but were ultimately aggregated into three classes: forest, non-forest, and no data (consisting of clouds, water bodies, and cloud shadows) (Forest and non-forest cover maps for each landscape in Supplementary Appendix 1). At a 30-m resolution, it can be challenging to distinguish smallholder agroforestry from the surrounding forests. Thus, for practical reasons, as well as the particular context of each site, we implemented slightly different definitions of “forest” to help account for the diversity of forest structural types and plantations seen across tropical regions (Table 1). Classification accuracy was assessed with a combination of field verification, iterative local knowledge pieces and high-resolution imagery (e.g., RapidEye, Google Earth). Accuracy assessments with local people’s knowledge about forest cover represent a best-case scenario for countries where local imagery and ground truthing were not feasible.

Our classified imagery was further used to characterize forest and land cover within zones and within a 2 km radius of each household to better understand the local context of forest resources. To do so, we used a suite of scale-independent forest change and fragmentation metrics, including % stable forest cover, % stable non-forest cover, % forest cover at time of third Landsat image, % forest loss, % forest gain, and % fluctuating forest over roughly 20–30 years (based on the three classified Landsat images for each site (Table 2; Figure 2). Stable forest was defined as pixels that remained forest throughout the three time periods, while similarly stable non-forest remained non-forest throughout the three time periods. Forest loss indicates a forest pixel that was deforested in the second or third period, and remained so. Fluctuating forest pixels switched land cover during time two but then reversed to original cover by time three. Stable non-forest and forest cover at third image were highly correlated with stable forest. As they did not produce as strong models, only stable forest was retained for statistical analysis. To explore the impacts of forest fragmentation, we also examined forest edge density using SDMtools. Edge density is quantified as the total length of forest edge divided by the total forested area, and in our case, we quantified this within a 2-km buffer of each household. While we also quantified additional landscape metrics (e.g., largest patch index, aggregation index), they were highly correlated. Thus, to reduce the impact of multi-collinearity on regression models, only edge density was used to characterize forest fragmentation in our analysis.

The preliminary estimate of forest dynamics showed a mean percentage of stable forest to be just under 34%, ranging from 0 to 86% (Table 3). Mean percentages of fluctuating forest and forest gain were under 10%, though forest gain was twice as variable as fluctuating forest (std. dev. 6.68% vs. 3.30%, respectively). On average, the percent of forest loss was 23%, yet ranged from < 1% to just < 60%. Forest dynamics differed by country and zone over the past several decades (Supplementary Appendix 2). The zones of land use intensification were distinguished by notable contrasts in percent coverage of stable forest (in six countries) as well as forest loss (in three countries). Stable forest was greatest in Zone 1 in six of the seven study countries, with the exception of Burkina Faso where it was nearly equivalent in Zone 1 (19.76%) and Zone 2 (20.84%). Forest loss was the most prominent form of change in forest cover, where loss surpassed both gains and fluctuating forest in four of seven countries (Burkina Faso, Cameroon, Indonesia, Nicaragua in all zones; Zambia in Zone 2; and Ethiopia in Zone 1). Forest loss was relatively consistent across zones in Bangladesh (ranging from 7.24% minimum to 9.24% maximum), and was greatest in Zone 1 in Ethiopia (7.0% in Zone 1, and ∼0% in Zones 2 and 3). Fluctuating forest cover was minimal across all study countries, but was greatest in Zone 2 for both Burkina Faso and Zambia. Similarly, few study sites displayed much forest gain, and those that did varied both in amount of gains and according to zone (Bangladesh witnessed a maximum of 19.7% in Zone 2 versus 12.2% in Zone 1; Burkina Faso hovered between 5.5 and 6.5% across all three zones; and Zambia had the greatest forest gain in Zone 1 of 21.9% and the least in Zone 2 of 2.3%). Both percent total forest cover at the third time period as well as stable non-forest were strongly correlated with stable forest cover, and were thus not included in the pursuant models.

The study applied a landscape-level approach through a three-stage nested hierarchical experimental design (i.e., land use zones-villages-households). Household surveys were undertaken in village/s or settlement/s (∼100 households in 1 village or combining several villages) within the dominant land use “zones”- representing a gradient/continuum of agricultural modification and decreasing tree cover (Details in Sunderland et al., 2017). Each landscape constitutes approximately 300 households surveyed (approximately 100 randomly selected households in each of the three zones, i.e., low intensification/Zone 1, moderate intensification/Zone 2 and high intensification/Zone 3. The sample households were selected using random numbers of the lists of households available within local census information or election records and arial photography. Across all countries, we surveyed 1,900 households (∼300 households from each landscape) between 2015 and 2017. The household response rates were 90% and we obtained 1,823 households with detailed information useful for this study. We obtained GPS coordinates for household locations to characterize forest conditions within 2 km radii around each household (described further below). The research collaborator in each of the seven countries recruited locally experienced field assistants for data collection. The household head or the respondents aged above 18 years old were considered for the survey.

The surveys covered information on household assets, the frequency of household respondent’s visits to different forest types and NTFP collection. Firstly, we gathered information on domestic, productive and non-productive and transport assets in each household to categorize them into wealth groups. We asked the respondents whether their family possessed (yes or no) any assets from a standard list provided to them. Domestic assets included roof materials, building walls, and toilets. Productive assets were different farming tools such as axes, machetes, chainsaws used, and the presence of livestock; transport assets included transportation modes, i.e., boats, bicycles and motorbikes.

Secondly, we asked the household respondents: Which forest types do you use and how often do you visit? By asking this question, we recorded the frequency of visits and the type of forests (in terms of ownership) visited at the household level. Participants selected from a range of forest ownership types under four broad categories: family/private land, common property (non-protected forests/small forest patches/riverine forest), reserve/protected forest, and community forest/managed forest (plantation forest/commercial forests). Participants also reported how frequently they visited different forest types from: daily, few times a week, once every 2 weeks, once a month, very occasionally, and never. We mainly targeted routine forest visits as visits that took place daily or weekly in the rest of the paper. Finally, in relation to forest product collection, participants were asked: What different forest products have been collected in the past 12 months? We provided a list of forest products including NTFPs relevant to the local contexts. Due to the rarity of any individual forest products gathered, we categorized them into forest foods, forest fuels, and forest construction materials for subsequent analysis. A fourth category called “other” was also created with miscellaneous items (bee products, dye or tanning, flowers or plants, gaharu—relevant for Indonesia, gums or resins, medicinal, reeds and papyrus, water, wildlife/pet trade) none of which were used in subsequent analyses.

We examined whether the forest visitations and use of forest products varied by forest cover types, agricultural intensification (zones) and wealth levels (Table 4). We conducted generalized linear mixed models (GLMM) to explain household forest visitation patterns and the collection of forest products across the forest types, intensification zones and wealth conditions with the countries as a random intercept. We used intensification zones (Zones 1–3), forest cover metrics (gain, loss, stable, fluctuating), forest configuration (edge density) and wealth conditions as predictors and the countries as a random effect. We developed a composite wealth index for each household comprising the asset variables to use as a proxy indicator of relative wealth conditions (relatively poor to relatively rich). Different domestic, productive and non-productive, and transport assets reported at individual household was used in the construction of wealth indices. Depending on country about 30–43 asset variables were found showing over 5–95% variations within the households that we considered for a Principal Components Analysis (PCA). The PC1 shows the maximum variation which was used to build a wealth index by dividing it into five quintiles representing different wealth groups. We categorized wealth index into five quintiles or wealth groups to apply wealth level in regression analysis to determine their effects on forest products use. Wealth quintiles calculated from the sample households represent the relative measure of how wealth is distributed within the population in three zones in each country (Supplementary Appendix 3).

In the first research question regarding forest visitation, we used a GLMM to determine which predictor variables best explained the likelihood of households making routine visits to forests and particular types of forests by fitting separate models for family, common, protected, and industrial forests. To address our second question regarding household forest product collection, we used a GLMM to examine a collection of three forest product categories (forest foods, forest fuels, and forest construction materials). We used logit models with a binary response for routine visits (visit/no visits) and forest product collection (product used/not used). In order to isolate the most significant factors affecting forest visitation and forest product collection, reduced models were fit via stepwise elimination based on Akaike’s information criteria (AICc) and log likelihood values using the buildmer package (Voeten, 2019) in R software. All numeric independent variables were standardized using Z-scores. In doing so, the difference in AICc and weight of AICc were calculated for each model. The reduced model was then compared to the full model using a Chi-square test (Supplementary Appendices 4, 5). We included a country random effect (random intercept) to account for unobserved cross-country differences while focusing inference on zone (between country variation, i.e., Interclass Correlation Coefficient [ICC] 0.270; Variance 1.215; Marginal R² = 0.145, Conditional R² = 0.376). All analyses were carried out in the software R-Core Package 3.4.2.

Almost two-thirds of the respondents surveyed (63% of households) reported visiting forests regularly across all seven landscapes, with relatively higher visitation in Zone 1 and Zone 3 (24 and 22%, respectively) compared to Zone 2 (18%) (Figure 3). Model of routine visits show a significant effect of zone, stable forest, fluctuating forest and forest gain on household’s regular visits (Figure 4). At zone level, household routine visits differed within the three intensification zones (P < 0.001). A significantly greater number of people visited in the low (Zone 1) and high intensification (Zone 3) zones compared to the moderately intensive land use Zone 2 (P < 0.001) (Figure 4). No significant difference in forest visitation was found between high and low intensification zones (P = 0.075). Increases in stable forests, fluctuations, and forest gains were positively associated with more households regularly visiting forests across the countries (Figure 4).

Further analysis of visitation by forest ownership showed higher regular visits to private/family-owned forests (30%) compared to common property (19%), reserve (16%), and community-managed forests (14%) among all households surveyed (Figure 3). Zone 1 recorded the highest household visits to forests across all ownership types except family-owned forests. No significant difference was observed in routine visits to family/private forests among intensification zones. Zone explained variation of routine visits in common forests (P < 0.001), managed forests (P < 0.001), and protected forests (P < 0.001), but not private forests (P = 0.344) (Figures 4, 5). However, only forest gains were positively associated with increased visits by private owners to their forests across countries (P < 0.001). With increases in stable forest areas, visits to reserved/protected areas are more likely, in contrast to a decline in common property forests. Additionally, only edge density showed a positive effect on community forest visitation (Figure 4).

Forest products commonly used by households include food, fuel wood and construction materials. The highest usage was recorded for fuel wood (83%), followed by forest-sourced food (71%) and construction materials (78%) (Figure 3). The use of construction materials varied significantly across intensification zones, with the highest usage in Zone 1 compared to Zone 2 (P < 0.05) and Zone 3 (P < 0.001). Households in areas with more intensive land use, i.e., Zone 3, were less likely to use forest-derived construction materials. No significant difference was found for forest sourced use of food and fuel wood by zones. Model estimates indicated a positive association between increased stable forest cover and fluctuating forest conditions with greater household use of forest-sourced food (Figure 6). However, the use of forest-based construction materials was less likely to increase in areas with forest loss or fluctuating forest cover across all countries. Overall, the analysis revealed that more frequent visits to forest areas are associated with higher household use of forest-sourced food and construction materials. Additionally, an increase in forest edge density was found to correlate with a higher number of households using construction materials.

Rather than deeply investigate the annual nuances and fluctuations of forest change, our study examined a few key points in time that corresponded well to the household surveys and oral histories of participants in the study region (who were intensively surveyed as part of the overall project, beyond what is discussed in this paper). While this approach does not create a classic quantitative accuracy assessment, it does reflect very realistic ground-based knowledge of the forest and non-forested conditions at each site which were examined in detail and revised in an iterative process. The studied landscapes include sparse dry forests and savanna ecosystems, which are often challenging to accurately classify land cover types using images with moderate spatial (30 m) and temporal (10 years) resolution (Dube et al., 2021; Gara et al., 2016; Alencar et al., 2020). Some of the sparser regions (e.g., Burkina Faso, Zambia) are less amenable to mapping with Landsat-derived land cover products due to their moderate resolution which we know (from our work on the ground) misses some key forest patches. The greater proportion of fluctuating landscape classified in Burkina Faso and Zambia study sites, both largely covered by savanna woodlands with relatively low forest cover (Zoungrana et al., 2018), may therefore be, at least in part, reflective of the challenges Landsat images present in classifying spatially diverse landscapes.