Modeling Cedrus atlantica potential distribution in North Africa across time: new putative glacial refugia and future range shifts under climate change

In recent years, species distribution models have been used to gain a better understanding of past and future range dynamics of species. Here, we focus on a keystone species of the North African forest ecosystem (Cedrus atlantica) by calculating a consensus model of the species current geographic potential distribution in North Africa, based on a weighted average method aiming to decrease uncertainty. The consensus model is obtained using seven species distribution model algorithms taking into account 24 environmental variables. The model is then applied to several past and future time slices. Past projections refer to the Middle-Holocene and the Last Glacial Maximum, whereas those of future are related to expect conditions around 2050 and 2070. We found that the current potential distribution of Cedrus atlantica is larger than its actual geographical distribution. For some explanatory variables, the analysis revealed their importance for the species current distribution. Among all obtained models, that for the Middle-Holocene showed the maximum expansion of the species potential distribution. The Last Glacial Maximum model provided new putative glacial refugia of Cedrus atlantica, not shown by other mechanistic models and palaeorecord localities. Future projections revealed a significant and fast contraction with shifting in altitude of the species range, showing more fragmented areas and even species disappearance in many North African localities. These findings can help to restore cedar forests and conserve them by ex situ strategies according to the future defined refugia in North Africa. Attention should be paid to the resolution of related output maps, the current biotic interactions, and those that may arise under climate change.


Introduction
Cedrus atlantica Manetti is endemic to North African mountains and constitutes the south-westernmost species of the genus Cedrus. It is considered the noble tree of Algerian and Moroccan forests due to its remarkable ecological, ornamental, and forest qualities (Benabid 1994;Derridj 1990;Quézel 1998).

Editor: Wolfgang Cramer
Electronic supplementary material The online version of this article (https://doi.org/10.1007/s10113-019-01503-w) contains supplementary material, which is available to authorized users.
In Algeria, C. atlantica occurs in the Tellian and Saharan Atlas Mountains, covering an area of 16,000 ha (DGF 2007). In Morocco, it occurs in the Rif Mountains, the peaks of Tazzeka, the Eastern Middle Atlas, the Central Middle Atlas, and the Eastern High Atlas Mountains, with a total area of 145,000 ha (M'hirit 1999). It occurs where a Mediterranean climate dominates and occupies semi-arid, sub-humid, humid, and per-humid bioclimatic floors with fresh to extremely cold thermal variants. It occupies regions with a mean annual rainfall of 400-2100 mm/year and an average of minimum temperature of the coldest month ranging between − 8.5 and 2.7°C (M'hirit 1999;Yahi et al. 2008). C. atlantica is indifferent to the substratum (Faurel and Laffite 1949;Maire 1924), and its fragmented distribution is primarily affected by the climate and orography of the Maghreb; it is restricted to mountain areas between 1300 and 2600 m (Benabid 1994;Derridj 1990;Emberger 1938a;Quézel 1998). It forms pure or mixed stands with other evergreen and/or deciduous tree species creating forests hosting remarkable biodiversity. Syntaxonomically, C. atlantica forests belong to the order of Querco-Cedretalia atlanticae Barbero, Loisel et Quézel, 1975(Benabid 1985Benabid 1994;Yahi and Djellouli 2010).
Given its ability to tolerate drought and cold, as well as the high quality of its wood, Atlas cedar has been introduced into several countries, especially in the Mediterranean region, where it is naturalized over large areas (Campredon 1934;Courbet 1991;El Azzouzi and Keller 1998;Pavari 1927;Slimani et al. 2014a;Toth 1973;Toth 2005;Touchan et al. 2010). However, it is infrequently considered in Algerian and Moroccan reforestation programs (Bensaid et al. 1998;Messaoudène et al. 2013;Quézel et al. 1990).
Cedar forests are declining dramatically in both Algeria and Morocco, and many stands are facing constraints on regeneration (Ezzahiri et al. 1994;Lepoutre and Pujos 1964) and destined to disappear in the coming decades (Abdessemed 1984;Quézel 1998). This is clearly seen in the Hodna forests in Algeria (Madoui and Gehu 1999;Slimani et al. 2014a). The lack of natural regeneration in some localities is coupled with other threats that strongly affect their existence, such as dieback, which is likely associated with climate change (Bentouati and Bariteau 2006;Ghailoule et al. 2012;Kherchouche et al. 2013), fires, and long-term browsing by domestic animals that all widespread in North Africa (Bensaid et al. 2006;Slimani et al. 2014b;Touchan et al. 2010). For example, the disappearance of Atlas cedar from Khroumire (Tunisia) is probably related to the Phoenician settlement (Ben Tiba and Reille 1982) and only place names attest to its past occurrence (Derridj 1990).
These factors have led to a remarkable contraction of Atlas cedar's natural range, with Algerian cedar forests decreasing from 37,910 ha (Combe 1889) to 16,000 ha in the last 130 years, a loss of almost 58%. Hence, it is not surprising that C. atlantica is in the red list of threatened species of the IUCN (Thomas 2013). The current state of North African cedar forests does not allow these forests to fulfill their ecological and socio-economical functions (Benabid 1994;Yahi et al. 2008).
Since the current distribution of Atlas cedar does not reflect its ecological range (Aussenac 1984;Derridj 1990), producing maps of its potential distribution and possible future refugia is necessary for effective conservation. The North African region is expected to experience drastic climate change (IPCC 2007;Klausmeyer and Shaw 2009), with an expected increase in temperature of 4-5°C and a reduction in precipitation of > 30% in the next few decades (Giorgi and Lionello 2008).
In view of this, we applied a correlative approach heavily based on bioclimatic predictors with the aims of obtaining modeling outputs of past, present, and future potential distributions and determining future putative refugia of Atlas cedar through the characterization of its ecological niche sensu Hutchinson (1957). We tested the assumptions that in the past the area of Atlas cedar had to be less fragmented and that the aridity of north-western Algeria and eastern Morocco was not a barrier to species colonization of the North African region (Emberger 1938a).
Ecological modeling techniques, especially species distribution models (SDMs), have undergone substantial development in recent decades (Guisan and Thuiller 2005;Peterson and Soberón 2012). Despite some limitations, SDMs are considered promising modeling tools for addressing conservation questions regarding species (Ferraz et al. 2012;Franklin 2013;Pearson 2010;Zhang et al. 2012) and terrestrial Johnson and Gillingham 2005;Thorn et al. 2009;Urbina-Cardona and Flores-Villela 2010) or marine ecosystems (Bentlage et al. 2009;Coro et al. 2013;Pike 2013). The development of SDMs relies on the geographic distribution range concept at both biogeographical and ecological levels (Peterson and Soberón 2012), as well as to new powerful statistical techniques and GIS tools (Austin 2002;Guisan et al. 2002;Guisan and Zimmermann 2000). Depending on the type of data required, SDMs use statistical methods (e.g., general linear models (GLMs), general additive models (GAMs)) or machine-learning methods (e.g., support vector machines (SVMs), artificial neural networks (ANNs), maximum entropy methods, and classification and regression trees (CARTs)) (Franklin 2010;Guisan et al. 2002;Guisan and Thuiller 2005;Peterson et al. 2011;Phillips et al. 2006;Stockwell 2006;Thuiller et al. 2003). They are intended to summarize relationships between the occurrence data of a given species that can be Bpresence-only,^Bpresenceabsence, or Babundance observations,^and environmental predictor variables with direct or indirect effects on the species' distribution (Austin 2002;Elith and Leathwick 2009;Franklin 1995;Guisan and Thuiller 2005). Data preparation before modeling is performed in a GIS environment, and the final result is a georeferenced map of occurrence probabilities (Franklin 2010;Phillips et al. 2006).
To our knowledge, such techniques are less well used in the study of tree species (examples are found in Attorre et al. 2008;Hidalgo et al. 2008;Vessella and Schirone 2013). In North African region, the present study constitutes one of the first attempts using the correlative approach (Arar et al. 2019;Tabet et al. 2018) while the mechanistic approach was already applied (Cheddadi et al. 2009;Demarteau et al. 2007).
Here, we used presence-only data and 24 environmental variables, including 19 bioclimatic, three topographic and two calculated variables. MaxEnt (Phillips et al. 2006) and six other algorithms run on openModeller Desktop (De Souza Muñoz et al. 2009) were used to produce a consensus model (CM) which could be projected in different time slices. In this study our objectives were: (i) to predict the current potential geographic distribution of Atlas cedar as well as future refugial areas and to identify the environmental variables mostly influencing its geographical distribution, (ii) to provide a better understanding of the potential response of Atlas cedar to past and future climate changes; and (iii) to help forest managers select appropriate areas to be restored and preserved given the urgent need to conserve such an iconic species.

Preparation of occurrence data
The current distribution of the North-African cedar forests was mapped, producing the most accurate picture to our knowledge, and which can be updated with future information (Fig. 1). To obtain this, GPS data for certain Algerian cedar forests were collected. Other data for C. atlantica forests in Algeria and Morocco were obtained from the scientific literature (Abdessemed 1981;Bouahmed 2012;Boudy 1955;Chbouki 1994;Cheddadi et al. 2009;Damnati et al. 2014;De Smet and Bouaza 1984;Derridj 1990;Faurel and Laffite 1949;Kaabeche 1996;M'hirit 1999;M'hirit and Benzyane 2006;Mouna 2009;Negre 1952;Slimani et al. 2014a;Slimani et al. 2014b;Terrab et al. 2006;Till and Guiot 1990;Touchan et al. 2010). We used the open source solution BSAS.Planet^(www. sasgis.org/) connected to several free online mapping services; 1423 polygons saved in BKML^format were defined using photointerpretation techniques used in forestry (Gougeon and Leckie 2003;Paine and Kiser 2012;Provencher and Dubois 2007). The BKML^file was converted to shapefile, then to raster, and finally to occurrence points (4585 entries) by using the Bconversion tools^of ArcGIS Desktop 10.2 (www.esri. com) with a resolution of 0.0083 decimal degrees (30 arcsec) and expressed in the same geographic coordinate system as that used for environmental variables (GCS_WGS_1984).

Current period
Twenty-four environmental layers with a resolution of 30 arcsec (~1km) were used (Table 1): (a) 19 bioclimatic layers for the period of 1960-1990 and elevation, which were downloaded from WorldClim database (Hijmans et al. 2005;www.worldclim.org/) (with elevation originating from SRTM (www2.jpl.nasa.gov/srtm/)), (b) two other topographic variables (slope and aspect) derived from elevation by using the 'Surface' extension of Bspatial analyst tools^in ArcGIS 10.2, and (c) two calculated variables referred to as the Emberger index BQ^ (Emberger 1930)  A priori, we considered the bioclimatic and topographic variables as most important in the distribution of C. atlantica and did not consider edaphic and geological variables, in accordance with many other authors (Benabid 1994;Boudy 1955;Derridj 1990;Quézel 1998;Yahi et al. 2008). Moreover, we aimed to find the most explanatory variables toward C. atlantica potential distribution. We explored first the prospective linear links that may exist among the quantitative environmental variables used to perform the model of the current period (present) clipped according to the presence points of the species. This was achieved by calculating a standard Pearson correlation coefficient Br^using the corrplot package (Wei and Simko 2016) of R v3.0.0 (R Development Core Team 2013) after extracting and exporting the values of the clipped environmental rasters.

Past periods
For the projections in past periods, we used 19 bioclimatic variables of WorldClim and the two calculated indices (Q and CI). Two periods were investigated: (i) the Middle-Holocene (hereafter, Mid-Holocene; ca. 6.5 kya BP), for which the variables have a resolution of 30 arcsec (for technical considerations, we accessed only to the data of the CCSM4 model); and (ii) the Last Glacial Maximum (LGM, ca. 21 kya BP), for which the variables, originally available at 2.5-min resolution, were downscaled to 30 arcsec using the nearestneighbor method. For LGM, the global climate models (GCMs) used were CCSM4, Miroc-ESM, and ESM-MPI-P (see Table S1). An arithmetic mean was then applied for each variable taking into account these three GCMs.
The topographic variables (elevation, slope, and aspect) were not considered when projecting in past periods, but they were included in statistical analysis (to plot the graphs of probability density). For the Mid-Holocene, we retained the same layers as the current period, whereas for LGM, the elevation was retrieved from ETOPO1 Global Relief Model (www.ngdc.noaa.gov/).

Future periods
For projections into the future, the same bioclimatic variables from WorldClim (in 30 arcsec) were used as for the past periods. We considered two scenarios of representative concentration pathways (RCPs 4.5 and 8.5) and two time slices (2040-2060 centred on 2050 and 2060-2080 centered on 2070). All GCMs available from Bwww.worldclim.org/ cmip5_30s^were used, i.e., 19 for RCP 4.5 and 17 for RCP 8.5 (see Table S2). For each of the four combinations of Bscenario-period,^an arithmetic mean was applied for each bioclimatic variable, taking into account all the GCMs. Q and CI were calculated from the averaged variables, while for topographic variables, we used those of the current period.
Potential distribution of C. atlantica First, we constructed C. atlantica niche models for a single time period (the present) using seven different algorithms with presence-only data to obtain a plausible CM (Marmion et al. 2009;Thuiller et al. 2009). We, then transferred the resulting CM to past and future periods. Two software applications were used to build current single models (SMs): openModeller Desktop 1.1.0 (De Souza Muñoz et al. 2009) and MaxEnt v3.3.3k (Phillips et al. 2006). The main settings used are summarized in Table 2.
Calculating a consensus model in SDMs involved combining probability values of occurrence of different models in order to decrease the predictive uncertainty (Araújo et al. 2005;Thuiller 2004). All SDMs had a continuous or semicontinuous scale of occurrence probabilities (De Souza Muñoz et al. 2009;Phillips et al. 2006;Phillips et al. 2004;Sutton et al. 2007). For Bioclim, 10 standard deviation cutoff were applied (see Table 2), with each model producing three ). An arithmetic mean was applied to obtain a single model using this algorithm, with a semi-continuous scale of occurrence probabilities. For the current period (the present), the CM was calculated using a weighted average (WA) method. This is considered more robust than other methods and requires a pre-evaluation of the predictive performance of the SMs (Marmion et al. 2009). Following these authors, we multiplied each model by its corresponding area under the curve (AUC) of the Receiver Operating Characteristic (ROC) curve. (see Evaluation of models).
The transferability or generality of a calibrated model regarding SDMs is defined by its spatial or temporal crossapplicability (Phillips 2008;Randin et al. 2006). In our study, the modeled niche conditions were transferred into past and future periods, assuming niche conservatism (Martínez-Meyer and Peterson 2006;Martínez-Meyer et al. 2004;Peterson et al. 1999), by using the Braster calculator^of ArcGIS. This was done using a script based on the minimum and maximum values of each bioclimatic variables (including Q and CI) extracted depending on the CM. Since the past periods included the LGM and the Mid-Holocene and the future projections comprised two periods (2050 and 2070), each represented by two scenarios (RCPs 4.5 and 8.5), this gave two maps for the past periods and four maps for future projections.
The final result was 7 georeferenced maps (including one of the present time) with the same resolution as the input data (30 arcsec) and showing the potential geographic distributions of C. atlantica across time based on its probability of occurrence (%).

Future refugial areas
The continuous probability maps of likely future potential distributions generated above were converted to binary maps (presence-absence) using a Blowest predicted value^method sensu Pearson (2010). A probability threshold of 50% was applied, defined as the lowest predicted value that corresponded to the observed occurrence records. The superposition of the presence cells of the different binary maps allowed us to extract putative future refugial areas of C. atlantica by 2070, including suitable areas for species persistence from the LGM to 2070 (Vessella et al. 2015 and references therein). This option appeared more realistic than deducing, globally, the future refuges according to future potential distribution maps of the species.
To explore the species response to the quantitative environmental variables, we plotted the graphs of probability density for these variables using the 50% threshold above for all the investigated periods (LGM, Mid-Holocene, current period, 2050, and 2070) taking into account two future scenarios of RCPs 4.5 and 8.5. This was performed using the ggplot2 R package (Wickham 2009).

Evaluation of models
Several methods are available to test the predictive performance of distribution models (Bahn and McGill 2013). To test the accuracies of the generated SMs, we used the area under the curve (AUC) of the receiver operating characteristic (ROC) method (Fielding and Bell 1997). This is a twodimensional depiction of classifier performance (Hand 2009) and Bequivalent to the probability that the classifier will rank a randomly chosen positive instance higher than a randomly chosen negative instance^ (Fawcett 2006). The AUC is widely used to assess SDMs, with values above 0.5 indicating a performance better than random (Franklin 2010;. We considered AUC values of 0.5-0.7 to indicate low accuracy, 0.7-0.9 indicative of moderate (good) accuracy, and > 0.9 indicative of high accuracy (Manel et al. 2001;Swets 1988).

Current potential distribution of C. atlantica (CM)
The obtained CM clearly shows the potential geographic distribution as much larger than the current distribution of C. atlantica, even if we focus on high probabilities of occurrence (Fig. 2a).
The predicted potential distribution shows that the southern end of the current distribution of the species is particularly suitable. Our findings also indicate the crucial importance of North African orography for the species distribution. Figure 2a shows  Table 3 shows the potential areas of C. atlantica based on threshold values applied on the CM map. The extent of these areas is perhaps optimistic, especially for thresholds of 0.5 and 0.6, as they do not take into account other limiting biotic (e.g., competition) and abiotic (e.g., soil, land cover) factors. Figure 3 shows the linear correlations that may exist between quantitative environmental variables used in modeling the current period. Corresponding values were extracted according to the presence of the species. This figure shows the existence of nuanced linear relationships between almost all studied variables at α = 0.05, with Pearson correlation coefficients ranging from 0.03 to 0.99. However, slope, annual mean temperature, isothermality, min temperature of coldest month, and precipitation of wettest month did not show correlations with the other bioclimatic variables.

Evaluation of the SMs and importance of environmental predictors
Evaluation of the current SMs (corresponding maps in Fig. S1) revealed AUC values exceeding 0.90, indicating high accuracy in all models except MaxEnt (0.87) (Fig. 4). The highest AUC was obtained for the environmental distance algorithm (AUC = 1.0). These high values indicate the quality of the resulting CM since we used the AUC values as a weighting factor.
To evaluate the importance of the studied variables in influencing the potential distribution of the species, we focused on the MaxEnt results. A jackknife test (Table 4) revealed the importance of almost all the variables in influencing the distribution of C. atlantica as generated by MaxEnt. The environmental variable with the highest gain, when considered in isolation, was annual mean temperature (Bio1), which appeared to be the most informative by itself. Mean temperature of coldest quarter (Bio11), mean temperature of warmest quarter (Bio10), minimum temperature of coldest month (Bio6), elevation, mean temperature of wettest quarter (Bio8), and precipitation of driest quarter (Bio17) are considered to be likely explanatory predictors as well, showing almost as high gain values to that of annual mean temperature. Some precipitation variables, mainly annual precipitation (Bio12) and precipitation of driest month (Bio14), also appeared to have some influence the species distribution, while aspect had a marginal contribution. Figure 2b, c shows that during past periods (LGM and Mid-Holocene), C. atlantica was potentially able to colonize more areas than it does today. Potential distribution reached a maximum during the Mid-Holocene. Potential regions of occupation outwards from the current geographic zones are notable, particularly east of Tell, Ouled Nail Mountains, Tlemcen-Debdou Mountains block, the Western High Atlas, and a restricted part of the Anti-Atlas Mountains.

Past projections
The LGM model (Fig. 2b) suggested that Zeghouane, Djebel Serj, the easternmost zone of Algeria, Ouled Nail Mountains, the mountains of Tiaret, Tlemcen, Beni-Snassen, Debdou, and the Western High Atlas might have constituted putative glacial refugia. During the LGM, several regions that correspond to current Atlas cedar forest were only weakly favorable or even unfavorable for the specie. These regions encompassed Theniet Elhad, Chrea, Djurdjura, Hodna, and Belezma in Algeria; High Atlas Mountains, Occidental Rif cedar forest, and the majority of Central Middle Atlas forests in Morocco. The model instead suggested the presence of C. atlantica in lower-altitude refugia, close to localities currently occupied by Atlas cedar forest.
During the Mid-Holocene, potential areas of Tlemcen and Debdou mountains showed their maximum expansion (Fig. 2c), attesting to the highly suitable climatic conditions for C. atlantica, whereas the potential areas in the regions of Ouled Nail and Western High Atlas were much smaller. Contraction over the entire potential area was evident, leading to the current potential distribution pattern. Low-lying refugia (of the Mid-Holocene) exhibited marked shifts up an altitudinal gradient, developing gradually into the current distribution. This is the case of Djurdjura and Babors in Algeria, much of the Western Middle Atlas Mountains and almost all of the Eastern High Atlas Mountains in Morocco.
In Tunisia, the models had low probabilities of presence of Atlas cedar in both the LGM and the Mid-Holocene, except in the Zeghouene, Djebel Serj, and Djebel Mghila regions, where the species was still very restricted. Indeed, its presence at that time still remains to be verified. It is important to note, however, the very low probabilities attesting the eventual

Future projections
The natural area occupied by Atlas cedar showed a significant contraction in the last century, mainly due to intense anthropogenic actions and dieback. However, the future projections of the species' ecological niche ( Fig. 2d-g)

Putative future refugia
Putative future refugia of C. atlantica, defined as 50% of the modeled potential areas throughout all the investigated periods from the LGM to 2070 (see BFuture refugial areas^), are shown in Fig. 5a, b based on two future scenarios of RCPs 4.5 and 8.5. The RCP 4.5 (Fig. 5a) shows potential persistence of the species in the Tlemcen Mountains in Tellian Atlas and Belezma and Aures in the Saharan Atlas (Algeria), the Rif, Debdou Mountains, Middle Atlas and the Western High Atlas Mountains (Morocco). However, the RCP 8.5 (Fig. 5b) shows that in Algeria, Atlas cedar would persist only at high altitudes in the Aures region, whereas in Morocco, future refugia include just three localities, the Rif, a restricted area of the Eastern Middle Atlas, and the Western High Atlas Mountains. Density probability graphs plotted using 50% of the modeled potential areas throughout all the investigated periods were constructed to summarize the potential response of the species to the studied quantitative environmental variables across time (Figs. S4 and S5). By ignoring low probability classes, our findings clearly showed that the strongest explanatory variables (see BEvaluation of the SMs and importance of environmental predictors^) changed slightly over time, whereas the less-explanatory variables showed an expected morepronounced shift in optima. However, the minimum and maximum values of species tolerance for the environmental variables remained more or less unchanged across time.

Discussion and conclusion
Focusing on past, current, and future ecological niche modeling of C. atlantica, we confirmed the usefulness of a consensus model based on the WA method after implementing many SMs (Marmion et al. 2009). Our findings provide new insight into the spatiotemporal dynamics of the species' potential range. The high AUC values from the SMs indicates high accuracy, comparative with other studies using the same algorithms (Benito et al. 2009;Giannini et al. 2010;Leidenberger et al. 2015;Pouteau et al. 2011;Tarkesh and Jetschke 2012).
The exploration of environmental variables extracted according to the presence points of C. atlantica revealed new information about the interactions between these variables within the current species distribution. It has become quite common to verify multicollinearity and spatial autocorrelation, particularly for statistical models, in selecting input predictor variables for SDMs (Franklin 2010). However, we did not do this due to the different kinds of algorithms we used to calculate a consensus model, as in similar cases (Marmion et al. 2009) and due to  our aim to explore, for the first time, the contribution of some variables to the potential distribution of C. atlantica. However, to aid future modelers applying the same geographical mask, we calculated the standard Pearson correlation coefficients for the correlations between all quantitative environmental predictors using the BCorrelations and Summary Stats^tool of the SDM Toolbox v1.1c (Brown 2014). Strong correlations have been recorded between some variables (see Table S3), if applying a threshold of |0.9| or |0.8| for Br^as adopted by many authors (Jaryan et al. 2013;Yang et al. 2013 Fig. 4 Area under the curve (AUC) values of the single models (SMs) used in this study

Current potential distribution of C. atlantica
To our knowledge, the implemented CM is the most geographically detailed model of the potential distribution of C. atlantica to date (Cheddadi et al. 2009;Demarteau et al. 2007). We consider the model to be plausible given that predicted favorable areas, particularly those with high values of probability of occurrence, were not found in unexpected locations, such as desert or very low altitudes. Such locations are unexpected given that the current geographic distribution of C. atlantica starts from 1300 m above sea level. It is noteworthy that at such a scale (~1 km), the model provides a valuable tool for reforestation efforts (Vessella and Schirone 2013) in North Africa aimed at species conservation (Pearson 2010) given the high degree of spatial fragmentation of its populations. Many authors agree that the current environmental potentialities for C. atlantica (Aussenac 1984;Boudy 1950;Derridj 1990) lie beyond its current geographical distribution pattern, which is consistent with our CM results. However, the potential areas we obtained based on probability thresholds of occurrence ≥ 50 and ≥ 60% appear to be somewhat overestimated (Table 3). This may be due to the spatial resolution adopted or to the over-predicted SMs used to calculate the CM. However, the potential areas are close to estimations of Boudy (1950) i.e., 456,000 and 128,000 ha for Morocco and Algeria, respectively. These areas correspond approximately to areas obtained with the CM at a probability threshold of 80% or possibly 70% after the application of other restriction parameters to the model. The potential distribution obtained in this study does not consider other variables influencing the dimensions of the species niche, such as biotic interactions (such as competition or parasites) and recent human activity (Guisan and Thuiller 2005;Guisan and Zimmermann 2000). This is the subject of ongoing work on reforestation in North Algerian populations; the loss of 58% of Algerian cedar forest over the last 130 years is striking.
C. atlantica will potentially be able to colonize some areas in the southern part of its range. However, we interpret this result with caution given the high rate of dieback reported in many locations of this region, most likely related to climate change (Bentouati and Bariteau 2006;Kherchouche et al. 2012). At the local scale, we suggest that the existence of suitable areas in such regions might be due to the output resolution of the CM or, probably, to the specific ecological conditions of those regions. These possibilities warrant further investigation. In addition, post-processing evaluations should be performed for some Tunisian localities that are predicted to be currently suitable (i.e., Zeghouene and Djebel Serj). For Tunisia, the CM confirms the inevitable extinction of C. atlantica (as happened in the past) due to climate change (Pons and Reille 1984;Stambouli-Essassi et al. 2007).

Variables explaining the distribution of C. atlantica
The MaxEnt results revealed that the temperature-related variables made a strong contribution to the current distribution of C. atlantica. In contrast, published data suggest that low winter temperature coupled with W/NW humid winds are of primary importance in affecting the species' distribution (Emberger 1938b;Faurel and Laffite 1949). The identified importance of elevation from our models is corroborated by the variation in altitudinal limits of Atlas cedar forests (Achhal et al. 1980;Quézel 1998). Precipitation variables had similar gain values to those of temperature; these predictors are affected by North African orography and are much more spatially variable than temperature (Meharzi 1994;Seltzer 1946). This does not limit the importance of precipitation for the distribution of Atlas cedar; rather, it is considered the most determining factor in the species' distribution (Faurel and Laffite 1949). The very low contribution of aspect might be explained by the resolution adopted; this variable plays an important role in some localities, particularly in the southern limits of the species' range (Abdessemed 1981). The present study revealed, for the first time, the contribution of some bioclimatic variables to the current distribution of Atlas cedar, such as mean temperature of warmest quarter and precipitation of driest quarter.

Insights from past projections
Our projections provided information on the past potential areas of the species in the past. Such projections can help clarify species range shifts (Benito Garzón et al. 2007) and identify past refugia (Buckley et al. 2009;Wilson and Pitts 2012) by comparing model simulations with palaeopalynological studies (Alba-Sánchez et al. 2010;Pearman et al. 2008). However, in North Africa, palaeopalynological studies are spatially limited, and some are not tied to date (Cheddadi et al. 2009). Two good syntheses summarizing the palaeorecords of C. atlantica were performed by Pons (1998) and Cheddadi et al. (2009); they and other studies have been very useful for performing comparisons and evaluating our findings. The models for past periods revealed that Atlas cedar had a potentially wider distribution than it does at present, with its maximum spatial extent recorded during the Mid-Holocene. The changes in potential area from the LGM are largely attributed to cold and aridity, which marked the climate of North Africa at that time (Rognon 1987;Stambouli-Essassi et al. 2007). We identified four putative glacial refugia of Atlas cedar that overlap putative refugia within the Mediterranean region (Médail and Diadema 2009), though Atlas cedar was not investigated by the authors: Petite kabylie/de collo, Tlemcen Mountains, Rif Mountains, and Anti-Atlas Mountains. By comparing the model of LGM with palaeopalynological data, we were able to divide the glacial refugia zones of C. atlantica into three categories: (i) Putative refugia predicted exclusively by our model (potentially useful for paleontologists). These refugia are, from East to West, Zeghouane and Djebel Serj in Tunisia, the greater part of the Algerian East (including lowlying areas), Ouled Nail Mountains, the mountains of Tiaret and Tlemcen in Algeria, Beni-Snassen Mountains, Debdou Mountains, and the Western High Atlas Mountains in Morocco. (ii) Refugia revealed by palynological investigations but not predicted by our model. This includes only the Dar Fatma locality in Tunisia (Ben Tiba and Reille 1982;Stambouli-Essassi et al. 2007), which may be due to imperfections of the SDMs regarding model transferability in time (Elith and Leathwick 2009) or the downscaling of the bioclimatic data of this period by interpolation, which originally had a resolution of 2.5 min. (iii) Refugia predicted by the model that overlap palaeorecords: Garaat El-Ouez (at El-Kala) in Algeria (Benslama et al. 2010), indicating that Atlas cedar could occur at very low altitudes (less than 100 m), and the Lac Ifrah at 1610 m in Moroccan Middle Atlas (Rhoujjati et al. 2010).
We observed good agreement between the predictions of the Mid-Holocene model and the palaeorecord localities of C. atlantica during this period. These regions are the Rif , Ifrah (Rhoujjati et al. 2010), Dayet Hachlaf , Ait Ichou (Tabel et al. 2016), Tigalmamine (Lamb and vanderKaars 1995), and Ksabi (Ballouche and Damblon 1988) with a low probability (30%). The only one unpredicted suitable site of the Mid-Holocene model was Lac Sidi Ali (Lamb et al. 1999), probably due to transferability of the CM in time.
Palaeodata indicate that the evolution of Atlas cedar geographic distribution over time from the LGM was characterized by species extinction at low and moderate altitudes in Dar Fatma (Stambouli-Essassi et al. 2007) and Djebel El Ghorra in Tunisia (Ben Tiba 1995), at Garaat-El-Ouez (Benslama et al. 2010), at Châtaigneraie in Algeria (Salamani 1993), and at Lac Ifrah (Rhoujjati et al. 2010) in Morocco. Nevertheless, C. atlantica appeared at high altitude for the first time during the Mid-Holocene at Tigalmamine (Lamb et al. 1989), Lac Sidi Ali (Lamb et al. 1999), Dayet Hachlaf , and Ait Ichou (Tabel et al. 2016).
The Mid-Holocene model showed, for example, that the Djurdjura was favorable for C. atlantica in only its eastern part while the Babors, Chelia, and Ouled Yakoub were not favorable at all during this time. These findings suggest that Atlas cedar made its first appearance at these localities after 6000 BP. Similarly, the Middle Atlas and the Eastern High Atlas Mountains became entirely favorable only in recent time, as corroborated by the predictions of the current period (CM). This interpretation is supported by the density probability graphs plotted for the quantitative environmental variables that were extracted according to current presence points of the species in all investigated time slices (see Figs. S2 and S3). Hence, we can conclude that the current Atlas cedar forests (at least in part) are the culmination of a long migration from glacial refugia along an altitudinal gradient (range shifts), leaving moderate altitudes to deciduous oaks, especially Quercus canariensis (Cheddadi et al. 2009;Salamani 1993;Stambouli-Essassi et al. 2007). Benito Garzón et al. (2007) focused on the same periods of time as our study, integrating SDMs for tree species in the Mediterranean region and showed similar range shifts.
Ongoing shifts and putative future refugia Future projections and elevation probability density graphs (see Figs. S4 and S5) indicate a continuing migration of the species towards higher altitudes (in some localities) and an important contraction of its future potential area associated with anthropogenic climate change. This is similar to future shifts in response to climate change identified for other tree species (Al-Qaddi et al. 2017;Morin et al. 2008;Thuiller 2003). Moreover, Allen et al. (2010) reported the sensitivity of tree species to drought during recent decades. Tree-ring based studies of the response of C. atlantica to recent climate variations indicate that the species is not immune to such variations (Slimani et al. 2014a;Touchan et al. 2010) which likely explains the dieback seen in many North African cedar forests, for example, in the Aures (Kherchouche et al. 2013;Kherchouche et al. 2012) and Middle Atlas Mountains (Linares et al. 2011).
Modeling outputs for future periods revealed an acceleration of range contraction between 2050 and 2070 compared with the present to 2050. Such rapid contraction may, counterintuitively, preserve high genetic diversity and induce low differentiation between Atlas cedar populations as an evolutionary strategy of the species (Arenas et al. 2012).
Future climate change might affect the fecundity of C. atlantica given that climatic conditions play important roles in the temporal variation of fecundity (Krouchi et al. 2004). The pace of this change will lead to North African cedar populations being the high vulnerability and facing likely extirpation. This is a particular concern in the Algerian Tell, where we recorded very low chances for species survival, given that the adaptation of conifers would require more than ten generations as reported by Rehfeldt et al. (2002).
The results showing putative zones of future refugia indicated an alarming spatial configuration characterized by the disappearance of Algerian Tell cedar forests and those of the Moroccan High Atlas. However, these results appear more optimistic than those from mechanistic approaches (Cheddadi et al. 2009;Demarteau et al. 2007), which predicted the Western High Atlas as the only future refugium of the species. This region and the Tlemcen Mountains will likely experience an expansion of the Mediterranean Climate Extent (MCE) sensu Klausmeyer and Shaw (2009). In view of this observation, they might constitute the most suitable areas for the ex situ conservation of the species. If RCP 4.5 is realized, other refugial areas close to the current Atlas cedar forests (Aures, Rif, Middle Atlas) will increase the possibility of natural persistence and adaptation of C. atlantica in these regions.
Overall, our contribution has shown the ability of SDMs to increase our understanding of potential distribution dynamics of C. atlantica across time. Potential palaeodistributions provided a complement to palaeopalynological studies and helped reduce the subjectivity of past geographic patterns of this species by offering a spatially explicit hypothesis at a broad scale. The current and future predicted potential distributions may greatly improve decision making in conservation planning and reforestation programmes. Special attention should be paid by conservation practitioners to (i) the resolution of the present and future output maps and (ii) the current biotic interactions and those that may arise under climate change.