Quantitative estimates of glacial refugia for chimpanzees (Pan troglodytes) since the Last Interglacial (120,000 BP)

Paleoclimate reconstructions have enhanced our understanding of how past climates have shaped present‐day biodiversity. We hypothesize that the geographic extent of Pleistocene forest refugia and suitable habitat fluctuated significantly in time during the late Quaternary for chimpanzees (Pan troglodytes). Using bioclimatic variables representing monthly temperature and precipitation estimates, past human population density data, and an extensive database of georeferenced presence points, we built a model of changing habitat suitability for chimpanzees at fine spatio‐temporal scales dating back to the Last Interglacial (120,000 BP). Our models cover a spatial resolution of 0.0467° (approximately 5.19 km2 grid cells) and a temporal resolution of between 1000 and 4000 years. Using our model, we mapped habitat stability over time using three approaches, comparing our modeled stability estimates to existing knowledge of Afrotropical refugia, as well as contemporary patterns of major keystone tropical food resources used by chimpanzees, figs (Moraceae), and palms (Arecacae). Results show habitat stability congruent with known glacial refugia across Africa, suggesting their extents may have been underestimated for chimpanzees, with potentially up to approximately 60,000 km2 of previously unrecognized glacial refugia. The refugia we highlight coincide with higher species richness for figs and palms. Our results provide spatio‐temporally explicit insights into the role of refugia across the chimpanzee range, forming the empirical foundation for developing and testing hypotheses about behavioral, ecological, and genetic diversity with additional data. This methodology can be applied to other species and geographic areas when sufficient data are available.


Abstract
Paleoclimate reconstructions have enhanced our understanding of how past climates have shaped present-day biodiversity. We hypothesize that the geographic extent of Pleistocene forest refugia and suitable habitat fluctuated significantly in time during the late Quaternary for chimpanzees (Pan troglodytes). Using bioclimatic variables representing monthly temperature and precipitation estimates, past human population density data, and an extensive database of georeferenced presence points, we built a model of changing habitat suitability for chimpanzees at fine spatio-temporal scales dating back to the Last Interglacial (120,000 BP). Our models cover a spatial resolution of 0.0467°(approximately 5.19 km 2 grid cells) and a temporal resolution of between 1000 and 4000 years. Using our model, we mapped habitat stability over time using three approaches, comparing our modeled stability estimates to existing knowledge of Afrotropical refugia, as well as contemporary patterns of major keystone tropical food resources used by chimpanzees, figs (Moraceae), and palms (Arecacae). Results show habitat stability congruent with known glacial refugia across Africa, suggesting their extents may have been underestimated for chimpanzees, with potentially up to approximately 60,000 km 2 of previously unrecognized glacial refugia. The refugia we highlight coincide with higher species richness for figs and palms. Our results provide spatio-temporally explicit insights into the role of refugia across the chimpanzee range, forming the empirical foundation for developing and testing hypotheses about behavioral, ecological, and genetic diversity with additional data. This methodology can be applied to other species and geographic areas when sufficient data are available.

K E Y W O R D S
Africa, diversification, ensemble, paleoclimate, species distribution modeling 1 | INTRODUCTION Paleoclimate reconstructions have greatly enhanced our understanding of how biodiversity has been shaped globally since the late Quaternary, including inferences on range shifts, extinctions, and the evolution of distinct lineages (Hewitt, 2000(Hewitt, , 2004Sandel et al., 2011;Svenning et al., 2015). For example, climatic variability through time has played a major role in determining the uneven distribution of biodiversity across the Afrotropics (Demenocal, 1995;Maslin et al., 2014;Sepulchre et al., 2006;Trauth et al., 2005). Although a high proportion of African biodiversity is concentrated in a handful of forested hotspots and centers of endemism (Kingdon, 1990;Myers et al., 2000), these areas are geographically diffuse and historically complex, with idiosyncratic characteristics including heterogeneous topography, hydrological features, and highly dynamic and unique forest histories. Many Afrotropical forest lineages across different taxonomic groups are assumed to have tracked available habitat as it shifted throughout the Quaternary, contracting into refugia during glaciation, and often expanding during postglacial periods to colonize (or recolonize) new geographic areas. However, spatio-temporal reconstructions of Pleistocene forest refugia across the Afrotropics are often limited to broad-scale maps without detailed information at local scales (Maley, 1996;Mayr & O'Hara, 1986). A more detailed quantification of the spatio-temporal distribution of forest refugia has been hampered by the lack of high-resolution paleoecological data with most previous reconstructions being limited to coarse spatial grains (e.g., 2.5 × 3.75°, Singarayer & Valdes, 2010) or to only a handful of temporal snapshots during the Quaternary (e.g., last interglacial, last glacial maximum, and late holocene; Hijmans et al., 2005). Here, we address this limitation by modeling the historical habitat suitability of chimpanzees (Pan troglodytes) across the Afrotropics using a comprehensive database of recently collected georeferenced presence points spanning their entire range, paleoclimate reconstructions, and human density data since the Last Interglacial (120,000 BP). Chimpanzees are an appropriate focal species for our approach due to their wide geographical distribution and high detectability. This means that their presence in an area is unlikely to be overlooked even in places that are only sporadically surveyed, as is the case for many poorly known and remote field sites across the Afrotropics. Chimpanzees occur across 21 African countries, from dense moist forest to arid savannah, from sea level up to an altitude of 2800 m and varying levels of anthropogenic pressure (Humle et al., 2016). Furthermore, their wide distribution encompasses high levels of ecological, behavioral, and genetic diversity amongst populations, providing a model system to study the availability of glacial refugia over time across diverse habitats and geographic regions. Four subspecies are currently described: P. troglodytes verus from West Africa, P. troglodytes ellioti from Nigeria and Cameroon, P. troglodytes troglodytes from Central Africa, and P. troglodytes schweinfurthii from Central and Eastern Africa. Chimpanzee diversity is also apparent in their social organization (community size ranges from 12 to 200, Langergraber et al., 2017), feeding ecology, and diet (proportions and species of vertebrates, insects and fruits consumed varies among populations (Basabose, 2002;McGrew et al., 1988;Nishida & Uehara, 1983;Wrangham, 1977), and range sizes (from around 3 km 2 to over 70 km 2 (Herbinger et al., 2001;Hunt & McGrew, 2002;Pruetz & Bertolani, 2009;Wessling et al., 2018). Behavioral complexity in chimpanzees is highly variable, with diverse repertoires of tool use (e.g. sticks, stones and leaves to access insects, honey, meat, seeds and algae, and water filtering) evident across forest and savannah populations (Galat-Luong et al., 2009;Kühl et al., 2019;Whiten et al., 1999). Additionally, genetic differentiation is also profound across the range, with four currently recognized subspecies and markedly variable population substructure and genetic diversity within each of these (De Manuel et al., 2016;Fünfstück et al., 2015;Lester et al., 2021;Mitchell et al., 2015;Prado-Martinez et al., 2013).
It is suspected that historical climatic changes and associated population size changes and migrations have been a major driver of chimpanzee genetic and cultural differentiation since at least the Last Interglacial (Prado-Martinez et al., 2013). To date, our understanding of the spatio-temporal climatic changes which could have influenced the mechanisms driving chimpanzee ecological, behavioral, and genetic diversity has been limited by the absence of information about historical climate change at fine spatial scales across the Afrotropics, a shortfall which this study aims to address. Specifically, we hypothesize that previously identified glacial refugia played a major role in the habitat suitability and persistence of chimpanzees (P. troglodytes) during late Quaternary climate fluctuations. Beyond chimpanzees, the data set created in this study could be an important tool used to study other taxonomic groups and biological communities within Afrotropical forest ecosystems.

| METHODS
The research described herein complied with protocols and laws approved by the relevant authorities in each country where data was collected.

| General workflow
To summarize our general workflow (Figure 1), we modeled contemporary chimpanzee habitat suitability using a database of presence points and six noncorrelated climatic and anthropogenic variables that could be expected to influence their distribution. As separate species distribution modeling algorithms may provide biased results because of their sensitivity to the data and the mathematical functions utilized, we built ensemble models which combine multiple replicates of several different modeling algorithms together (Araújo & New, 2007). These models based on contemporary chimpanzee habitat suitability were then projected into the past at paleoclimate timescales (back to the Last Interglacial, 120,000 years BP) to assess how habitat suitability has changed through time. We then calculated various estimates of stability per pixel to represent the relative variation in habitat suitability over time, identifying areas of long-term habitat suitability.
We provide methodological details for each of the steps in our workflow in the following subsections, with extended technical information for Reconstructing African paleoclimatic conditions back to the Last Interglacial and Species distribution modeling of chimpanzees made provided in Online Supplementary Material (Text S1). Unless otherwise specified, all analyses were implemented in R (version 4.1.0, R Core Team, 2019).

| Reconstructing African paleoclimatic conditions back to the Last Interglacial
We obtained data on paleoclimatic conditions in Africa since the Last Interglacial from Bell et al. (2017). These data are based on the Hadley Centre Coupled Model version 3 general circulation model (HadCM3), a climate change prediction model routinely utilized by the Intergovernmental Panel on Climate Change (IPCC, 2007(IPCC, , 2014. These data represent 62 snapshots of global climatic conditions at 1000 year intervals from the present back to the Last Glacial Maximum (LGM,~22,000 BP), at 2000 year intervals from the LGM to 80,000 BP, and at 4000 year intervals from 80,000 BP back to the Last Interglacial (120,000 BP). These data comprised of eight climate variables: mean annual temperature (bioclim_01), temperature seasonality (bioclim_04), mean temperature of the warmest (bioclim_10) and coldest quarters (bioclim_11), mean annual precipitation (bio-clim_12), precipitation seasonality (bioclim_15), and precipitation of the wettest (bioclim_16) and driest quarters (bioclim_17) for each snapshot. We set the geographic extent of our models to encompass both the contemporary and potential historical habitat suitability of chimpanzees (McBrearty & Jablonski, 2005), from −18°to 32°l ongitude and −10°to 16°latitude. This extent includes the current distribution range of bonobos (P. paniscus), the sister species of chimpanzees, spatially separated from the latter by the Congo River.
We included the bonobo distribution in our models due to previous detection of historic admixture between chimpanzees and bonobos based on genomic data (De Manuel et al., 2016;Kuhlwilm et al., 2019), suggesting that chimpanzees may have occasionally been found in these areas.
To account for anthropogenic effects on chimpanzee habitat suitability, we also included a spatial layer available from the HyDE database (Klein Goldewijk et al., 2011), which provides information on modeled past human population density based on population and agricultural data. The calculation of this spatial layer is based on historical population, cropland, and pasture statistics combined with satellite information and specific allocation algorithms (which change over time) to create spatially explicit maps covering the last 12,000 years BP. We linearly extrapolated the human population densities back in time from 12,000 BP (i.e., the last available data from the HyDe database) uniformly to 120,000 BP and downscaled it so that it matched the temporal and spatial scale of our paleoclimate reconstructions. We selected these climatic and anthropogenic variables to reflect biologically informative conditions likely to have influenced chimpanzee habitat suitability, and because these variables were available at paleoclimatic timescales. However, we acknowledge that additional variables may be important to influence chimpanzee distributions (see Caveats heading of Discussion).

| Species distribution modeling of chimpanzees
We compiled a database of 139,902 georeferenced chimpanzee presence point records (nests, sightings, feces, and footprints) from our own fieldwork, from the IUCN SSC A.P.E.S. database, a collaborative initiative to centralize great ape population surveillance data (http://apes.eva.mpg.de/) and from the Pan African Programme: The Cultured Chimpanzee (http://panafrican.eva.mpg.de/). For each occurrence record, subspecies were determined by their geographical location. In total, this data, collected between 1984 and 2006, represent over fifty collaborating wildlife research institutions across the entire distributional range of chimpanzees.
We followed recommended guidelines for constructing SDMs (Araujo et al., 2019;Merow et al., 2013, Merow et al., 2014 to minimize bias in our models. This included spatially rarefying the presence data so that occurrences were not overclustered together (minimum distance of 10 km between each point), resulting in 1677 unique presence points remaining that were used to build SDMs. We selected background points (pseudo-absences) from a 0.5°buffer radius around presence points, to emphasize factors locally relevant in distinguishing suitable from unsuitable habitat, while adequately sampling the range of climatic conditions for chimpanzees (VanDerWal et al., 2009). Furthermore, to minimize spatial autocorrelation in our models we reduced our nine variables to six due to high co-correlations (Pearson's r > 0.6). Our retained variables were: annual mean temperature, temperature seasonality, annual precipitation, precipitation of wettest quarter, precipitation of driest quarter, and human density.
As individual modeling algorithm approaches may introduce additional biases, we built ensemble models combining multiple replicates of several different modeling algorithms to represent alternate possible states of the system being modeled (Araújo & New, 2007). Due to their combined power, ensemble models are widely accepted to provide more accurate results than single models F I G U R E 1 General workflow for our modeling approach, using georeferenced presence points, paleoclimatic reconstructions, and human density data to the Last Interglacial period as input data. Ensemble Species Distribution models for chimpanzees are built and evaluated, projected overall temporal snapshots (n = 62, until 120,000 years ago) and then stability over time for each pixel is calculate (three metrics, dynamic, static stability, and co-efficient of variation). These three output stability maps are then quantitatively compared with assumed forest refugia (Maley, 1996), as well as species richness of major keystone food resources for primates (figs, Kissling et al., 2007) and palms (Blach-Overgaard et al., 2013) (Forester et al., 2013). We constructed an ensemble species distribution model (SDM) representing the contemporary habitat suitability of chimpanzees (i.e., all combined subspecies) using the "sdm" R package version 1.0-46 (Naimi & Araújo, 2016). We used five crossvalidated replicates each of 14 different modeling algorithms available in the sdm R package, evaluating the performance of each modeling algorithm based on the Area Under the Curve of a Receiver Operating Characteristics plot (AUC, Fielding & Bell, 1997) and True Skill Statistic (TSS, Allouche et al., 2006). We retained only modeling algorithms that performed adequately (AUC > 0.8, and TSS > 0.5, Bell et al., 2017) in an ensemble model prediction using five crossvalidated replicates of each algorithm, weighting the contribution of each modeling algorithm to the ensemble model by their AUC. Given that we do not expect the variables to be equally important in influencing chimpanzee presence, for each modeling algorithm and the ensemble model we also measured the permutation importance of predictor variables across each model iteration using the "getVarImp" function in the "sdm" R package.

| Paleoclimatic modeling and habitat stability
To identify geographic areas where high habitat suitability for chimpanzees has remained stable since the Last Interglacial (i.e., to identify refugia where the effects of climate change have been less pronounced than surrounding areas), we projected our SDM back in time to obtain a habitat suitability estimate for each pixel at each paleoclimate snapshot. We projected our contemporary ensemble model onto the 62 snapshots of downscaled paleoclimate reconstructions and human density data since the last interglacial using the "predict" function of the "sdm" R package. In doing so, we assumed that the realized ecological niche of chimpanzees represented by our predictor variables remained relatively constant over this period. Though there is some controversy with this assumption, especially for species with narrow ecological niches (see Veloz et al., 2012), the large geographical extent and wide range of climatic conditions that our spatial and climatic data encompass (i.e., most of the Afrotropics), and the ability of our study species to track climate change minimize our concerns for a geographically widespread species such as the chimpanzee. However, we acknowledge there may be exceptions at local scales over such long temporal intervals where ecological niches have shifted over time, for example in subspecies.
Based on these per pixel estimates of habitat suitability over time, we then evaluated suitability changes through space and time using three approaches: static stability (Barratt et al., 2017, Hugall et al., 2002, and the coefficient of variation % (CV), which do not consider dispersal between pixels, and dynamic stability (Graham et al., 2010;Rosauer et al., 2015), which account for how chimpanzees may have tracked suitable climate conditions across time.
To calculate the CV stability in each pixel across all of the temporal snapshots of habitat suitability, we created a rasterstack of model outputs for each temporal snapshot and used the cv function of the "raster" package version 3.1-5 (Hijmans et al., 2011). For the static stability estimate, we summed the negative log of suitability through time for each pixel and took the exponent of this value to give a value between 0 and 1 to represent the degree to which the pixel has continuously provided suitable habitat. Similarly, dynamic stability reads the suitability at each time period but permits dispersal across time periods (i.e. paleoclimate time snapshots) using a cost distance (a function of habitat suitability), based on a user-defined maximum dispersal distance. To calculate dynamic stability per pixel we allowed a conservative 5 m/year dispersal rate, which amounts to a maximum dispersal distance of 600 km over the 120,000 year time period for chimpanzee populations since the last interglacial, approximately matching the historical distribution inferred by McBrearty and Jablonski (2005) shown in Figure 2. Using the dynamic stability approach, a pixel is considered as stable as long as pixels within the defined dispersal distance have suitable climate in adjacent time steps (i.e., so an area of suitable habitat for chimpanzees does not necessarily need to be consistently stable if they can disperse there from adjacent pixels to track climate change). All static and dynamic stability surfaces were calculated using modified R scripts from Graham et al. (2010) provided by Jeremy VanDerWal.

| Validating our inferred refugia using additional data
To supplement the qualitative assessment of stability maps, we quantitatively compared our modeled areas of long-term habitat stability (i.e., refugia) against the most detailed available estimates of Afrotropical refugia by Maley (1996), based on a combination of paleopalynological (fossil pollen) data and species richness and endemism patterns. Though Maley's refugia are widely used, we feel that they may not accurately reflect microrefugia at smaller spatial scales due to a lack of comprehensive surveying effort across the Afrotropics. To do this we first defined the threshold for pixels in our model outputs to be considered as refugia as having a dynamic stability value of 0.97 over the 62 temporal snapshots covering 120,000 years. We then calculated the number of pixels present in each of Maley's refugia and expressed the proportion of each of then known refugia that were captured by our analyses as a percentage. Moreover, we report additional areas inferred as refugia which are not inferred by Maley (1996) in terms of their size and location.
As an additional validation of our putative refugia, we tested our obtained habitat stability estimates for chimpanzees against species To this end we utilized data from Onstein et al. (2020), which incorporated distribution data for figs from Kissling et al. (2007), and for palms from Blach-Overgaard et al. (2013) . Fig species richness was calculated from digitized range map GIS shapefiles which were converted to binary presence-absence raster maps and then summed in R to give species richness in each pixel. Palm species richness was calculated in a similar way, but directly from the species distribution models, so instead of binary presence-absence per pixel the values represent probabilistic values of a species occurring in that pixel.
Both the fig and palm data were clipped to the exact same size as our study region so that direct comparisons could be made with our stability estimates. To assess correlations between the data sets, we stacked the raster layers (dynamic stability, static stability, CV, fig species richness, palm species richness in R using the "raster" package version 3.1-5 and conducted Pearson's correlation tests between them. We visualized these correlations usin the "corrplot" version 0.84 R package (Wei & Simko, 2017).

| Sensitivity analyses
Our final contemporary chimpanzee habitat suitability model ( Figure 3) represents all four currently recognized subspecies combined (Humle et al., 2016), but as the four subspecies diverged from one another within the last one million years, we also repeated the SDM, paleoclimate projections and stability estimations described above for each subspecies separately to ensure that our model shown in Figure 3 did not underpredict habitat suitability at local scales that may be ecologically important for each subspecies. The number of presence points (and background points) used for the subspecies, P. troglodytes verus, P. troglodytes ellioti, P. troglodytes troglodytes, and P. troglodytes schweinfurthii were 519, 134, 663, and 451, respectively). As a further sensitivity analysis, we also repeated all analyses (i.e., combined subspecies together and for each subspecies separately) using presence data spatially rarefied so that no points were within 25 km of each other to assess whether spatial bias in the presence data was influencing the output models (number of presence points and background points for combined subspecies and for P. troglodytes verus, P. troglodytes ellioti, P. troglodytes troglodytes, and P. troglodytes schweinfurthii separately were 658, 225, 57, 212, and 164, respectively) ( Figure S1). To statistically test how well our contemporary habitat suitability models and stability models (static stability, dynamic stability, and CV stability) matched one another across sensitivity analyses (with different spatial rarefaction in the input presence points) we conducted Pearson's correlation tests in R using the "raster" package version 3.1-5, comparing both matching sets of output rasters. We also visually compared our model outputs against other published contemporary habitat suitability estimates available for chimpanzees which are available at different spatial resolutions and extents than our models (Heinicke et al., 2019;Jantz et al., 2016;Junker et al., 2012;Sesink Clee et al., 2015;Strindberg et al., 2018).

| Contemporary chimpanzee habitat suitability
Our contemporary habitat suitability model performed generally well at the species level (i.e., all four subspecies combined), with AUC > F I G U R E 2 Current chimpanzee distribution and presence points from our own fieldwork, from the IUCN SSC A.P.E.S. database (http://apes. eva.mpg.de/) and from the Pan African Programme: The Cultured Chimpanzee (http://panafrican.eva.mpg.de/). Currently recognized subspecies ranges are depicted by colored polygons (Humle et al., 2016), with black dots representing sampling localities (n = 139,902 georeferenced presence points) used for species distribution modeling, historical distribution range (ca. 500kya, McBrearty & Jablonski, 2005) represented by diagonal shading. Country borders, topography (lowlands = green, highlands = brown), and major hydrological features are also shown 0.8 and TSS > 0.5 for 12 of the 14 tested modeling algorithms (Table 1). This was similar for sensitivity analyses on subspecies models P. troglodytes schweinfurthii (n = 13), P. troglodytes troglodytes (n = 13), and P. troglodytes verus (n = 11) (Table S1). Bioclim and bioclim.dismo algorithms performed poorly for most models (i.e., low AUC < 0.8, TSS < 0.5) and were therefore excluded from the ensemble models. For P. troglodytes elliotti, all 14 modeling algorithms passed our AUC and TSS thresholds and were retained in the ensemble. The best performing modeling algorithms (i.e., higher AUC and TSS) were most often Maxent, Generalized Additive Models, and Random Forest (Tables 1 and S1). We believe these differences between model performance to be due to the nature of the input data, where some model assumptions simply match the characteristics of the predictors, presences, and background better (for a detailed overview, see Araújo et al., 2019). These results were also consistent when performing the sensitivity analyses using only presence points rarefied to a minimum distance of 25 km, indicating that our results are robust against spatial bias in the input data (Table S1), albeit with generally lower AUC and TSS, most probably due to the reduced number of presence points used to train the models.
The contemporary habitat suitability model for the full species ( Figure 3a) showed more extensive areas of high suitability compared with sensitivity analyses when subspecies were modeled separately, and suitability was also lower using the 25 km rarefied input presence points ( Figure S2 and Table S2), most probably in both cases due to the lower sample sizes (presences) to capture the true ecological niche. Our contemporary habitat suitability models for chimpanzees ( Figure 3a) approximately matched those previously published (Heinicke et al., 2019;Jantz et al., 2016;Junker et al., 2012;Sesink Clee et al., 2015;Strindberg et al., 2018) but there were several differences, likely due to the different variables used to build models that are unavailable or inappropriate at the timescales we projected (e.g., forest cover and change, distance to roads/rivers, and disease dynamics), and the different spatial resolution and extent of the models. For example, our models predicted much higher habitat suitability in west Africa for P. troglodytes verus across parts of Sierra Leone, Liberia, and Guinea than some previous work (Jantz et al., 2016;Junker et al., 2012), though they are similar to the modeled chimpanzee density patterns reported by Heinicke et al. Summarizing variable permutation importance indicated that precipitation-related predictors (bioclim_16precipitation of wettest quarter, bioclim_17precipitation of driest quarter, and in particular bioclim_12annual rainfall) were more important than temperaturerelated predictors (bioclim_01annual mean temperature and bioclim_04temperature seasonality) and human population density in explaining chimpanzee habitat suitability (Figure 3b,c). These results were consistent across all sensitivity analyses including models using the 25 km rarefied input presence points and for separate subspecies models except for P. troglodytes ellioti (bioclim_16 with the highest permutation importance rather than bioclim_12) and P. troglodytes schweinfurthii, which had similar variable importance between all temperature and precipitation predictor variables, and low importance for human density ( Figure S3).  (Maley, 1996;Mayr & O'Hara, 1986), with our data providing estimates at finer spatial scales than the highest resolution previous work available (i.e., the most detailed polygons, Maley, 1996). Dynamic stability and CV stability estimates showed slightly larger areas of suitability than static stability estimates (Figure 4, S5-S7, Table 2), though all were highly congruent, highlighting the same general regions as areas of long-term habitat stability. However, when thresholding models to define refugia (>0.97 stability for dynamic and static stability, less than 3% variation per pixel in CV stability), dynamic stability and CV stability estimates captured higher proportions of each of Maley's (1996) refugia than static stability ( Figure S8). Dynamic stability and CV stability estimates corroborated that refugia suggested in previous studies have remained more stable than surrounding areas, capturing major parts of several known refugia (numbered in  (Maley, 1996;Mayr & O'Hara, 1986). However our models did not capture parts of these refugia in some cases, and also failed to capture several others including those in West Africa (#3) and Central Africa (#7, #8, and #9), as well as some of the central African microrefugia described by Leal (2004), though the CV stability estimate did also capture a proportion of refuge #7.

| Habitat suitability through time and identifying potential refugia
Our results suggest that chimpanzees may have occupied stable habitat across a wider range than indicated by previously identified rainforest refugia for equatorial Africa (Figures 4a and S8; Table 3 Based on the other two stability metrics (static stability and CV stability), additional pixels of previously unrecognized refugia are fewer than dynamic stability estimates but still represent over approximately 12,500 km 2 and~55,000 km 2 , respectively. Habitat suitability and stability estimates across our sensitivity analyses using different levels of presence points spatial rarefaction (10 km vs. 25 km) were highly consistent (Pearson's r > 0.971 for all contemporary habitat suitability models and dynamic and static stability estimates, and Pearson's r > 0.677 for all CV stability estimates, Table S2, Figures S5-S7). Our dynamic and static habitat stability estimates for the 10 km rarefied data set (Figure 4d Figure 4b,c), with Pearson's r ranging between 0.535 and 0.862. Our F I G U R E 3 (a) Contemporary habitat suitability model for chimpanzees (ensemble model of 5 cross-validated replicates of 12 modeling algorithms) and predictor variable permutation importance b) averaged across all modeling algorithms, and c) per individual modeling algorithm. Axes represents predictor variable permutation importance (y axis), for each modeling algorithm (x axis), bioclim_01 = Mean annual temperature, bioclim_04 = temperature seasonality, bioclim_12 = mean annual precipitation, bioclim_16 = precipitation of the wettest quarter, bioclim_17 = precipitation of the driest quarter, human = human density. Country borders (gray lines) and chimpanzee subspecies ranges (green lines, Humle et al., 2016)

| DISCUSSION
We hypothesized that the geographic extent of Pleistocene forest refugia and suitable habitat fluctuated significantly in time during the late Quaternary for chimpanzees (P. troglodytes), testing this by building spatio-temporally explicit models of habitat suitability across Africa over the last 120,000 years, since the Last Interglacial. We found several areas of high habitat suitability that have remained stable over this time period, largely matching estimates of putative refugia from previous studies (Maley, 1996;Mayr & O'Hara, 1986), and areas of higher species richness for figs and palms, both keystone tropical plant families known to be important for primates generally (Onstein et al., 2020). Our refugia also approximately match patterns of higher effective diversity across the species range inferred from genetic data (microsatellites) by Lester et al. (2021). Our results suggest that for chimpanzees, most of these previous refugia are underestimated in size, implying that additional geographic areas (potentially representing over 60,000 km 2 in total) may have supported suitable habitats during glacial climate fluctuations. This is further supported by fig and palm species richness patterns, which are not solely restricted to known refugia (Maley, 1996) showing T A B L E 1 Performance of each individual modeling algorithm based on five cross-validated replicates for full species (Pan troglodytes) and each of the four currently recognized subspecies using presence points rarefied to minimum 10 km distance from one another

| Understanding historical and recent factors shaping contemporary biodiversity patterns
Many previous studies have shown that climate change throughout the late Quaternary had a profound influence on the contemporary distribution and structure of species and populations globally (Davis & Shaw, 2001;Ellegren & Galtier, 2016;Hewitt, 2000Hewitt, , 2004Sandel et al., 2011;Svenning et al., 2015). In Africa it has been hypothesized that forest refugia, often in mountainous regions, are important areas for buffering the effects of climate change on biodiversity over paleoclimatic scales, supporting rare species and unique genetic variants (Fjeldsaå & Lovett, 1997;Mayr & O'Hara, 1986). Although initially these refugia were identified using species distributional patterns for F I G U R E 4 Stability of chimpanzee habitat suitability over 62 snapshots of paleoclimate reconstructions representing the past 120,000 years using the full species data set. (a) Our refugia (in blue) are defined as all pixels with a dynamic stability value > 0.97 in habitat suitability maps over the 120,000 year time period, compared with refugia defined by Maley (1996) (in gray, numbered and referred to in Table 2 (c), more yellow colors represent higher species richness (previous esimates of forest refugia (Maley, 1996) drawn in red, and green, respectively. In maps (d)-(f), warmer colors represent areas of higher suitability over time that have remained stable compared with previous estimates of forest refugia (dotted black lines, Maley, 1996). Country borders (gray lines) and chimpanzee subspecies ranges (green lines, Humle et al., 2016) are also shown Mikula, Patzenhauerová, et al., 2014;Gaubert et al., 2016;Mizerovská et al., 2019;Nicolas et al., 2011), including other primates, (Anthony et al., 2014;Clifford et al., 2004;Gonder et al., 2011;Pozzi, 2016;Telfer et al., 2003), amphibians (Charles et al., 2018;Leaché et al., 2019;Portik et al., 2017), and plants (Faye et al., 2016;Hardy et al., 2013;Piñeiro et al., 2017Piñeiro et al., , 2019, albeit with some differences between species due to different ecological characteristics and idiosyncratic responses to climatic changes (Lowe et al., 2010).
Our results are concordant with most of the broadly defined glacial refugia for the Afrotropics (Maley, 1996;Mayr & O'Hara, 1986) but suggest that potentially stable areas of chimpanzee habitat may have been larger. This may well be due to the ecological and behavioral flexibility exhibited by chimpanzees Wessling et al., 2018;Whiten et al., 1999), with their ability to tolerate drier savannah-like conditions, for example at the edges of their current distribution range where forest interdigitates with savannah, resulting in a broad ecological niche. This flexibility may also account for our estimates of chimpanzee habitat suitability and stability over time being generally higher when modeling all four subspecies combined compared with modeling each subspecies separately, as the ecological niche of each subspecies is narrower than that of the whole species and there may well have been adaptations to local climatic conditions for each subspecies.
Since the turn of the 20th century, the higher deforestation rates around the edges of chimpanzee range (e.g., in West Africa) than those at the core of their range (Aleman et al., 2018), may have caused the local extirpation of some populations, and potentially stimulated ecological and behavioral diversification in others. Historical climate changes combined with these more recent anthropogenic effects may help to explain some of the complex patterns of behavioral and cultural diversity described in chimpanzees across different parts of Africa Whiten et al., 1999) and particularly why some behaviors are geographically disjunct (e.g., accumulative stone-throwing, Kühl et al. (2016), algae fishing (Boesch et al., 2017). However, these complex patterns of diversity are further complicated by multiple factors, including disease dynamics (Leroy et al., 2004), the bushmeat trade (Bennett et al., 2007), and local genetic adaptations (Schmidt et al., 2019).

| Caveats
Although our reconstructions of paleoclimatic refugia align with those of Maley (1996) and Mayr and Hara (1986), and with more recent molecular evidence in vertebrates and plants, we advocate Note: Each of the stability estimates are thresholded so that only highly stable areas over time are classified as refugia (dynamic and static stability threshold = 0.97, CV threshold = maximum 3% coefficient of variation). These areas are then quantified in terms of their size in pixels (each pixel represents~5 km 2 ), and compared with Maley's (1996) refugia expressed as the % of each refuge recovered by our models.
Abbreviation: CV, coefficient of variation %; CIV, ivory coast; DRC, democratic republic of the Congo.
T A B L E 3 Additional refugia identified by our models that are not captured by Maley (1996)

| Future directions for understanding diversification mechanisms
By modeling habitat suitability and stability over the last 120,000 years at high spatial and temporal resolution, we provide a foundation to investigate the diversification mechanisms underlying patterns of ecological, behavioral, and genetic diversity in chimpanzees. We suggest that future research should exploit our data set to generate testable hypotheses about chimpanzee population diversity by combining our spatio-temporally explicit models with ecological, behavioral, and molecular data. For example, habitat suitability fluc- Beyond correlative approaches such as those described above, mechanistic approaches would enable a deeper exploration of biogeographical patterns and processes that have affected chimpanzees, especially with the recent availability of comprehensive behavioral  and molecular data (Lester et al., 2021) for a large number of wild chimpanzee populations. Rapid developments in the generation and analysis of genome-wide molecular data over the past decade have revealed detailed demographic histories, enabling the identification of diversification mechanisms due to forest refugia, which are characterized by divergence, isolation, and secondary contact as refugial habitats fragment and reconnect with each other during glacial cycles (Barratt et al., 2018;Charles et al., 2018;Feng et al., in press;Leaché et al., 2019;Portik et al., 2017). The ability to distinguish signals of forest refugia from other diversification mechanisms such as landscape barriers, ecological gradients, and BARRATT ET AL | 13 of 18 anthropogenic habitat fragmentation, would represent a powerful approach for gaining a more mechanistic understanding of population diversification.

| CONCLUSIONS
Deeper insights into historical climate change in highly diverse tropical regions are essential to understand their contribution towards shaping current biodiversity patterns, and how biodiversity loss might be predicted given future projections of anthropogenic and climate change. By increasing the spatial resolution and number of time periods used to project SDMs, the data set we present here is able to improve on existing paleoclimate data that are temporally limited to build fine-scale models of changing habitat suitability for species through time accounting for their dispersal ability. Integrating this data set with other data types in the future (e.g., ecological, behavioral, and molecular) will help to increase our general understanding of how climate change impacts biodiversity, and how we may mitigate against predicted biodiversity loss in the future.

ACKNOWLEDGMENTS
We are indebted to several organizations including the Rwanda De-