Microrefugia, Climate Change, and Conservation of Cedrus atlantica in the Rif Mountains, Morocco

Introduction

The genus Cedrus has been present in the Eastern Mediterranean for more than 23 Ma (Biltekin et al., 2015). An ancestral species of the Himalayan cedar, Cedrus deodara, initially diverged into lineages that exist today as Cedrus libani and Cedrus brevifolia (Bou Dagher-Kharrat et al., 2007), before deriving a North African species, Cedrus atlantica around 8 Ma (Qiao et al., 2007). Pollen data confirm the occurrence of the genus Cedrus in Morocco since the Messinian period (ca. 7–5 Ma) and throughout the Pliocene and Pleistocene (Feddi et al., 2011; Magri, 2012). Although species-level identification is not possible from the pollen, the fossil cedar species is presumed to be C. atlantica in NW Africa.

North African plant species are highly restricted in their migratory space when climate switches between cold and warm periods. The narrow strip (<5° latitude) between the Mediterranean Sea and the Saharan Desert currently supports rather arid climates, though this has not always been the case (Tierney et al., 2017). In this setting, mountains provide habitat heterogeneity and a variety of microclimates. Plant species have survived past climate changes by adjusting their ranges locally and/or adapted physiologically/genetically to the new conditions. Although the Sahara was wetter than today between 11,000 and 5,000 BP up to 31°N (Tierney et al., 2017), it did not contribute to the long-term persistence of the Mediterranean species. Under this more restrictive set of conditions species either went extinct or adapted. Recently, UNESCO declared the mountain landscape of Morocco to be a critical habitat for conservation, justifying the recent inclusion of the Atlas cedar Biosphere Reserve. In addition to its natural value, these are “cultural landscapes” that have been used by people for millennia, adding another reason to foster their protection (Sarmiento, 2011). Thus, understanding how a species such as, Atlas cedar persisted in the Rif Mountains, despite profound past climate change, is essential for its long term conservation.

Mountain areas are of particular interest to conservation biologists as they have played a major role in the survival of tree species during the Ice Ages by providing habitats with more ecological stability, i.e., availability of moisture, varied topography, and short dispersals to maintain climatic equivalence (Tzedakis et al., 2002; Hewitt, 2004). High habitat heterogeneity increases the probability that a viable set of growing conditions exists for any given species within a realistic migratory distance. Whether those conditions persist through time, however, is another matter. Conservation strategies for mountain tree species in a warmer climate are either to establish migration corridors and/or establish “stepping stones” populations that could link current and future ranges (Hannah et al., 2007) or to provide assisted migration in which individuals are moved to new habitats of predicted suitability (Aitken et al., 2008; Rehfeldt and Jaquish, 2010). These two options require that we know with confidence the climatically optimal areas for the survival of species. An alternative option is to help species populations persist in areas where growing conditions may remain locally suitable.

Model simulations could help to depict potentially suitable areas and evaluate the risk of extinction based on a variety of climate scenarios (Thuiller et al., 2005). These models, however, are based on the relationship at time t between the range of a species and its contemporaneous climate. Individuals of long-lived tree species may have undergone significant climate change during a lifetime that could encompass several centuries. For example, some individuals of C. atlantica (Manetti), which can live for 1,500 years, withstood strong precipitation variability over the last nine centuries (Till and Guiot, 1990; Esper et al., 2007). Indeed, conservation policies of such a long-lived organism, would benefit from a temporal perspective of the relationship between the species and its past environments (Willis and Birks, 2006; Willis et al., 2007; Willis and Bhagwat, 2010) including a quantification of past climate variability.

Palaeoecological studies showed that during glacial periods in the temperate latitudes, the geographical distribution of all tree species shrank to scattered populations that persisted in areas known as glacial refugia (Bennett et al., 1991). These refugial areas, which were not necessarily spatially delineated, were often located in mountainous terrain (Hewitt, 2000; Taberlet and Cheddadi, 2002; Bennett and Provan, 2008; Hughes et al., 2011; Keppel et al., 2012). Over the late Quaternary, these refugia played a major role in maintaining and shaping modern biodiversity (Davis and Shaw, 2001; Tzedakis et al., 2002; Rull, 2011). At a more reduced geographical scale, there were also microrefugial habitats that may have harbored small or reduced populations of a species during periods when the broader population collapsed or migrated. Local microclimates buffered these isolated populations relative to the surrounding regional or global climate (Rull, 2009; Dobrowski, 2011). Microrefugial habitats may have existed either during cold or warm periods. The restricted range and the time span over which populations remained isolated in a microrefugium may have promoted its local adaptation/evolution compared with their original wider-ranging population.

In the present study, we aim to identify microrefugial areas in the Rif Mountains in Morocco where the Atlas cedar persisted over the last millennia. Identifying such areas may be important for preserving this endangered tree in the future. Our approach is based on the analysis of the relationship between climate and the species distribution using (1) a set of fossil pollen records to reconstruct its range changes over the Holocene; (2) a quantitative reconstruction of several climate variables to identify and compare the local climate changes with more regional or global trends and (3) a set of simulations of the species distribution over the last decades using instrumental climate data as inputs into a process-based vegetation model. One of the rationales for this study is to explore conservation actions to preserve Atlas cedars in the identified modern microrefugia.

Materials and Methods

Fossil Records

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).

Chronological Frame

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.

Modern Range and Basic Climate Requirements of Atlas Cedar

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.

Climate Reconstruction Approach

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.

Vegetation Model Simulations

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).

Results

The variation of pollen percentages through time in a fossil record may provide information about the first known occurrence of species around the site, their expansion, regression or local extinction. In general, the closer the species is to the studied site the higher is its pollen percentage in the fossil sample. Obviously, such pollen/plant relationships are modulated by the population density, pollen production of the species, and landscape topography as well as the wind velocity and direction. One should keep in mind that one fossil sample may integrate the pollen grains deposited in a decade or more. Atlas cedars are relatively poorly dispersed pollen grains with Wright (1952) estimating the great majority of grains falling within ca. 800 m of their source. Beyond this distance, the pollen percentages of Atlas cedars fall to <1% in sediment samples (Hajar et al., 2008). Magri and Parra (2002) have shown that Atlas cedar pollen grains may disperse to much greater distances from the originating population due to atmospheric circulation. However, these occurrences did not reach 1% of the total pollen sum in their samples.

Based on these general considerations, the three pollen records from the western part of the Rif Mountains indicate that populations of Atlas cedars occurred during the Holocene period (Figure 2A) at a much lower altitude than today. The record from MHD (at 754 masl) shows that this species grew close to the site between 5.6 and 2.5 cal ka BP (Figure 2B, pollen percentages >1%). Afterwards, its occurrence remained intermittent until ca. 1.25 cal ka BP before becoming extinct from around the MHD site. The BEK record (at 1,178 masl) shows consistently high pollen percentages (15–80%) between 9.2 (base of the record) and 1.6 cal ka BP; prior to it also becoming locally extinct. Today, the BEK site is surrounded by Quercus pyrenaica (Willd.) and the closest Atlas cedars are located in TNP, ca. 8 km from the marsh (Figure 5). The overall pattern of range contraction of the Atlas cedar in the Rif Mountains is confirmed in the ANS pollen record. This record, collected at 1,342 masl, also shows a reduction of the cedar representation after 1.6 cal ka BP. However, Atlas cedar continues to represent about 1% of the pollen assemblage until today. This finding is consistent with isolated populations of cedar occurring in Jbel Tizirene <1 km from the ANS site (Figures 1, 5).

Climate reconstruction from the BEK fossil pollen record (Figure 4) shows that Tjan fell between 7.5 and 5.7 cal ka BP, with a minimum value around 6 cal ka BP, where average temperatures were about 2°C lower than the long-term pattern of the Holocene. The period from 7.5 to 5.7 cal ka BP was also the wettest within the last 9 cal ka BP with a marked increase of spring rather than winter precipitation. This coolest and wettest period between 7.5 and 5.7 cal ka BP corresponded to the maximum expansion of C. atlantica. In this interval, the pollen representation of C. atlantica was continuously higher than 20%, with a maximum of 80% recorded at ~6.2 cal ka BP (Figure 2A).

Our model simulation (Figure 6) shows the decadal change in the distribution of C. atlantica over the last 50 years (1960–2010), together with the minimum of available soil water calculated with the hydrological module of the CARAIB model in relative units (i.e., in terms of the variable (SW-WP)/(FC-WP) where SW, WP, and FC are the soil water content, the wilting point and field capacity in mm, respectively). This improved bucket model (Hubert et al., 1998) evaluates the average soil water content in the root zone, the snow cover and all related water fluxes. The inputs of the soil water reservoir are the rainfall (rainfall reaching the ground) and the snowmelt. The output fluxes of the reservoir are the evapotranspiration, the drainage and the surface runoff. The soil water (SW) is allowed to vary between the wilting point (WP) and saturation.

The simulated distribution of Atlas cedar mirrors the observed range at altitudes higher than 1,400 m over the western part of the Rif Mountains (Figure 7). The presence of C. atlantica is also simulated in the eastern part of the range, where it is not observed in the field. These grid cells represent areas where climatic conditions would theoretically allow the occurrence of the plant species. The absence of Atlas cedars in these simulated potentially suitable areas may be linked to human disturbances or other factors not taken into account in the model simulations (e.g., migration and competition with other species).

The transient model simulation shows a gradual reduction of the potential distribution between 1960 and 2010, particularly in the eastern part of the Rif Mountains (Figure 6). The area of the modeled distribution decreases by about 75% over 50 years to a mean area of 10,278 ha in the last decade (Figure 7), which corresponds quite well to the modern coverage. A marked reduction is observed between 1980 and 1990, followed by a slight regeneration between 1990 and 2000. The evolution of the distribution is mainly driven by soil water availability, which has been reduced significantly over the last 50 years. A decrease in the minimum available soil water to levels below the water stress threshold for the species (set to 9.5% for the relative soil water in the model simulation) should indicate a loss of suitable habitat. In the simulation, this pattern of drought is evident in many grid points in the eastern part of the distribution after 1980; explaining the disappearance of the species in these areas.

Discussion

Today, C. atlantica is considered to be an endangered species, and its population continues to fall (A2cd ver 3.1; IUCN, 2017). As with many tree species characteristic of montane forests, Atlas cedars are being threatened by both abrupt climate change and anthropogenic activities. Although this species may withstand some drought (Aussenac and Valette, 1982; Aussenac and Finkelstein, 1983), the ongoing climate aridification and the more frequent occurrence of extreme soil drought (Ladjal et al., 2005) is already impacting its range. C. atlantica is dying out in the drier portions of its range leading to a net upslope retreat along its lower elevation limit in the Middle Atlas (Rhanem, 2011). Tree-ring data evidence the negative influence on its growth of recurrent droughts and the steep temperature rise since the 1970s (Linares et al., 2011), underlining the negative impacts of ongoing global warming.

Simulated Range Changes over the Past Decades

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.

Climate Change and Range Shift over the Holocene

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).

Persistence of the Atlas Cedar in a Few Modern Microrefugia

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.

Conclusions

Paleoenvironmental data from the Rif Mountains and model simulations show clearly that Atlas cedars had a more extensive range during the Holocene and even in the last 50 years than today. Both paleoecological data and model simulations show that Atlas cedar populations are declining in Northern Morocco due to a combined effect of climate and human impact. The eastern populations of the Atlas cedars seem to be clearly more threatened than those in the western part of its range. The modern remaining populations in Jbel Kelti, Talassemtane, Jbel Tizirene, Issaguen, Oursane, and Tidighine represent microrefugial areas within which the temperature and precipitation ranges seem to be still within the species climatic range. The persistence of the Atlas cedar in these microrefugia in the Rif Mountains may be considered as a variant on traditional views of metapopulation dynamics where the long-term stability of the local climate is a temporary substitute for gene flow between populations in permitting species survival (Mosblech et al., 2011). Such genetic isolation cannot continue indefinitely without substantial harm from genetic drift and in-breeding. Hence, as observed for other species' metapopulations (Hanski, 1998, 2005), the habitat fragmentation of Atlas cedar, either natural or human-induced and the progressive reduction of some populations will lead to a reduction of gene flow between populations and therefore may increase the chances of its extinction.

Our simulation underscores the threat posed by the expected 4–6°C temperature increase and more than 20% decrease of the annual amount of precipitation in the Mediterranean predicted by 2100. To maximize the probability of long-term survival of Atlas cedar we suggest (1) that all the modern remaining microrefugia should be protected from human activities, and in particular to protect seedlings at the upper treeline limit; (2) the climate stability and suitability in all microrefugia should be evaluated (Hylander et al., 2015) through continuous climate monitoring, using climate data loggers (Ashcroft et al., 2012) and collecting climate time series from sites as close as possible to microrefugia; (3) perform an exhaustive DNA survey, using different genetic markers (Petit et al., 1998) for evaluating the genetic diversity, the risk of extinction, and to further prioritize populations for conservation, and (4) sites identified as potentially favorable microrefugia under conditions of 2,100, but presently unoccupied by Atlas cedar, should become recipients of seeds or seedlings. Ultimately, we should consider some artificial gene introgressions between the microrefugia so that we can improve heterozygosity and the capacity of that hybridized generation to adapt in situ (Wolf et al., 2001; Fady et al., 2003; Soliani et al., 2012).

Author Contributions

RC designed the study, performed the past climate reconstructions, made the figures and has written the original manuscript. AJH and LF performed the model simulations. RC, AJH, and LF have integrated the interpretations of the model simulations into the text. MB has corrected and improved the last version of the manuscript. All co-authors have contributed to improving the manuscript in the frame of the VULPES project funded by the Belmont-Forum consortium.

Conflict of Interest Statement

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.

Acknowledgments

This work is a contribution to the Belmont Forum funded project VULPES (Project ID: ANR-15-MASC-0003). We acknowledge support for this research from F.R.S.-FNRS under research grant FRS-FNRS X.3041.17—VULPES-ULg. We thank Juanma Rubiales, Henry Lamb and William Fletcher whose comments helped improve and clarify this manuscript.

References