This study focuses on North and South Kivu provinces in eastern DRC, which straddle the equator and extend across the Albertine Rift between ~27° and 30°E, with a combined surface area of over 120,000 km² as illustrated in Fig. 1 (Ngalamulume et al. 2025). These provinces form part of the Albertine Rift, a core component of the Eastern Afromontane Biodiversity Hotspot and one of the most species-rich regions in continental Africa (Ayebare et al. 2018; Plumptre et al. 2021). Elevation ranges from 770 m in lowland valleys to over 5,000 m on volcanic and montane massifs. The region experiences a tropical humid climate, with mean annual temperatures between 17 and 24°C and distinct dry and wet seasons. Protected areas like Kahuzi-Biega National Park hold significant portions of its ecological diversity and carbon stock (Cirezi et al. 2025).

Figure 1. 

Study area showing North and South Kivu provinces in eastern Democratic Republic of the Congo, with Cinchona occurrence records overlaid on land-cover classes.

A total of 456 unique presence points for Cinchona spp. were assembled from multiple sources, including field surveys conducted in 2024 across North and South Kivu, and validated datasets provided by local partners such as Pharmakina S.A. All occurrence coordinates were projected in WGS 84 (EPSG:4326) and duplicates and spatially imprecise records were removed to minimize sampling bias (Palacio et al. 2021; Davis et al. 2024). After processing, 125 high-confidence, georeferenced occurrence records were retained for integration with environmental predictor layers in the species distribution modelling framework (Suppl. material 1).

Environmental variables were selected based on their documented ecological relevance to the distribution of Cinchona spp. along elevational climatic gradients typical of humid tropical montane forests (Ayebare et al. 2018; Coronel-Castro et al. 2024).

A total of 32 environmental variables were initially considered to model the potential distribution of Cinchona spp. in eastern DRC (Suppl. material 2). This dataset included 19 bioclimatic variables and one environmental variable (solar radiation), all obtained from WorldClim v.2.1 at a spatial resolution of 30 arc-seconds (~1 km) (Fick and Hijmans 2017; accessed on 21 Jun. 2025). Three topographic variables (elevation, slope, and aspect) were derived from the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM), downloaded from the CGIAR Consortium for Spatial Information (https://srtm.csi.cgiar.org; accessed on 21 Jun. 2025). Nine edaphic variables were extracted from SoilGrids (https://soilgrids.org; accessed on 21 Jun. 2025) at 250 m resolution. Future climate effects were evaluated using WorldClim v.2.1 CMIP6 ensemble projections for two mid-century scenarios in eastern DRC (SSP2-4.5 and SSP5-8.5) at 2050, defined as 2041–2060 (https://www.worldclim.org/data/cmip6/cmip6_clim30s.html; accessed on 21 Jun. 2025). All environmental predictors were resampled to 250 m with bilinear interpolation and aligned to EPSG:4326. Raster processing and data management were performed using QGIS v.3.34 (QGIS Development Team 2023) and R v.4.3.0 (R Core Team 2025).

To reduce multicollinearity among environmental predictors and enhance model reliability, a two-step variable selection procedure was implemented. First, pixel values for the 32 initial environmental variables were extracted at all 125 occurrence locations using the Point Sampling Tool in QGIS. Pairwise Pearson correlation analysis was then conducted to assess inter-variable relationships. Variables with correlation coefficients greater than 0.7 (|r| > 0.7) (Smith and Santos 2020) were considered redundant and excluded to minimize collinearity (Bradie and Leung 2017; Schnase et al. 2021).

In the second step, Variance Inflation Factors (VIFs) were computed on the reduced dataset using the vifstep function in the R package usdm v.2.1-7 (Naimi et al. 2014). Variables with VIF values exceeding 10 were iteratively removed until all remaining predictors met the threshold of VIF < 10 (Ab Lah et al. 2021; Kim et al. 2025; Lorenț et al. 2025). This process resulted in a final set of 10 predictors retained for MaxEnt modelling: bio2 (mean diurnal range), bio3 (isothermality), bio4 (temperature seasonality), bio16 (precipitation of wettest quarter), bio18 (precipitation of warmest quarter), elevation, slope, aspect, pH water, and coarse fragments. All the ten spatial layers were processed and formatted as ASCII grids compatible with MaxEnt software.

Species distribution modelling was performed using the Maximum Entropy algorithm implemented in MaxEnt v.3.4.4 (Phillips et al. 2025). The 10 environmental predictors retained after multicollinearity filtering (bio2, bio3, bio4, bio16, bio18, elevation, slope, aspect, pH water, and coarse fragments) were used as input. The modelling extent was restricted to the North and South Kivu provinces to avoid extrapolation into environmentally dissimilar areas. Model calibration used 125 georeferenced occurrences together with the ten predictors. Predictive performance was assessed by 10-fold cross-validation (ten replicates), with each fold trained on 90% of the occurrences and validated on the remaining 10%; replicates were generated independently (Ali et al. 2023; Kim et al. 2025). The background was defined within the study area, and the maximum iterations were set to 5,000 to ensure model convergence (Boussouf et al. 2023). Output was generated in logistic format on the [0,1] scale, producing habitat suitability scores ranging from 0 (unsuitable) to 1 (highly suitable) (Xiao et al. 2024). For cartographic interpretation, the continuous logistic surface was reclassified into four suitability classes, following prior work adapted to our context (Vergara et al. 2023): very low (≤ 0.100), low (0.100–0.300), moderate (0.300–0.500), and high (> 0.500).

Predictive performance was quantified with the threshold-independent Area Under the ROC Curve (AUC), a standard metric in ecological niche modelling (Fielding and Bell 1997; Jiménez-Valverde 2012). AUC values below 0.60 indicate very poor discrimination, values between 0.60 and 0.69 are considered poor, 0.70–0.79 satisfactory, 0.80–0.89 good, and ≥ 0.90 excellent (Lissovsky and Dudov 2021).

The contribution of each environmental variable to the MaxEnt model was assessed using jackknife tests and response curves. Jackknife tests quantified the gain reduction when each variable was omitted, identifying those with the highest unique explanatory power (Li et al. 2020; Vergara et al. 2023). Response curves were examined to characterize the relationship between habitat suitability and key predictors (Ab Lah et al. 2021; García et al. 2022). These analyses helped to interpret ecological drivers of Cinchona distribution and to compare our findings with those of earlier studies.

Predictive performance was evaluated with 10-fold cross-validation (ten independent replicates; 90/10 splits). The model achieved high discrimination, with Test AUC = 0.8808 (SD = 0.0437) and Training AUC = 0.9036 (Suppl. material 3). For comparability, we report two threshold-dependent statistics from the average model: the 10th percentile training presence threshold (0.1552) yielded a test omission = 0.1352, while the maximum test sensitivity + specificity threshold (0.3196) yielded test omission = 0.1596. The ROC summarizing cross-validated performance is shown in Suppl. material 4; the omission-predicted area curve across cumulative thresholds is provided in Suppl. material 5.

Together, relative contributions and jackknife diagnostics indicate that a small set of predictors governs Cinchona distribution. Mean diurnal range (bio2) accounts for 37.6% of training gain but has low permutation importance (3.7%), which is consistent with overlap among temperature metrics. In contrast, elevation (17.1%) and temperature seasonality, bio4 (16.5%), show high permutation importance (20.2% and 23.2%, respectively), indicating substantial performance loss when perturbed in permutation tests (Suppl. material 6). Jackknife curves confirm that elevation carries the most unique information; it yields the highest “with-only” gain and its omission produces the largest decline (Suppl. material 7). Isothermality (bio3) contributes 8.2% yet attains the highest permutation importance (24.2%), revealing a strong non-redundant signal after accounting for collinearity. In 2050 projections under SSP2-4.5 and SSP5-8.5, this ranking (temperature regime: bio2/bio3/bio4 plus elevation) remains broadly stable, reinforcing the robustness of these drivers.

Response curves further clarify the climatic, topographic, and edaphic controls on suitability. In the marginal curves (red = mean across 10 replicates; blue = ±1 SD; other predictors held at their means), suitability rises steeply with elevation and levels off at mid to high elevations; it increases monotonically with mean diurnal range (bio2) across the sampled range (to ~13°C) and declines with temperature seasonality (bio4). For edaphic predictors, both coarse fragments and pH water show sustained negative responses: suitability decreases as coarse fragment content increases, and declines with increasing pH. Taken together, these patterns support a composite thermal regime (bio2/bio3/bio4) coupled with an altitudinal constraint, rather than reliance on a single temperature metric (Fig. 2).

Figure 2. 

Response curves key predictors. A. Elevation: suitability rises quickly from ~540 m, reaches a near-plateau around ~2000–2500 m, and remains high up to > 4000 m. B. Mean diurnal range/bio2: suitability increases monotonically across the full range, highest near 12.909°C. C. Temperature seasonality/bio4 (8.774–71.280): suitability declines steadily as seasonality increases, high around ~9–20. D. Isothermality/bio3 (69.162–94.228%): suitability is highest at lower–moderate values (~69–80), lowest near ~94.

The modelling results provide a clear spatial baseline for anticipating climate change and land-use impacts on Cinchona habitats in eastern DRC (Fig. 3). High suitability (> 0.500) is concentrated along the eastern escarpment from Lubero in the north to Walungu and Fizi in the south, forming an elongated corridor where plantations coincide with remnant montane forests, thereby delineating landscapes where land-use choices directly intersect with suitable Cinchona habitats. This spatial configuration highlights a close overlap between optimal Cinchona habitat and montane landscapes shaped by long-term human land use, indicating that present-day land-use decisions strongly influence habitat persistence. Moderate suitability (0.300–0.500) surrounds these hotspots and extends into adjacent mid-elevation zones. Low suitability (0.100–0.300) covers much of the western foothills, reflecting transitional landscapes increasingly shaped by agriculture and settlement. Unsuitable areas (≤ 0.100) dominate the lowlands and interior plateaus.

Figure 3. 

Current habitat suitability for Cinchona spp. in North and South Kivu derived from MaxEnt (logistic output reclassified into four classes: ≤ 0.100, 0.100–0.300, 0.300–0.500, > 0.500). Warmer colours indicate higher suitability. District boundaries are shown for orientation.

Area statistics (Suppl. material 8) indicate that 4,799.31 km² (4.13%) of the study area falls in the high-suitability class, 7,533.39 km² (6.49%) in the moderate class, 32,142.64 km² (27.67%) in the low class, and 71,671.57 km² (61.71%) is unsuitable. High suitability is most extensive in Lubero (1,509.58 km²), Walungu (811.17 km²), Oicha (718.58 km²), and Kabare (584.47 km²), all located on montane ridges between roughly 1,400 m and 2,300 m above sea level, a narrow elevational band that constrains the spatial extent of optimal habitat. Moderate suitability is more widespread but largely confined to mid-elevation slopes, and only trace amounts of high suitability occur in lowland territories such as Shabunda and Fizi. These statistics confirm the restricted distribution of optimal Cinchona habitat under present conditions.

The composite climatic niche of Cinchona in eastern DRC is governed by a cool, thermally stable regime combined with a clear altitudinal constraint. Our MaxEnt model highlighted diurnal temperature range (bio2), isothermality (bio3), and temperature seasonality (bio4) as key climatic predictors, with elevation acting as a major non-climatic driver. Diurnal range (bio2) accounted for 37.6% of training contribution but only 3.7% permutation importance and was not supported by the jackknife test (Suppl. material 7). This discrepancy is consistent with collinearity among temperature variables, since percent contribution can be inflated when predictors are correlated (De Marco and Corrêa Nóbrega 2018; Feng et al. 2019). Permutation importance, by contrast, assesses how much predictive power is lost when a variable is permuted, and is generally less biased in the presence of correlated variables (Hooker et al. 2021). The high contribution but low permutation of bio2 therefore indicates that mean diurnal range strongly co‑varies with other temperature variables; it helps fit the training data but carries little unique information. The jackknife results confirm this: removing bio2 caused only a small drop in gain, whereas the model built with elevation alone had the highest gain and the largest drop when excluded. We therefore interpret bio2 as part of a composite thermal regime (bio2/bio3/bio4) rather than a uniquely driving factor. Elevation (17.1% contribution; 20.2% permutation) and the two indices of temperature variation (bio3 and bio4) provided non-redundant signals. Bio3 (isothermality) had moderate contribution (8.2%) but the highest permutation importance (24.2%), indicating that the ratio of diurnal to annual temperature range is critical for predicting suitable sites; bio4 (seasonality) followed closely (16.5% contribution; 23.2% permutation). This coincides with the findings of García et al. (2022), who also demonstrated that elevation, isothermality (bio3), and temperature seasonality (bio4) were among the most influential predictors of Cinchona distribution in Peru. Their study similarly emphasized the ecological relevance of thermal stability and altitudinal gradients in shaping the species’ niche, reinforcing the robustness of our model across distinct biogeographic contexts. Comparable results were reported by Coronel-Castro et al. (2024) in north-eastern Peru, where bio3 and bio4, alongside edaphic factors such as cation exchange capacity (CEC), emerged as key contributors to the distribution of multiple Cinchona species.

Marginal response curves help translate these statistics into ecological interpretations. Elevation showed a rapid increase in suitability from ~540 m, with a plateau between 2,000 and 2,500 m above sea level. This plateau coincides with the crest zones of the Albertine Rift, where montane forests and subalpine grasslands provide cool, moist microclimates (Martin and Burgess 2023). Although bio2 contributed substantially to model fit, its low permutation importance and lack of jackknife support suggest that it does not act independently. Nevertheless, the marginal response shows monotonic increases within the observed window (8.173–12.909°C). The response to bio4 declined steeply from 8.77 to 71.28, indicating that high temperature seasonality reduces suitability. Similarly, suitability declined with increasing isothermality (bio3) beyond 69.16%, suggesting preference for climates where daily variability does not dominate the annual range. Suitability remains high across acidic soils (pH ≲ 5.5–6) and declines steadily towards neutral–slightly alkaline conditions (up to ~7.8), indicating a preference for moderately acidic substrates (Nair 2010). Suitability declines sharply with increasing coarse-fragment content. Suitability is highest at very low contents (≲ 2–4%), drops to moderate levels around ~8–10%, and approaches zero above ~25%, consistent with a preference for fine-textured, stone-poor substrates. Together, these patterns depict an ecological niche defined by cool, seasonally stable montane climates at mid-to-high elevations, with moderate day–night thermal variation, relatively low isothermality, and fine-textured, moderately acidic, stone-poor soils.

The present suitability map shows a continuous corridor of high to moderate suitability along the eastern escarpment, notably in the territories of Lubero, Walungu, Oïcha, Kabare, and around Butembo. Climate projections under SSP2‑4.5 and SSP5‑8.5 suggest profound habitat contraction and upslope fragmentation by 2050. This shift from moderate and low suitability to unsuitability reflects the combined effects of warming, increased temperature seasonality and edaphic stress, pushing the niche upslope and constricting its extent.

Spatially, the continuous corridor breaks into isolated patches on crest zones. High suitability persists in parts of Lubero, Oïcha, and Kabare, but these areas contract and shift upslope. Walungu exhibits a local persistence signal: under SSP5-8.5, a few crest-linked cells remain above the 0.5 threshold, consistent with topographic buffering and potential microrefugia at high elevations. In contrast, lowland territories such as Shabunda and Fizi become almost completely unsuitable, with Shabunda being ~99.8% unsuitable under both scenarios and the Fizi-Mwenga highlands retaining only residual suitable patches. These shifts mirror broader ecological predictions that montane species must track cooler conditions upslope, losing area as mountain surfaces narrow (Conniff 2018). Afrotropical montane birds provide a parallel: upslope shifts of lower elevational range limits combined with restricted dispersal lead to range contractions in fragmented forests (Neate-Clegg et al. 2021). Comparable patterns have been documented for indigenous montane flora of the Albertine Rift, where plant species distributions are strongly constrained by elevation and increasingly fragmented by land-use change, as shown by regional niche-based conservation planning analyses (Ayebare et al. 2018; Plumptre et al. 2021). The fragmentation observed in our projections implies that Cinchona may experience a similar “escalator to extinction” effect (Conniff 2018), if dispersal and regeneration cannot keep pace with climatic changes.

The spatial reconfiguration of suitable habitat has direct consequences for conservation and land‑use planning in the eastern DRC. Our projections indicate that future climatic suitability for Cinchona will be increasingly restricted to crest-linked montane areas, reinforcing the importance of identifying and safeguarding climatic refugia. Such refugia are widely recognised as critical components of climate-adaptation strategies, as they provide relatively stable microclimatic conditions that can buffer species against regional warming (Hannah et al. 2014; Morelli et al. 2016).

By 2050, most high-suitability areas lie on ridge crests around Lubero, Oïcha, Kabare, and parts of Walungu. These small refugia should be prioritised for strict protection. These areas are also subject to strong land-use pressure, as montane zones in the Albertine Rift are increasingly targeted for agriculture and settlement, often overlapping with climatically resilient habitats. Prioritising crest-linked refugia for conservation therefore requires their explicit integration into provincial land-use plans and forest zoning frameworks. We recommend establishing a network of 10–50 ha micro-reserves around the remaining high-suitability patches, either as strictly protected nature reserves or as legally recognised community-managed conservation areas, depending on local land tenure and governance contexts, with formal legal recognition and community-based co-management, particularly in fragmented landscapes (Neate-Clegg et al. 2021). Management should also preserve montane forest canopy to maintain microclimate stability, as canopy cover lowers understorey temperatures and limits soil evaporation (Watts et al. 2022).

As the climate warms, suitable habitat shifts uphill. Without connections, high-elevation patches become isolated. We recommend elevational corridors that link low, mid, and high zones, especially between 1,400 and 2,600 m identified by our response curves. Maintaining such altitudinal connectivity is increasingly recognised as a core element of climate-wise conservation planning, as it facilitates dispersal, gene flow and adaptive range shifts in mountainous landscapes where available habitat narrows with elevation (Hilty et al. 2020). This reduces dispersal barriers and supports gene flow. Evidence from Afrotropical montane birds shows that fragmented forests block upslope shifts and cause range losses (Neate-Clegg et al. 2021). In the Albertine Rift, corridor-based conservation would also complement existing transboundary initiatives and regional biodiversity strategies.

Many suitable areas overlap with settled highlands where Cinchona is already cultivated. Agroforestry can buffer heat, stabilise soils, and diversify livelihoods. In such human-dominated landscapes, strict protection alone may be neither feasible nor socially acceptable, making climate-smart agroforestry a necessary complement to reserve-based conservation. Shade trees can lower air temperature by ~4°C and soil temperature by 6–10°C, while regulating humidity and improving soil moisture (Watts et al. 2022). We recommend integrating diverse canopy species such as Ficus spp. and native Albizia spp. with Cinchona plantations, as well as using mulching and organic amendments to improve soil structure. Avoid coarse-fragment and high-pH soils, and conserve moist valley bottoms where microclimates are more stable.

More broadly, our results demonstrate the relevance of species distribution models as decision-support tools for land-use policy and climate adaptation. Climate-informed suitability maps can guide environmental impact assessments, restoration planning, and the spatial targeting of national climate frameworks such as REDD+ and National Adaptation Plans. Embedding such spatially explicit ecological information into planning and governance processes would strengthen anticipatory decision-making and enhance the long-term effectiveness of conservation and land-use policies in eastern DRC (Guisan et al. 2013; Pavón-Jordán et al. 2015; Babanezhad and Naqinezhad 2025 and references therein).

Several sources of uncertainty qualify our conclusions. First, the occurrence data may be biased towards accessible areas such as roads and villages, a well‑known accessibility bias in species distribution modelling. This well-documented accessibility bias (Meyer et al. 2015) may lead to overrepresentation of human-modified habitats and under-sampling of remote regions. As a result, model predictions may not fully capture the species’ ecological breadth. Future work could address this by applying bias correction techniques such as spatial filtering or background manipulation (Kramer-Schadt et al. 2013).

Second, 10‑fold cross‑validation used here assumes that presence records are independent and identically distributed. Random cross‑validation can overestimate model performance because spatial autocorrelation violates this assumption, leading to optimistic statistics (Tziachris et al. 2023). Spatially blocking data into non‑overlapping subsets is recommended to reduce this bias (Koldasbayeva and Zaytsev 2025). Future studies should evaluate alternative partitioning strategies (spatial blocking, environmental clustering) and report sensitivity of AUC and omission rates to these methods.

Third, collinearity among temperature variables complicates the interpretation of variable importance. Metrics such as percent contribution can be skewed by correlated predictors, whereas permutation importance more accurately reflects each variable’s unique influence. In this study, we reduced redundancy by applying Variance Inflation Factor (VIF) and Pearson correlation analysis to retain only five relatively uncorrelated predictors (bio2, bio3, bio4, bio16, and bio18). Nevertheless, residual collinearity may persist, and its potential influence cannot be entirely excluded. Future research could explicitly compare the effectiveness of different dimensionality-reduction approaches, such as VIF/Pearson filtering versus Principal Component Analysis (PCA) (De Marco and Corrêa Nóbrega 2018), to assess which better minimizes multicollinearity while preserving ecological interpretability.

Finally, beyond technical limitations, it is crucial to consider the socio-ecological context shaping species distributions. Research should explore socio‑economic drivers of land‑use change and the potential for incentive schemes (e.g. payment for ecosystem services) to support conservation. By addressing these interconnected ecological, methodological, and socio-economic uncertainties, future studies will enhance the reliability and policy relevance of GeoAI and climate-informed species distribution models in the eastern DRC.