Predicting range shifts of African apes under global change scenarios

Aim Modelling African great ape distribution has until now focused on current or past conditions, whilst future scenarios remain scarcely explored. Using an ensemble forecasting approach, we predicted changes in taxon-specific distribution under future scenarios of climate, land-use and human population changes. Location Sub-Saharan Africa Methods We compiled occurrence data on African ape populations from the IUCN A.P.E.S. database and extracted relevant human-, climate- and habitat-related predictors representing current and future (2050) conditions to predict taxon-specific distribution under a best- and a worst-case scenario, using ensemble forecasting. Given the large effect on model predictions, we further tested algorithm sensitivity by considering default and non-default modelling options. The latter included interactions between predictors and polynomial terms in correlative algorithms. Results The future distributions of gorilla and bonobo populations are likely to be directly determined by climate-related variables. In contrast, future chimpanzee distribution is influenced mostly by anthropogenic variables. Both our modelling approaches produced similar model accuracy, although a slight difference in the magnitude of range change was found for Gorilla beringei beringei, G. gorilla diehli, and Pan troglodytes schweinfurthii. On average, a decline of 50% of the geographic range (non-default; or 55% default) is expected under the best scenario if no dispersal occurs (57% non-default or 58% default in worst scenario). However, new areas of suitable habitat are predicted to become available for most taxa if dispersal occurs (81% or 103% best, 93% or 91% worst, non-default and default, respectively), except for G. b. beringei. Main Conclusions Despite the uncertainty in predicting the precise proportion of suitable habitat by 2050, both modelling approaches predict large range losses for all African apes. Thus, conservation planners urgently need to integrate land-use planning and simultaneously support conservation and climate change mitigation measures at all decision-making levels both in range countries and abroad.


| INTRODUC TI ON
Currently, a major conservation challenge is to assess the potential future effects of climate and land use changes on species distributions, typically through the use of species distribution models (SDMs) and usually under a range of future environmental scenarios. SDMs are widely used to predict and map species' ecological niches through time and space Guillera-Arroita et al., 2015;Hao et al., 2019). Importantly, SDMs can inform spatial prioritization decisions for conservation and management actions, such as identification of strategic locations for new conservation or survey sites, and predicting future distributions (Araújo & New, 2007;Guillera-Arroita et al., 2015).
Changes in climate and land use are among the main global threats to biodiversity, and therefore, how the synergistic interactions between these drivers impact species are an important area of research (Oliver & Morecroft, 2014). Newbold (2018) assessed the single and combined effects of future climate and land use change on local vertebrate biodiversity. They found that climate change is likely to be the principal driver of species range change in coming decades, equalling or surpassing the potential effects of land use change by 2070. Similar results were reported for orangutans (Struebig et al., 2015). Because human population growth is already an extinction threat to many species (McKee et al., 2013), it is also important to determine how human distribution will impact future species presence (Jones & O'Neill, 2016).
Many primates are facing imminent extinction, due to the direct impact of extensive habitat loss and fragmentation, land use change and hunting, and indirect effects linked to global commodity growth and trade (Estrada et al., 2018). Climate change is a delocalized, multi-faceted driver to add to the list. It exposes many species, especially forest-dwelling primates, to climatically unsuitable conditions . Primates have relatively limited dispersal abilities for their slow reproduction, low population densities, dietary requirements and poor thermoregulation, and a predicted reduction of up to 86% of Neotropical primate ranges with >3°C warming is likely to constrain their dispersal, resulting in elevated risks of extinction .  (Junker et al., 2012) as well as large population losses (Kuehl et al., 2017;Plumptre et al., 2016;Strindberg et al., 2018) caused by human activities and/or infectious epidemics (Walsh et al., 2003).
Many African apes live in areas that are suitable for agricultural expansion and 58.7% of oil palm concessions currently overlap with African ape ranges (Wich, Garcia-Ulloa, Kühl, et al., 2014). Moreover, massive development corridors (Heinicke et al., 2019) and mining activities (Howard, 2019) in their geographic ranges are projected to expand considerably and to disrupt ape habitat connectivity and accelerate habitat loss.
Most African apes occur outside protected areas (PAs) (Freeman et al., 2018;Heinicke et al., 2019;Strindberg et al., 2018;Wich, Garcia-Ulloa, Kühl, et al., 2014). Importantly, PAs will not be exempt from climate change (Araújo et al., 2011), and shifts in species range as predicted by future scenarios would certainly determine the degree of species representation inside and outside PAs. Improving the effectiveness of conservation efforts inside and outside PAs as well as habitat connectivity would allow apes to disperse to new climatically suitable areas, and favour ape population survival in the long term.
The influence of the combined effects of current climate conditions and anthropogenic disturbances on African ape distribution have been widely explored (Clee et al., 2015;Hickey et al., 2013;Junker et al., 2012;Plumptre et al., 2016;Strindberg et al., 2018). In contrast, few studies have only examined future effects of climate change (Clee et al., 2015;Lehmann et al., 2010;Thorne et al., 2013) or human disturbances (Wich, Garcia-Ulloa, Kühl, et al., 2014), but how future synergistic interactions among climate, land use and human population changes will affect African apes and their habitat has been largely unexplored.
Here, we combine data on projected climate, land use and human population changes to model taxon-specific distribution of African apes for the year 2050. We use the most comprehensive database on ape populations available, the IUCN SSC Ape Populations, Environments and Surveys database (A.P.E.S.) to predict the distribution of great apes on the African continent under best-and worstcase scenarios. We subsequently employ an ensemble forecasting approach to reduce the uncertainty among different models and future scenarios (Araújo & New, 2007;Thuiller, 2004) and estimate the proportional change in range size in 2050 relative to current estimated range sizes for African apes by considering (1) only areas outside PAs by assuming complete effectiveness of PA management and consequently complete range stability within PAs, (2) the entire study region and (3) interspecies range overlap. Specifically, we addressed the following questions: (a) What is the extent of range loss by 2050? and (b) What is the proportion of new range predicted under future scenarios? Given that range loss and gain occur at different time scales, we predict that massive range loss is likely to occur in the next 30 years, but range gain is more uncertain given that African apes will not be able to occupy these new areas immediately due to their limited dispersal capacity (Schloss et al., 2012), migration lag and ecological constraints.
We obtained a total of 62,469 presence records across all (sub)species (hereafter taxon) (Appendix S1 in Supporting Information, Table   S1.1). We first checked the spatial autocorrelation of these records using Ripley's K-function (spatstat package; Baddeley et al., 2015) in R version 4.0.2 (R Development Core Team, 2020) and then rarefied the presence data by removing those points within a certain distance of one another (ecospat package;Di Cola et al., 2017), resulting in 5,203 presence records ( Fig. S1.1, Table S1.1).
Although these data may be spatially biased as sampling effort is unevenly spread over the ape ranges, presence-only data are commonly the most available and hence most frequently used in SDMs (Phillips et al., 2009). The taxon occurrence data we used were collected during systematic site-based wildlife and human impact surveys, often in or close to PAs, Forest Stewardship Council (FSC)certified and other logging concessions, or from habituated populations. Those surveys were generally based on some prior knowledge of occurrence which can distort an SDM (Phillips et al., 2009). Different approaches have been applied to account for biased datasets: random background, bias background, geographic thinning/ filtering and environmental filtering (Aiello-Lammens et al., 2015;Fourcade et al., 2014;Phillips et al., 2009;Varela et al., 2014). Thus, we considered all approaches, and we included distances to roads, villages and PAs for the bias background as they are known to influence the distribution of African apes (Carvalho et al., 2013). We extracted data on roads and villages (from http://sedac.ciesin.colum bia.edu/) and PAs (from https://www.prote ctedp lanet.net/) within each taxon's range. For each taxon, we selected the approach with the best performance by visually inspecting the greatest overlap between taxon occurrence and each sampling bias layer ( Fig. S1.2).
Given that the geographic thinning approach performed best for all taxa, we integrated it into the SDMs for sampling bias correction  Table S1.1). We also checked the spatial autocorrelation of this bias layer using Ripley's K-function ( Fig. S1.1).
We delineated taxon-specific study regions to avoid unrealistic geographical predictions (Anderson & Gonzalez, 2011). For this, we created buffers bounding IUCN range polygons (IUCN (2020)) and including all occurrence data for each taxon (Table S1.1) (Jantz et al., 2016;Junker et al., 2012;Thorne et al., 2013). We defined the size of each buffer according to the range size of each taxon (Table   S1.1). Whenever the buffer overlapped with a known geographic barrier to ape dispersal (e.g. major rivers), we disregarded that area.
Model algorithms require presence and absence data, so we randomly generated a set of 10,000 pseudo-absence occurrences (Barbet-Massin et al., 2012;Guillera-Arroita et al., 2015;Phillips et al., 2009) within the study region of each taxon, except for G. b. beringei. Only 1,000 background occurrences were created for mountain gorillas due to their small range.

| Predictor variables
We selected predictor variables based on their importance for African ape ecology (Clee et al., 2015; Lehmann et al., 2010), while guaranteeing data availability for current and future (2050) conditions under best-and worst-case scenarios and minimizing correlation between variables. We compiled the following climatic variables for present and future conditions from Worldclim (periods of 1950-2000and 2050Hijmans et al., 2005): annual mean temperature (bio1), maximum temperature (bio5), annual precipitation (bio12) and precipitation seasonality (bio15). For future predictions, we chose a best scenario (i.e. high mitigation scenario, CCSM4 RCP 4.5) and a worst scenario (i.e. low mitigation scenario, HadGEM-ES RCP 8.5; for more details see Carvalho et al., 2019).
Land use/cover data for current conditions and 2050 projections were compiled from the Land use Harmonization Project (period of 1500-2100; Table S1.1; Chini et al., 2014;Hurtt et al., 2011). This dataset represents a set of land use change and emission scenarios for studies of human impact on the past and future global carbonclimate system. Again, we considered a best scenario (MiniCam RCP 4.5) and a worst scenario (MESSAGE RCP 8.5) .
We focused on the land use states that best represent land cover types where great apes can be found: primary land (i.e. natural vegetation, either forest or non-forest, undisturbed by humans), secondary land (i.e. natural vegetation previously disturbed by agriculture or wood harvesting) and cropland.
We based human population scenarios on a new set of future societal development scenarios, namely Shared Socioeconomic Pathways (SSP) (Jones & O'Neill, 2016). These future scenarios are based on both qualitative narratives of future development and quantitative projections of key elements such as human population growth at the national level, educational composition, urbanization and economic growth. These data are available from 2010 to 2100 for urban and rural populations. We used two future scenarios, SSP1 and SSP3, given that they represent best and worst scenarios, respectively.
Firstly, we extracted all variables for the extent of the study region of each taxon, resampled onto a 5km x 5km equal-area grid and transformed them into the WGS 1984 geographic coordinate system. Secondly, for each taxon, we computed and visualized a Pearson correlation matrix to assess the collinearity among variables  Marco & Nóbrega, 2018). For that, we considered both present and future conditions to perform a PCA for each taxon and then used PCA loadings for each scenario to create PCA-derived variables. Finally, we selected the loadings of the first four PCA axes because they explained more than 79% of variance (Table S1.2). We performed data analyses using the software R and ArcMap version 10.7.1 (ESRI, 2011).

| SDM performance and ensemble forecasting
We predicted future African ape distributions using an ensemble forecasting approach (i.e. combining predictions from individual models into an ensemble as implemented in the biomod2 package; Thuiller et al., 2016). We selected two correlative algorithms, generalized linear model (GLM) and generalized additive model (GAM), and three machine-learning techniques, Maxent, random forest (RF) and artificial neural networks (ANN) to build predictive SDMs for each species. These algorithms have been shown to perform well in previous SDMs (Elith et al., 2006;Thuiller et al., 2009). We decided to keep the default settings of "biomod2" for each algorithm to avoid an overwhelming complexity of the study outcomes and for ease of comparison between taxa.
For the present time period only, we assessed the predictive performance of each model through cross-validation using a bootstrap approach, that is partitioning of the presence data, using 80% of presences, randomly selected, for model calibration and 20% for evaluation, and repeating this procedure five times (Thuiller et al., 2009). We evaluated the performance of each model by the "true skill statistic" metric (TSS) (Allouche et al., 2006). TSS is an accuracy measure that accounts both for omission errors (i.e. the percentage of true presences predicted as absences are minimized) and commission errors (i.e. the percentage of true absences predicted as presences are minimized), is affected by prevalence (Leroy et al., 2018;Somodi et al., 2017) and ranges from −1 to 1, with a prediction accuracy considered similar to "random" when ≤0, "poor" in the range 0.2-0.5, "useful" in the range 0.6-0.8 and "good" to "excellent" when >0.8 (Allouche et al., 2006). Ensemble forecasting has been widely employed to reduce the uncertainties associated with using a single algorithm and is a useful method to account for uncertainties of extrapolation of speciesenvironment relationships outside the environments sampled by the species data (Araújo & New, 2007;Hao et al., 2019;Thuiller et al., ,2009Thuiller et al., , , 2019. We chose to apply the weighted mean ensemble method, which scales predictions of different models by weights based on some measure of predictive performance (Araújo & New, 2007;Thuiller et al., 2009). We included only individual models that reached at least "useful" predictive accuracies (TSS>0.7) in ensemble models whenever possible (except for G. b. graueri, G. g. gorilla and P. t. troglodytes TSS>0.5; P. paniscus and P. t. verus TSS>0.6), to map the current and future range predicted for each taxon . For each modelling approach, we repeated the modelling five times (cross-validation) and given the five modelling algorithms and the three repetitions for variable importance (see below), we obtained an ensemble of 75 predicted distributions for each taxon for each time period (present and 2050) and future scenarios (best and worst scenarios).

| Relative importance of PCA-derived variables
For each taxon, we calculated the importance of PCA-derived variables by correlating the fitted values of the full models with those from the model in which the values of the PCA-derived variables were randomly permuted. We repeated this procedure three times (default settings of "biomod2" were used) and used the average Pearson's correlation to measure variable importance. A high correlation between the values from the full and permuted models indicates that the PCA-derived variable has a low importance, contributing poorly to the model. We then ranked each PCA-derived variable value based on the correlation coefficients and reversed its relative importance and scaled from 0 to 1, the more influential PCAderived variables for the model representing those with a higher relative importance (Thuiller et al., 2009). By identifying the most influential PCA-derived variable for the model, we further selected the original variables with the strongest loading (i.e. >|0.4|) to better describe the most important variables influencing each taxon's distribution.

| Species range change
For each taxon, we estimated the proportional change in range size, in 2050 compared to the present, by subtracting the future prediction ensemble output from the SDMs for the best and worst scenarios from that under current conditions. Firstly, we considered as a cut-off the maximum value of TSS to create binary predictive outputs from ensemble models . Secondly, we identified areas of range loss (i.e. sites where the species is present at the moment but is likely to be absent in the future), gain (i.e. sites where the species is absent at the moment but is likely to be present in the future), stability (i.e. sites where the species is potentially present at the moment and is likely to be present in the future) and absence (i.e. sites where the species is absent at the moment and is likely to be absent in the future). For this, we considered range change under two contrasting dispersal scenarios: full dispersal, which assumes that the species can disperse to new suitable areas in the future; and no dispersal, which assumes that the species will be unable to disperse and only the overlap between present and future distributions will be the expected range for the species (Thomas et al., 2004).
Finally, we extracted this information for (1) areas outside PAs only, (2) the entire study region and (3) range overlap where areas of loss, gain, stability and absence were the same between taxa.

| RE SULTS
In general, predictive performance of the individual models based on TSS was "poor" to "excellent," depending on the algorithm and taxon (Appendix S2 in Supporting Information, Fig. S2.1a). On average, RF models performed best relative to GLMs which performed worst at predicting species distributions. Importantly, with TSS scores >0.5 ensemble models had good predictive performance and clearly out-  Table S1.2). Our ensemble models indicated that the distribution of all taxa is strongly influenced by all predictors, particularly by human-related variables (Figures 1b, 2; Table S1.2).
We do not show results regarding range change for mountain and Cross River gorillas given the extreme range loss obtained, that is complete loss of suitable habitat and no new suitable habitat predicted under both future scenarios. Under the assumption of complete range stability in PAs, both future scenarios agree that, on average, more than half of suitable range is likely to be lost if no dispersal occurs (50% best, 61% worst) (Figures 3, 4, Fig S2.2, F I G U R E 1 Results of the ensemble models for each African ape taxon. (a) Predictive performance (mean TSS values and respective standard deviation (SD)) and (b) PCA-derived variable importance (mean and SD of the correlation values). Taxon name abbreviations: gbb -Gorilla beringei beringei, gbg -G. b. graueri, ggd -Gorilla gorilla diehli, ggg -G. g. gorilla, ppan -Pan paniscus, pte -Pan troglodytes ellioti, pts -P. t. schweinfurthii, ptt -P. t troglodytes, ptv -P. t. verus. Background colours in plot (a) corresponds to TSS performance: red -"poor," yellow -"useful," and green -"good to excellent." For details about the PCA loadings in plot (b), see Figure 2  Table 1). Most suitable range is predicted to be lost if the entire study region is considered (85% best, 94% worst). However, if dispersal occurs, most suitable range gain is predicted outside PAs under both future scenarios (52% best or 21% worst versus 66% best or 24% worst for the entire study region) (Figures 3, 4, Table 1). The range of G. b. graueri overlaps completely with that of P. t. schweinfurthii, but only 16% of the latter overlaps with the former ( Figure 5).
In contrast, both ranges of G. g. gorilla and P. t. troglodytes fully overlap (99% and 95%, respectively). Under both future scenarios, no reduction in range overlap between G. b. graueri and P. t. schweinfurthii is predicted if no dispersal occurs, but half of gains are expected outside PAs if dispersal occurs ( Figure 5, Table 1). In contrast, more than half of losses and gains are predicted where the ranges of G. g. gorilla and P. t. troglodytes overlap, particularly outside of PAs. Moreover, range stability is expected outside PAs (G. g. gorilla: 27% best, 4% worst; P. t. troglodytes: 18% best, 4% worst).

| Gorilla beringei beringei (mountain gorilla)
All model algorithms performed equally well at predicting mountain gorilla distribution (Figure 1a, Fig S2.1a). Cropland, human population and secondary land were important predictors in the majority of individual and ensemble models, whereas annual precipitation and secondary land were the strongest determinants of mountain gorilla distribution in ANN models (Figures 1b, 2, Fig S2.1b; Table S1.2). This taxon is confined to fragmented habitat remnants in a sea of agriculture, within which human population (2-10 people km −2 ) is low, and secondary land (>60%) and annual precipitation (1,200-1,600 mm) are high (Fig. S1.4). All predictors will increase by 2050 under both future scenarios, except for annual precipitation and secondary land which are predicted to decrease under the worst scenario. The representativeness of PAs within the study region of mountain gorillas is the highest (36%) among all taxa.

| Gorilla beringei graueri (Grauer's gorilla)
On average, RF, Maxent and GAM models performed best in predicting the distribution of Grauer's gorillas (Fig S2.1a), in which annual and maximum temperatures were the most important explanatory variables (Figures 1b, 2, Fig S2.1b; Table S1.2). In contrast, cropland, human population and primary land were the most influential variables in those models that performed more poorly such as GLM and ANN models. The study region of this taxon is characterized by mean annual temperatures of 14-20°C, maximum temperatures of 20-30°C, low human population (5-10 people km -2 ), high primary land cover (>80%) and very low cropland cover (<8%) (Fig. S1.4).
Both climatic variables, proportion of cropland and human population density are expected to increase and primary land to decrease under both future scenarios.
Only one quarter of the study region is in PAs. Under the assumption of complete range stability in PAs, Grauer's gorillas are predicted to lose three quarters of their range under both future scenarios if no dispersal occurs, with most range predicted to be lost if the entire study region is considered (Figure 4, Fig S2.3; Table 1).
No new suitable areas are expected outside PAs, but range gain is predicted inside PAs under the best scenario if dispersal occurs.
Under both future scenarios, range loss is expected outside PAs where losses are also predicted for eastern chimpanzees ( Figure 5, Table 1). Under the best scenario, almost half of the gains overlap with range gains predicted for eastern chimpanzees.

| Gorilla gorilla diehli (Cross River gorilla)
All model algorithms performed equally well at predicting the distribution of Cross River gorillas (Figure 1a, Fig S2.1a), and all predictors ranked equally in importance in both individual and ensemble models ( Figure 1b, Fig S2.1b). Annual temperature, maximum temperature and primary land were strongly associated with PC1, and annual precipitation, seasonal variation in precipitation (precipitation seasonality), human population and secondary land with PC2 ( Figure 2, Table   S1.2). Annual temperatures of 20-26°C, maximum temperatures of

| Gorilla gorilla gorilla (western lowland gorilla)
On average, RF and ANN models performed best at predicting western lowland gorilla distribution (Fig. S2.1a). All individual and ensemble models showed that seasonal variation of precipitation and the proportion of secondary land are important predictors of western lowland gorilla distribution (Figures 1b, 2, S2.1b; Table S1.2). Areas characterized by high seasonal variation in precipitation (30-80 mm) and presence of secondary land (20%-80%) provide suitable conditions for the persistence of this taxon (Fig. S1.4). According to both future scenarios, seasonal variation in precipitation will not change, but secondary land cover is predicted to increase under the best scenario.
Only 17% of the study region of western lowland gorillas is in PAs. Assuming no dispersal, more than half of predicted range loss will occur outside PAs under both future scenarios (Figure 4, Fig.   S2.3; Table 1). If the entire study region is considered, a loss of more than three quarters of their range is predicted under the best scenario, and most of the taxon's range is likely to disappear under the worst scenario. With dispersal, however, substantial range increases outside PAs are predicted under both future scenarios, with a slight increase if the whole study region is considered (Figure 4, Fig. S2.3; Table 1). Under both future scenarios, most losses are expected outside PAs where losses for central chimpanzees were also predicted ( Figure 5, Table 1). Slightly higher values were found for the entire study region. Gains in range overlap are predicted outside PAs ( Figure 5, Table 1). Importantly, 27% or 4% of range stability is expected where central chimpanzees are also predicted to be present under the best and worst scenarios, respectively.

| Pan paniscus (bonobo)
On average, RF and ANN models performed best in predicting bonobo distribution (Fig. S2.1a). Annual temperature, maximum  Table 1). The taxon's range is predicted to expand into new areas and, if bonobos disperse, substantial range gains are predicted outside PAs under future scenarios, with a slight increase expected for the entire study region (Table 1).

| Pan troglodytes ellioti (Nigeria-Cameroon chimpanzee)
All model algorithms performed equally well at predicting Nigeria-Cameroon chimpanzee distribution (Figure 1a, Fig S2.1a). Annual F I G U R E 4 Predicted percentage change outside protected areas (filled bars) in African ape ranges by 2050 under the best-and the worst-case scenario, assuming either no dispersal (loss) and dispersal (gain). Shaded and filled bars together represent the results for the entire study region. Taxon name abbreviations: gbg -G. b. graueri, ggg -G. g. gorilla, ppan -Pan paniscus, pte -Pan troglodytes ellioti, pts -P. t. schweinfurthii, ptt -P. t troglodytes, ptv -P. t. verus precipitation and seasonal variation of precipitation were the best predictors in both individual and ensemble models (Figures 1b, 2 Table S1.2). Areas with high annual precipitation (2,000-3,500 mm) and pronounced seasonal variation of precipitation (50-90 mm) offer suitable conditions for Nigeria-Cameroon chimpanzees ( Fig. S1.4). Under the worst scenario, both annual precipitation and seasonal variation of precipitation are predicted to decrease.
Only one tenth of the study region is covered by PAs. If no dispersal occurs, most of the taxon's range is predicted to be lost outside PAs under both future scenarios (Figure 4, Fig. S2.3; Table 1).
Greater losses can be expected when the entire study region is considered. In contrast, if dispersal occurs, substantial range gains are predicted outside PAs, particularly under the best scenario, and a slight increase is expected for the entire study region (Table 1).

| Pan troglodytes schweinfurthii (eastern chimpanzee)
On average, RF and Maxent models performed best in explaining eastern chimpanzee distribution (Fig. S2.1a). Annual and maximum temperatures were important predictors in most individual and ensemble models, except for RF models, where annual precipitation, seasonal variation in precipitation, primary land and secondary land performed best (Figures 1b, 2, Fig. S2.1b; Table S1.2). Eastern chimpanzees encounter suitable conditions where the annual temperature is low (10-20°C), and maximum temperature is between 15 and 30 ⁰C (Fig. S1.4). Under the worst scenario, both temperature variables are predicted to increase.
Only one fifth of the study region of eastern chimpanzees is in PAs. According to both future scenarios, one third of the taxon's range is expected to be lost outside PAs if no dispersal occurs, with most range predicted to be lost if the entire study region is considered (Figure 4, Fig. S2.3; Table 1). In contrast, if dispersal occurs, a slight gain is expected outside PAs under both future scenarios, with a slight range expansion into new suitable areas expected for the entire study region (Figure 4, Fig. S2.3; Table 1). None of the predicted losses are expected where their range overlaps with that of Grauer's gorilla, but all gains are expected where both ranges overlap ( Figure 5, Table 1).

| Pan troglodytes troglodytes (central chimpanzee)
On average, RF, ANN and Maxent models performed best in predicting central chimpanzee distribution (Fig. S2.1a). Secondary land and seasonal variation of precipitation were the predictors of greatest importance in individual and ensemble models, except for GAM and Maxent models, where cropland and human population were slightly better predictors (Figures 1b, 2, Fig. S2.1b; Table   S1.2). The study region of central chimpanzees is characterized by a relatively high percentage of secondary land (>40%), a human population density between 5 and 15 people km −2 , seasonal variation of precipitation between 30 and 80 mm and low percentage of cropland (<15%) (Fig. S1.4). The best scenario predicts secondary land expansion, an increase in human population and a reduction in cropland area. Large increases in all variables are predicted under the worst scenario.
As was found for eastern chimpanzees, only one fifth of the study region of central chimpanzees is covered by PAs. A reduction of three quarters of range is expected outside PAs under both future scenarios if no dispersal occurs, with most range expected to be lost if the entire study region is included (Figure 4, Fig. S2.3; Table 1).
Predictions of range gains for central chimpanzees suggest that substantial suitable habitat will become available outside PAs under the best scenario, with a slight increase if assuming the entire study region (Table 1). Under both future scenarios, most losses were predicted outside PAs in the same geographic areas where losses for western lowland gorillas are also expected, with a slight increase expected for the entire study region ( Figure 5, Table 1). The same trend was predicted for gains (Table 1). *range overlap between gbg and pts; **range overlap between ggg and ptt TA B L E 1 Results of the predicted change (%) in African ape ranges, assuming either no dispersal (loss) and dispersal (gain), for areas outside protected areas (PAs) only (assuming complete management effectiveness of PAs) and for the entire study region, by 2050 under the best-and the worst-case scenario

| Pan troglodytes verus (western chimpanzee)
On average, RF and ANN models performed best in predicting the distribution of western chimpanzees (Fig. S2.1a). Annual precipitation, primary land and human population were the most important variables in individual and ensemble models (Figures 1b, 2,   Fig. S2.1b; Table S1.2). Current environmental conditions found in the study region of western chimpanzees are annual precipitation below 2,000 mm, a very high presence of human population (>10 people km −2 ) and high primary land cover (>60%) (Fig. S1.4) Table 1). On the other hand, range gains are mostly anticipated outside PAs, particularly under the best scenario, if there is dispersal (Table 1).  (Schloss et al., 2012), migration lag and ecological constraints. The 30-year time frame considered here represents a bit more than an ape generation (Kuehl et al., 2017;Plumptre et al., 2016) and it is unlikely that migration into new areas during this time occurs to any greater extent. It is therefore very important that these results are not interpreted as indicating that range gain will definitely occur as effective protection of new suitable areas will need to be ensured for a great ape population to shift to such habitat. Importantly, massive range loss can be anticipated in the next 30 years given the 2%-7% of annual population decline previously estimated for great apes (Kuehl et al., 2017;Plumptre et al., 2016;Strindberg et al., 2018;Wich et al., 2016). correlative and mechanistic models to address the future effects of F I G U R E 5 Predicted change in range overlap where areas of loss, gain, stability and absence were the same between (a) and (b) Gorilla beringei graueri (gbg) and Pan troglodytes schweinfurthii (pts), and between (c) and (d) Gorilla gorilla gorilla (ggg) and Pan troglodytes troglodytes (ptt), by 2050 under the best-and the worstcase scenario. The category "overlap" represents areas where taxa overlap but future range conditions are unlikely to be the same between them. The black lines represent protected areas climate change on mountain gorillas and found that correlative models predict more dramatic changes in range change than mechanistic models. In contrast, mechanistic models revealed that, if plant productivity remains unaffected, gorillas will be able to adapt to warming temperatures (Thorne et al., 2013). Lehmann et al. (2010) employed a physiological/behavioural mechanistic approach to investigate how climate change under a worst scenario would influence African ape survival and reported that chimpanzees might lose 10% of current range and gorillas 75%

| D ISCUSS I ON
by the year 2100. Our study concurs with these results for most gorilla taxa, but more than three quarters of chimpanzee range is predicted to be lost under both future scenarios if no dispersal occurs. However, our full dispersal scenario predicts range gains in new areas under both future scenarios for most taxa, in agreement with Lehmann's study. The different results regarding range loss predicted for chimpanzees can be explained as follows: (1) community size rather than animal density was considered in Lehmann's model and the minimum party size requirements were not allowed to vary in response to vegetation cover, otherwise losses up to 39% may be expected for chimpanzees; (2) the effect of climate change on great ape distribution is critically dependent on the minimum viable group size that apes require for survival, and a conservative value for minimum viable group sizes (i.e. 10 individuals) based on minimum observed group sizes (which is not exactly the same) was used in Lehmann et al. (2010). In contrast, if a greater number of individuals had been considered (i.e. 45 individuals, Lehmann et al., 2007) up to 50% of range loss would be expected for chimpanzees; (3) by including not only climatic variables as in Lehmann's study, but also relevant human-related variables known to have a strong effect on current great ape distribution, our models provide a better understanding of the combined effects of these global change drivers and imply that negative effects of climate change on African apes can be reduced through appropriate land use planning and management.
Nevertheless, both approaches agree that future climate is predicted to dramatically change across African ape ranges.
Among the important climatic variables in determining resource availability and species distributions, and consequently, their effects on African ape time budgets, gorillas and chimpanzees are more sensitive to variations in temperature than in precipitation and they persist better in habitats with lower monthly temperature variation (Lehmann et al., 2010). Moreover, gorillas are predicted to be affected more than chimpanzees given the more restricted behavioural flexibility of gorillas to cope with temperature variation (Lehmann et al., 2010). Our study suggests that annual and maximum temperatures influence the distribution of most gorillas, bonobos and eastern chimpanzees, but not that of the other three subspecies of chimpanzee. Additionally, annual precipitation, and particularly its distribution over the wet and dry seasons, affects the distribution of most gorillas and most chimpanzees. These findings were also identified in shaping the distribution of Nigeria-Cameroon and central chimpanzees (Clee et al., 2015) and are indirect evidence of the marked influence of both temperature and precipitation on species niche with regard to dehydration and thermoregulation (Wessling et al., 2018).  investigated the representativeness of oil palm concessions within African ape ranges and found that more than half of the oil palm concessions are located within their current ranges. Moreover, potential future oil palm development is widely expected across their ranges, particularly in land outside PAs (Wich, Garcia-Ulloa, Kühl, et al., 2014). Thus, we can expect climate change to exacerbate range loss for African apes and consequently pose serious threats to species persistence, as they are anticipated to impact orangutans (Struebig et al., 2015). By integrating future climate and land use changes as well as human population scenarios, our predictions provide strong evidence for synergistic interactions among these global drivers constraining species distributions. We suggest that future studies assess how much of the new predicted suitable areas will African apes be able to colonize by considering a mechanistic approach that integrates population dynamics (as in Lehmann et al. 2010) and dispersal abilities.
Using two future scenarios, two dispersal scenarios, an ensemble forecasting approach and including only a few but highly important predictors of the distribution of African apes, should have addressed potential sources of uncertainty in our distribution models (Brun et al., 2019;Thorne et al., 2013). A recent study proposed that SDMs include historical records to produce better predictions of range shifts rather than relying on contemporary records alone (Faurby & Araújo, 2018). This is important for large vertebrates given the direct effects of anthropogenic disturbances on their distribution, and many ranges being far from equilibrium under current environmental conditions (Faurby & Araújo, 2018 To simplify the interpretation of results and for a better comparison between taxa, we decided to keep the default settings of the algorithms, which may have resulted in uncertainty regarding model performance and resulting species distribution due to the differential sensitivity of each algorithm to species modelled, sampling effort and evaluation metric (Hallgren et al., 2019). We recommend investigating the sensitivity of the algorithms to their configuration settings on the resulting projected species distribution (Hallgren et al., 2019). In addition, we relied on pseudo-absences instead of true absences and considered a greater number of randomly selected pseudo-absence points for those taxa with a larger study region, which may have introduced uncertainty given the differential sensitivity of each algorithm to the number of pseudo-absences (Barbet-Massin et al., 2012;Guillera-Arroita et al., 2015;Phillips et al., 2009). We thus recommend selecting the number of pseudoabsences and choosing the method to generate them depending on each algorithm (Barbet-Massin et al., 2012) and using true absences whenever available.
Model performance was evaluated through cross-validation using a bootstrap approach, which is a common procedure in SDMs (Thuiller et al., 2009). However, it can be problematic and lead to over-fitting models and inflated performance metrics, particularly in the presence of spatial autocorrelation (Telford & Birks, 2009;Wenger & Olden, 2012). Despite rarefying species points and creating sampling bias layers from these points, we could not reduce spatial autocorrelation completely and potentially it inflated the performance metric (i.e. TSS). To better deal with the effect of spatial autocorrelation, it would be useful if non-traditional crossvalidation schemes such as block spatial CV and h-block CV would be implemented in future versions of the biomod2 package (Telford & Birks, 2009;Wenger & Olden, 2012).
Model performance also depends on spatial resolution of the environmental variables, particularly when variables are originally available at coarse resolutions (Wenger & Olden, 2012). The original resolution of the spatial layers of land use (at 50 km) and human population changes (at 15 km) potentially influenced the performance of our models, particularly in those species with narrow distributions such as mountain gorillas (Thorne et al., 2013) and Cross River gorillas, given the extreme trends in range loss predicted (Wenger & Olden, 2012). Therefore, it would be important to consider future scenarios for climate, land use and human population changes at finer scales once such spatial layers become available. Additionally, mechanistic approaches may be more appropriate for species with narrow distributions as large differences in range change can be obtained with different modelling approaches (Thorne et al., 2013).
Mining concessions and granted mining claims are increasing dramatically across Africa, particularly threatening large ape populations in Guinea, Gabon and Liberia (Howard, 2019). It will be important to model the influence of this threat on future African ape distributions once appropriate spatial datasets become available.
Nevertheless, despite these potential sources of uncertainty in our distribution models, the predicted ranges under current conditions in this study (Fig. S2.2) are well in line with results from previous studies, importantly also for those taxa for which model algorithms had "poor" or "useful" performance (e.g. G. g. gorilla and P. t. troglodytes (Strindberg et al., 2018); G. b. graueri (Plumptre et al., 2016); and all African apes (Junker et al., 2012)). Strindberg et al. (2018) found that western lowland gorillas and central chimpanzees, two sympatric taxa with 97% range overlap, occur mostly outside PAs, and argued for "reinforcement of anti- analysis to optimize conservation land use planning and identify priority areas for these species (Freeman et al., 2018;Jones et al., 2018). This is extremely important given that many African

| Conservation implications
PAs in ape ranges are separated from each other -although there is often transboundary connectivity (Santini et al., 2016). This is also of particular concern because great apes occur mostly outside PAs and have a low dispersal capacity due to their small population sizes, low population densities, dietary requirements and poor thermoregulation. It will be important to ensure objective assessments of human pressures and habitat conditions in potential PAs to avert species extinctions in the long term (Jones et al., 2018).
Taxon-specific frameworks of environmental and socioeconomic trends (Estrada et al., 2018;Strindberg et al., 2018;Tranquilli et al., 2014) should be considered at all major decisionmaking levels in range countries and abroad to (1)  where an NGO or other organization takes on management responsibility for a given site over one or more decades (Scholte et al., 2018). should also be used to guide the prioritization of conservation efforts for these flagship species to avoid irreversible losses. Our study goes beyond previous work, which focused only on modelling the effects of climate change alone. Here, we also consider land use and human population changes, which advances our understanding of the joint effects of these key drivers on African ape distribution. Importantly, our findings suggest that some of the negative effects of climate change on African apes can be mitigated if appropriate land use planning and management action is taken. Given that all great ape taxa will find most suitable areas outside PAs, and that the existing network of PAs is inadequate for ensuring the long-term conservation of African apes (Strindberg et al., 2018), we support the argument that effective conservation strategies require taxon-specific conservation planning that focuses on existing and proposed PAs, the creation and/or management of which can be informed by our habitat suitability models.

| CON CLUS IONS
Additionally, efforts to maintain connectivity between the habitats predicted to be suitable in the future will be crucial for the survival of African apes. As an example, a country-wide approach has been undertaken in Gabon, where planning for the development of agriculture, road and rail links, and mineral extraction has been informed by wildlife and vegetation data in order to locate these activities in areas that are already degraded, and to avoid closed-canopy oldgrowth and remote forests (Government of Gabon, 2012;Strindberg et al., 2018). This will be an effective way of promoting habitat connectivity to maintain African ape populations and sympatric wildlife.

ACK N OWLED G EM ENTS
We are grateful to the governments and national authorities for research permissions, and field staff for logistical support and guidance during data collection. We would like to thank all organizations for Open Access funding enabled and organized by Projekt DEAL.

CO N FLI C T S O F I NTE R E S T
The authors declare there is no conflict of interest.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13358.

DATA AVA I L A B I L I T Y S TAT E M E N T S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section.