Bioclimatic habitat limitations for argan trees (Argania spinosa (L.) Skeels) in Northern Africa and Spain

Argania spinosa L. Skeels is an Algerian-Moroccan endemic tree. This species is part of various plant communities consisting of Mediterranean, Macaronesian and Saharan floristic elements. It has been introduced and perhaps sometimes naturalized in various regions of the Mediterranean basin. Due to its role in combating desertification, high socio-economic value, and traditional use as fodder and food, the southwestern Moroccan argan grove (Arganeraie) was declared Biosphere Reserve. It had already been subject to conservation and reforestation programs a century earlier. Its cultivation for oil production could be, besides an economic objective, an effective method to conserve its genetic diversity. Therefore, this study aims to estimate its potential distribution and establish efficient breeding programs by determining its ecological requirements, identifying its different habitats, and predicting habitat suitability models for Morocco, Algeria, Tunisia, and Spain. Using 53 occurrence points, wind speed and direction data, and 29 bioclimatic variables, multivariate methods were applied to describe the ecological profiles and characterize the heterogeneity of its habitat to subsequently, train a Maxent model that establishes, besides Morocco and Algeria, suitable cultivation areas in Tunisia and Spain. The North African potential area is limited to the western Mediterranean coast of Algeria and flat and coastal areas of eastern Tunisia. The increased likelihood of suitability remains in the southeast Iberian Peninsula. A high probability of argan cultivation is also evident in the Canary Islands. These results provide possibilities for future expansion of argan crop and a window of opportunity to improve its genetic diversity and conservation.


Introduction
Argania spinosa L. Skeels (hereafter argan) is an Algerian-Moroccan endemic tree. It has been introduced and locally naturalized mainly in Tunisia, Israel, and Spain (Nerd et al. 1994;Mezghenni et al. 2014;Martínez-Gómez et al. 2018Louati et al. 2019a). The species is well known for the high socio-economic impact that its traditional exploitation has in large areas of Morocco by the extraction and use of oil from its seeds, wood, and other purposes (Benaouf et al. 2020). These knowledge, techniques and practices related to argan were recognized as Intangible Cultural Heritage of Humanity on 27 November 2014 by UNESCO. In some parts of Morocco, due to its ethnobotanical uses, argan takes the role of the olive tree as a multipurpose specie, since it is used as fodder, wood, and even human food. For example, argan is the fuel of the Berber society, especially near Essaouira and, when food is scarce, goats climb the tree to forage. In 1998, UNESCO declared an area of 2 568 780 ha of the argan grove (Arganeraie) as a Biosphere Reserve situated in the Souss-Massa National Park, in southwestern Morocco near Agadir (Ait Aabd et al. 2019).

Communicated by José Valentin Roces-Diaz.
This native tree is more widely spread in southwestern Morocco (Fig. 1), from the Saharan littoral plain to the valleys formed by oueds (wadis) draining the Great Atlas and the Anti-Atlas and disjoints with small populations growing more than 700 km away from it in northern Morocco and Algeria (Benaouf et al. 2014;Faouzi et al. 2015). This dispersion pattern can be interpreted because of a habitat reduction or losing process. There are also ideas based on the molecular polymorphism of the species, explained by the human action that would have dispersed the species for its economic and ethnobotanical interest (El Mousadik and Petit 1996). Despite the open debate about the possible fragmentation of its natural distribution area, different sources agree on the subspontaneous nature of these populations (Msanda et al. 2005;López-Sáez and Alba-Sánchez 2009;Kechairi and Benmahioul 2019). Some localities may even respond to ancient cultivation attempts. In fact, there is evidence of this dispersion of argan in the Iberian Peninsula where it was cultivated during the historical period of Al-Andalus (VII-XV centuries) (Carabaza-Bravo et al. 2004). There are other extensive regions in Argelia (northwest NW of the Tindouf region) where the argan grove is found in a spontaneous state (Kechairi et al. 2018).
The analysis of the plant communities (Boudy 1950;Tregubov 1963;Benabid 1985) in which argan appears evidence firstly its paleo-relict character from the late Tertiary vegetation, currently very fragmented and preserved in the Macaronesian flora, and secondly, its paleo-tropical character like other species sheltered in the Mediterranean littoral of Tyrrhenian area.
On the other hand, argan was probably part of the woodlands that condensed the oceanic fog as nowadays the Laurel Forest does, and it still preserves that fog-harvesting ability, like other species that appear in argan formations such as Pistacia lentiscus, Olea maroccana, Whitania frutescens or Buxus balearica according to Ritter et al. (2008) and Navarro-Cerrillo et al. (2011). Fog and dew regimes have been studied as a potential source of water in areas where argan is distributed (Font-Tullot 1949;Marzol et al. 2011;Lekouch et al. 2011), demonstrating the importance of the contributions of humidity through dew and fog as a complement to the scarce rainfall. This happens especially in the southwestern arid and semiarid Morocco, where these two humidity sources have a significant role in different periods.
Since the end of the 1990s, numerous studies have attempted to explain the A. spinosa biological and ecological aspects to improve the conservation of its genetic diversity. The allogamous nature of argan with a small percentage of self-compatibility justifies a geographical model of distribution of genetic diversity that shows greater values of intrapopulation than inter-population variability (El Mousadik and Petit 1996;Bani-Aameur 2004;El Bahloul et al. 2014;Mouhaddab et al. 2015Mouhaddab et al. , 2016Yatrib et al. 2017;Ait Aabd et al. 2019) and some independence from environmental or physiographic factors such as distance between populations or proximity to the coastline. At the same time, this variability level is also a complex problem when establishing propagation, cultivation, and forest management mechanisms using sexual propagation (Mouhaddab et al. 2017).
The argan tree grows in areas of varied orography; however, very little attention has been paid to its ecological range and the heterogeneity of its habitat (Kechairi and Abdoun 2016;Ait Aabd et al. 2019). This lack of information on the bioclimatic and physiographic factors that could determine argan presence throughout its natural range is extended to the suitability assessments aimed to determine the optimal areas to either afforestation or crop introduction. This is an intention already pointed out by López-Sáez and Alba-Sánchez (2009), who demanded the development of habitat suitability models. They also highlight the need to better understand the ecophysiological behaviour of this plant (but see Díaz-Barradas et al. 2010), especially when assessing its mechanisms against drought. Some of them suggested by Peltier et al. (1990) and Msanda et al. (2005) when referring to argan as an arid manifestation of "fog forests". Potential growing areas have been attempted to investigate in Argentina (Falasca et al. 2018) managing environment variables but have not dealt with the application of predictive suitability models.
A key aspect in Conservation Biology programs is the habitat heterogeneity studies (Holtmann et al. 2019) that help to answer how variable the habitat of a species is, the diversity of its elements and the complexity of its relationships. During decades and before the possibility of studying the genetic diversity of plant species at molecular level, forest science has used the analysis of the habitat diversity as an estimating parameter of the species genetic diversity. Nowadays, the habitat heterogeneity is incorporated to Geographical Information Systems (GIS), when establishing forest management models aimed to conservation, economic exploitation or sustainable use of species and landscapes applying ecogeographic variables. These models are very useful when choosing climatic refuges and establishing the areas and localities where its cultivation or introduction can be more efficient or carried out with better guarantees of success (Navarro-Cerrillo et al. 2011;Sharma et al. 2018;Pratumchart et al. 2019). Thus, the main purpose of this study is to analyse the habitat diversity of argan to estimate its potential distribution and establish efficiency levels of potential breeding programs, pursuing three specific aims: (1) to determine its ecological requirements, (2) to assess habitat heterogeneity, and finally (3) to predict its habitat suitability by using the species distribution model (SDM) algorithm MaxEnt for its native countries (i.e. Morocco and Algeria), plus Tunisia and Spain, where cultivation for farming purposes may expand its genetic diversity and improve its conservation status.

Considered area
The area considered in this study extends over Morocco, Algeria, Tunisia, and Spain, between 15° N and 45° N of latitude and 20° W and 15° E of longitude ( Fig. 1). There are arid and semi-arid ecosystems with Atlantic influence between 29º and 32º N in southwestern Morocco. Inland, the larger populations occur along the Oued Souss on the S-side of the High Atlas, and on the N and S of the Anti-Atlas where it reaches 1 300 to 1 500 m, while the disjoint populations are influenced by the Mediterranean climate. In western Argelia, it sometimes appears in microhabitats of certain humidity (usually at the bottom of the oueds).
Considering this evidence and to establish the potential area of argan cultivation, the occurrence points in its native area were used to explore the ecological requirements and habitat heterogeneity. Moreover, there are other locations where it has also been introduced and naturalized but never only cultivated (López-González 2001; Ruiz de la Torre 2006). For this reason, samples from Spain (Murcia and Alicante provinces), Algeria (Mostaganem, and Mascara provinces) and Morocco (Berkane, Agadir, and Taroudant provinces) were included to make the predictive suitability analysis more robust.

Presence data
Twenty-seven records were downloaded from the Global Biodiversity Information Facility (GBIF), an open-access network that includes occurrence records of any species around the world (https:// www. gbif. org/ es/). This database was completed by adding coordinates obtained from a bibliographic search in Morocco (Díaz-Barradas et al., 2010;Faouzi et al., 2015;Sjölund et al., 2019), Algeria (Benaouf et al. 2015;Kechairi and Benmahioul 2019), and Spain (Justamante et al. 2017). All the occurrence site records represent the natural and naturalized distribution of A. spinosa in The WGS84 coordinate system (see Online Resource 1).

Selecting environmental data
Nineteen bioclimatic variables and the elevation model (SRTM) were obtained from WorldClim data portal (Fick and Hijmans 2017). The SRTM was used to calculate three topographic variables: slopes, hillshade and aspect. In addition, two indicators of aridity were included: The Global Aridity Index and the Potential Evapotranspiration (Trabucco and Zomer 2019). The NDVI is a vegetation index that serves as an indicator of the contrast between vegetation and bare soil. To include the state of the earth's surface where argan is found, four satellite images of NDVI and Land Surface Temperature during April and July 2002 were extracted from the standard MODIS Land products that combine data records from Terra and Aqua satellites (https:// daac. ornl. gov/). Daily wind speed and direction data from the Global Surface Summary of the Day (GSOD) stations were download from the 1973-1999 period using the Weather Data Client of GSOD in R (GSODR) package (Sparks et al. 2017). In total, the 3 stations located at the closest airports to the study area were used: Safi, Essaouira, and Sidi Ifni (Fig. 1).
The twenty-seven climatic and environmental variables were compiled in raster (grid) format with 1-km 2 spatial resolution. The data file was generated using the point sampling tool QGIS plugin.
Additionally, the argan preference for basic soils is known so (https:// pfaf. org/), although soil pH has not been considered for statistical analysis, the global gridded of pH soil information at 0-5 cm depth has been included as a restrictive layer to refine the results of the suitability model.
The complete environmental dataset was classified into seven categories (Table 1) and was managed and overlapped using the Quantum Geographic Information System (QGIS) and MaxEnt (Phillips et al. 2006).

Ecological requirements
The frequency distribution of the sampling points (n = 51) was studied in the environmental variables of the blocks of temperature, precipitation, topography, and dryness, to characterize the values that represent the ecological optimum for the species and its ecological amplitude concerning each environmental variable (Gounot 1969).
The ecological profiles display the distribution of the values of each variable and reveal their mean values, variance (ecological amplitude), and possible adjustments to gaussian, bimodal, or exponential distribution curves. To compute and illustrate the histograms of the observations, the generic R function hist was used. Thus, the height of each rectangle is proportional to the number of points with the same value (Dalgaard 2008).
The role of the wind is essential in determining the presence of fog, and it is necessary to consider both direction and speed. Due to the nature of the wind data, the time Average function of the openair R package (Carslaw and Ropkins 2012) was applied to calculate the monthly and annual average value of wind speed (u i ) and orientation ( i ) at each station. This function calculates the east-west (Vni= ui * Cos i) and north-south ( Vei = ui * Sin i) components of all observations (N), then gets an average of each and direction ( tan = Ve Vn ). To convert the vector direction (towards which the wind is blowing) into the meteorological direction (from which the wind is blowing), add 180º to the direction if Ɵ < 180º or subtract 180º if Ɵ > 180º.

Habitat heterogeneity
In order to be able to apply the different multivariate analyses using data that have different units of measurement, a Z-score method (the difference between the input value and the mean is divided by the standard deviation of each variable) has been used.
Kmeans Clustering was applied as a non-hierarchical clustering technique to establish groups from the partition of the set of n = 51 observations in k groups in which each observation belongs to the group whose mean value is closer. It was calculated using the kmeans function in the stats R package (Hartigan and Wong 1979). It began by creating k random clusters, then the mean of each group and the distances of each individual. The individuals were placed within new resulting groups with the minimum intraspecific variance and maximum interspecific variance (members most similar to each other and the most different from those of other classes).
Then, a principal component analysis was applied to summarize the large number of raster layers or descriptors and easily understand their role in the formation of the groups. The pre-selection of the descriptors followed two criteria: (a) More contribution of each one to the new dimension and (b) less correlation between them.
The Partial Least Squares Discriminant Analysis (PLS-DA) is used to determine heterogeneity in environmental systems (Kwon et al. 2014). The function plsda of the R package mixOmics v2.7-1 (Gromski et al. 2015) was applied to obtain the Variable Importance in Projection (VIP) score and select the variables with the higher discriminating power (VIP > 1 criterion).
In addition to this, the R 2 (measures the degree of fit of data) and Q 2 (measures the predictive power) values describe the quality of the PLS-DA model. Q 2 is an estimation of the predictive ability of the model and is calculated via cross-validation (CV). In each CV, the predicted data are compared with the original data, and the sum of squared errors is calculated. The prediction error is then summed over all samples (Predicted Residual Sum of Squares or PRESS). For convenience, the PRESS is divided by the initial sum of squares and subtracted from 1 to resemble the scale of the R 2 . Good predictions will have low PRESS or high Q 2 (Szymańska et al. 2012).

Modelling habitat suitability
The MaxEnt algorithm was used to construct the potential distribution model of A. spinosa to North Africa and the Iberian Peninsula (Phillips et al. 2006;Elith et al. 2011). This algorithm represents the distribution of the species by the probability function P over a set X of locations in the studied area. P assigns a positive value to every place x so that the Wind direction Restrictive Rph pH soil information at 0-5 cm depth  sum of P(x) is the unit. In addition to the occurrence points used in the statistical analysis, we have added records of the presence of specimens introduced in Spain and Algeria, using the seven bioclimatic variables selected in the previous section and the soil pH layer. The output obtained predicts the suitability of the argan habitat satisfying all the associated bioclimatic limits and maintaining the distribution to its maximum entropy. The lowest suitability areas are symbolized by 0, while 1 symbolizes areas with maximum suitability (Navarro-Cerrillo et al. 2011;Sharma et al. 2018;Pratumchart et al. 2019). The ROC curve (receiver of characteristics) graphically represents the discriminative capacity of the model from all possible cut-off points and how adjusted this distribution is expressed by the value of the area under the curve (Poor 0.75 AUC > 0.95 Excellent) (Kogo et al. 2019). The function ENMevaluate from ENMevalin R package (version 1.3.1093) helps in selecting model setting that balance model goodness-of-fit and complexity (Muscarella et al. 2014).
The analyses carried out to evaluate the effect of the bioclimatic variables extracted from the values associated with the points with the presence of the species in their spatial distribution are summarized in Fig. 2.

Potential limitations
The scope of such studies depends on the quality of the input data. In our case, free access data (both occurrence points and explanatory variables) have been used, for this reason a one-by-one coordinates check and filter was made and a biasfile has been calculated ( Barbet-Massin et al. 2012). This type of data has resolution and precision limitations that could create biases in the results. For example, it does not identify topographic or climatic features occurring in specific areas smaller than 1 km 2 .
On the other hand, the wind station data belong to the global airport network and the ones closest to the study area were selected with data available in the period 1973-2000.

Results
The ecological amplitude of the species (Fig. 3, Online Resource 2,3) showed clear evidence of the argan adaptation to xeric climates, by observing the low annual precipitation (200-400 mm) (Fig. 3i), the high land surface temperature (up to 40 ºC) (Fig. 3r), and the aridity index ranging between 0.01 and 0.24 (Fig. 3s), typical from arid and semi-arid ecosystems in the sub-Saharan regions of Morocco and Algeria.
The thermometric variables show a clear trend of having a binomial distribution (Fig. 3). Figure 3 a, c, e, and g revealed very centred distributions with just 4 ºC between the minimum and maximum values. At the same time, the driest quarter of the year (Fig. 3e) is strictly centred between 22 and 24 ºC while the coldest quarter (Fig. 3h) has a wider range between 8 and 16 ºC. Argan tolerates minimum temperatures in a wide range (Fig. 3d) from 1 to 10 ºC and a narrow diurnal thermal range (Fig. 3b) between 10 and 12 ºC.
A Gaussian distribution with 40-60 mm (precipitation of wettest month) is showed in Fig. 3 j, while the precipitation of driest month (Fig. 3k) concentrates between 0 and 2 mm. The precipitation of wettest quarter has a maximum between 120 and 160 mm, which seems to indicate that the rain is concentrated in that period. In summary, 50% of the rain may occur during the wettest and coldest quarter in the argan habitats, but it is close to zero in the warmest and driest quarter. These findings further support the idea of the argan capacity to acquire humidity from wet wind condensation and fog-harvesting.
Evapotranspiration shows a gaussian distribution with values between 1500 and 2200 mm/year (Fig. 3t). The  (Fig. 3n) shows a clear multimodal distribution, recognizing the existence of three populations: a coastal is distributed between 0 and 400 m, an inland one located in the mountains of the Saharan and Grand Atlas, with altitudes concentrated to the 900 m, and another one with higher altitudes between the 1 200 and 1300 m, being able to rise in the Saharan Atlas above the 1500 m. The hillshade graph indicates values from 180 to 190º, and Fig. 3 p shows a clear preference for the N-NW orientations. The logistic distribution of the slopes (Fig. 3q) shows the preference of argan for flat areas (slope 15%), littoral plains, open plateaus, and oueds.
At the Sidi Ifni station, prevailing winds with a WNW-NNW direction are observed in the months of April, May and June and speeds between 1.5 and 3 m/s (Fig. 4a-c). However, in the northern stations of Safi (Fig. 4d-f) and Essaouira (Fig. 4g-i), the dominant humid winds occur during the summer months (July, August, September) in this case with NNW-NE direction and speeds between 2 and 6 m/s. The non-hierarchical clustering analysis shows four main groups applying all the 29 bioclimatic variables (see Online Resource 7) and the PCA analysis shows the distribution of the variables (Fig. 5).
The first three dimensions explain 74.59% of the total variability. Principal component 1 (PC1) explains 37.93% of the variability and is related to precipitation, temperature, and elevation variables (Temp, Tncm, Tmar, Prec, Pmdm, Pmdq, Elev). On the other hand, PC2 describes 27.46% with the more extreme values such like high temperatures, aridity index and evapotranspiration (Tmdr, Tmwm, Pmwm, Dai, Devt). The weighting loadings indicate the importance of each variable in the construction of the principal components (see Online Resource 8). The variable distribution reveals a significant difference in the dimension 1 between a group with higher elevation and annual precipitation during the driest month and quarter, besides a wide temperature range in contrast with another group that is affected by higher annual mean temperatures and lower temperatures during the coldest month. The dimension 2 separates the groups considering stress variables that get down localities with high aridity index, precipitation during the wettest month and quarter and get up localities with higher diurnal range, maximum temperature of the warmest month and evapotranspiration (Fig. 5).

Desert group (1)
In this group, the annual and mean diurnal range of the temperature varies between 12 to 14 ºC and 28 to 34 ºC, respectively. The localities of Tindouf (01) and Assa Zag (02) have the higher temperature during the warmest month which causes more evapotranspiration. Along with the Laayoune (03), Tiznit (51) and Guelmim (49,50) localities, all together comprise a geographical area influenced by the Sahara, with a hyperarid climate and annual rainfall above 70 mm and less precipitation during the wettest month and the driest quarter. It is important to mention that these three last samples have less evapotranspiration and a lower temperature range, that is the reason why they are classified in the sublittoral group.

Littoral group (2)
Due to the Atlantic influence, the Littoral group of 18 localities has the least range temperature, this means, the lowest temperature during the warmest month and the highest temperature during the coldest month. This area that also include the Mediterranean coast localities, has less evapotranspiration, and is considered as Arid climate. In this class, the elevation ranges between 26 and 300 m except for two more inland localities (14, 15) that reach 500 m. The recorded annual rainfall is between 200 and 328 mm and becomes lower during the wettest month in a range of 47 to 62 mm. Besides, the precipitation does not exceed 8 mm during the driest quarter.

Mountain group (3)
This group is made up of 18 localities characterized by having lower temperature values both in the annual mean (14-18 ºC) and in the minimum of the coldest month (1.30-5.70 ºC). Likewise, annual precipitation varies in a range of 200 to 427 mm with lower precipitation during the wettest month (35-69 mm) and the driest quarter (5-16 mm). The altitudinal range of this group is the highest with elevations from 750 to 1500 m. This group also includes the two isolated localities from Oued Grou (39) and Chlef (35) that have lower elevation and higher temperatures than the others.

Sublittoral group (4)
The sublittoral group consisting of 13 localities corresponds to a more inland and meridional area located between 88 and 488 m, not reach higher altitude zones. This group is characterized by having intermediate climatic conditions (Tmdr = 11 ºC, Tmar between 21.2 and 24.8 ºC) higher than in the Mountain group, but lower than in the Littoral group. The precipitation of the driest quarter varies in a range of 2-5 mm. It is the group with the lowest aridity index after the desert group.
The most important variables (VIP > 1) identified by the PLS-DA method were Tmdr, Pmdm, Tncm, Elev, Tmar, Dai, and Pmdq and the boxplots show the distribution of the VIP variables in each group (see Online Resource 10). The nonuniform concentration of the desert group samples coincides with the similarities between the southernmost Fig. 6 Clustering results using k-means in the native area. The desert group is represented by orange rhombuses, littoral by magenta circles, mountain by blue triangles, and sublittoral by green squares. Darker areas show mountain ranges above 1 000 m sublittoral localities (Laayoune, Tiznit, and Guelmim), even so, for most variables the remaining three groups are well differentiated.
In total, 66 occurrence records of A. spinosa and the 7 most important environmental layers were used in this suitability analysis, including the Soil pH in H 2 O (Dph) as a restricting variable (see Online Resource 11).
Ten MaxEnt runs were executed to determine the distribution model resulting in a robust and stable model (see Online Resource 12). The mean AUC value of 0.98 indicates that the distribution explains very well from the environmental variables.
The MaxEnt model's jackknife test of variable importance indicated that the aridity index (Dai) is the most effective single variable for predicting the distribution of argan in the considered area. However, the elevation (Elev) is not (by itself) useful for estimating the distribution of the species. Omitting the minimum temperature of the coldest month (Tncm), temperature annual range (Tmar) or Dai lightly decreases the predictive performance, which therefore appears to have the most information that was not present in the other variables (see Online Resource 13).
The response curves (see Online Resource 14) plot the dependence of the predicted suitability model both on the selected variable and on dependencies induced by correlations between the selected variable and other variables. Taking as reference the values for P = 0 (absence) and P ≥ 0.50 (presence), these graphs substantiate the analysis of ecological profiles.
Areas with a high probability of occurrence are shown in red, gradually decreasing to low probability in blue (Fig. 7). Figure 7 (a, c, e) shows the coherent zones that the model considers ideal for the cultivation of argan in all the studied territories. In any case, the introduction of the pH of the soil considerably reduces the potential area of argan cultivation, giving the model even more coherence Fig. 7 (b, d, f).
In Morocco, the suitability model is clearly more restrictive when including the pH layers especially as regards the level of probability of success, but there is little change in the area given the extensive presence of non-acid soils with pH ≥ 7 in most of the Moroccan territory and the extensive area of distribution of argan in the SW of the country.
The model in Algeria shows the limited potential area for argan cultivation, reduced to the western Mediterranean coast. Further south, some level of suitability reappears, coinciding with the Hoggar and Tassili massifs.
In the Peninsular Spain, there is a cultivable area in the SW in the province of Huelva and the ideal zone is mainly spread over the Iberian S, SE, and E always looking for the most thermal areas at the same time as close to the sea as possible and without a special preference for mountainous areas, with greater probability in the province of Almeria. This is also observed in the Balearic Islands.
In the Canary Islands, it is not surprising the high probability of success for the cultivation of argan, considering the presence of Sideroxylon canariense, a Macaronesian endemic species close to A. spinosa, which lives in phytosociological communities, equivalent to argan in Morocco, in shrubs of succulent plants dominated by species of Euphorbia (tabaibos and cardones), Kleinia, Ceropegia, Dracaena and other species of Macaronesian distribution.
Finally, regarding Tunisia, the model recognizes the suitability of the flat coastal areas that appear on the eastern slopes of the small mountains and Tunisian reliefs. Towards the interior, despite the presence of relief of the Keirouan mountains, the model is not wrong in recognizing their suitability as a hyperarid orographic system, where only Phoenix dactylifera palms developed at the bottom of the drainage systems. The suitability model that includes pH reduces the area of success in the northernmost part of the country by removing acidic soils. For more detailed description, see Online Resources 15, 16, 17.

Discussion
The multimodal distributions of the ecological profiles demonstrate the different ecological optimum of argan in its distribution area and confirm the existence of diverse habitats, something already mentioned in studies on argan plant communities (Msanda et al. 2021). Thus, the elevation and mean precipitation of the driest month are trimodal and the minimum temperature of the coldest month is quadrimodal, suggesting that the argan tree lives in frost-free areas and adopts strategies to compensate for extreme climatic conditions, for example, being located at the bottom of oueds and taking advantage of the moisture supply from the fog and dew regimes in the more desert-like areas.
It is well-known that the Atlantic influence in the argan area located between 27º and 32ºN partly compensates its xericity at the littoral habitat (López-Sáez and Alba-Sánchez 2009; Faouzi et al. 2015;Mouhaddab et al. 2015). According to Hammou et al. (2018) and Nerd et al. (1998), the warm winds of the Sahara (named Chergui) and the amount of water available affect the flowering of argan trees in southwestern Morocco. This area benefits from the oceanic influence whose atmospheric circulation is dominated by the Azores anticyclone with N-NE winds that drag summer fog inland.
The argan drought-avoiding strategy has been studied in the areas of its ecological optimum, according to Zunzunegui et al. (2018), deep soil layers are its main source of water. However, 130 km to the southwest of there (sublittoral zone) where more extreme conditions occur, fog and dew events become more relevant. This dependence on fog as water supply has been studied with other xerophytic species in hyperarid coastal zones like Namib (Welwitschia mirabilis) and Atacama (Tillandsia spp.) deserts (Wolf et al. 2016;Mitchell et al. 2020).
The advection fog appears in southwestern Morocco when moist air from the Atlantic Ocean flows inland, where it meets the relief of the region. The wind speed (1.2-2.4 m/s) and orientation (WNW-NNW) form the ideal conditions for generating dense fog layers (Marzol and Sánchez-Megía 2008;Lekouch et al. 2012). These conditions coincide with those observed in the months of April, May, June in Sidi Ifni (SW), and July, August, and September in both more to the north Safi and Essaouira.
According to Marzol and Sánchez-Megía (2008), although there are more fog days on the coast than inland, the collected amount of water could be higher in the mountains (perhaps due to adiabatic cooling, a phenomenon typical of mountain fogs). A correlation between the fog events, an efficient condensation, and the presence of phenological phases sensitive to water deficit could be suspected of this outcome. In fact, the months with most fog events inland coincide with the flowering period of some argan populations (April, May, and June) (Nerd et al. 1998).
These argan abilities to adapt to the different habitats are achieved thanks to its polymorphism, a character that has Fig. 7 Potential distribution of A. spinosa generated by MaxEnt ecological niche modeling using 7 (excluding pH) and 8 (including pH) environmental variables: Mean diurnal range (Tmdr), minimum temperature of the coldest month (Tncm), temperature annual range (Tmar), precipitation of the driest month and quarter (Pmdm, Pmdq), elevation (Elev), aridity index (Dai) and Soil pH in H 2 O (Dph). Morocco and Algeria ((a) excluding pH, (b) including pH), Tunisia ((c) excluding pH, (d) including pH), and Spain ((e) excluding pH, (f) including pH). Bluer shades represent less suitability and redder ones represent more suitability been demonstrated by numerous authors at molecular level (El Mousadik and Petit 1996;Majourhat et al. 2008;Yatrib et al. 2015Yatrib et al. , 2017Mouhaddab et al. 2017;Ait Aabd et al. 2019). At the same time, the differences observed in their phenology (Zahidi et al. 2013;Hammou et al. 2018) are also correlated with this polymorphism. Surprisingly, the model does not recognize suitability in Tindouf (west Algeria) where the presence of argan is naturalized or may be even spontaneous, and in any case a cultivated species (Kechebar et al. 2013;Kechairi et al. 2018). The unusualness of the environmental values of this area of southwestern Algeria is not very representative with respect to other larger areas, making these characteristics outliers and do not effectively train the suitability model. As in many other species of economic interest, the presence and human footprint are so ancient that often confounds interpretation whether they are palaeophytes or relict species in certain locations.
However, although samples from Tunisia were not used, the model recognizes suitability in the provinces of Tunis, Nabeul, Sfax, Sousse and Medenine named in recent studies (Hanana et al. 2018;Louati et al. 2019b). The same occurs in Canary Islands where there is evidence of argan introduction in some locations (https:// www. atlas desem illas decan arias. org).
Finally, the identification of the ecological euroicity of argan justifies the comprehensiveness of the suitability models obtained in this study and makes it possible to find solutions to mitigate the effects of climate change on natural argan populations.

Conclusions
The simultaneous application of three complementary quantitative analysis methods (ecological profiles, multivariate analysis, and suitability models) has allowed identifying the ecological requirements of A. spinosa, analysing the heterogeneity of its habitat and establishing the successful potential areas for cultivation of argan in the western Mediterranean Basin.
Outcomes from this ecogeographic study further support the polymorphic character of argan and provide insights into the specific adaptation strategies that it has allowed growing in littoral, sublittoral, mountainous, and desert habitats.
In addition to its natural range in Morocco and Algeria, potential argan cultivation areas have been characterized as more suitable in Spain than in Tunisia. The Mediterranean coast of south-eastern Spain represents the area that best will resemble the ecological requirements of the argan tree and can be considered an ideal choice as a cultivation area and climatic refuges to improve the conservation of this multipurpose species.
Funding Open Access funding provided thanks to the CRUE-CSIC (Universidad de Córdoba/CBUA) agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.