Benin is located in West Africa between 6°25′ N to 12°30′N and 0°45′E to 4°00′E. Generally, Benin is a low lying country, with altitudes varying from sea level to 400 m a.s.l., although it can go up to 650 m in the northwest of the country, where temperatures can also be exceptionally high (35–40°C) (Adomou et al., 2006). Agriculture is the main source of livelihood in rural areas. Like most developing countries, Benin is experiencing rapid population growth, with the recent human population estimated at 12.5 million people (United Nations, 2021) and associated with increasing urbanisation and land-use change (Guidigan et al., 2019).

Benin is divided into three contrasting climatic zones (Adomou et al., 2006): Guinean (in the south), Sudanian (in the north), and Sudano-Guinean (a transition zone in the centre). Basically, the Guinean zone is characterised by a sub-humid climate, bimodal rainfall (April–June & September–November) averaging 1,200 mm/year, a temperature range of 18–33°C, relative humidity ranging between 30 and 98%, ferrallitic and hydromorphic soils without concretions, and woodlands and fallows as the dominant habitats and/or vegetation types. While, sub-humid climate, unimodal rainfall (May–October) ranging between 900 and 1,100 mm/year, an annual temperature range of 20–36°C, relative humidity of 31–98%, ferruginous to somewhat tropical ferruginous and ferrallitic soils on crystalline basement characterise the Sudano-Guinean zone. This zone is dominated by a mosaic of savanna woodlands, dense shrub and tree forests, and gallery forests. On the other hand, the Sudanian zone is characterised by Sudano-Sahelian climate, unimodal rainfall (700–900 mm/year), a relatively wide temperature range (17–42°C), and a wide relative humidity (18–99%). Furthermore, this zone is dominated by ferruginous soils on crystalline basement and mostly dry shrubby forests and savannas (Adomou et al., 2006; Assogbadjo et al., 2011; Hounkpatin et al., 2022). Like in most parts of the sub-Saharan Africa, the soils in Benin are mostly sandy and infertile (Hounkpatin et al., 2022).

Vigna CWRs are herbaceous, annual or perennial, climbing, scrambling or prostrate plants. Perennial species generally have large, woody rootstocks which often die back in cooler months only to grow again in warm weather or following burning e.g., Vigna frutescens A. Rich (Maxted et al., 2004). Their sizes generally range from <1 m for species like Vigna laurentii De Wild. to over 7 m e.g., Vigna racemosa (G. Don) Hutch. & Dalziel. Their stems may be glabrous or with various levels of pubescence. They are found in a wide range of habitats such as savannas, grasslands, open woodlands, and shrublands usually at low altitudes, with some species such as Vigna luteola (Jacq.) Benth. and V. laurentii often associated with wet areas. Most Vigna CWR taxa in Benin flower and reproduce between August/September and December. Like other leguminous plants, wild Vigna species have a ballochory (explosive dehiscence) seed dispersal as their primary seed dispersal mechanism (Lush et al., 1980; Trzeciak-Limeira et al., 2013); thus, they are likely to be mostly short-distance dispersed plants, probably dispersing their seeds over a 1–5 m radius (Vittoz and Engler, 2007; Parker et al., 2021). A few species like V. luteola are known to be dispersed over long distances by sea-drifting (Miryeganeh et al., 2014). While those that are used as forage such as V. racemosa and Vigna reticulata Hook. f. (Catarino et al., 2021), are likely to be dispersed over long distances by secondary agents such as herbivores (Ramírez-Rodríguez et al., 2021; Wang and Hou, 2021). Vigna CWR taxa favour warm temperatures, but higher temperatures (>36°C) are detrimental to their metabolic processes such as photosynthesis (Farooq et al., 2017).

Currently, Benin has about 31 Vigna CWR taxa that have been described (Supplementary Table 1), but this study used only 15 of these; the rest were left out as they had ≤ 10 occurrence records (Yesuf et al., 2021). The full list of the studied taxa were Vigna comosa Baker, Vigna filicaulis Hepper, V. frutescens, Vigna gracilis (Guill. & Perr.) Hook. f., Vigna heterophylla A. Rich., V. laurentii, Vigna longifolia (Benth.) Verdc., V. luteola, Vigna multinervis Hutch. & Dalziel, Vigna nigritia Hook. f., Vigna oblongifolia A. Rich., V. racemosa, V. reticulata, Vigna unguiculata subsp. baoulensis (A. Chev.) Pasquet, and Vigna unguiculata var. spontanea (Schweinf.) Pasquet. Of these 15 taxa, only three (V. heterophylla, V. unguiculata subsp. baoulensis, and V. unguiculata var. spontanea) were at the time of this study not yet assessed using the IUCN Red List Categories and Criteria (Supplementary Table 1). The spatial distributions of the 15 studied taxa are presented in Figure 1.

This study is a follow-up on the national inventory and prioritisation of CWRs for Benin (Idohou et al., 2013). The Global Biodiversity Information Facility (GBIF, www.gbif.org) was the main source of species occurrence records (1,952) [accessed on July 14, 2020]. The RAINBIO (https://gdauby.github.io/rainbio) [accessed on July 16, 2020], and the Genesys Global Portal on Plant Genetic Resources (https://www.genesys-pgr.org) [accessed on July 15, 2020], provided additional records, 128 and 24, respectively.

In addition, random field visits were made to verify location data related to the distribution of some species and partly minimise the sampling bias in the database (Meng et al., 2021). To ensure data quality, visual inspection was used to identify outliers and these were clipped out. Records collected earlier than 1990 were removed from the data set to reduce the effects of temporary bias (Idohou et al., 2017b). To reduce clumping bias, duplicate records were removed from the data set, and occurrence records were spatially thinned to a geographic distance of 1 × 1 km² (Idohou et al., 2017b) using Environmental Niche Modelling (ENM) tools (www.ENMTools.com) (Warren et al., 2010) performed in QGIS version 3.8.1 (QGIS Project https://qgis.org). Spatial thinning not only reduces model overfitting but also improves the performance of models better than background manipulation (Ratnayake et al., 2021). The Flora of Benin (Akoegninou et al., 2006) and experts were consulted to verify the adequacy of the observed range distribution of the target species.

The study used a combination of bioclimatic and soil variables. Bioclimatic data represent annual trends in climate conditions, seasonality and climate extremes, which may impact reproduction and survival of species over broad extents (Aguirre-Gutiérrez et al., 2017; Idohou et al., 2017a). On the other hand, soil variables may directly constrain the establishment and development of species (Aguirre-Gutiérrez et al., 2017), and their incorporation in SDM, especially at the local scale, which was the case in the current study, appears to improve the predictive capacity of models (Hageer et al., 2017; Zuquim et al., 2020).

Twenty-one bioclimatic variables for the present and future scenarios were downloaded from the AfriClim database version 3.0 at a spatial resolution of 30 arc seconds (~1 km²) (https://www.york.ac.uk/environment/research/kite/resources/) (Platts et al., 2015) [accessed on March 17, 2020]. The database spans 10 general circulation models (GCMs), downscaled using five bias-corrected regional climate models (RCMs) and four contemporary baselines, under two representative concentration pathways of the IPCC-AR5 (RCP4.5 and RCP8.5) (Platts et al., 2015). For the current conditions, the WorldClim v1.4 option with the baseline year of 1975 (1950–2000) was used from the database. For the future climatic conditions—horizon 2055 (2041–2070), the AfriClim Ensembles 3.0 with the WorldClim as the baseline was used, and two scenarios, RCP 4.5 (optimistic scenario) and RCP 8.5 (pessimistic scenario) (Platts et al., 2015) were considered. The two scenarios seem to be plausible for Africa (Platts et al., 2015), and have been widely used in SDM studies across Africa (Lompo et al., 2021; Zuza et al., 2021; Assogba et al., 2022). The mid-Century horizon (the 2050s) was chosen to align with the United Nations framework of global challenges in agriculture and food security (Zuza et al., 2021), which also dovetails well with the Agenda 2063 for Africa.

Soil data were obtained at 250 m resolution from the Africa Soil Profiles Database (http://www.isric.org) (Hengl et al., 2017). Eleven sets of bio-physiochemical soil characteristics were downloaded. These were bulk density (t/m³), soil organic carbon (g/kg), pH in water, clay content (%), sand content (%), silt content (%), cation exchange capacity (cmolc/kg), exchangeable acidity (cmolc/kg), exchangeable Ca (cmolc/kg), exchangeable K (cmolc /kg), and exchangeable Mg (cmolc/kg). These soil variables have been used in previous studies in Benin (e.g., Idohou et al., 2017b; Vihotogbé et al., 2021). Since soil-plant interactions seem to be critical within the 0–16 cm soil depth (Goebes et al., 2019) and given that Vigna species are herbaceous small-statured plants, only a maximum of three soil depth horizons (0–5, 5–15, and 15–30 cm) were considered for this study. A total of 30 soil layers were thus downloaded. The soil data were resampled at 1 km resolution to match the resolution of the bioclimatic variables.

Maximum Entropy (Maxent ver. 3.4.1) algorithm (Phillips et al., 2006) was used. As with other correlative models, the MaxEnt procedure establishes the relationship between species occurrence records at sites and the environmental variables and/or spatial characteristics of such sites (Phillips et al., 2006; Elith et al., 2011). Despite its limitations (Lissovsky and Dudov, 2021), MaxEnt is among the many SDM tools (see Elith and Graham, 2009) with increasing use in conservation-oriented studies owing to its high predictive accuracy, stability and reliability even with presence-only data and small data sets (Elith et al., 2011; Phillips et al., 2017; Vincent et al., 2019; Çoban et al., 2020; Mponya et al., 2021; Ratnayake et al., 2021). Moreover, it produces spatially open habitat suitability maps and evaluates the significance level of individual environmental variables using the built-in Jackknife test (Çoban et al., 2020), which were among the core objectives of the current study.

Dealing with variable collinearity in SDM is evolving just like the algorithms themselves, with the use of the Pearson's correlation test in ENMTools (Warren et al., 2010) and priori selection of variables based on the ecological system (Scridel et al., 2021) as some of the common approaches to reducing potential multicollinearity. However, Feng et al. (2019) showed that the exclusion of highly correlated variables does not significantly influence model performance, especially those built by MaxEnt, as the algorithm accounts for redundancy in variables. According to Mod et al. (2016), neglecting eco-physiological meaningful predictors could result in incomplete niche quantifications, thereby limiting the predictive power of SDMs. The authors thus suggested that the selection of climatic-related variables should be determined by the environmental conditions of the study site and the requirements of the target species. Therefore, this study followed the approach used by Singh et al. (2021) of eliminating variables that consistently contributed less or nothing to the model for three successive runs. Firstly, two variables, minimum temperature of the coldest month (BIO 6), and mean temperature of the coldest quarter (BIO 11), were deemed unsuitable and were thus removed based on expert knowledge (Vihotogbé et al., 2021). Thereafter, all the remaining 49 variables were used in the model pre-assessment (but see Dormann et al., 2013) from which variables that contributed <5 % to each of the model explanatory power after every run for three successive runs were excluded following a Jackknife test (Singh et al., 2021). Finally, five variables that contributed the most to each of the 15 models were retained for modelling (Table 1).

Since the area included in the background can influence MaxEnt model fit (Merow et al., 2013; Pang et al., 2021), and given that models are more reliable when built at a larger scale (Barve et al., 2011), as this reduces the risk of niche truncations (Pang et al., 2021), species occurrence records were thus first filtered for West Africa, and models trained by projecting the present-day variables over West Africa. To further improve the model performance, the maximum number of iterations in MaxEnt was adjusted to 5,000, the cross-validation method was used, and the number of replications was increased to 15 (Abrha et al., 2018). Increasing the number of iterations and replications provides, respectively, ample time for the model to converge and run multiple times to develop superlative averaged results (Merow et al., 2013; Abrha et al., 2018). The “clogclog” was used as an output format (Favi et al., 2021; Scridel et al., 2021), since it appears to provide the model output with a stronger theoretical justification than the logistic transform (Qu et al., 2018) and seems to give realistic binary predictions of species distributions (Scridel et al., 2021). The jackknife test was used to determine the variable contribution to the models.

The model performance was assessed through two widely used metrics, the threshold-independent Area Under the Curve (AUC) of the Receiver Operating Characteristic (ROC) curve and the threshold-dependent True Skill Statistic (TSS) (Fielding and Bell, 1997; Allouche et al., 2006; Favi et al., 2021; Ratnayake et al., 2021). The AUC measures the model's ability to distinguish between random and background points (AUC = 0.5) with values ranging from 0 to 1 (Fielding and Bell, 1997; Favi et al., 2021). AUC values closer to 1 indicate good-performance models, and a stronger correlation between predictor variables and the distribution of the target species (Fielding and Bell, 1997). As such, models were considered acceptable if 0.7 ≤ AUC < 0.8, good if 0.8 ≤ AUC < 0.9 and excellent if AUC ≥ 0.9 (Ratnayake et al., 2021). The TSS is a measure of accuracy i.e., the capacity of the model to detect the true presence (sensitivity) and true absences (specificity), expressed as the sensitivity plus specificity-1 (Allouche et al., 2006). Its values range from −1 to +1, and like AUC, TSS values closer to 1 indicate good-performance models (Allouche et al., 2006). Therefore, models were described as poor (TSS <0.4), acceptable (0.4 ≤ TSS < 0.8), and very good (TSS > 0.8) (Ratnayake et al., 2021). Response curves were used to further evaluate and quantify the biological plausibility of the models since they show the predicted relative occurrence rate (ROR) of species against the value of a predictor variable (Merow et al., 2013). The variables were finally clipped to Benin and converted to Ascii format using SDM Tools for use in modelling.

MaxEnt was used to predict the habitat suitability distribution for each of the 15 species under the current and future climate change scenarios using 1,432 occurrence records (range: 14 to 333) and a combination of bioclimatic and edaphic variables as inputs (Table 1). The same settings as in the model calibration were used. Habitat suitability maps were developed following Ramírez-Rodríguez et al. (2021) with slight modifications. Briefly, output rasters from MaxEnt were converted into binary layers and exported to QGIS to develop habitat maps. Binarization of continuous habitats, in spite of its shortcomings (Santini et al., 2021), is still one of the popular approaches for delimiting habitats in SDM, quantifying range changes and building species richness overtime in conservation studies (Politi et al., 2021; Ramírez-Rodríguez et al., 2021; Singh et al., 2021; Vihotogbé et al., 2021; Lima et al., 2022). It is believed that binairization avoids the effects of model over-fitting (Vihotogbé et al., 2021), and that it makes interpretation of distribution maps much easier compared with the more liberal interpretation of a continuous habitat output (Singh et al., 2021). Using the 10th percentile presence threshold to separate suitable from unsuitable habitats Ramírez-Rodríguez et al. (2021), all values above and below this threshold were considered suitable and unsuitable habitats, respectively. According to Politi et al. (2021), compared with other percentile thresholds, the 10th percentile maximises the correct prediction of the percentage of presences and absences, thereby providing conservative species distributional ranges.

Subsequently, to calculate range changes between the current and future habitat distributions, binary-thresholded future rasters were subtracted from the current rasters (Ramírez-Rodríguez et al., 2021) using the range shift tool in the SDM Toolbox. This tool classifies output layers into four classes: expansion in range (absence in current, presence in future), no occupancy (absence in both current and future), occupancy or stable (presence in current and future) and contraction in range (presence in current, absence in future). Accordingly, three range changes in habitat suitability (stability, expansion and contraction) were mapped and areal extents were calculated.

To prioritise areas for in situ conservation, binary raster layers for the range change of all the 15 species and each of the three distributions (current, 2055 RCP 4.5 and 2055 RCP 8.5) were summed up to show overlaps of potential habitat distributions (Ramírez-Rodríguez et al., 2021). Each suitable distribution was thereafter reclassified into three classes (low suitability, moderate suitability, and high suitability) using the SDM Toolbox. Finally, these were overlaid onto the PA network of Benin to further estimate the extent to which the PAs would conserve the 15 Vigna CWR taxa. The area of each distribution was calculated using the extract by mask tool in the SDM Toolbox.

The AUC and TSS values ranged from 0.914 to 0.997 (median = 0.956; mean = 0.957) and 0.617 to 0.876 (median = 0.771; mean = 0.774), respectively (Figure 2). Likewise, the ROC curves were away from the random distribution (not shown here), and there were lower differences between AUC values for training and test samples for the model corresponding to all the 15 taxa (range = 0.001–0.025; mean = 0.008), indicating good performance and high accuracy of the model for generalisation. The models were thus considered excellent and highly informative to describe the distribution patterns in the 15 modelled Vigna CWR species in Benin.

The modelled species demonstrated differential responses to different environmental variables. Based on the frequency of the power of predictive contribution of each variable across all the 15 species, three variables, precipitation of the driest quarter (BIO 17), isothermality (BIO 3), and mean diurnal range in temperature (BIO 2), in that order, were found to contribute the most in explaining the majority of the models (Table 2). For instance, precipitation of the driest quarter was consistently found among the dominant three variables for 10 of the 15 models (V. frutescens, V. gracilis, V. heterophylla, V. longifolia, V. luteola, V. nigritia, V. oblongifolia, V. racemosa, V. reticulata, and V. unguiculata subsp. baoulensis). This was closely followed by isothermality (seven models: V. heterophylla, V. longifolia, V. multinervis, V. racemosa, V. reticulata, V. unguiculata subsp. baoulensis, and V. unguiculata var. spontanea); whilst mean diurnal range in temperature predicted five models (V. filicaulis, V. frutescens, V. heterophylla, V. laurentii and V. racemosa). Also worth noting was the influence of soil variables on some models. For instance, exchangeable acidity (KCl) for 10 cm depth (EACKCL_T_M _sd2) was important for V. filicaulis and V. unguiculata subsp. baoulensis,; while soil texture fraction sand for 2.5 cm depth (SNDPPT_T_M_sd1), was among the key variables for V. laurentii, V. nigritia and V. oblongifolia. Another important result was the influence of the length of the driest month (LLDS) on V. comosa, where it was the top-most predictor variable (53.7 %).

However, the Jackknife test of variable importance (not included here) flagged the mean diurnal range in temperature as containing the most important information by itself and, therefore, explaining the gain in six models (V. filicaulis, V. frutescens, V. heterophylla, V. laurentii, V. oblongifolia, and V. racemosa). This was followed by isothermality for four models (V. longifolia, V. multinervis, V. unguiculata subsp. Baoulensis, and V. unguiculata var. spontanea) and precipitation of the driest quarter for three models (V. gracilis, V. luteola, and V. nigritia). On the other hand, isothermality appeared to have contained the most unique information that could not be found in all other variables for nine models (V. frutescens, V. heterophylla, V. longifolia, V. multinervis, V. oblongifolia, V. racemosa, V. reticulata, V. unguiculata subsp. Baoulensis, and V. unguiculata var. spontanea).

Suitable ranges of the three most important environmental variables for each model are represented by their respective response curves (Supplementary Figure 3). Generally, suitable habitats for the majority of the models were suggested highly likely to be in areas characterised by precipitation of the driest quarter (BIO 17) ranging between 5.5 and 133.5 mm; isothermality (BIO 3), with an average range of 64.3 to 80.4; and mean diurnal range in temperature (BIO 2), with an average range of 7.4–12.2°C. For the edaphic-related factors, suitable habitats were likely to be found in sites having a wide range of sandy soils, with an average soil texture fraction sand for 2.5 cm depth (SNDPPT_T_M_sd1) ranging between 17 and 90%, and somewhat acid soils, with an average exchangeable acidity (KCl) for 10 cm depth (EACKCL_T_M_sd2) below 0.6 cmolc/kg.

Figure 5 represents the results from the predicted potential distribution of the suitable habitat richness of the 15 modelled taxa under both the current and future climatic conditions, taking into account the PA network in Benin. Under the current conditions, the largest portion of the highly suitable area was predicted to be in the Sudano-Guinean Zone, with a few patches in the Sudanian and Guinean Zones. The current potential suitable habitat accounted for 39,647.88 km², representing ~33.69% of the land surface area of Benin. Out of this area, an estimated 3,931.72 km² (9.92%) fell in the existing PA network.

As for the future climate change scenarios, only about 27,107.87 km², representing about 23.03% of the land surface area was predicted potentially suitable under the moderate emission scenario (RCP 4.5). This represented a substantial loss (~31.63%) relative to the present conditions, mostly in the northerly located patches. Out of this suitable area, 2,718.18 km² (10.03% of the land surface area), was predicted to fall in the existing PA network. A somewhat different situation was projected under the severe climate change scenario (RCP 8.5). Here, the suitable habitat appeared to expand northwards, while the southern portion became less favourable. An estimated 48,933.01 km² (~41.57%) of the total land surface area was predicted to become potentially highly suitable, representing an increase of about 23.42% relative to the present conditions. From this area, about 5,322.57 km² (~10.88%) of the land surface area of Benin, was predicted to fall under the existing PA network.

The results from this study showed variations in taxa response to different environmental factors, thus supporting the first hypothesis that different Vigna CWR species would respond differentially to a different suite of environmental variables. Jarvis et al. (2008) reported a similar tendency among Vigna species. Several other studies have also reported similar generic tendencies such as in Cucurbita in Mexico (Lira et al., 2009), Piper and Oryza in Sri Lanka (Ratnayake et al., 2021), Vaccinium in the Netherlands (van Treuren et al., 2020) and Adansonia in Madagascar (Tagliari et al., 2021). Differential response of species to environmental factors underscores the need for species-specific studies on climate change effects (Jarvis et al., 2008).

While climate-related factors are generally the key environmental factors influencing species distribution (Jarvis et al., 2008; Amissah et al., 2014; Lompo et al., 2021), several studies have also shown that habitat distribution of species, especially at a local scale, is influenced by a combination of climatic and edaphic factors (Hageer et al., 2017; Idohou et al., 2017b; Zuquim et al., 2020; Assogba et al., 2022; Liu et al., 2022). The results from this study were consistent with these observations. The influence of edaphic factors appears to be particularly critical for narrowly distributed or understorey species (e.g., Hageer et al., 2017; Idohou et al., 2017b; Wang et al., 2019; Roe, 2020), and may be related to specialised habitat requirements by species that constrain their distribution (Corlett and Tomlinson, 2020). For instance, Idohou et al. (2017b) observed that the potential cultivable areas for the relatively localised and understorey palms such as Raphia hookeri and R. vinifera in Benin were more characterised by soil factors compared with overstorey palm species. Similarly, soil factors were the most discriminating factors in the distribution of an endemic orchid Spiranthes parksii (Navasota ladies' tresses) in the USA (Wang et al., 2019). Indeed, disregarding edaphic factors in SDM may overestimate future habitat adaptability of many plant species (Bertrand et al., 2012; Zuquim et al., 2020, but see Feng et al., 2020). Soil texture, for instance, is important for plant root development, especially for the relatively high-biomass rooted plants like most Vigna CWR species (Iseki et al., 2018), while exchangeable acidity (KCL) is crucial for soil nutrition and texture balance (Liebenberg et al., 2020).

The importance of precipitation of the driest quarter (BIO 17) to Vigna species as found in this study may relate to their reproductive fitness and adaptation to arid- and semi-arid areas. Given that reproduction and maturity in the modelled species largely coincide with the driest quarter, precipitation during this period might be important for gamete formation and viability, pod-set, and pod-filling (Nadeem et al., 2019). It may also play a role in seed dormancy (Smýkal et al., 2014). Precipitation of the driest quarter was also found to be one of the most important predictor variables for four legume species (Adenocarpus mannii, Afzelia bella, Afzelia bipindensis, and Baphia nitida) from the Nigeria—Cameroon border in West Africa (Salako et al., 2021) and for A. digitata in Benin (Assogba et al., 2022). Except for species that largely thrive in wet areas including V. laurentii, V. luteola, and V. multinervis, most Vigna CWR species are renowned for their persistence in arid areas, partly owing to their absorptive root systems or tuberous rootstocks (Maxted et al., 2004; Iseki et al., 2018). This might explain why the majority of species in this study demonstrated less demand for heavy precipitation. In Tibet in China, Xin et al. (2021) attributed the absence of excessive demand for precipitation in Sophora moorcroftiana and the distribution of this species in drought-prone areas to its strong absorptive root system.

On the other hand, large isothermality against the affinity for low temperatures as found in this study may, according to Zhang et al. (2018), suggest that the species use the relatively high temperatures during the day for photosynthesis while reserving energy at night through decreased respiration when temperatures are relatively low. Isothermality was suggested to be the second most important predictor factor for barbed goatgrass (Aegilops triuncialis) in Iran (Mousavi Kouhi and Erfanian, 2020), and its importance in shaping the distribution of plant taxa in the tropics has been widely reported (Amissah et al., 2014; Xin et al., 2021; Zuza et al., 2021).

These results provide an understanding of the key environmental factors that are shaping the distribution of wild Vigna species in Benin, and how changes in these factors as a result of future climate might affect the distribution of the studied taxa. This knowledge is critical in effective planning for adaptive management approaches of wild Vigna taxa in the face of climate change (Iriondo et al., 2021).

Several studies have demonstrated that climate change is causing species to shift their climatically suitable habitats towards higher latitudes and elevations (Aguirre-Gutiérrez et al., 2017; Lenoir et al., 2020). The observed northerly expansion in suitable habitats found in this study was thus consistent with these global findings. This partly confirmed the second hypothesis that there would be an increase in the northern habitat. However, it is also becoming apparent that some species may shift their distributions towards other directions (Tagliari et al., 2021; Balima et al., 2022). This was also the case in this study with some models like those of V. frutescens, V. heterophylla, V. racemosa, V. reticulata, and V. unguiculata var. spontanea that appeared to shift their suitable habitats eastwards (Figure 3).

The predicted patterns in the expansions and contractions in the current study might relate to an increase in the number of suitable habitat patches and a reduction in the sizes of those patches, respectively, as was reported by Jarvis et al. (2008). The discrepancies between this study and Jarvis et al. (2008), including the higher contraction rate of over 50% as observed by Jarvis et al. (2008), might be due to the differences in the sources of environmental data, types of species studied, the scale of the study, and the methods used. For example, Jarvis et al. (2008) did not incorporate edaphic variables in their study. Bertrand et al. (2012) suggested that edaphic factors may increase the tolerance of a species in confronting climate constraints, which could have been the case in the current study.

Climatic models have predicted that global climate change will lead to increased night temperatures and prolonged droughts, with West Africa being one of the most affected regions (IPCC, 2021). As a result, lower isothermality values and reduced precipitation of the driest quarter around the 2055s are anticipated than currently observed (Platts et al., 2015). The predicted contractions and expansions in the suitable habitats in the mid-Century observed in the studied species might partly reflect these changes.

The suggested minimum negative effects of climate change on the habitat suitability distribution of the seven species (V. frutescens, V. heterophylla, V. laurentii, V. nigritia, V. oblongifolia, V. racemosa, and V. reticulata) (Figures 3, 4) may be explained by several factors. First, the stability in the key predictor variables, for species like V. frutescens which appeared to be spatially confined to the north of Benin (Figure 1). Second, the relatively higher habitat heterogeneity that may provide a wider range of microhabitat options (Jarvis et al., 2008; Foden et al., 2019) for species such as V. heterophylla, V. racemosa, and V. reticulata. Indeed, these three species were among the most spatially distributed, spanning the different eco-geographical zones of Benin (Figure 1) and were predicted to have comparably wider habitat distributions under the present and future conditions (Figure 3). According to Hirst et al. (2017), common taxa as these within a clade are expected to perform relatively well across a wider range of novel environmental conditions than their rarer relatives. Lastly, it might be due the constraining influence of edaphic factors (Bertrand et al., 2012), especially on the narrowly-distributed species like V. laurentii, V. nigritia, and V. oblongifolia. These three species were among the least spatially distributed (Figure 1) and appear to have specialised habitats, with a predilection for marshy or seasonally inundated areas that have poor soils (Supplementary Table 1).

The habitat distribution maps generated in this study provide insights into potentially suitable habitats for the 15 Vigna taxa which may help in planning for taxa-targeted conservation measures including in situ and ex situ approaches. Further, they may aid in the promotion of cultivation of populations of some taxa to enhance their direct use by local communities so as to diversify the quantity and quality of the food basket, while acting as an incentive for the conservation of wild populations.

The area that was predicted to remain stable under future climatic conditions may be considered for conservation (Iriondo et al., 2021), since stability suggests the existence of favourable climatic conditions that may provide refugia for genetic diversity of species (Cobben et al., 2011). Given that the creation of new PAs solely for the conservation of CWRs may have huge cost implications (Maxted et al., 2016), and may also escalate the existing human-conservation conflicts over land (Tranquilli et al., 2014), the existing PA network together with its bordering landscape in the predicted stable sites is thus recommended for the conservation of Vigna species.

The suggested PAs in this regard include (not in the order of importance) Oueme Superior, Wari-Maro, Mont Kouffe, Agoua, Savalou, Oueme Boukou, and Dogo Forests [in the Sudano-Guinean Zone]; Pendjari National Park [Sudanian Zone]; and Lama, Ahozou (Pahou) and Drabo Gbo Forests [Guinean Zone]. However, it has been suggested that aligning these sites with hotspots of other plant species would help optimise conservation resources (Maxted et al., 2016; Vincent et al., 2022). In this respect, the relative high species richness and endemism in Pendjari National Park (Akoegninou et al., 2006; Neuenschwander et al., 2011) and its high legal protection status (Neuenschwander et al., 2011) would probably make conservation efforts within and around Pendjari National Park more capturing, less expensive, and practically manageable.

For effective conservation efforts in the suggested conservation sites, the following are recommended: species population monitoring, floristic inventories, habitat characterisation, genetic diversity, and ethnobotanical studies (Iriondo et al., 2008; Jarvis et al., 2015; FAO, 2017). Equally important are studies on the reproductive ecology, seedling recruitment, seed longevity, dispersal mechanisms, and responses to abiotic and biotic stresses of Vigna CWR species. Given that nearly 40 of the estimated 63 Africa Vigna species (54 out of an estimated 105 Vigna spp. globally) have presently been assessed using the IUCN Red List Categories and Criteria (https://www.iucnredlist.org) (accessed 02 June 2021), a risk assessment of the remaining 23 African Vigna spp. (51 spp. globally) is also recommended.

The favourable conditions in disturbed areas outside PAs where most CWRs have often persisted for long periods (Jarvis et al., 2015), coupled with slight increases in potentially suitable habitats in non-PAs (Figure 3) appeared to have rendered the PA network cover a comparably small habitat suitability area for the modelled species. This finding supported the third hypothesis that the existing PA network in Benin would be less effective in conserving Vigna CWR diversity than the non-PAs. This poses a threat to the conservation of these species, considering the increasing anthropogenic pressure outside the PAs (Guidigan et al., 2019), which, if taken into account, would substantially reduce the predicted suitable habitat (Riordan and Rundel, 2014). These results also corroborate previous reports of increased mismatches between climatically suitable habitats for CWRs and existing PAs, such as in Ethiopia (Davis et al., 2019), Mexico (Lira et al., 2009), the Netherlands (Aguirre-Gutiérrez et al., 2017), and Sri Lanka (Ratnayake et al., 2021). These results, therefore, buttress the calls for effective conservation of CWRs both within and outside PAs, using an integration of multiple approaches (Riordan and Nabhan, 2019; Iriondo et al., 2021) that takes into account “other effective area-based conservation measures” (OECMs) (IUCN, 2019), since PAs may not be the only best option (Goettsch et al., 2021). To achieve this, and for long-term monitoring and active management that involves participation of various stakeholders including local communities, a wide range of guidelines have long been made available (e.g., Iriondo et al., 2008, 2021; Maxted and Kell, 2009; Hunter and Heywood, 2011; FAO, 2017; IUCN, 2019).

The main caveat of this study is that the models represent only potentially suitable habitats for the modelled species. Distribution, access and persistence may be controlled by many other factors including anthropogenic such as land use, biotic interactions such as pollination, parasitism and diseases, and dispersal (Feng et al., 2020; Spicer et al., 2021). This study did not consider these factors, and, therefore, caution should be exercised when interpreting these results.

Although MaxEnt has proved to be a robust SDM tool in delineating habitat suitability maps for many taxa under the changing climate (e.g., Aguirre-Gutiérrez et al., 2017; Phillips et al., 2017; Ratnayake et al., 2021), the use of a single SDM algorithm does not provide for comparisons, and therefore accuracy in predictions may be put into question (Zuza et al., 2021). Future studies may consider a combination of models.

Another weakness of this study could be the use of publicly available occurrence data from the herbarium and/or online databases. It has been posited that such data may not randomly sample the true occurrences of species since such data is often biassed towards easily accessible areas such as roads and human settlements (Barlow et al., 2021). According to Barlow et al. (2021), such a spatial bias may present an over-representation of environmental conditions associated with regions of higher sampling effort. While this shortcoming is recognised, it should also be pointed out that scientific evidence is growing indicating that CWRs are often associated with such anthropogenically disturbed areas (Jarvis et al., 2015; Iriondo et al., 2021). As such, given that sites that have been accessible to species for a long time are ideal for species distribution modelling (Barve et al., 2011), it might as well be argued that data from such disturbed sites may provide a true reflection of the potential habitat distribution ranges of these taxa.

Nevertheless, this study has identified potential areas for prioritising conservation efforts for the 15 Vigna taxa in Benin. Further, it has raised a red flag for species that may need more attention, considering the predicted vulnerability of their potentially suitable habitats. The taxa include V. comosa, V. filicaulis, V. multinervis, V. unguiculata subsp. baoulensis, and V. unguiculata var. spontanea. Likewise, close attention ought to be paid to the predicted locally narrowly-distributed species, V. laurentii and V. oblongifolia. Although the suitable habitats for these two species were suggested to be negligibly affected by climate change, the species may fail to move and thus become predisposed to intense anthropogenic pressure (Qu et al., 2018; Corlett and Tomlinson, 2020; Spicer et al., 2021). Especially worrying about the two models of V. laurentii and V. oblongifolia is that their distribution ranges are localised at the southern tip of Benin, where land-use change is the greatest due to increasing urbanisation and agriculture (Guidigan et al., 2019). Moreover, V. laurentii is classified as Endangered on the IUCN Red List of species with a declining population (McFarlane and Maxted, 2019). Therefore, population monitoring of the seven species ought to be considered to inform the development of appropriate area-and species-based measures. Equally important are complimentary ex situ collections for long-term conservation of these taxa.