We focussed on the four endemic Commiphora species occurring in Socotra Island (C. ornifolia, C. parvifolia, C. planifrons, C. socotrana; see Figure 1). Location points of individual trees were taken in the field by hand using a GPS through a custom-designed application via mobile phone using ArcGIS Field Maps, which was also used for the field inventory of Dracaena in Oman (Vahalík et al., 2020). The island was investigated during field surveys in 2022 in all known areas of the distribution of the four species to cover the actual distribution, but also to cover new areas with the intention to improve our knowledge of the geographic distribution. During the surveys we encountered many individuals exhibiting morphological features intermediate between different species of Commiphora, including probable hybrids. These individuals were not included in the database to not introduce any bias in the model processing. Location points were then filtered by a grid of 50 m to avoid overprediction of the species distribution models (Mendes et al., 2020).

Since it is impossible to take absence points on the island (i.e., for example, we do not know if the recent extinction events such as extreme climate and overgrazing in Socotra have destroyed trees that were recently present), a number of pseudo-absence points equal to the number of available presence points were randomly generated in RStudio platform (R Core Team, 2022), to have balanced datasets (Barbet-Massin et al., 2012).

Different environmental layers were used for the modelling of the species distribution in RStudio platform (R Core Team, 2022). Topographic, bioclimatic, and edaphic variables were used as predictors of the potential distribution (Table 1). Altitude, slope and aspect have been assessed using a DTM at 30 m resolution from the SRTM (Farr et al., 2007) and for the climate, we used WorldClim V.2 bioclimatic variables at a resolution of 1 km (Fick and Hijmans, 2017). Finally, SoilGrid variables (Hengl et al., 2017) were used for the edaphic conditions. While we know that these are global data layers and local differences are possible (e.g., local microclimates), these are currently the best variables available (e.g., there is no detailed soil map for Socotra). Several of these layers were used in recent studies to update land classes, understand tree cover, potential distribution, and productivity of (tree) vegetation (Rezende et al., 2020, 2022; Riccardi et al., 2020).

We harmonised the different resolutions of the environmental layers with a bicubic interpolation neighbours using the software Google Earth Engine (Gorelick et al., 2017). It creates high-resolution layers and allows the selection of the final resolution. We chose a 50 m resolution corresponding to the cell size used to filter GPS points of Commiphora species.

Multicollinearity is a problem for the species distribution models since the environmental variables should be independent in order not to weaken the model’s statistical power. For this reason, to detect any correlation between the variables, a Variance Inflation Factor (VIF) test was conducted for each species on RStudio (R Core Team, 2022) using the “usdm” package (Naimi et al., 2014). This step excludes the highly correlated variables through stepwise processing. All the environmental variables with a threshold ≤5 were selected for the elaboration of the model (Table 1).

The CMCC-CESM2 from CMIP6 (Coupled Model Intercomparison Project) was downloaded from WorldClim (Fick and Hijmans, 2017) and used as a global climate model since it works more accurately than others in this geographic area, considering the effects of the winter and summer monsoons (Mohan and Bhaskaran, 2019). The same bioclimatic, topographic and edaphic variables (Table 1) were employed for each species to train the model, jointly with the layers harmonisation conducted in Google Earth Engine software mentioned before.

We compare two different standard scenarios, the SSP126 (Sustainable and “green” pathway with 2.6 W/m² radiative force by the year 2100) and SSP585 (Fossil-fuelled Development with 8.5 W/m² radiative force by the year 2100). The selected time range is 2081–2100 (mean 2090).

All the following analyses were performed on RStudio (R Core Team, 2022) using the package “sdm” (Naimi and Araújo, 2016). The datasets were split into 70% for training and 30% for testing the species distribution model. For each dataset of species occurrences, a distribution model was produced. Three modelling approaches–GAM (Hastie and Tibshirani, 1990), Random Forest (Breiman, 2001), and MaxEnt (Phillips et al., 2006)–were used, and the results were averaged in an ensemble model (Araujo and New, 2007). Sometimes using single models based on different assumptions, makes it difficult to choose which algorithm performs best. Therefore, ensemble models have been used to combine and average predictions in different ways (Thuiller et al., 2009) but also to provide information about the overall trend (Guisan et al., 2017; Thuiller et al., 2019). The models were run with 25 replicates for each method (for a total of 75 replicates), using the “subsampling” partitioning method. Before doing it, a seed (9505) was set to ensure the analysis’s replicability. The variable importance and response curves were elaborated with the same R package.

The Area Under the Receiver Operating Characteristic Curve (AUC) was used as a model performance indicator. When the AUC values were <0.6, 0.6–0.7, 0.7–0.8, 0.8–0.9, or 0.9–1.0, the predictions were considered as invalid, poor, fair, good, or excellent, respectively (Swets, 1988).

Instead of using the arithmetic mean or median, the models were ensembled with a weighted average, in particular applying the predictive performance of the True Skill Statistic (TSS), a measure of omission and commission errors and among the other accuracy measures is the more realistic and practical (Shabani et al., 2018). Finally, we selected an occurrence threshold through the sensitivity-specificity sum maximisation on the ensemble models (Manel et al., 2001) to create the Habitat Suitability maps. To have an idea of the potential distribution of the species in the future and under different climate change scenarios, the ensemble model was run with the global climate model chosen using the same parameters for the ensemble model as before.

The results achieved through this work, together with additional threats assessments during the surveys (unpublished data), will allow us to reassess the IUCN conservation status of the four endemic Socotran Commiphora species according to the Guidelines for Using the IUCN Red List Categories and Criteria (IUCN Standards and Petitions Committee, 2022). Since we only discuss distributional data in this study, the species are assessed for now by criterion B (yet only partly); this is only a first step and preliminary because we do not include the exact population estimates and fragmentation of the species in this study (which may alter the final IUCN Red List reassessment). The calculation of Extent of Occupancy (EOO) and Area of Occurrence (AOO) was made with the “redlistR” package in RStudio (Lee et al., 2019).

The data collected during the fieldwork, allowed us to explore the actual (current) distribution of the four endemic species of Commiphora on the island. We collected 345 occurrences of C. ornifolia, and, as shown in Figure 2, it is widely distributed with a preference for the central and eastern parts. Similarly, C. parvifolia seems to prefer the central and north-eastern parts with few isolated populations on the west and for a total of 279 individuals recorded. C. planifrons records (208 occurrences) are more frequent in mountainous areas of the island, especially in the Haggeher Mountains, however, there are many records of the species on the high plains in the western part (Ma’alah and Qatariyah). The species has no geographical preference, but it is present at higher altitudes, with some exceptions in the northwest part where it is present between 200 and 300 m a.s.l. Finally, C. socotrana is concentrated in the northern part of the island, with a relatively lower number of occurrences (179 records), and one southern population, in the Zahr Basin.

Averaging all the results from each model (Figure 3), it turned out that the most important variable explaining the distribution potentiality is bio03 (isothermality; Table 1) which is shared between three species (C. ornifolia, C. parvifolia, and C. socotrana) with a relative importance of 19.7, 42.9, and 45.1%, respectively. The second shared variable is bio13 (precipitation of wettest month) which turned out to be important for C. ornifolia (23.1%), C. parvifolia (11.6%), and C. planifrons (13%).

Other important predictors of the actual species distribution are: bio09 (mean temperature driest quarter) for C. planifrons (46%), bio07 (temperature annual range) for C. socotrana (15.7%) and C. parvifolia (14%). Regarding soil variables, clay plays a role for three species, C. ornifolia (13.2%), C. parvifolia (11.6%) and C. planifrons (11.7%).

The potential optimal ecological conditions for each species are shown in Figure 4. The occurrence of C. ornifolia is related to clay, bio03 (isothermality) and bio13 (precipitation of wettest month). The occurrence probability increases when the amount of clay in the soil is higher than 200 g/kg, the isothermality is lower than 60% and the precipitation of the wettest month is between 25 and 30 mm (Figure 4A). Regarding C. parvifolia (Figure 4B), we have the same response to the clay variable, with a higher occurrence when the clay is higher than 200 g/kg. The occurrence probability sharply increases where bio13 (precipitation of wettest month) is 22 mm and when bio07 (temperature annual range) reaches values between 10.5°C and 12.5°C. In C. planifrons (Figure 4C) the optimum conditions occur when bio09 (mean temperature of driest quarter) is 21°C and bio13 (precipitation of wettest month) is higher than 50 mm, as well as the amount of clay in the soil which increases the occurrence when it is 300 g/kg. Finally, C. socotrana (Figure 4D) shows an increase in occurrence when bio03 (isothermality) and bio07 (temperature annual range) are 55% and 13°C, respectively. On the other hand, the optimum for bio15 (precipitation seasonality) is between 60 and 70%.

The potential distribution (niche modelling according to the current environmental variables), i.e., all suitable areas for the actual, current occurrence of the four species, is shown in Figure 5. From the actual records and potential distribution of C. ornifolia, we found that this species occupies mostly wadi/canyon areas (Figure 5). However, is not unusual to find the latter species in elevated zones, on high slopes. There is only an exception in the southern part of the island where it is recorded on the coastal plains of Noged (Figure 2). In total the potential calculated distribution of the species is 751.2 km² and it currently occupies 23% of the suitable areas.

Commiphora planifrons can be found in different areas (Figure 2). The model extends the potential distribution in the eastern and western areas (Figure 5), following the topographic conditions needed by the species, reaching a potential coverage on the island of 611.8 km² of which 28% is currently occupied.

The last two species share a similar distribution, mostly concentrated on the northern side of the island. Nevertheless, there are differences since records of C. parvifolia can be found in the south and C. socotrana in the central plain of the island (Figure 2). Concerning the extent of potential distribution, they cover 494.1 and 271.1 km² respectively (Table 3), while the actual distribution cover only 17% for C. parvifolia and 33% for C. socotrana which has the highest percentage.

Relating to the performance of the species distribution models, they have been classified as very good for C. planifrons and C. socotrana and excellent for C. ornifolia and C. parvifolia with a mean testing AUCs of 0.859, 0.870, 0.912, and 0.959, respectively (Table 2 and Supplementary Table).

Among the four species, only two show a clear reduction of their distribution in the future in both climate scenarios (SSP126, SSP585; Table 3). Commiphora socotrana is the only species that could potentially disappear in 2090 in a SSP585 scenario, and in the same case, C. planifrons will encounter a remarkably high reduction of its potential distribution (−97.9%).

The situation for these two species does not change in the “better” scenario (SSP126) since habitat suitability is forecasted to decrease extremely, by −97.2% for C. socotrana and −92.6% for C. planifrons (Table 3).

The species least affected by the climate change scenarios seems to be C. parvifolia. In both scenarios, the species shows an actual increase in its potential distribution, more in the worst-case scenario where it could potentially reach an increase of 108.4% by 2090 (Table 3).

Commiphora ornifolia shows a different behaviour in the two scenarios (Table 3). In the “better” one (SSP126), the potential distribution will increase by 45.4% by 2090 whilst in the worst case (SSP585) it will slightly decrease by −7.8%.

The actual and expected Areas Of Occupancy and Extent Of Occurrence (AOO and EOO) are shown in Table 3, while the maps of EOO are in the Supplementary Figures 1–11. Notes on the assessment of the species (only using criterion B) suggests a change in the IUCN Red List status for at least two species (Table 4). Commiphora planifrons and C. socotrana could be listed as “Vulnerable” while C. ornifolia and C. parvifolia would remain in the same category (Near Threatened). However, a future reassessment is still likely to change because all these species have, besides limited occurrence, strongly fragmented subpopulations with severe risks of decline due to lack of regeneration (because of overgrazing). A more detailed threats assessment and detailed population size results are part of ongoing studies.

In our work, we explored the current and potential distribution of four endemic Commiphora species of Socotra Island using an ensemble of three modelling approaches. This work confirms the predictive power of the method in investigating the potential distribution of species. Ensemble models are shown to have a better predictive performance compared to single models, reducing the uncertainty and the inevitable problems of over-fitting and under-fitting (Marmion et al., 2009; Battini et al., 2019). The potential distribution of this genus has been analysed only by Kulloli and Kumar (2014) for C. wightii in Rajasthan (India), however, no one has covered the potential future impact of climate change for other species of the genus. Despite this, the impact has been studied and highlighted in other similar works with focus on the trees of Socotra. Many of these studies are limited to the application of a single SDM method, however, yielding useful results (Attorre et al., 2007; Rezende et al., 2022).

Highlighting the ecological factors which shape the potential distribution, the species in our study have some similarities with the Indian C. wightii only in terms of seasonal factors, such as the precipitation seasonality, which it has in common with C. socotrana. The Socotran and Indian species likely have a different phytogeographical history, even though the section to which they belong originated at the beginning of the Miocene (Gostel et al., 2016). The island of Socotra was in contact with the Arabian Peninsula before rifting started in Oligo-Miocene times, and it is possible that Commiphora species colonised before this (Brown and Mies, 2013). Commiphora wightii probably arrived in India through the Arabian Peninsula being a potential divergence area (Lakhanpal, 1970). However, the ecological characteristics of the endemic Commiphora of Socotra can hardly be compared with species in other continents, because of the specific local conditions and potential local radiation in Socotra, which is often the case for island endemics (unusual morphologies and/or ecological characteristics), and because of the wide range of niches in which Commiphora may occur worldwide.

Comparing the current and potential distribution of the Socotran Commiphora according to their ecology, we found that these four species occur in many of the areas predicted as currently suitable, shown by the high values of the sensitivity index. Another interesting aspect is how the prediction for C. socotrana conserves a prevalent northern distribution, whilst C. ornifolia does not seem to prefer lowland areas, although in the south there is a significant population of this species, as shown by the current distribution. This population could be the remnant of an ancient lowland forest, yet there is also a possibility that some trees may have been planted by local communities in previous times because of its ethnobotanical value. Population genetic analyses could help to understand if the distribution is entirely natural. Thinking about the differences between the actual and potential distribution, the concept of realised and fundamental niche comes to mind (Peterson, 2006). The actual distribution is how the species interacted with biotic interactions and/or geographical/historical constraints. Consequently, the extent is inferior to the potential distribution, which corresponds to the fundamental niche.

Predictions of actual vs. potential niche occupancy have been made for other endemic trees of Socotra, for example the Dragon’s Blood Tree Dracaena cinnabari (Attorre et al., 2007; Rezende et al., 2022). However, there is a clear difference between the wider actual occupancy of Commiphora vs. Dracaena, the latter covering only 5.5% of its potential distribution in Socotra Island (Attorre et al., 2007), despite the fact that both genera are facing similar local pressures of overgrazing and climate impacts, as well as being used for resin extraction (Maděra et al., 2019a). This is likely primarily a result of a difference in biology and age-structure of the populations of these endemic trees. Dracaena cinnabari is extremely slow-growing (Maděra et al., 2019a) and the populations of Commiphora could be relatively younger. In general, suitable areas were studied for trees as well as reptiles in Socotra using similar ecological niche modelling (Attorre et al., 2013; Fasola et al., 2020).

Our study reveals key details about the Socotran Commiphora species, as we identified their optimal ecological conditions for the first time. Whereas the climatic variables provide us with information about the optimal temperature and precipitation ranges for a given species, the edaphic ones tell us about species ecology. We found that relatively higher levels of clay in the soil correspond to a higher probability of occurrence of C. ornifolia, C. parvifolia, and C. planifrons. Clay is an important soil component because it helps the hydric retention and nutrients release through the cation exchange for the plants (Rakshit et al., 2021). Unfortunately, the ongoing soil degradation in Socotra due to overgrazing does not allow the formation of stable vegetation communities (Pietsch and Morris, 2010; Van Damme and Banfield, 2011) and leads to more erosion and a loss of potentially suitable areas for these species in the future, threatening them further with extinction. Suitable soil substrates are not only crucial for Commiphora, but also for other endemics living in Socotra’s terrestrial ecosystems, such as the Boswellia elongata woodlands, which are undergoing a degradation phase as shown in some areas (Lvončík et al., 2020), and the reptiles, where types of substrates and vegetation play a crucial role in their ecology (Fasola et al., 2020). Such detailed ecological information and identifying suitable areas has been used in other species to identify proper management, for example in reptiles where a large part of the endemic endangered species and vulnerable populations were not covered by the areas with highest protection in Socotra (nature sanctuaries) (Vasconcelos et al., 2018). With the current information, this can now be also evaluated for the Socotran Burseraceae in the future (reserve design).

Insular ecosystems are more vulnerable to the effects of climate change (Duvat et al., 2017; Cartwright, 2019) even more if other pressures such as invasive species and land use change are added (Reaser et al., 2007; Shao et al., 2017). This is why more studies are urgently needed to understand the potential effects of climate change on the endangered tree species endemic to Socotra. The two different (SSP) scenarios which we examined give us some idea about the reaction of the species in the next 60 years of climate effects, taking into account the current ecological needs of these taxa to observe how potential distribution changes in the future. A similar study was conducted in the New Caledonia Archipelago, a biodiversity hotspot rich in endemics, where changes in potential distribution of tree species were assessed (Pouteau and Birnbaum, 2016).

Commiphora ornifolia seems to be relatively less vulnerable to predicted climate change impacts since the future potential habitat area will actually increase by over 45.4% until 2090 according to SSP126, whilst according to SSP585 the species’ habitat will decrease by just 8%. In any case, the actual distribution of the species on the island is relatively large since the local communities still use it as a food source (e.g., due to fruit properties to lessen thirst) and several local communities have replanted it in their villages (Miller and Morris, 2004). On the contrary, C. parvifolia seems to be favoured by in a general climate change scenario with an increase in the future potential habitat suitability by 2090 of 55 and 108% for SSP126 and SSP585, respectively. This enlargement of the distribution may be due to the species’ tolerance for low precipitation regimes (optimum precipitation during the wettest month is 22 mm), driving it to resist and thrive under changing climate effects. The last two species, C. planifrons and C. socotrana share a similar fate, with a decrease in the potential distribution in both SSP scenarios, making them the relatively more vulnerable Commiphora species on the island. The first one clearly shows a future likely upward shift to the higher parts of the island, in search of higher precipitation and colder temperature during the dry season, disappearing in many areas of the island with a decrease of its distribution of 98% in the worst scenario (Figure 6).

The movement of species to higher altitudes and latitudes as refuge areas is recorded worldwide. It is widespread in mountainous areas (Parmesan, 2006; Engler et al., 2011; Zu et al., 2021) and in insular ecosystems (Pouteau and Birnbaum, 2016), where it could lead the species to extinction since they cannot find other suitable mountainous areas due to the topographic constraints. Even though fluctuations are normal, our climate is now changing faster than before (Jansen et al., 2007). Species moving in newly suitable areas sometimes may create “climate refuges” (Hegazy and Lovett-Doust, 2016), places where isolated populations of a migrating species remain stuck in small areas with suitable conditions but surrounded by unsuitable ones. In fact, the mountains of Socotra are already “wet refugia” for many plants, which are thought to have had a wider distribution in the past (Miller and Morris, 2004). This leads us to think that the populations of C. planifrons in Socotra’s central granite mountains (the Haggeher mountains) could further become a relict population (and even more endangered than now).

On the other hand, C. socotrana could probably face a different path. On Socotra Island, the species is not so widespread at present, but it forms populations with a high density of individuals. Despite this, the species potentially will decrease by 2090 of 100% of its potential distribution (Figure 6H), indicating that the species is at a high level of extinction risk. It seems to be affected by seasonal factors, while a major role is played by the winter and summer monsoons that hit the island, causing different rainfall and temperature regimes between the northern and the southern part (Scholte and De Geest, 2010; Brown and Mies, 2013). The possible extinction of C. socotrana on the island by 2090 could be related to the strengthening of the winter monsoon due to the effects of the Arabian Sea warming on the atmospheric circulation, which will increase extreme events, e.g., cloudbursts, high-speed winds, etc. (Mishra et al., 2022).

However, all the modelled potential distributions under different climate scenarios are based only on climatic, edaphic and topographic parameters, which are absolutely fundamental for the ecological niche of a species (Franklin, 2010; Wan et al., 2021), but they do not take into account the availability of suitable habitats and the current and future anthropogenic pressure on these species. On Socotra, these two aspects are strictly interdependent. In fact, thousands of years of humans living there, with introduced livestock, have definitely changed the pristine habitats (Scholte and De Geest, 2010; Van Damme and Banfield, 2011). What we are seeing now is just the result of a prolonged nature disturbance activity (Van Damme and Banfield, 2011). At present, climate, soil, temperature, and altitude are stronger predictions of forest landscapes in Socotra, relatively more than anthropogenic effects (Riccardi et al., 2020), however, from suitable niche modelling it is clear that several tree species do not cover their full potential and are relictual populations already, e.g., Dracaena (Rezende et al., 2022). In addition, our applied climate prediction models do not consider the rapid degradation of woodlands and forests in Socotra due to cyclones (Lvončík et al., 2020), the effects of overgrazing, increase of landslides and erosion in the highlands and slopes (Rezende et al., 2020), and economic-driven threats such as wood collection (Maděra and Van Damme, 2022). Therefore, our predictions with the current models could be under-estimations: considering the likelihood of more cyclones, and ongoing trends of soil degradation, overgrazing and woodland degradation, as well as the lack of regeneration, the rate of potential decline of some very vulnerable species (e.g., Commiphora socotrana) is likely much faster. Such species should be not seen separately from their plant communities, but as local indicator species in Socotra to understand the threats for entire habitats.

Many areas are still in good conservation condition but are threatened by the grazing activity of goats (Van Damme and Banfield, 2011; Brown and Mies, 2013). Overgrazing poses under pressure the plant communities in several ways. One is the change in plant composition since goats eat only the palatable plants, with no chance of survival, highly stressed and with development problems (Trlica and Rittenhouse, 1993; Herrero-Juregui and Oesterheld, 2018), whilst the others are totally untouched. Another issue is the lack of natural regeneration from seeds, due to the grazing pressure which leads the goats to eat whatever plants they encounter in particular small seedlings, except for those situated in inaccessible areas or in a protected location (e.g., rocky slopes, surrounded by stones), as seen for Boswellia spp. (Lvončík et al., 2020; Hamdiah et al., 2022) or Dracaena cinnabari (Habrová and Pavlis, 2017; Maděra et al., 2019a). Ultimately, an indirect effect of overgrazing is land degradation and soil erosion. The lack of plant regeneration due to the high number of goats on the island, together with drought periods, exacerbate the soil erosion through strong wind and surface runoff, leading to a loss of suitable habitats for Commiphora spp. and other taxa (Brown and Mies, 2013).

In general, the new insights of current and potential distribution of the endemic Commiphora species of Socotra Island provided in our study, as well as a prediction of changes under climate change scenarios, can be used to design specific species management plans for the future and strategize direct conservation actions, prioritising the most vulnerable species. The ethnobotanical uses of these species locally (Miller and Morris, 2004) is a starting point to work with indigenous communities for their conservation. With the current information on specific ecology and sensitivity of Socotra’s endemic Commiphora, and suitable areas for the future, well-informed conservation actions can be taken to protect these culturally relevant insular species.