Three fossil records were collected in the Rif Mountains at various altitudes and distances to the Atlas cedar forests (Figure 1). A sediment core of 8.5 m length was obtained from the Bab El Karn (BEK) wetland [35.022524 N, 5.206978 W, 1,178 masl (meters above sea level)]. A detailed pollen record has already been published (Cheddadi et al., 2016). BEK is located in the Western part of the Rif Mountains ca. 8 km SW of the closest mixed Atlas cedar/Moroccan fir [Abies pinsapo var. marocana (Trab.) Ceballos and Martín Bol.] forest in the Talassemtane National Park (TNP). The lowest-elevation Atlas cedar stands in the TNP lie at ca. 1,650 masl. Another 6.5 m sediment core was obtained from the Aanasser (ANS) peat bog (35.017629 N, 4.992503 W, 1,342 masl); located at the foot of Jbel Tizirène. Preliminary data of ANS are published (Cheddadi et al., 2015), but here we provide a more detailed pollen record. The ANS site harbors today a small isolated Atlas cedar forest (ca. 250 hectares) between 1,700 and 2,200 masl. A third 5 m core was collected from M'Had marsh (MHD) (35.128395 N, 5.438668 W, 754 masl), which was located SW of the Bou Hachem natural Park (BNP). MHD lay ca. 30 km west of the Atlas cedars in the TNP. Preliminary data of MHD have also been published (Cheddadi et al., 2015).

All three cores were collected using a Russian corer. We extracted the pollen grains from these three records using a standard chemical procedure (HCL, KOH, ZnCL₂, acetolysis). The pollen percentages were computed on a sum that excludes aquatic plants (Figures 2A–C).

The three sedimentary sequences were dated using AMS ¹⁴C dates (Cheddadi et al., 2015, 2016) and calibrated using CALIB 7.0 software (Stuiver et al., 2013) using the IntCal13 calibration dataset (Reimer et al., 2013). To establish a chronology for the three records we interpolated the calibrated ages over the sample depths using a linear fit between neighboring calibrated ages. The time spans covered by BEK, MHD, and ANS were 9.2 cal ka BP (thousands of calendar years Before Present; hereafter cal ka BP), 5.7 and 3.5 cal ka BP, respectively.

C. atlantica is an endemic species in Algeria and Morocco with its densest populations in the Middle Atlas, Morocco. The overall range of Atlas cedar in the Rif Mountains covers ca. 12,000 ha. We have refined the mapping of the range of the species via GPS data that we collected in the field. Atlas cedars in the Rif Mountains (Figure 1) occur in six main scattered populations (Jbel Kelti, Talassemtane, Jbel Tizirene, Issaguen, Oursane, and Jbel Tidighine) between ca. 1,400 and 2,300 masl, where the annual precipitation ranges roughly between ca. 500 and 900 mm (Figure 3). The minimum and the mean temperature of the coldest month (January) varies between ca. −1 and 4°C, and ca. 4 and 9°C, respectively (Figure 3).

In a study of seedling regeneration in natural forest areas at high altitudes (over 2,000 masl), field observations show that the germination process is often inhibited by both the water deficit and temperatures below −5°C (Ezzahiri and Belghazi, 2000). At lower altitudes the snow cover does not persist for long, which provides more moisture when melting. Seedlings are vulnerable to direct light and may be harmed during the summer dry season if they grow outside the canopy (Ezzahiri and Belghazi, 2000).

In the Middle Atlas, there is a clear positive effect of the west and northwest exposures on Atlas cedar germination due to prevailing westerlies that bring moisture from the Atlantic Ocean. This humidity lowers evapotranspiration, which facilitates germination and allows young seedlings to survive the summer heat (Ezzahiri and Belghazi, 2000). Conversely, southern and eastern exposures, which are generally the hottest, have a negative influence on the young seedlings.

We reconstructed January temperature (Tjan), winter (Pdjf), spring (Pmam), summer (Pjja), and annual (Pann) precipitation from the BEK pollen record, which covers the longest time span (Figure 4). The climate reconstruction is based on the assignment of pollen taxa, identified in each fossil sample, to a corresponding modern plant species. Once the fossil pollen/modern plant species assignment is established, we combine the median value of the climate range encompassed by each selected modern species and infer a climate value and its standard deviations. Similar statistical approaches have already been used to reconstruct past climate variables based on the mutual climatic range of insects (Elias, 1997), plant fossil remains (Mosbrugger and Utescher, 1997; Pross et al., 2000), mollusks (Moine et al., 2002), and ostracods (Horne, 2007) or to derive climate probability density functions (Kühl et al., 2002) from fossil pollen data.

In the present study, the reconstructed median value of each climate variable (Tjan, Pdjf, Pmam, Pjja, and Pann) is obtained using a leave-one-out approach. For each fossil sample we removed one known taxon and computed the weighted median value, using the pollen percentages as a weight, of all the remaining species. This calculation is iterated as many times as there are taxa assigned to a modern plant species in each fossil sample. The final reconstructed value for each sample corresponds to the median value of all iterations. The standard deviations correspond to the median value of the standard deviations of all iterations. This method allows us to minimize the effect of some specific taxa that are either over- or under-represented or may have a strong variation throughout the record. The method was written using R software version 3.4.0 (2017-04-21) (R Core Team, 2014) with the following libraries: akima (Akima and Gebhardt, 2016), RMySQL (Ooms et al., 2016), and stats which is part of R.

Our modern plant database includes species distributions from Flora Europaea (Jalas and Suominen, 1973, 1979, 1980) and additional data from GBIF (data.gbif.org). The modern climate variables used to define the climate range of the modern plant species, were obtained from the WORLDCLIM database (Hijmans et al., 2005) and interpolated onto the species georeferenced occurrences.

A transient simulation of the geographical distribution of C. atlantica was performed using, as Cheddadi et al. (2009), the CARAIB (CARbon Assimilation In the Biosphere) dynamic vegetation model. CARAIB is a physically-based model developed to study the role of vegetation dynamics in the carbon cycle at the global, continental or regional scale (Gérard et al., 1999; Dury et al., 2011; Fontaine et al., 2014). CARAIB includes coupled hydrological, biogeochemical and biogeographical modules, respectively describing the exchange of water between the atmosphere, the vegetation and the soil, the evolution of biospheric carbon stocks and fluxes, and the establishment, growth, competition, mortality, seed production, and regeneration of plants groups or individual species (Laurent et al., 2008; Raghunathan et al., 2015). In this study, only one tree species (C. atlantica) is considered in the overstorey. The understorey is assumed to consist of C3 herbs.

CARAIB requires the input of several climatic variables: mean air temperature, precipitation, diurnal temperature range, air relative humidity, cloud cover (converted into percentage of sunshine hours), and surface horizontal wind speed. We used CRU (Climate Research Unit, http://www.cru.uea.ac.uk/) monthly climatic anomalies (0.5° dataset, Harris et al., 2014) combined with WORLDCLIM 1951–2000 climatology (30 arc sec dataset, Hijmans et al., 2005) to obtain high-resolution (~1 km) transient climatic inputs over the period 1901–2010. Atmospheric CO₂ concentrations are prescribed from time series of global values (Meinshausen et al., 2011; Dlugokencky et al., 2014). The soil texture was obtained from the global HWSD v1.21 (30 arc sec) database [FAO/IIASA/ISRIC/ISS-CAS/JRC, Harmonized World Soil Database (version 1.2). FAO, 2012]. In our simulation we have assumed a homogeneous soil layer 1.2 m in depth. This value is not based on field data, but was obtained by model calibration. It is slightly lower than the standard rooting depth of Mediterranean woodlands in CARAIB (1.70 m, following Schenk and Jackson, 2002). The CARAIB simulations are performed at 30 arc sec resolution, which represents a significant improvement compared with the work performed with the same model at 10 arc min resolution by Cheddadi et al. (2009). The improved spatial scale more realistically captures the distribution of C. atlantica relative to topographic variability than prior models. Other differences with this previous study are that the current simulations are transient model runs over the period 1901–2010 (compared with steady-state simulations using mean climatological values in Cheddadi et al., 2009) and that the new version of the model contains a fire module (Dury et al., 2011) which impacts tree mortality.

CARAIB requires a set of parameters describing climate tolerances of the simulated species, which are controlling plant stress and germination. These factors were obtained by superimposing the above CRU/ WORLDCLIM gridded climatology interpolated at a 30 arc sec resolution with the observed distribution of C. atlantica (based on GPS data collected on the field), and then selecting a given quantile in the distribution of the climatic variables.

The potential distribution of C. atlantica simulated with CARAIB is discussed and compared with the modern geographical distribution of the species. C. atlantica is considered as present in a grid-cell when the biomass is >3 kg C m⁻². This threshold is based on the range and variability of tree productivity in the model simulation and would be among the lowest values for Mediterranean biomes according to the data over Eurasia (Schepaschenko et al., 2017).

Our transient model simulation suggests that the soil water availability has a major effect on the Atlas cedar distribution in the Rif Mountains (Figure 6). Populations at lower elevations seem to be particularly impacted (Figures 6, 7), consistent with the observed recent 200 m uphill shift of the lower treeline limit in the Middle Atlas (Rhanem, 2011). The climate data sets used show that the water stress increased in recent decades. This increase is stronger on the eastern side of the Rif Mountains (Figure 6). Consistent with this model outcome, our field observations document desiccation of Atlas cedar leaves in many populations located in Oursane (Figure 1).

The simulated potential range has decreased by 75% between 1960 and 2010, with a marked reduction between 1980 and 1990 (Figures 6, 7). The latter decade corresponds to a period of severe drought in Northern Africa (Tucker et al., 1991; Esper et al., 2007); an event that was even more marked on the eastern than the western sides of the Rif Mountains (Figure 6). The species staged a recovery between 1990 and 2000. However, it is worrying to observe that despite the fact that the last decade (2000–2010) was not as dry as that of 1980–1990, the Atlas cedar forests have gone into a rapid decline, with modern coverage as low as that of 1980–1990. Low water availability, varying according to altitude and longitude, accounts for complete absence of Atlas cedars in the eastern part of the Rif Mountains and the pronounced fragmentation of its range.

In addition to climatic disruption of populations, Atlas cedar seedlings have been undergoing strong anthropo-zoogenic impacts, especially over the last 60 years (Medail and Quezel, 1997). In the Middle Atlas, the combined human and climatic impacts have altered both the structure and dynamics of the species (Navarro-Cerrillo et al., 2013). The model simulation overpredicts probable range in the east, but it is clear that the disproportionately high loss of the eastern population tracks a drying climate (Figure 6). The discrepancy between the observed and simulated ranges may be linked to the additional effect of human activities which had a major impact on all mountain forests in Morocco (Reille, 1977; Lamb et al., 1991; Lamb and van der Kaars, 1995; Ajbilou et al., 2006), accelerating the Atlas cedar extinction along the lower edge of its range (Cheddadi et al., 2015). Such human influences are thought to underlie several periods of deforestation over the last three millennia in many North African forest ecosystems (Mercuri et al., 2011; Alba-Sánchez et al., 2015).

Empirical data and our simulations clearly show that the soil water content has a direct impact on the overall distribution of Atlas cedar. Atlas cedar seedlings may fail to establish either due to the lack of moisture (Epron, 1997) or direct exposure to the sun during the summer dry season (Ezzahiri and Belghazi, 2000). Takos and Merou (2001) have shown that the germination of Cedrus deodora seeds is less effective when temperature is higher than +5°C. In the Rif Mountains, the minimum temperature of the coldest month (January) today is below +5°C and the mean temperature is slightly higher than that value (Figure 3). If the germination of C. atlantica seeds is controlled by similar temperature thresholds as those of C. deodora, then a further increase could represent a threat to its germination. Both these winter and summer conditions appear to become more frequent as a result of anthropogenic climate destabilization. While we recognize that a multitude of other factors will influence treeline changes, such as, fire frequency, carbon and nitrogen availability in the soil, tree phenology, topography, exposure, competition with other species, and soil composition, there is not time to investigate all parameters fully before conservation measures are enacted.

Although it is the climate that drives the vegetation changes, one should note that in this study the past climate variables have been inferred from the pollen data. To avoid circularity we cannot interpret them as drivers of any vegetation change.

The persistence of Atlas cedars in the western part of the Rif Mountains (Jbel Kelti, Talassemtane, and Jbel Tizirene) is favored by the existence of microclimates modulated by the influence of moisture from the Atlantic Ocean and the Mediterranean Sea. The reconstructed climate variables from BEK record (Figure 4) show that when Atlas cedars dominated the landscape (pollen percentages between 20 and 80%) between 7.5 and 5.5 cal ka BP climate was cooler and wetter. Between 8 and 6 cal ka BP the upwelling along the western African coast, including offshore Morocco, was stronger than today. The strengthened upwelling led to a transition to cooler sea surface temperatures along the Western African coast (DeMenocal et al., 2000) including that of Morocco. These oceanic conditions induced a cooling effect over northern Africa (McGregor et al., 2007). Simultaneously evaporation over the warmer-than-modern Alboran Sea (Cacho et al., 2001) between ca. 7.5 and 6 cal ka BP (Figure 4) increased atmospheric moisture transport. This additional moisture input may explain the higher water availability over the Rif Mountains in the mid-Holocene.

Pollen records from the Rif Mountains as well as from the Middle Atlas (Reille, 1976, 1977; Lamb et al., 1995, 1999; Cheddadi et al., 1998, 2009, 2016; Nour El Bait et al., 2014; Tabel et al., 2016; Zielhofer et al., 2017) show that Atlas cedar forests have experienced considerable variations since the last glacial and during the Holocene period. The reconstructed climate from the Tigalmamine pollen record (Figure 4) shows that an expansion of Atlas cedars took place under a cooler and wetter climate than today in the Middle Atlas after 7.4 cal ka BP (Cheddadi et al., 1998). In the Rif Mountains, populations of Atlas cedars were more widespread than today between 9 and 2 cal ka BP, extending to lower elevations in BEK (1178 masl, Figure 2A), and close to MHD (754 masl) between 5.5 and 2.5 cal ka BP (Figure 2B). Atlas cedar does not occur today at these elevations. However, the interpretation of the low pollen occurrences in the MHD record requires some caution because of the potential input from remote populations. The BEK site and the TNP are geologically connected through a corridor that lies at ca. 800 masl (Figure 5). Thus, a lowering of the tree line limit may have allowed the Atlas cedar to reach the BEK site and therefore cover a wider range than today. Another pollen record collected in the BNP (1100 masl) north of MHD (Reille, 1977), indicates that Atlas cedars formed extended populations in the area until 810 ± 50 BP. Two high-resolution and well-dated pollen records from Ifri Oudadane (Zapata et al., 2013) and Ifri n'Etsedda (Linstädter et al., 2016) suggest that Atlas cedar expanded its range between ca. 10 and 6 cal ka BP. The expansion to the north-east and to lower elevations than the 800 masl isoline (Figure 1), penetrated areas that are too warm and dry to support populations today. The 600 m upslope erosion of lower montane populations between 800 masl and the modern lower treeline limit at ca. 1,400 masl seems to have taken place progressively between 6 and 2 cal ka BP.

Our climate reconstruction shows that there was a decrease of ca. 2°C between ca. 7.5 and 6 cal ka BP (Figure 4). Taking into account the low dispersal capacity of the Atlas cedar pollen grains (Wright, 1952), the high pollen percentages during that time span indicated that the species substantially expanded its geographical range in comparison with its modern one (Figure 1). After 6 ka, we observe a decreasing trend of pollen occurrence, which translates into a substantial range contraction of the Atlas cedar toward higher elevations. The range contraction between 6 and 2 cal ka BP corresponds to changes of both the winter temperature and the seasonal distribution of precipitation. Earlier simulations using a vegetation model suggested that an increase of ca. 2°C in the annual temperature over the Mediterranean would not have had a major impact on the temperate conifer forest biome (which includes Atlas cedars) unless the annual amount of precipitation decreased drastically (Cheddadi et al., 2001).

The temperature gradient in mountainous landscapes is a major factor controlling species dynamics (Körner, 2016) through speciation and extinction processes (Graham et al., 2014). Temperature is also considered as the primary climate variable controlling the treeline (Hessl and Baker, 1997). More precisely, winter temperatures have a direct influence on soil temperature, which is an important factor in determining treeline globally (Körner and Paulsen, 2004). For example, winter temperature increase caused treeline to move upslope by 115 m in the European Alps between 1901 and 2000 (Leonelli et al., 2011), and it is expected to cause further impacts over the next century (Lenoir et al., 2008; Engler et al., 2011; Moritz and Agudo, 2013).

The modern microrefugial areas in the western Rif Mountains are still cooler and wetter (with snow occurrence during winter) than the North African climate average. The future persistence of Atlas cedar will rely on the local climate stability (and the management of the direct human impacts) in these microrefugia relative to the Mediterranean climate that is expected to increase by 4–6°C (IPCC, 2014). Note, this increase is two to three times higher than the reconstructed temperature shift of the last 9,000 years reconstructed in this study. Such climate scenarios are expected to have a major impact on Mediterranean biodiversity hotspots (Malcolm et al., 2006). One approach to evaluate the pertinence of the modern microrefugia for the conservation of Atlas cedars is to quantify the local topoclimate in each microrefugium, using data loggers, and compare them with the regional climate (Ashcroft et al., 2009, 2012). The spatial resolution of the climate data used in the present work does not allow us to evaluate the impact of a local climate on a specific microrefugium. Micro-climate monitoring is a key issue for a more accurate evaluation of the local climate stability and suitability for the persistence of the Atlas cedar in their modern microrefugia.

Besides the instrumental climate monitoring, we need a thorough genetic survey of populations within all the modern microrefugia. Presently, the DNA surveys of the modern Atlas cedars in Morocco have focused mainly on its overall genetic diversity (Renau-Morata et al., 2005; Terrab et al., 2006; Cheddadi et al., 2009) which makes it difficult to identify areas with high adaptive capacities within the modern range. Genetic studies suggest that long-term forest exploitation leads to a loss of tree genetic diversity, possibly including rare alleles which are important for population adaptation to environmental change (Pautasso, 2009). Generating novel genetic variation in populations through artificial hybridization and introgression is one of the ways to increase genetic diversity of Atlas cedars (Fady et al., 2003). A phylogeographic study of two widespread Nothofagus species in Patagonia has shown that hybridization may have prevented the species from going extinct over the Quaternary climate cycles and that hybridization should be considered as a conservation strategy (Wolf et al., 2001; Soliani et al., 2012). Another approach based on the analysis of gene markers such as, allelic richness (Petit et al., 1998) should help identifying populations that may be better adapted to environmental changes and therefore prioritized for conservation.

A genetic survey of the modern populations suggests that three populations found on different mountain ranges within Morocco were genetically distinct. Within each mountain range, however, there was no clear geographical diversity structure (Terrab et al., 2006) and that each montane population may have been isolated from one another for a considerable period of time (Renau-Morata et al., 2005; Terrab et al., 2006). The lack of a clear spatial genetic structure within each mountain range could indicate moderate genetic flow between sub-populations. Alternatively, it could indicate the relative success of one genotype compared with others that went extinct during adverse climatic conditions. A combined study of pollen records and the modern genetic diversity of Atlas cedar indicated that it may have persisted in several areas in the Moroccan mountains through local altitudinal migration (Cheddadi et al., 2009). The lack of a clear spatial genetic structure may be related to the persistence of the species in a discrete range where some populations may have gone extinct under less suitable global climate while others have either adapted locally into microrefugia or migrated into climatically more suitable areas.

Metapopulation theory (Hanski, 1994) predicts that species persist in fragmented landscapes as discrete sub-populations occupying only part of the suitable habitats (Hanski and Ovaskainen, 2003). A comparison of the simulated range of Atlas cedars for the most recent time span (2000–2010) with the observed occurrences, shows that the range is highly fragmented and that the species does not occur in all simulated suitable habitats (Figure 6). Each sub-population within a metapopulation has a finite probability of chance extinction. Thus, populations can appear and disappear from the overall landscape without necessarily threatening the long-term persistence of the species itself. The metapopulation concept was originally designed for animal organisms (Hanski, 1994) with a strong turnover that is related to high dispersal and migration capacities over a short life-time span (Hanski, 1991), which allows the observation and study of several generations. Long-lived plant species (such as, Atlas cedars) are static with low turnover and their dispersal capacities are much lower both in time and space than those of animals (Ouborg and Eriksson, 2004). Deriving a full metpopulation-based approach to considering microrefugia (Mosblech et al., 2011), would require more data on the spatial structure, migration capacity, seed dispersal, extinction rate, turnover, gene flow, than is currently available for Atlas cedars. However, what can be stated is that based on metapopulation theory, the six current populations of Atlas cedar in the Rif Mountains cannot be assumed to be permanent features of the landscape even in the absence of climate change. Equally, from a theoretical perspective, there is no a priori reason why Atlas cedar should not colonize an unoccupied suitable patch. In practice, management to establish populations on unoccupied, but suitable sites would be a safeguard against potential loss.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.