Introduction

The change in land cover and its impacts on biodiversity are present all over the globe (Foley et al. 2005; Newbold et al. 2015). Land cover changes are often due to deforestation and the decline of natural forest habitats (Foley et al. 2005; Curtis et al. 2018). Main drivers behind these processes are forestry, urbanization, commercial agriculture and shifting cultivation under subsistence farming (Seto et al. 2012; Curtis et al. 2018). Although the probability of deforestation and its extent is highly variable across landscapes (Mertens et al. 2004; Busch and Ferretti-Gallon 2017), some consistent predictors have been identified. Income competition with alternative land uses (e.g., agriculture), proximity to roads, populous areas (e.g., urban centers) and already converted areas accelerate deforestation, whereas inaccessibility of forests (e.g., slope, elevation, status of protection) and law enforcement reduce deforestation rates (Busch and Ferretti-Gallon 2017).

Land use patterns in many low-income countries of Sub-Saharan Africa are characterized by shifting cultivation to sustain rural livelihoods (Curtis et al. 2018). Although this land use system is predicted to be on the decline globally (Heinimann et al. 2017), it is still widespread in the biodiversity hotspot of Madagascar (Myers et al. 2000; UNDP 2016; Curtis et al. 2018). Here, forest cover is steadily declining since decades due to an increased demand for arable lands (Harper et al. 2007; Styger et al. 2007; Vieilledent et al. 2018) and a strong dependence on forest resources (Kiefer et al. 2010; Urech et al. 2015). An impressive example for this is the eastern fraction of the island that was almost entirely covered with rain forest about a century ago, but has lost 61% of its forest cover until 2014 (from 112,000 km2 (Green and Sussman 1990) to 44,100 km2 (Vieilledent et al. 2018)). Forests in south- and central-eastern Madagascar are mainly restricted to higher elevations while lowland areas are almost completely deforested. In contrast to that, northeastern Madagascar is still characterized by a high percentage of forested areas and large continuous forest tracts (Vieilledent et al. 2018). A possible explanation for this might be the remoteness of the region (Weiss et al. 2018). However, the vast majority of the local population lives in a rural setting depending on local resources and shifting cultivation systems to sustain their livelihoods (Styger et al. 2007; Urech et al. 2015; Zähringer et al. 2015, 2016). Recent studies have rigorously documented major land use types and most recent transitions: In accordance to global predictions (Heinimann et al 2017), a shift from traditional land use systems mainly based on shifting cultivation towards land use intensification and the establishment of permanent cultivation can be observed (Zähringer et al. 2015, 2016; Llopis et al. 2019). However, these shifts are severely impeded by poverty, environmental events and the lack of infrastructure (Pfund et al. 2011). Shifting cultivation is therefore still widespread (Zähringer et al. 2015; Llopis et al. 2019) but also due to its high cultural value (e.g., Desbureaux and Brimont 2015) and the fact, that it represents a major strategy for local households to cope with external threats like price volatility of cash crops or environmental events like cyclones (Harvey et al. 2014; Hänke et al. 2018; Llopis et al. 2019).

Apart from this growing body of scientific literature reporting socio-economic and land use trajectories, northeastern Madagascar is particularly rich in biodiversity: About a quarter of all endemic primate species of the island occur in an area covering only 9% of Madagascar’s surface (Mittermeier et al. 2010; Brown and Yoder 2015). These lemurs are particularly useful as an indicator for overall species richness (Muldoon and Goodman 2015) as they receive the highest attention during comprehensive biodiversity inventories. Still, recent descriptions of several new lemur species from the region indicate that its actual biodiversity is not yet fully known (e.g., Radespiel et al. 2008; Ramaromilanto et al. 2009). Same is true for amphibians and reptiles with the recent discovery of new species (Gehring et al. 2010; Vences et al. 2010). Reliable data on other taxa is sparse, but the above given studies indicate high biodiversity on the one hand but also a very high degree of endemism with many species occurring only in very limited areas.

Here we (1) document recent patterns of forest cover and settlement development to (2) determine predictors for the observed changes. Using these data we aim to finally (3) model future deforestation probabilities to delimitate crucial areas for the conservation of biodiversity and locally endemic species.

Material and methods

Study area

Northeastern Madagascar is characterized by a humid climate with a mean annual precipitation of 2086 mm (± 530 mm SD) and rain forests as the natural vegetation (Green and Sussman 1990; Kottek et al. 2006; Fick and Hijmans 2017). The study area covers about 21,919 km2 and represents the entire humid bioclimatic zone between the Maningory river in the South and Mahalevona river in the Northeast (Fig. 1). A total of twelve rivers subdivide the region into 13 inter-river-systems (IRS). In as late as 1953, the study area was covered by a continuous rain forest belt that stretched along the entire east coast of Madagascar (Fig. 1; Vieilledent et al. 2018). Several protected areas are located in the region of which Makira Natural Park (NatP) with connections to Masoala National Park (NP) and Anjanaharibe-Sud Special Reserve (SR) represent the largest complex of protected landscapes in Madagascar (Fig. 1). The majority of the region lacks basic infrastructure and most villages are only accessible by motorbike, boat or on foot (Kremen et al. 1999; Golden 2009; Schüßler et al. 2018b; Andriamihaja et al. 2019). Coastal villages can be also reached via the national road RN5. The terrain is hilly with an underlying altitudinal gradient ranging from the coast (0 m a.s.l.) to about 1353 m a.s.l. at the western border of the study area.

Fig. 1
figure 1

Overview of the study area depicting the 1953 forest extent (derived from Harper et al. 2007), major rivers and protected areas. The red rectangle indicates the area used for the detailed investigation of settlement development

Full size image

Field sampling

We conducted field work in August to October 2017 and 2019 to acquire independent vegetation plots to be used as training and validation data for the land cover classification (see below). Additionally, we used survey walks to verify the location of villages from our settlement map and to delimitate important pathways (see below) using key informant interviews. In 2017, we sampled 99 georeferenced 30 × 30 m vegetation plots and conducted extensive survey walks (115.8 km) covering all land cover types. This sampling was extended in 2019 with a further 277 vegetation plots and 204.1 km of survey walks for validation.

Acquisition and pre-processing of satellite images

We acquired satellite images from the USGS archive (https://earthexplorer.usgs.gov/, Table S1, Suppl. Mat.) from two different time stages. The early 1990s are covered by multiple Landsat 5 Thematic Mapper (TM) scenes while the latest images are Sentinel-2 scenes from 2017/2018 and additional Landsat 8 OLI/TIRS scenes from 2014 to fill remaining gaps from cloud masking (see below). Cloud coverage was a serious problem in the study area. To overcome this issue, we used a longer time period for the early 1990s time stage (Table S1) and relied on Sentinel-2 images for the latest time stage as they offer a higher temporal resolution. Landsat scenes are level-1-terrain corrected, georeferenced products (L1T) which were radiometrically corrected and converted to Top-of-Atmosphere reflectance (Chander et al. 2009; Schüßler et al. 2018a). Sentinel-2 images were already acquired as Top-of-Atmosphere reflectance (L1C products) but required rescaling to actual reflectance values (ESA 2015). As Landsat and Sentinel-2 scenes have different spatial resolutions of 30 × 30 m and 10 × 10 m, respectively, we resampled Sentinel-2 scenes to the coarser resolution using nearest neighbor sampling. Normalized difference vegetation index layers (NDVI) were calculated for each scene and stacked together with red, green, blue and infrared bands to a multi-layer raster dataset (Schüßler et al. 2018a). All calculations and analyses were conducted in ArcMap (ArcGIS Desktop 10.6.1, ESRI, Redlands, USA).

Classification

Pre-processed images were subsequently classified using a supervised maximum likelihood algorithm based on manually selected training areas comprising four different land cover classes: open water bodies (rivers, lakes, ocean), forested areas (primary rain forest of different degradation stages, canopy cover > 80%), non-forest areas (deforested, covered with secondary vegetation of varying heights or agroforestry usually with < 80% canopy cover) and open soils with only marginal plant coverage. We focus our analyses here on the transition between forest and non-forest which is used as a data basis for the subsequent modeling approach (see below). For a more detailed perspective on small-scale land cover transitions in our study area, please refer to Zähringer et al. (2015, 2016) and Llopis et al. (2019).

The classification of the latest time stage was done first and based on our field data (from 2017) to acquire experience on visual appearance of land cover classes in the satellite images. Based on this experience, the classification of the first time stage (early 1990s) was conducted. We therefore retrieved further 84 point locations from available high resolution google earth images for which we could unambiguously assign the land cover class. All images were mosaicked after classification to represent one raster dataset per time stage.

Areas which were covered with clouds were masked and clipped out of the multi-layer stack of bands prior to classification so that previously clouded areas were set to “NoData”. By mosaicking several classified images from the same location (e.g., same path and row for Landsat), these previously masked areas were filled with values. Remaining NoData-areas covered only 2.4% of the study area in the early 1990s and 0.0% in 2017/2018.

The quality of the classified land cover maps was assessed using a three-step procedure (see Schüßler et al. 2018a for a workflow). First, we visually compared each classified land cover map with the underlying satellite images to search for large-scale classification errors. In a second step, we used our field data (from 2017) to identify further misclassifications. When incorrect areas were detected, we refined the training areas accordingly in an iterative process (s. Figure 2 in Schüßler et al. 2018a). When land cover maps showed an acceptable degree of correctness, we used 258 stratified and randomly distributed sample points throughout the study area for the accuracy assessment. For each sample point, we compared the assigned land cover class with the underlying satellite picture for its correct assignment, and user’s, producer’s and overall accuracies were calculated for each time stage (Olofsson et al. 2014). As an additional validation step, we used our field data from 2019 (277 vegetation plots) to independently assess the accuracy of the 2018 land cover classification following the same framework. Classified and validated maps were then used for change detection calculating class-by-class transitions between the two time stages.

Fig. 2
figure 2

Forest cover development from 1990–2018 and change map depicting deforestation during this period. Red arrows indicate the two major forest disruptions of the once continuous rain forest belt of eastern Madagascar

Full size image

Modeling the determinants of deforestation and settlement expansion

Preparation of the input dataset

A variety of approaches has been used in previous studies to model the predictors of deforestation with machine learning algorithms offering the highest potential (e.g., Park et al. 2017). We applied three different modeling approaches to account for this variety and to explore the drivers of deforestation in northeastern Madagascar: a rather classical logistic regression approach (LR, Ludeke et al. 1990; Lin et al. 2011), multivariate adaptive regression splines (MARS) and artificial neural networks (ANNs). Recent studies have demonstrated that the latter two approaches (MARS and ANN) are well suited to capture the complexity of forest cover change (Lin et al. 2011; Tayyebi and Pijanowski 2014; Park et al. 2017).

Every model was supplied with nine explanatory variables to test for their respective contribution in explaining the observed pattern of deforestation and to produce predictions of deforestation probability. These variables were elevation and inclination (both derived from a digital elevation model with a 30 × 30 m resolution, JAXA 2015), distance to the nearest village and pathway (derived from field sampling and google earth images, see below), distance to forest edge and density of forest per km2 in 1990 (derived from the land cover maps), distance to major river and distance to valley bottom (calculated from the digital elevation model, JAXA 2015), as well as the status of forest protection (Suppl. Table S2). Explanatory variables were prepared as raster data layers (30 × 30 m resolution) in ArcMap, and distances were calculated as Euclidean distances in meters. The response variable was binary with the states “forested in 2018” or “deforested since 1990”.

Due to the remoteness of the study area, no reliable settlement map or population estimates are available. We therefore manually extracted the position of villages from recent Google Earth imagery and categorized villages in size according to the number of houses visible. Settlements with 5 to 50 houses were defined as category 1, 51 to 100 houses as category 2, 101 to 150 houses as category 3 etc., up to > 350 houses which were assigned to category 8. Key informant interviews were conducted to identify major pathways of the region and find out which villages or cities they connect. Based on this knowledge, we manually digitalized the exact course of these pathways also from recent Google Earth images. The course of the key informant interviews was approved by the ethics committee of the University of Hildesheim (Germany) and adhered to international and national guidelines (e.g., Wilmé et al. 2016).

We distributed a total of 5000 random points across the forested areas in the early 1990s and values for all explanatory variables (data layers as described above) and the binary response variable were extracted for each point. The points had a minimum distance of 1 km between each other to minimize spatial auto-correlation (Aguiar et al. 2007). This dataset was then exported and analyzed in R (R Core Development Team 2019).

To counteract modeling bias resulting from the random sample of input data points, we calculated ten different models for each of the three modeling approaches (LR, MARS, ANN), each based on a random subset of 4000 points (out of the whole dataset of 5000 points). This subset was then split randomly into a training and a testing dataset (50/50 split).

Modeling approaches

For each of the three different modeling approaches, we calculated ten models, each based on a different subset of data points and a training dataset derived from this subset. The same training dataset was used once for each modeling approach to achieve comparability. For logistic regression, a global model containing all explanatory variables was fitted to the randomly generated training dataset. All possible combinations of the explanatory variables were then computed as candidate models (“dredge” command from “MuMIn” package; Barton 2018) and ranked according to their AIC-values (Akaike Information Criterion). The best model (lowest AIC value) was selected and used for validation and determination of variable importance (see below).

We trained our MARS models accordingly using the “earth” function (“earth” R package, Milborrow 2019) by allowing two interaction terms (determined after tuning). Pruning of models was achieved using the backward method (Milborrow 2019).

ANNs were trained using “neuralnet” function (“neuralnet” R package, Fritsch et al. 2019) with cross-entropy error function, resilient back propagation with weight backtracking (rprop +), and with 20 repetitions. Several ANN architectures were tested and an ANN composed of one hidden layer with four nodes returned the best fit to our data (lowest error). Out of the 20 repetition computed for each randomly generated subset of the data, we selected the ANN returning the lowest error for subsequent validation and evaluation. Explanatory variables were min–max-normalized prior to modeling of ANNs.

Model validation and evaluation of relative variable importance

Trained models were validated with the testing dataset using the area under the Receiver Operating Characteristics curve (AUC, Swets 1988.; “pROC” package, Robin et al. 2011) and a suite of frequently used evaluation metrics including absolute maximum error (AME), mean average error (MAE), root mean squared error (RMSE), Akaike information criterion (AIC), coefficient of determination (R2) and the index of agreement (IoAd; see Dawson et al. (2007) for definitions). These metrics were calculated using the “validann” R package (Humphrey et al. 2017) and used to select the best model from each of the three different approaches for prediction (see below).

We assessed relative variable importance of the explanatory variables using partial dependence plots (Friedman 2001) implemented in the “pdp” R package (Greenwell 2017). This method allows to judge about the relative variable importance in so-called “black-box-models” like ANNs, but also provides comparability across various machine learning approaches (Greenwell et al. 2018). To further promote comparability between the three different modeling approaches, variable importance values were expressed as the percentage of the variable receiving the highest importance in the respective model (i.e., highest importance = 100%).

Settlement development

Settlement development in previously forested areas was modeled in correspondence to the approach described above. The response variable was also binary and defined as villages that were established in areas still covered with forest in 1990 and 627 random points distributed over the forested areas in 1990. We used the proportion of forest per km2 in 1990, elevation, inclination as well as distances to paths, major rivers, minor valleys and villages in 1990 as explanatory variables (see above; Suppl. Table S2).

To further investigate changes in village development, we used high resolution Google Earth images which were available for two different dates (21.05.2000, 06.05.2018) for a southeastern portion of the study area (about 3% of study area, red frame in Fig. 1). We randomly selected 66 villages out of the 301 villages located in the area covered by both images to estimate their size class and the proportion of houses with corrugated iron roofs instead of leaf roofs at both time stages. We used this as an indicator for change as purchasing corrugated iron requires access to the local market system in which produced goods (e.g., from agroforestry) are sold to buy material for building houses. Furthermore, having a proper house with a solid roof is perceived as a status symbol (Diebold 2018). Altogether, we assume that acquisition of corrugated iron shows at least some degree of market access and the shift from pure subsistence farming.

Prediction of deforestation probability

Based on the best models returned from each of the three different modeling approaches concerning forest development, we calculated deforestation probabilities using three new data layers (distance to forest edge and villages in 2018 (in meters) as well as proportion of forest per km2 in 2018) and the set of layers used for training the model algorithms. We did not integrate a temporal dimension in our models as this would require more time slices in between the two time stages we used, which was complicated by frequent cloud coverage. Deforestation probabilities were calculated under a “business-as-usual scenario”. Although change in land use patterns are already documented to some degree from our study region (Zähringer et al. 2015; Llopis et al. 2019), using such a scenario is appropriate for the aim of this study, since forest frontier areas are rather remote and may show the slowest change in land use strategies. Furthermore, projected population growth rates (e.g., UNDP 2016) and the need for arable lands for cash crop cultivation (e.g., Llopis et al. 2019) may still require clearing of currently forested areas and therefore continuing deforestation.

Results

Change in forest cover 1990–2018

Overall accuracies of the 1990s and 2018 land cover classifications were 94.9% and 96.5% while user’s and producer’s accuracies by class did not fall below 93.3% and 81.8%, respectively (Tables S3 and S4, Suppl. Mat.). The independent land cover dataset from 2019 returned comparable overall, user’s and producer’s accuracies of 97.8%, > 89.4% and > 94.4%, respectively (Table S5, Suppl. Mat.).

The forest cover of the entire study area declined from 12,039 km2 (54.9%) in the early 1990s to 7501 km2 (34.2%) in 2018 totaling to a loss of 4537 km2 (20.7%, Table. 1). This forest cover change corresponds to a mean annual deforestation rate of 0.83%. Previously forested areas were converted to non-forest all over the region, and no deforestation hotspot was identified (Fig. 2, change map). About 52% of the remaining forest cover in 2018 is located inside protected areas and afforestation was low with only 394.  km2 or 1.8% of the region.

Table 1 Extent of major land cover classes in 1990s and 2018 with percentage of change and class-to-class changes
Full size table

We detected two major interruptions of the formerly continuous rain forest belt (red arrows in Fig. 2). One of these interruptions is located south of Makira NatP (Fig. 2) between the Fahambahy and Mananara rivers, while the second covers the southern edge of the study area along the Maningory river. Additionally, Zahamena NP lost its direct connection to the wider Ankeniheny-Zahemena Corridor (CAZ), and Betampona Special Nature Reserve (SNR) is separated from the continuous forests of the CAZ by at least 20 km (aerial distance) of deforested land. The area between Mananara-Nord NP, Ambatovaky SR and Marotandrano SR was still forested in 1990 but is also characterized by a matrix of small forest fragments within deforested lands in 2018.

In the northern half of the study area, deforestation around the borders of Makira NatP was most striking with only small-scale signs of forest loss inside the protected area (Fig. 2). In fact, the larger forest tracts in this region were still rather continuous, despite recent deforestation. Makira NatP, Masoala NP, Anjanaharibe-Sud SR and Marojejy NP were still connected by a continuous rain forest belt in 2018.

Determinants of deforestation

All selected models (one per modeling approach) returned similar and relatively high AUC values of 0.827, 0.815 and 0.841 for LR, MARS and ANN, respectively, indicating accurate model performance (Swets 1988). However, ANNs noticeably outperformed LR and MARS models based on lower AME, MAE and RMSE scores as well as higher R2 and IoAd values (Tables 2, S6).

Table 2 Validation metrics for the best models in explaining the observed pattern of forest and settlement development returned from each modeling approach
Full size table

Based on the relative variable importance values returned from the ANNs (highest accuracy scores), status of forest protection was most important followed by distance to the forest edge (Fig. 3a). The remaining seven variables were of lower importance but all contributed to explaining the observed patterns of deforestation. Although this observation is consistent for most variables across the different modeling approaches, some discrepancies should be highlighted: LR models assigned the highest importance to distance from the forest edge, intermediate values to status of forest protection and the lowest to elevation, while in MARS models this pattern was reversed giving highest weight to elevation and the status of forest protection (Fig. 3a).

Fig. 3
figure 3

Comparison of relative variable importance returned from ten different calculations of artificial neural networks (ANN, red), logistic regression (LR, green) and multivariate adaptive regression splines models (MARS, blue) based on ten randomly generated training and testing datasets. a Models concerning forest cover development and b models concerning settlement development. Abbreviations: dfor density of forest per km2, elev elevation, edge distance to the forest edge, incl inclination, path distance to major pathway, prot status of forest protection, rive distance to major river, vall distance to minor valley, vill distance to nearest neighboring village

Full size image

Settlement development

A total of 3,136 villages was detected in the study area of which 627 (20.0%) were found in areas that were still forested in the early 1990s. Out of these 627 newly established villages, 557 (88.8%) were small (< 50 houses) and 9.9% contained between 50 and 100 houses. In general, larger villages can be mainly found in areas that were already deforested before 1990 and closer to the coast (Fig. 4). The highest number of houses per km2 can be found south of the city of Soanierana Ivongo (south east corner of study area), in the middle of the study area around Mananara Avaratra (north of Mananara-Nord NP) and around the city of Maroantsetra (east of Makira NatP; Fig. 4). In general, all deforested areas in 2018 (Fig. 4, in white) are dotted with villages.

Fig. 4
figure 4

Distribution of villages throughout the study area, estimated population densities (approximated as the number of houses per km2) and major pathways along the coast and towards the highlands

Full size image

The detailed analysis of the south-eastern, densely populated portion of the study area revealed, that a majority of the selected villages (72.7%, n = 66) increased in size, and the proportion of houses with corrugated iron roofs strongly increased from on average 9.1% (± 10.8%) in 2000 to 81.7% (± 15.6%) in 2018. We also noticed, that some villages were still newly established (after 2000), even though the area was already deforested in 1990. However, this development was not quantified.

Determinants of settlement expansion since the 1990s (in areas that were still forested in 1990) were also modeled using the three different approaches. ANNs again outperformed LR and MARS models indicated by lower AME, MAE, RMSE scores and higher R2, IoAd and AUC values (Table 2). According to the best ANN model, settlement expansion is best explained by (ordered according to importance): inclination, distance to nearest neighboring village, elevation, distance from major pathways, density of forest in 1990, proximity to the valley bottom and the nearest major river (Fig. 3b). The results from the ANNs were mostly congruent with those obtained from LR (except for elevation). Major discrepancies can only be found between ANN and LR compared with MARS models: here, density of forest in 1990, proximity to pathways and major rivers received a noticeable lower weight while the relative importance of inclination, elevation and distance to neighboring villages was confirmed (Fig. 3b).

Prediction of deforestation probabilities

Although the three modeling approaches differed in weighting relative variable importance, all models predicted the same trends concerning deforestation probabilities. Lower probabilities (green areas in Fig. 5) are predicted for the protected areas and the high altitude forests (e.g., cells A9, A10, B8, C2, C3 in Fig. 5). However, some parts of the protected areas are likely to become deforested in the future and this projection is particularly strong in the MARS model (Fig. 5).

Fig. 5
figure 5

Forest extent in 2018 and deforestation probabilities predicted under three different modeling approaches (LR logistic regression, MARS multivariate adaptive regression splines, ANN artificial neural network)

Full size image

Almost all remaining lowland forests outside the protected areas (cells D4, D5, D7, D8 in Fig. 5) are predicted to face a high probability of deforestation in the future. Moreover, most remaining forested areas in between larger forest blocks and protected areas in 2018 are already highly fragmented (e.g., cells B6, B7, C7, D6, D7) and are also predicted to face high deforestation probabilities in the future.

Discussion

Recent forest cover developments and determinants

Decreasing forest cover has been documented all over Madagascar by various authors (e.g., Harper et al. 2007; Vieilledent et al. 2018) and our study area is no exception to this development. Although the mean annual deforestation rate did not deviate from country-wide estimates (e.g., Vieilledent et al. 2018), more localized studies already reported significantly higher rates from the study region (e.g., Schüßler et al. 2018b; Llopis et al. 2019). This discrepancy is probably due to geographic restrictions towards the forest frontier compared to the broader scope of larger-scale studies and inter-annual variations that can be triggered by political instability, the establishment of protected areas and/or price fluctuations of major cash crops (Llopis et al. 2019).

Shifting cultivation was already identified as the main driver of forest cover loss in northeastern Madagascar (Zähringer et al. 2015). Our results illustrate that despite deforestation happening on small and local scales, it is scattered all over the region and has major effects on regional forest cover. Major determinants of forest conversion are the status of forest protection and proximity to the forest edge. Protected areas effectively prevented deforestation inside their boundaries, although some illegal selective logging is nevertheless reported from national parks such as Masoala NP (Burivalova et al. 2015) and Mananara-Nord NP (pers. comm. Madagascar NP staff member). In contrast, deforestation is omnipresent outside of the protected areas and increasingly along their boundaries as it can be seen around Zahamena NP and along the eastern border of Makira NatP (Fig. 2). Areas within the proximity of the forest edge are primary subjects to deforestation with smaller fragments being particularly vulnerable. Their ongoing loss from the landscape mosaic leads to landscape homogenization and the loss of habitat complexity and connectivity (Urech et al. 2015; Zähringer et al. 2015).

We detected two major interruptions of the former continuous rain forest belt of eastern Madagascar in which only few smaller forest fragments remained in 2018 (red arrows in Fig. 2). We want to highlight that these landscapes are distinctively different from those in other ecoregions of the island. In northern and northwestern Madagascar, for example, at least some lemur species (used here as indicators, Muldoon & Goodman 2015) appear to be able to maintain population connectivity among forest fragments that are embedded in a rather open Savannah matrix (e.g., Quéméré et al. 2010; Salmona et al. 2015; Aleixo-Pais et al. 2018; Steffens and Lehman 2018). Although less studied in the humid ecoregion, it is much more difficult for lemurs and other forest-dwelling species to migrate through the degraded areas of secondary vegetation (Nunziata et al. 2016; Schüßler et al. 2018b; Galán-Acedo et al. 2019). These areas are usually characterized by a mosaic of dense stands of the perennial Aframomum spp. intermixed with smaller tree species like Harungana madagascariensis and Psiadia altissima or, when severely degraded, dense and almost impassable fern- and grasslands with Sticherus flagellatis and the exotic Psidium cattleianum as typical species (Styger et al. 2007; Burivalova et al. 2015; Miller et al. 2018; Schüßler et al. 2018b). Homogenization of landscapes and the disruption of previously continuous forest tracts have major implications for landscape connectivity (Haddad et al. 2015). Increasing fragmentation lowers the resilience of species communities to adapt to environmental changes (Brown et al. 2015; Brown and Yoder 2015). Most alarmingly, landscape connectivity and thus the potential for range shifts along the rainforest belt was already largely restricted and definitely precluded in at least two areas of forest interruption in 2018.

The remaining explanatory variables, despite having lower relative importance values, nevertheless contributed to explaining the observed pattern of deforestation. For example, the relevance of the proximity to villages and to major pathways indicates that deforestation is mainly the result of expanding the size of arable lands within daily walking distance from villages (Kremen et al. 1999; Allnutt et al. 2013). These distances are usually between 1.17 h and 2.5 h of one-way walking time from the village (Fedele et al. 2011; Allnutt et al. 2013) but can also reach almost 3 h corresponding to 7–8 km (Schüßler unpubl. data). Indeed, the degree of forest disturbance and degradation was shown to be inversely correlated with distance to the nearest village (Schüßler et al. 2018b).

Lowest relative importance values were returned for inclination and the distance to the next major river. Shifting cultivation in this area is mainly focused on hill rice, beans or cassava, for which even steep slopes can be suitable for cultivation if located within daily walking distance from the village. Major rivers can be a significant way of transportation (e.g., Kremen et al. 1999), but appear to have had less pronounced effects on deforestation in the region within the last three decades.

We documented a relatively low natural forest regrowth. This is in accordance with other studies (Styger et al. 2007; Manjaribe et al. 2013; Zwartendijk et al. 2017). Natural succession is assumed to be particularly slow when soils are severely degraded (Klanderud et al. 2010; Randriamalala et al. 2012). Limited forest regrowth further increases the pressure on the remaining forests concerning the extraction of wood for construction or cooking as well as for non-timber-forest products (NTFP) like medical plants or fibers. Although fallow lands can provide some of these products (Pfund et al. 2011), forest resources play a key role in local livelihoods (Ratsimbazafy et al. 2012). To counteract the progressive isolation of forest fragments and the loss of crucial forest products for local communities, active reforestation measures are urgently needed. Land users in the region are particularly experienced in raising and planting tree saplings, for example to plant clove, coffee and cacao trees. Cash crops are obviously beneficial to invest in, since income opportunities from crop sales are likely. Plantation of non-cash crop trees in terms of reforestation, however, will require incentives to facilitate the involvement of local communities. Frameworks still need to be developed, evaluated and communicated scientifically, but payments for certain numbers of planted trees with controlling for their survival after some time may be conceivable. Only few restoration projects have been scientifically evaluated to date, but promising approaches are based on raising fast growing early successional native tree species in tree nurseries in order to restore soil conditions and to promote natural succession (see Manjaribe et al. 2013; Rajaonarimamy et al. 2017; Zwartendijk et al. 2017). In turn, we highly discourage the use of non-native tree species, since unpredictable side-effects can have massive negative effects on Madagascar’s fragile endemic biodiversity.

Settlement development and determinants

Many low-income countries are characterized by strong population growth (UNDP 2016). Madagascar’s population rose from about 7.9 Mio. in 1975 (UNDP 2005) to 24.2 Mio. in 2015 and may further increase (3.0% annual growth rate) to 36.0 Mio. people by 2030 (UNDP 2016). It is difficult to predict how human population density will develop in our study area. Apart from birth rates, domestic immigration triggered by recent price booms in major cash crops (vanilla and cloves) may accelerate population growth in the region (Hänke et al. 2018; Zhu 2018). So far, Jones et al. (2018) documented only rather local migratory trends in the southern part of the study area (< 50 km movements reported by about 90% of interviewed households from the Ankeniheny-Zahamena Corridor). These short distance movements were mainly driven by the search for accessible lands (Jones et al. 2018). At least one fifth of all villages in our study area were newly established since 1990 in areas that were still forested at that time. We identified inclination and proximity to the nearest neighboring village as the major predictors for places in which new villages were built. Additionally, the altitude of the area and the proximity to the nearest major pathway also contributed to explaining the observed settlement expansion. This resembles a housing sprawl along the forest edge, which is in accordance with short distance movements to access and acquire new lands for agricultural use (Nawrotzki et al. 2012; Jones et al. 2018). Although there are first tendencies towards land use intensification and the establishment of permanent land use systems (Zähringer et al. 2015; Llopis et al. 2019), this documented housing sprawl until 2018 is alarming, since it concerns almost all forests outside protected areas and steep mountains, irrespective of their potential importance as landscape corridors and biodiversity reservoirs.

Although we could not quantify this development in areas that were already deforested prior to 1990, it occurred to some degree. Most villages in these places, the southeastern part of the study area, grew in size in accordance with the documented population growth (UNDP 2016). Houses are traditionally built from natural material extracted from the forests around the villages. Interestingly, the proportion of houses with natural leave roofs significantly declined and most houses were equipped with iron roofs in 2018. This development can be explained in two ways: first, we investigated an area that has only restricted access to forest resources and intact forest tracts are far away from the villages. Under these circumstances, the availability of leaves to tile roofs may be low, and people may be forced to use other materials. On the other hand, this observation implies that households may be shifting away from pure subsistence farming towards the production and sale of products in regional and even global markets. In fact, recent price booms of vanilla and cloves have resulted in significant income increases for smallholder farmers and in remarkable improvements of the socio-economic situation in some parts of northeastern Madagascar (Hänke et al. 2018; Zhu 2018). Furthermore, our study shows that most of the larger villages lie in areas close to the coastline that were already almost entirely deforested in 1990. These places are still well suited for land use intensification (Zähringer et al. 2016) and may therefore allow higher population densities than the more rugged terrain along the remaining forest frontier.

Future development in forest cover

Forest cover loss in northeastern Madagascar is mainly driven by smallholder farmers practicing shifting cultivation (Zähringer et al. 2015). This is similar to most of the world’s tropics, though, the extent of shifting cultivation is expected to decline within the next few decades (Heinimann et al. 2017). First indications for this development were found in our study area with land use intensification and the establishment of permanent cultivation systems (Zähringer et al. 2015; Llopis et al. 2019). However, the transition from traditional shifting cultivation towards permanent land use practices and thus a significant decrease in deforestation is unlikely to happen within a short time interval in those remote areas that are at the frontier of deforestation activities. We therefore aimed at identifying areas that are particularly threatened by further deforestation under a “business-as-usual” scenario.

Although our models weighted input variables differently, the results of the predictive modeling are mostly consistent across modeling approaches. Protected areas and highland forests face the lowest deforestation probabilities due to their inaccessibility and the relative efficiency of forest protection inside protected areas compared to those outside their boundaries. However, our predictions also illustrate that the pressure on at least some parts of the protected areas may also be increasing in the future.

The risks of landscape homogenization, i.e., the decrease of habitat complexity towards more uniform land use types, were already noted by other authors (Urech et al. 2015; Zähringer et al. 2015) as well as by our own analysis that predicts further loss of smaller forest fragments from the landscape mosaic. Despite presumed negative effects on biodiversity due to decreasing habitat size and diversity (Haddad et al. 2015), the complete loss of forests in some areas will have adverse effects on local climates and water availability. Land users in northeastern Madagascar already noticed problems in rice paddy irrigation due to decreasing forest cover and water availability around their fields in two study sites (Llopis et al. 2019). Our survey walks revealed that this problem is rather widespread across the region. Diminishing forest resources will therefore most likely decrease the productivity of intensive land use systems that produce the most important staple crop in Madagascar. Unfortunately, awareness about the scarcity of forest resources and their ecosystem services seems to be highest in areas with limited or absent forest cover, whereas people in currently forest-rich areas appear to be less aware about decreasing forest coverage (Urech et al 2015). This mismatch in perceptions needs to be targeted in forest frontier areas to illustrate the finiteness of forests and the crucial services they provide. Promising approaches using environmental education are described elsewhere (e.g., Schüßler et al. 2019).

Almost all of the already restricted lowland forests of the study area are predicted to face high deforestation probabilities (except Mananara-Nord NP). These findings highlight that lowland forests require increased conservation attention to protect locally endemic lowland forest species, but also that effective conservation management is crucial for Mananara-Nord NP as the only officially protected tract of lowland rain forest in the region. Why do comparably small forests matter as long as there are larger forested areas with low deforestation probabilities available in the highlands of the southern part of the region and in Makira NatP in the northern part? The study area can be subdivided into 13 inter-river-systems (IRSs) that may differ in their species inventories. Even single IRSs can be of particular importance for micro-endemic species like mouse lemurs (Radespiel et al. 2008, 2012), amphibians (Gehring et al 2010; Vences et al. 2010) or reptiles (Gehring et al. 2010). Furthermore, the lowland rain forests of eastern Madagascar have been stated the key for lemur conservation in the long term (Campera et al. 2019). At least two of these IRSs are predicted to have almost no areas of low deforestation probability (fields C7 and D8 as well as field D6 in Fig. 5), but some forests are already protected by the local community (see Miller et al. 2018; Schüßler et al. 2018b). These forests still harbor a high lemur diversity and undisturbed habitats despite being surrounded by a deforested matrix. Although these community-managed areas are so far effectively protected against deforestation (Miller et al. 2018; Schüßler et al. 2018b), it appears that shrinking resources around these forests will most likely increase the pressure on the currently protected forest patches (A. Miller pers. com.; Schüßler et al. 2018b).

Four other IRSs contain rather limited areas of low deforestation probability or areas predominantly restricted to higher elevations (fields C9, C10, D5 in Fig. 5) without any current forest protection scheme. Deforestation of whole IRS will inevitably increase extinction risk of micro-endemic species in an area that is still not well covered with biodiversity surveys and even new species to be found. Certain areas should be therefore highlighted in the following to be considered for the expansion of the protected area network of the region. These forests could act as stepping stones offering habitats for dispersing species or refuges for species retreating from human-dominated landscapes. We focus on areas that are predicted to have low deforestation probabilities as they may have lower value for communities in terms of expanding shifting cultivation or agroforestry systems due to their remoteness and/or inaccessibility to minimize land use conflicts. However, our village map illustrates that all forests are surrounded by villages. Protection schemes should therefore be promoted that encourage communities to set aside a certain portion of their forest area for biodiversity conservation. As described above, community-managed forests are already established in some places away from the major protected areas as well as in the buffer zones around Makira NatP (Brimont et al. 2015). This scheme of community-managed forest should also be promoted in other areas by focusing particularly on land users most vulnerable to environmental or other external changes (Harvey et al. 2014; Brimont et al. 2015). Recent evaluations of community-based approaches have illustrated that the danger of marginalization of poorer households, non-targeting of remote but most affected households and elite capture are severe issues to tackle (Brimont et al. 2015; Cullman 2015; Ward et al. 2017). The growing body of studies illustrating the pitfalls of previous approaches should be considered carefully when establishing new community-managed forests.

One of these potential stepping stone forests is located in the highlands west of Mananara-Nord NP (field C6 in Fig. 5) spanning two IRS in the source region of the Anove river. As of 2018, there is still a large tract of forest left in this area and a zonation scheme integrating different land use types may be beneficial to achieve protection of a core forest zone. Further forested areas that should receive more attention are the lowland forests of Manompana (see above) that are currently protected by the local community, as well as the lowland and mid-altitude forest south of Makira NatP (fields C5 and D5 in Fig. 5). The latter area represents the only forests with low deforestation probabilities for this IRS and could therefore become an important stepping stone forest. Another important stepping stone forest is located in the highlands of the southern part of the region (fields A8 and B8 in Fig. 5). Encouraging forest management based on a zonation system should conserve again important habitats spanning two IRS. We want to emphasize that we are not promoting to set apart extensive forest areas to establish new protected areas. The forests we highlighted here rather represent the minimal recommendation to protect most of the regional biodiversity at least in a shorter term. Community-based forest management approaches appear to be valuable in several parts of the study area and should receive further attention in protecting key forest areas.

Synthesis and conclusion

With this study we have illustrated the decrease of forest cover over the last three decades in northeastern Madagascar. The predominant driver for this development is small-scale but ubiquitous conversion of forests into agricultural fields with accessibility to forests as major predictors. A large proportion of villages were newly established since 1990 at the expense of previously forested areas. We pointed out, that this housing sprawl may have significant and negative effects on landscape connectivity for forest dwelling species.

We strongly advocate for integrative land use planning that combines aspects of sustainable landscape development that harness the benefits of land use intensification (i.e., reduction of the need for shifting cultivation) with the valuable role of agroforestry systems and widely grown cash crops as well as biodiversity conservation. In addition, we proposed several areas that could act as important biodiversity reservoirs and stepping stone habitats that may even facilitate future range shifts in one of the most diverse regions of Madagascar (Brown and Yoder 2015).

Changes in land use patterns are already taking place throughout the study area and are triggered by diminishing resources, establishment of protected areas and recent price booms in major cash crops (Urech et al. 2015; Llopis et al. 2019). However, the integration of human and wildlife needs on a regional scale requires external input by governmental or non-governmental stakeholders. In fact, in regions of weak governance and remoteness such as in northeastern Madagascar, we see the highest potential for the involvement of national or experienced international non-governmental organizations (NGOs). They could focus on important areas like the stepping stone forests that we identified with target-oriented interventions such as community-managed land use zonation schemes or reforestation projects (Schüßler et al. 2019; Nicli et al. 2019). Such non-governmental stakeholders may best encourage sustainable land use planning and integrate the livelihood needs of rural communities with the challenge of long-term conservation of the unique biodiversity in the fragile landscapes of northeastern Madagascar.