High genetic diversity and low future habitat suitability: will Cupressus atlantica, endemic to the High Atlas, survive under climate change?

Cupressus atlantica is a narrow endemic species of semi-arid and sub-humid habitats in the western High Atlas, Morocco. We explored the possible dynamics of the species’ range under climatic changes using species distribution modelling (SDM) to identify populations vulnerable to range changes. Additionally, we investigated the spatial genetic structure (SGS), the effective population size and genetic connectivity in natural populations, which may provide important data on demo-genetic processes and support the conservation management of this critically endangered species. The SDM results showed that the current species range constitutes only 49% of the potential distribution. Under the most pessimistic scenarios (RCP6.0 and RCP8.5), a significant reduction in the species range was predicted. However, the projection based on RPC4.5 revealed possible extensions of the habitats suitable for C. atlantica. Potentially, these areas could serve as new habitats that could be used with the assisted migration approach aiming to mitigate the current fragmentation. In terms of the SGS, significant and positive aggregation of relatives was detected up to ca. 100 m. In comparison to other fragmented and endemic species, the detected SGS was weak (Sp = 0.0112). The estimated level of recent gene flow was considerable, which likely prevented a strong SGS and allowed diversity to accumulate (HE = 0.894). The most alarming results concern the effective population size, which was very low in the studied populations (< 53), suggesting a possible increase in inbreeding and loss of diversity in the near future. More effective conservation actions integrating in situ and ex situ measures should be undertaken to prevent extirpation of the species.


Introduction
Biodiversity boosts the productivity of natural ecosystems and is positively related to ecosystem stability and resilience (Oliver et al. 2015). Hence, one of the serious challenges for conservation biology is global biodiversity loss, which affects the functioning of natural ecosystems and brings about great risks for human wellbeing (Thomas et al. 2004;Bellard et al.

Communicated by Wolfgang Cramer
Supplementary Information The online version of this article (https:// doi.org/10.1007/s10113-020-01711-9) contains supplementary material, which is available to authorized users. 2012; Cardinale et al. 2012). Currently, global climate change is one of the major threats to biodiversity, leading to its rapid depletion, which may be difficult to reverse (Thomas et al. 2004;IPCC 2013). This threat has led to increasing interest in conservation studies focused on the development of effective approaches, concepts and practices aimed at the preservation of diversity at the ecosystem, species and genetic levels (Bellard et al. 2012). Proper identification of conservation priorities is a crucial stage in any conservation programme. For example, conservation of natural populations of an endangered species requires reliable data on its population abundance, dynamics and genetic resources. This information is then a basis for assessments regarding species vulnerability and adaptability that are considered during conservation planning (Frankham et al. 2002;Aravanopoulos 2016;Fady et al. 2016a, b).
The Mediterranean is one of the 36 biodiversity hotspots of global importance due to its high species diversity, high endemism and alarming rate of biodiversity loss (Médail and Baumel 2018;Myers et al. 2000). However, this region struggles with strong climate-driven environmental changes. Progressive aridification is the most critical threat in the Mediterranean (Lionello and Scarascia 2018). According to the Mediterranean Experts on Climate and Environmental Change network, the warming effect in this region is 20% faster than the global average (Klausmeyer and Shaw 2009;MedECC 2019). Thus, the region, especially the mountainous areas, which support the majority of the biological diversity, is of particular interest to conservation biologists (Médail and Quézel 1997). Mediterranean mountains served as major refugia during the Pleistocene glaciations for species of Tertiary and Quaternary floras (Médail and Diadema 2009) and acted as centres of diversification (Magri et al. 2007;Linares 2011;Molina-Venegas et al. 2017). These refugial areas critically contributed to the current patterns of diversity and differentiation of the Mediterranean and European woody species (Hewitt 2000;Dering et al. 2017).
However, mountain ecosystems are currently considered particularly susceptible to climate change, which entails a great risk for rare and narrow-ranged woody species growing in these regions (Parmesan 2006;Nogués-Bravo et al. 2007). Because of their geographic isolation, limited range size and unique adaptations to local conditions, endemic tree species are seriously threatened by environmental changes because they may exist in only stable climate refugia (Kaya and Işik 2002;Harrison and Noss 2017). It is highlighted that some species may not be able to keep pace with rapid climate changes (Parmesan 2006;Nogués-Bravo et al. 2007). Therefore, the preservation of rare and narrow-ranged endemic tree species raises the highest concern (Myers et al. 2000; Thomas et al. 2004). To successfully prevent extirpation of this specific group of trees, conservation efforts require detailed assessments of the factors contributing to their rarity and more robust predictions of future changes in distributions (Guisan et al. 2013).
Within the Mediterranean, a considerable contribution of plant species richness and endemism is concentrated in the mountains of North Africa (Maghreb), including the High Atlas in Morocco (Médail and Quézel 1997;Fennane and Ibn Tattou 1999). In Morocco, 22% of all plant species are endemic, and over half of them are endangered (Rankou et al. 2013;Fennane and Rejdali 2019). Nearly half of these species occur only in the High Atlas, where most of the 22 northern African Key Biodiversity Areas and the eight Important Plant Areas are situated (Taleb and Fennane 2011;Rankou et al. 2018). The region is of high conservation significance because of the exceptional degree of habitat loss and the high abundance of threatened species. According to climatic simulations, the trend of increasing aridification that has been observed for the second half of the twentieth century in North Africa is expected to continue. However, climate-driven environmental changes will likely have a stronger effect in Morocco than in other areas of the region since the decrease in precipitation may be as high as 53% by the end of this century (Schilling et al. 2012;Verner et al. 2018;Hssaisoune et al. 2020). Humid, sub-humid and semi-arid areas of North Africa, which have experienced shrinkage since the 1970s, are most susceptible to those changes (Verner et al. 2018). Apart from climate changes, ongoing overexploitation of the mountain Moroccan woodlands over the last few millennia also causes a serious problem for endangered species (Esteban et al. 2010;Cheddadi et al. 2017;Taib et al. 2020). Such a situation calls for urgent conservation action to preserve those unique species assemblages that represent the evolutionary heritage of the Mediterranean region.
In recent decades, conservation challenges have been successfully supported with various scientifically grounded methods and techniques. Among them, species distribution models (SDMs) are increasingly applied in conservation planning and management supporting decision making (Franklin 2010;Guisan et al. 2013). These tools are widely used to predict past and future species distributions and may help to identify the areas with the highest conservation priority and to evaluate the extinction risk of species (Guisan et al. 2013). Additionally, these methods have great potential to effectively support assisted migration, which is one of the possible methods of species conservation under climate change that seems to be especially appropriate for highly fragmented landscapes (Vitt et al. 2010;Gallagher et al. 2015). Assisted migration assumes deliberate movement of species populations according to the predicted direction of range shifts to mitigate the effects of climate change. SDMs are thus able to define climate-suitable habitats under current and future conditions; therefore, these methods support the decision of whether translocation is necessary and help to locate potential recipient sites (Gallagher et al. 2015). Therefore, the SDM approach is now a vital element of studies aiming at delivering reliable projections of climate-induced distributional changes for species that will allow tailored conservation strategies to be developed (Guisan et al. 2013).
Effective conservation strategies may benefit greatly from the knowledge of the level and distribution of genetic diversity of targeted species (Hoban et al. 2018). High intraspecific genetic variation is crucial for the ability of species to cope with changing environments in the future and establishes the fundamental criterion for conservation (Fady et al. 2016a, b;Aravanopoulos 2016). However, the maintenance and conservation of genetic resources also require comprehensive knowledge of additional genetic criteria that include parameters such as inbreeding coefficients, the effective population size, the rate of gene flow and the spatial genetic structure (SGS). Understanding the drivers of the spatial organization of genetic diversity is crucial for the effective targeting and sampling of populations for in situ and ex situ conservation. Thus, conservation genetics play a key role in developing a strategy for short-and long-term management and preservation of endangered species and their genetic resources (Pertoldi et al. 2007). Mapping genetic variability in the spatial dimension may efficiently support the implementation of conservation management, such as assisted migration or genetic rescue.
Cupressus atlantica Gaussen (Moroccan cypress) is a narrow endemic long-lived conifer that is a component of open forest communities in the semi-arid and locally sub-humid areas in the western High Atlas, Morocco (Fig. 1a). The range of this species is severely fragmented, and the major distribution area covers mainly the Oued N'Fiss, where the most representative stands are the Aghbar Forest, Ijoukak and Talat N'yaâkoub, with several small scattered stands or groups of individuals located in Allous (Assif Ougdemt), Taghzout (Assif N'aït Toumat) and Mzouzit (Oued N'Fiss). The species occurs at elevations between 900 and 2200 m a.s.l. (El Wahidi 2004;Gardner and Griffiths 2013). Moroccan cypress has undergone intensive decline during the twentieth century (Griffiths 1998). Currently, the total area of species occurrence is estimated as approx. 1960 ha, but it is still decreasing. According to the International Union for the Conservation of Nature (IUCN), C. atlantica is classified as critically endangered, and the Food and Agriculture Organization (FAO) classifies the species as one of the 17 world forest species whose genetic diversity is impoverished (FAO 1976;Gardner and Griffiths 2013). As a result of the negative climatic predictions anticipated for North Africa (IPCC 2013), the decline of the C. atlantica population will likely accelerate without any conservation intervention, finally leading to species extinction. The scenario of profound range reduction has recently been reported for other endangered North African conifers (Arar et al. 2019;Bouahmed et al. 2019;Taib et al. 2020). According to recent studies focused on the evolutionary history of Mediterranean cypress taxa, the genetic diversity of C. atlantica currently remains relatively high (Sękiewicz et al. 2018;Bagnoli et al. 2020). However, there is a justified concern for the loss of variability over the short term in the face of demographic fluctuation. The question of whether C. atlantica will cope with these projected changes or may share the same fate of being near extinction with C. dupreziana, for which the prospects are rather pessimistic (Abdoun et al. 2013).
Here, we combined species distribution modelling and population genetics methods to efficiently support further conservation planning and management of critically endangered C. atlantica. An SDM was applied to (1) predict the impact of future climate change on the species distribution and evaluate the potential risks involved and (2) estimate the suitable areas that may serve as potential new habitats for the species in the future in frame of the assisted migration concept. Finally, we characterized the genetic diversity and SGS within the natural populations of C. atlantica, which may provide important data on demo-genetic processes and support the conservation management of the species. Given that fragmentation intensifies the SGS, we expected to find a significant strong SGS within

Species distribution modelling
The maximum entropy approach implemented in MaxEnt 3.4.1 (Phillips et al. 2004) was applied to construct the models for the past, present and future potential habitat suitability for C. atlantica. The data on the occurrence of the species were accessed from the chorological study presented by Sękiewicz et al. (2014). In total, the datasets host 214 georeferenced records, but these occurrences were validated to reduce spatial autocorrelation, and only one occurrence point per grid cell (1 × 1 km) was considered, resulting in 88 unique occurrence points used for modelling. The 19 bioclimatic variables (Bio1-Bio19) for both current climatic conditions  and future projections (2061-2080) were obtained from CHELSA 1.2 (Karger et al. 2017) with a 30 arc-sec resolution, while for past projection including the late Holocene (Meghalayan; 4.2-0.3 ka) from PaleoClim (Fordham et al. 2017;Brown et al. 2018) with a 2.5-min resolution. The raster layer of soil data based on the World Reference Base (WRB) classification was downloaded via SoilGrids™ (Hengl et al. 2017) at 250-m 2 resolution and upscaled to match the resolution and extent of the bioclimatic variables. To reduce the multicollinearity effects of environmental variables on the species distribution projections, Pearson's correlation test (p < 0.05) was applied with the raster.cor.matrix function implemented in the ENMTools package in R (Warren et al. 2019). Thus, the highly correlated variables ((r) > 0.9) were removed.
To evaluate the possible future distributional changes in the C. atlantica range, climate projections obtained from the Coupled Model Intercomparison Project Phase 5 (CMIP5) were used (Collins et al. 2014). Three representative concentration pathway (RCP) scenarios were considered for the CCSM4 model (Gent et al. 2011;Collins et al. 2014) (Online Resource 1). For intermediate scenarios, the global mean surface temperatures and the atmospheric CO 2 concentration for 2081-2100 were projected to increase by 1.1 to 2.6°C and 650 ppm in RCP4.5, while in RCP6.0, they will increase by 1.4 to 3.1°C and 850 ppm. Under a pessimistic scenario, RCP8.5, the increase in the global mean surface temperature is likely to be 2.6 to 4.8°C, and the CO 2 concentration will be 1.350 ppm by the end of the twenty-first century (Collins et al. 2014).
The modelling was conducted with the 'random seed' option, where 20% of the occurrence points were randomly sampled and set aside as test data, the remaining points were used as training data, and a random test partition was used for each run. Bootstrapping analysis was run with 100 iterations, the maximum number of iterations was set at 10 4 and the convergence threshold was set at 10 −5 with the logistic output of the model prediction for suitability. To evaluate the accuracy of the model, the area under the curve (AUC) values of the receiving operator curve (ROC) were used as a thresholdindependent evaluation metric (Wang et al. 2007;Mas et al. 2013). The resulting distributions were projected with QUANTUM GIS 2.18.26 Lyon (QGIS Development Team 2012).

Sampling, DNA extraction and amplification
We examined three populations of C. atlantica located in natural cypress woodlands in the western High Atlas including the Idrarane and Tanmmart populations from the most representative stands, the Aghbar Forest, which are separated from each other by ca. 2.6 km and one population represents the Allous stand, which is separated from the Aghbar Forest populations by ca. 12 km (Table 1 Genomic DNA was extracted from dried leaves using the CTAB protocol (Dumolin et al. 1995). Nine nuclear microsatellite markers (nSSRs) developed for C. sempervirens L. were tested (Sebastiani et al. 2005); however, only six were selected for genotyping in C. atlantica: cyp84, cyp101, cyp174, cyp257, cyp258 and cyp293, as they gave clearly interpretable and high-quality amplification products (Online Resource 2). PCR was performed according to the procedure described in Sebastiani et al. (2005), and amplification was conducted following the touchdown PCR protocol presented in Sękiewicz et al. (2018). PCR products were run on an ABI PRISM 3130 genetic analyser with an internal size standard, GeneScan™ 600 LIZ® (Thermo Fisher Scientific, Waltham, USA), and fragment sizes were scored using GeneMapper 4.0 software (Thermo Fisher Scientific, Waltham, USA). For classifying observed SSR allele sizes into representative discrete alleles, Raw2Gen software (IJ Chybicki, unpubl. data) was used.

Genetic diversity and differentiation
GENEPOP v 4.2 (Raymond and Rousset 1995) was used to evaluate the linkage disequilibrium between each pair of loci at each site and to test departures from Hardy-Weinberg equilibrium for each population. To identify loci putatively under selection, the F ST outlier test was assessed based on the multinomial-Dirichlet model using BayeScan v 2.1 (Foll and Gaggiotti 2008). Twenty pilot runs were performed at a length of 5000 iterations to estimate model parameters. The sample size was set to 5000 with a 20 thinning interval and an initial burn-in of 50,000 steps, resulting in a total length of 150,000 MCMC iterations. To reduce the occurrence of false positives, the prior odds for the neutral model were set at 100 (García-Verdugo et al. 2015).
Genetic diversity parameters, including the average number of alleles (A), effective number of alleles (A E ), number of private alleles (P A ), observed (H O ) and expected (H E ) heterozygosity were computed at the population and locus level using INEst 2.1 (Chybicki and Burczyk 2009). Allelic richness (A R ) was calculated using the rarefaction method implemented in FSTAT 2.9.3 to avoid the bias of different sample sizes (Goudet 2001). An individual inbreeding model (IIM) implemented in INEst was applied to estimate the frequency of null alleles (Null) and multilocus inbreeding coefficient (F IS Null), including 'null alleles' correction (Chybicki and Burczyk 2009). Estimation was conducted with 5 × 10 5 MCMC iterations, storing every 200th value with burn-in steps of 5 × 10 4 . The deviance information criterion (DIC) was used to compare the tested models to assess the factors affecting the homozygosity level.
BayesAss 3.0 was used to infer recent migration rates among the studied populations using the Bayesian approach (Wilson and Rannala 2003). The migration rate was estimated based on the proportion of individuals in each population that were assigned to other populations with high probability. Initially, runs with the default values for mixing parameters (allelic frequency, migration rate and inbreeding) were used. Subsequent runs with different delta values were tested to obtain the best value for acceptance. However, even for maximum values of mixing parameters, the acceptance rates were high (0.72 for inbreed coefficients, 0.68 for allele frequencies and 0.50 for migration rates). To ensure convergence, three independent runs were performed with MCMC iterations set at 10 7 , a burn-in of 10 6 and a sampling frequency set at 10 3 . To choose the optimal run with the lowest value of the Bayesian deviance, the script of Meirmans (2014) was applied in R. Finally, a model-based Bayesian clustering approach implemented in Structure 2.3.4 (Pritchard et al. 2000) was used to infer population substructuring. The non-spatial admixture model with correlated allele frequencies was applied. The analysis was conducted with ten independent runs for four groups (K) with a burn-in period of 10 5 and 2 × 10 5 MCMC iterations for each run. CLUMPAK (Kopelman et al. 2015) was used to align replicated runs and to estimate the optimal number of clusters according to the mean log posterior probability per cluster (lnP[K]) method (Pritchard et al. 2000).

Population demography
The contemporary effective population size (N e ) was estimated using a single-sample method based on a bias-corrected version of the linkage disequilibrium approach, assuming a random mating model (Waples and Do 2010) implemented in NeESTIMATOR v2 (Do et al. 2014). Following the recommendations of Waples and Do (2010), low-frequency alleles P Crit ≤ 0.02 were excluded from the analysis to minimize potential bias caused by rare alleles.
To infer past demographic fluctuations in population size, the M-ratio approach (Garza and Williamson 2001) implemented in INEst was performed to estimate recent bottlenecks. The simulated analysis was run under the two-phase model (TPM) model with recommended default parameter sets. The Wilcoxon signed-rank test was used to verify the significance of the deficiency in the M-ratio based on coalescent simulations consisting of 10 6 iterations.

Spatial genetic structure
The SGS was assessed by spatial autocorrelation analyses of kinship between individuals within each population using SPAGeDi 1.5 (Hardy and Vekemans 2002). The average multilocus pairwise kinship coefficients (F ij ) were estimated according to Nason's formula (Loiselle et al. 1995) and averaged within each of the analysed distance classes. The upper distances of the intervals for the Idrarane and Tanmmart populations were 5, 10,15,20,25,30,35,40,45,50,60 and 100 m, while those for Alouss were 5, 10, 15, 20, 25, 30 and 40 m. To define the distance classes, we considered the recommendation of Hardy and Vekemans (2002) assumes that the number of pairwise comparisons belonging to each interval (% partic) should be > 50% and the coefficient of variation of the number of times each individual (CV partic) set as ≤ 1. Average values of F ij were regressed on the natural logarithm of the spatial distance between individuals ln(d ij ) to obtain the regression slope, b F , which quantifies the extent of the SGS. To assess the significance of the SGS, the 95% confidence intervals of F ij for each distance class and b F were obtained from 2 × 10 4 permutations of the spatial position of the individuals within locations. To visualize the SGS, average F ij values were averaged over a set of distance intervals and plotted against the geographic distance. The intensity of the SGS was estimated according to the Sp statistic developed by Vekemans and Hardy (2004) as follows, is the slope of a log-linear regression between observed kinship and the distance between individuals, and F (1) is the mean F ij for individuals in the first distance class.

Species distribution modelling
The multicollinearity among the analysed bioclimatic and soil predictors was prevented by discarding eleven correlated variables, keeping nine variables as predictors for the modelling process (Table 2; Online Resource 3). The MaxEnt model showed high levels of predictive performance for the predictions of the past and current C. atlantica distribution, as expressed by the accuracy value (AUC) reaching 0.991. Similarly, a high accuracy score (AUC = 0.991) was obtained for the future theoretical range including all tested climatic scenarios, indicating that the model is highly reliable for predicting the potential species distribution in Morocco at this time. The results of the variables contributions used in distribution modelling are shown in Table 2. The most important factors limiting the current potential C. atlantica range were isothermality (Bio3; 30.2%), the minimum temperature of the coldest month (Bio6; 25.0%), the precipitation of the driest quarter (Bio17; 19.6%) and the type of soil (WRB; 13.5%). For future predictions, the examination of the relative contribution of the predictor indicated that these above factors (Table 2) were also identified as important factors limiting the future distribution of C. atlantica. The set of factors with the highest contribution used to the past distribution model was similar to other models. However, the major difference relates to isothermality (Bio3). It was defined as the second most important variable (35.4%), while the minimum temperature of the coldest month (Bio6; 35.6%) had the highest impact on the past distribution model for C. atlantica ( Table 2).
The present theoretical range of C. atlantica is much larger than the current distribution (Fig. 3). The theoretical range covered mainly the area of the western part of the High Atlas with high habitat suitability (65%). However, the current area of occurrence of the species consists of only ca. 49% of the predicted potential areas with medium and high habitat suitability (Table 3). In addition, prediction indicated that the Anti-Atlas in some regions seems to possess climatically suitable habitats for the species (65%), but currently no observations confirm its occurrence there (Fig. 3).
Modelling the species distribution during the late Holocene revealed a profound reduction in areas with suitable habitat for C. atlantica ( Fig. 4a; Table 3). However, very limited areas with low habitat suitability (35%; ca. 43 ha) were found in the Anti-Atlas and northwestern parts of the High Atlas ( Fig. 4a; Table 3). For future climatic scenarios, the projection under the intermediate RCP4.5 scenario was more optimistic and showed the potential range expansion of C. atlantica by 2070 (Fig. 5). The regions with high habitat suitability for species (> 65%) encompassed the western High Atlas Mts. and most of the Anti-Atlas, with very limited areas in the Middle Atlas (< 5%). However, a dramatic contraction of the potential range of the species by 2070 was projected under more pessimistic scenarios (RCP6.0 and RCP8.5) and altitudinal shifts of suitable areas reached an average of 103 m a.s.l. and 407 m a.s.l, respectively (Table 3). Specifically, except for very limited areas with low habitat suitability (< 5%) in the western High Atlas Mts., the species practically disappeared (Table 3, Fig.  4b, c).

Genetic diversity and differentiation
All analysed nSSR loci were selectively neutral, and no evidence of linkage disequilibrium between each pair of loci in each population was detected (Online Resource 2). In total, 145 alleles were observed in all screened individuals, and the number of alleles (A) ranged from 14 to 34 with an average of 24.167 alleles per locus. Null alleles (Null) were detected in all loci, and the frequency of this parameter ranged from 0.014 (cyp174) to 0.228 (cyp257), resulting in an average of 0.065 (Online Resource 2). The detected high frequency of null alleles in cyp257 resulted from its compound microsatellite structure; this type of microsatellite is prone to failed amplifications, producing null alleles.
Genetic diversity parameters did not significantly differ between the analysed study sites (Table 1). All populations were characterized by relatively high genetic variation,  However, the estimated F IS including 'null alleles' correction was low at each site (F IS = 0.046; F IS = 0.046 and F IS = 0.055). Generally, null alleles were detected in each population with an average frequency of 0.073. According to the DIC, inbreeding was indicated as a likely causative factor of the observed deficiency of heterozygotes in populations, except for Idrarane (Table 1). BayesAss analysis indicated asymmetric migration, with a much stronger inflow from Idrarane to Tanmmart and Alouss (Aghbar Forest) populations than in the opposite direction, which was minimal (0.6%) (Online Resource 4). The average proportion of migrants was 0.1652 and 0.239, indicating that 16.5% and 23.9% of individuals in these populations are migrants from the Idrarane population. A moderate level of bidirectional gene flow between Tanmmart and Alouss was noted, corresponding to 8.3% and 4.8%, respectively (Online Resource 4).
The clustering analysis did not reveal a substructure in the studied population. All individuals could be pooled into four gene pools with a similar proportion of memberships and a high gene admixture rate (Online Resource 5).

Population demography
The contemporary effective population size (N e ) was comparable among populations (Online Resource 6) and reached N e = 49.8 (CI 95% 82.4-220.8) in Idrarane, N e = 49.9 (CI 95% 291.6-infinite) in Tanmmart and N e = 52.7 (CI 95% 58.4-134.2) in the Alouss population. No evidence of a recent bottleneck was detected in the analysed populations (Online Resource 6).

Spatial genetic structure
Details on the statistics that characterize each distance classes are given in Online Resource 7. A significant SGS depicted by a negative regression slope (b F ) was found for all study sites with the highest positive value of the kinship coefficient (F ij ) in the first distance class (0-5 m) and decreasing continuously with increasing spatial distance, especially in the Alouss population (Table 4; Fig. 6). Positive inter-individual F ij values were significant at only short distances up to 10 m (p < 0.05) in both the Idrarane and Alouss populations, while in Tanmmart was detected in only the first distance class (0-5 m; p < 0.05). A negative significant kinship coefficient was observed at 480 m in Idrarane and at 360 m and 480 m in the Tanmmart population (p < 0.05; Fig. 6). A weak spatial pattern of genetic structure was reported for all analysed cypress stands, as indicated by the low value of Sp (Table 4). In both the Idrarane and Alouss populations, the intensity of the SGS was at very similar levels (Sp = 0.0135 and Sp = 0.0112), and the lowest value of Sp was detected in the Tanmmart population (Sp = 0.0088). Due to detected high level of null alleles in locus cyp257, the additional analyses excluding this locus were performed  (Table 4).

Species distribution patterns
The estimation of the theoretical present distribution of C. atlantica (Fig. 3) indicated a significantly larger suitable habitat than the confirmed species range in the Oued N'Fiss. The current range constitutes only 49% of the estimated range, considering the area with the highest and medium suitability (Table 3). However, the current species range overlaps the modelled areas with the highest suitability. Therefore, high confidence can be placed on the predictions, as was also evidenced by the high accuracy score of the model. Although potentially suitable locations for species also appeared in the Anti-Atlas, they had relatively low habitat suitability (< 35%).
Three main factors limit the current potential C. atlantica range: isothermality, minimum temperature of the coldest month and precipitation of the driest quarter (Table 2). Among them, the isothermality had the highest contribution, which indicates the high sensitivity of species to oscillations Fig. 6 Correlograms of spatial autocorrelation coefficients (F ij ) (solid lines) plotted against maximum distance classes, estimated for Cupressus atlantica populations. Dashed lines represent the 95% confidence intervals Table 4 Spatial genetic structure (SGS) parameters for the studied populations of Cupressus atlantica (SGS parameters estimated without locus cyp257 are in italics) of diurnal and annual temperature ranges. This finding suggests specialization in the seasonality of the temperature regimes with hot summers and chilly winters. The significant impact of the precipitation of the driest quarter reflects the species adaptation to semi-arid and locally sub-humid conditions. Currently, C. atlantica occupies regions where the precipitation of the driest quarter ranges from 2.1 to 11.7 mm and the mean annual precipitation does not exceed 250-400 mm, except for some regions with rainfall of 600-700 mm/year (El Wahidi 2004). The low contribution of the mean temperature of the driest quarter might be explained by the species preferences regarding the wide temperature range. The species occurs in the area with mean summer temperatures reaching 12 to 30°C (El Wahidi 2004). However, the potential distribution of C. atlantica is limited by the minimum temperature of the coldest month, meaning that the species is strictly related to areas with relatively cold winters. The temperature during this period mainly ranges from c. 8.9 to − 1.4°C and does not exceed − 10°C (El Wahidi 2004). Considering the warming trend, warmer winters may affect the species distribution by allowing broadleaved species, primarily oaks, to shift upwards in altitude. According to climatic projections for oak species distribution, habitat suitability is likely to increase, and upward migration is predicted (Al-Qaddi et al. 2017;López-Tirado and Hidalgo 2018). This would be a limiting factor due to competition and would affect cypress woodlands, for which altitudinal shifts are also expected under future climate change. Apart from climatic variables, the type of soil positively affected the current species theoretical occurrence, especially in the Anti-Atlas. When soil type was excluded from the model (data not shown), the western High Atlas seems to be the only location that may highly support the occurrence of C. atlantica; in the Anti-Atlas, low suitability (< 5%) was reached in some areas. The most suitable habitat for the species included areas with leptosol, that is, the most common soil type in the mountainous regions. Nevertheless, the fact that C. atlantica currently does not occupy all of its potential range suggests that non-climatic factors are involved in species distribution. In our opinion, the most important factor would be strong human activity. In the mountains of Maghreb, the conifers have been subjected to strong human pressure for several millennia, which has resulted in the overexploitation of many of these species, which are now under considerable conservation concern (Barbero et al. 1990;Esteban et al. 2010;Taib et al. 2020). Griffiths (1998) reported a ca. 73% contraction of the C. atlantica range over a period of 36 years. This significant reduction in a relatively short time is a consequence of the strong anthropogenic pressure, mostly illegal logging and branch cutting (El Wahidi 2004;Arjouni et al. 2011). In 1998, Griffiths (1998 estimated that 84% of cypress trees were damaged by humans in several stands in the Oued N'Fiss, which was also confirmed in a more recent survey   (Fig. 1b-d).
There are no available most recent data about population abundance or demographic trends. However, the factors that contributed to the decline of the Moroccan cypress in the recent past still operate. For example, the maintenance of small-scale agriculture near human settlements, livestock and traditional pastoral activities are still frequently practised in the Maghreb (El Wahidi 2004;Hammi et al. 2010;Sękiewicz et al. 2014;Taib et al. 2020). Overgrazing negatively affects the natural regeneration of the species. Natural regeneration is also limited due to the high degree of soil erosion and physical dormancy caused by a hard seed coat (Sfairi et al. 2012). Additionally, in most of the populations, senile individuals prevail, which probably also decreases the species reproductive abilities. For example, most trees in the Aghbar Forest are 250-500 years old (Bechir et al. 2004). Some archaeological evidence (Ruas et al. 2011) may suggest a wider range of C. atlantica in the past, which additionally supports our hypothesis related to human-induced changes in species occurrence. In the Anti-Atlas, the findings of wood charcoal identified as Cupressus/Juniperus dated to ca. the eleventh-twelfth century found at the Iĝîlîz site may suggest local species presence (Ruas et al. 2011). According to modelling data, it seems that climatic conditions for C. atlantica during the late Holocene (4.2-0.3 ka) were less suitable than the current, and the species may survive in some refugial areas located in the Anti-Atlas and northwestern parts of the High Atlas, supporting the hypothesis that C. atlantica may have occurred in the Anti-Atlas in the past. However, verification of the hypothesis of a wider range of species is rather difficult because of the sparse information about Cupressus fossils due to their poor preservation in arid and semi-arid regions and uncertainties concerning taxonomic identity (Willis and McElwain 2002).
The future climate predictions for Maghreb highlight that the region will be exposed to intensified aridification, which will entail the risk of reduction/loss of unique plant assemblages (Rankou 2017). Considering the demographic trends of C. atlantica populations (Griffiths 1998;Sękiewicz et al. 2014), the reduction in species range seems to be unavoidable. However, future modelling of species distribution revealed a puzzling pattern of the species response to climate change. The most pessimistic future projections of the C. atlantica distribution under RCP6.0 and RCP8.5 (Fig. 4b, c) assume the serious retraction of suitable habitats for the species. By the end of the twenty-first century, the projected reduction would cover ca. 65% of the current theoretical range. However, the reduction mostly affects sites with the highest suitability, which mirrors progressive deterioration of climatic conditions for C. atlantica in the High Atlas (Fig. 4b, c, Table 3). The contraction would be especially severe in the Oued N'Fiss, probably reflecting a marked aridity trend in the region (Schilling et al. 2012). The increase in the temperature regimes, especially isothermality and minimum temperature of the coldest month, and the reduction in precipitation of the driest quarter, likely by 15% compared to the current conditions, had a primary influence on limiting the species' distribution (Online Resource 1 and 8). The model shows the possible persistence of the species in the westernmost part of the High Atlas beyond its core range but with very low suitability (Fig. 4b, c). Recently, such a pessimistic scenario has also been reported for other endemic trees in North Africa, such as Cedrus atlantica, Juniperus thurifera subsp. A definitely more optimistic projection for C. atlantica was obtained under RCP4.5, indicating the extension of the potentially favourable areas for the species (Fig. 5). Apart from the theoretically suitable areas (suitability > 65%) that overlapped with the core of the current species' range, the optimal conditions for the species appeared beyond the Oued N'Fiss. The mountainous area between Imi n'Tanout and Tizi-n-Test, near Argana, appeared likely to be suitable, as well as the areas situated in the Anti-Atlas and areas located southeast of the Oued N'Fiss among Askaoun, Anezal and Agouim (Fig. 5). However, even if suitable conditions would appear for the Moroccan cypress beyond its current range, the colonization ability of those new habitats might be restricted by dispersal limitations due to a lack of habitat connectivity related to topographic complexity, competition and low regeneration ability (El Wahidi 2004;Franklin 2010;Sękiewicz et al. 2018).

Genetic diversity and spatial structure
Given the narrow distribution of C. atlantica and notable reductions in its population size, the low level of genetic diversity and bottleneck effect was anticipated (Jump and Peñuelas 2006). Nevertheless, the results obtained here underline the persistence of high genetic diversity in C. atlantica (H E = 0.894) and provide no evidence of a recent bottleneck. This result is in line with previous genetic studies that reported a considerably high level of diversity for C. atlantica (Bechir et al. 2004;Sękiewicz et al. 2018;Bagnoli et al. 2020) and other narrow endemic conifers (Eliades et al. 2011;Cobo-Simón et al. 2020;Dering et al. 2014). In fact, this pattern of diversity observed for narrow endemic plants from the Mediterranean is not an exception. According to the recent review, half of the studied endemic species have moderate to high genetic diversity (Médail and Baumel 2018). In all analysed populations, a moderate level of homozygosity excess was observed, and in Tanmmart and Alouss, it likely resulted mostly from inbreeding. However, there is justified concern for the loss of genetic variability because in all analysed populations, a similarly low effective population size was detected (N e < 53). According to recent recommendations of the 100/1000 rule as the minimum viable population size (Frankham et al. 2014), the estimated N e for the C. atlantica population is too low to avoid the consequence of random genetic drift and thus prevent genetic loss.
The high level of diversity as well as lack of evidence for the bottleneck in C. atlantica populations could be explained by the relatively recent fragmentation, which has not yet been reflected in the current genetic structure of this extremely long-lived species (Lowe et al. 2015;Sękiewicz et al. 2018). Other explanations consider large extensive gene exchange among populations, which generally favours high variability in tree populations (Petit and Hampe 2006). Moreover, the high migration rate may result in false negative bottlenecks because gene flow can mitigate the genetic effects associated with small population sizes by replenishing alleles lost through genetic drift (Jump and Peñuelas 2006). In the studied populations, a significant rate of migration and high genetic admixture were inferred. Interestingly, a higher rate of gene flow and its unidirectional manner was revealed between the spatially separated Alouss and Aghbar populations (ca. 10.5 km and 12.5 km) than between the closely located populations in the Aghbar Forest (Idrarane to Tanmmart; ca. 2.6 km; Fig. 2). Cupressus atlantica forms an open woodland community with low-density populations that can facilitate effective pollen and seed dispersal. The species is wind pollinated and primary anemochoric, and this vector of dispersion is potentially very effective in gene dispersal in trees. Longdistance dispersal generally requires pollen movement to higher altitudes. However, the topographic complexity of the habitats, on the one hand, and the species occurrence at low altitudes in valley bottoms, on the other hand, may interrupt the effectiveness of spread and the model of its directionality (Di-Giovanni et al. 1996;Szczepanek et al. 2017).
Generally, fragmentation increases SGS due to nonrandom mating, low population density and potential aggregation of reproductive individuals, but the risk of its negative outcomes depends particularly on dispersal capabilities (Vekemans and Hardy 2004;De-Lucas et al. 2009). In all analysed populations, a significant but weak SGS was detected (Table 4, Fig. 6), indicating that current gene flow is partially restricted, which has led to the kinship structure. In all populations, the kinship coefficients among pairs of individuals decreased with the increasing logarithm of geographical distance between individuals. However, in Idrarane and Tanmmart, this decrease was not uniformly linear over the whole distance range (Fig. 6). A significantly positive SGS was detected up to 90 m in the Idrarane and Alouss populations and up to 40 m in the Tanmmart (Fig. 6). Interestingly, in populations from the Aghbar Forest after being negative for several distance classes, the kinship coefficients reached positive values between individuals that were over 500 m apart, suggesting non-random positive spatial aggregation of related genotypes at a larger scale. In the Aghbar populations, a significant negative SGS was observed over distances reaching 360 m and 480 m (Fig. 6). The low level of relatedness in these distance ranges suggests the large effective dispersal distance in this species. The weak spatial patterns of genetic relatedness detected in this study confirm the general assumption of a weak SGS in wind-pollinated and wind-dispersed tree species with an outcrossing mating system (Vekemans and Hardy 2004;Dering et al. 2015). The SGS was also lower than those reported for other endemic conifers with highly fragmented populations (Vekemans and Hardy 2004;Eliades et al. 2018). Thus, our results suggest that fragmentation has not yet brought about fully harmful outcomes for the spatial aspects of genetic diversity in C. atlantica. Currently, the species occurs in low-density populations, which generally leads to a reduced overlapping of seed shadows, which increases the probability of kinship within individuals and hence may contribute to strengthened SGS patterns in the future (Vekemans and Hardy 2004).

Conservation implication
Considering the projection of the future changes of the C. atlantica distribution and the demographic trends of the populations (Griffiths 1998;Sękiewicz et al. 2014), there are serious challenges to the sustainability of the species. Primary passive protection is not a sufficient strategy to prevent species loss, particularly in the face of climate change, habitat loss and human pressure. For that reason, there is a pressing need for the development of integrated conservation approaches, including in situ and ex situ measures. This approach should be oriented to the protection of both the species and its habitats and effective conservation of the species genetic resources in response to detected reduced habitat suitability (Koskela et al. 2013;Volis 2019).
The results obtained here and previously (Sękiewicz et al. 2018;Bagnoli et al. 2020) underline the persistence of high genetic diversity of C. atlantica, which is a very positive factor in terms of long-term species performance. However, there is a justified concern for the loss of variability over the shortterm related to external factors. Ongoing habitat fragmentation and overexploitation reduce the effective population size and increase the SGS, exposing populations to genetic stochasticity processes (Jump and Peñuelas 2006). The additional threat for C. atlantica is hybridization with non-native species, such as C. sempervirens and C. arizonica, which were introduced in close proximity to the natural sites of the species in the Oued N'Fiss, for instance, to the Aghbar Forest or Amizmiz sites (Bechir et al. 2004;Maerki and Lamant 2014). These factors pose a risk for C. atlantica genetic integrity. The loss of genetic diversity may become a limiting factor that may reduce the long-term stability and adaptability of populations to environmental changes, increasing the risks of local extinction. Therefore, conservation strategies should focus on monitoring temporal changes in genetic variability and structure due to demographic threats and promoting genetic diversity and gene flow in populations to maintain the evolutionary potential of the species (Koskela et al. 2013). Management should prioritize the translocation between populations in the frame of the concept of assisted migration/gene flow or genetic rescue. These conservation strategies are feasible for increasing the population density and size in cypress stands, which could reduce breeding among relatives and boost diversity. The elimination of plantations of non-native cypress species would also be crucial. According to pessimistic future scenarios, assuming the serious retraction of highly suitable habitats for C. atlantica, populations occurring in the main core of the species range in the Oued N'Fiss should be included as a priority in situ conservation measures, and regarded to require ex situ conservation methods.
Most of the areas today occupied by C. atlantica woodlands were predicted to match those indicated under the optimistic RPC4.5 scenario, which involved future species expansion. Consequently, at first, conservation efforts should be oriented towards promoting the in situ persistence of this species (Koskela et al. 2013), particularly by reducing intensive human pressure, especially uncontrolled pastoral activity, which affects natural regeneration of the species. In terms of this threat, management should aim at monitoring and protecting natural regeneration. Such activities were effective for other endemic tree Abies nebrodensis, which is currently represented by only 30 mature individuals in the Madonie Mountains (Sicily). As a result of in situ conservation, an increase in natural regeneration by approx. 94% was observed (Raimondo and Schicchi 2005). Additionally, the establishment of an in situ collection of the species would be highly required. Several breeding and planting programmes for C. atlantica have been undertaken; however, initial observations have shown that success rates remain low (El Alaoui El Fels et al. 2017). For that reason, some efforts have been made to improve the natural regeneration and germination of C. atlantica (Sfairi et al. 2012;El Alaoui El Fels et al. 2017). For example, the use of mycorrhizal inoculation of seedlings and nurse plants of C. atlantica has been considered for the reforestation of degraded areas (Birouk et al. 1996;Ouahmane et al. 2007;Sfairi et al. 2012;Arjouni et al. 2013;Hafidi et al. 2013). Recently, restoration guidelines have been developed to restore heavily degraded forest ecosystems in the Oued Ourika (High Atlas). The authors recommend using native plants, among others, C. atlantica during the reforestations. These activities were developed as part of the research project addressing the management of water resources and adaptation to climate change in Morocco (Etienne 2017). Accordingly, our SDM results support the possibility of the species growing in the upper the Oued Ourika ( Fig. 5; habitat suitability 35%).
In addition, the future extension of potentially suitable habitats for C. atlantica beyond the Oued N'Fiss invites the inclusion of assisted migration in conservation management. In particular, the introduction of the species into the sites in the Anti-Atlas and in areas southeast of the main core of the species range (Fig. 5) could be considered. However, effective conservation planning for C. atlantica should not neglect the knowledge of local adaptation, which has the potential to mitigate the effects of climate deterioration on the long-term persistence and range shifts of the populations (Razgour et al. 2019).
Despite the lack of a negative effect of fragmentation on the current level of diversity in C. atlantica, the following questions remain: which populations should serve as main seed sources and what sampling strategy should be applied to capture the widest spectrum of diversity? The proper ex situ conservation strategy, including assisted migration, breeding programmes and gene banks, and the establishment of in situ living collections and habitat restoration should minimize the recognized genetic and demographic threats associated with the creation of non-natural populations, including avoidance of inbreeding and a loss of genetic diversity due to founder effects (Ottewell et al. 2016). Here, the SGS characterized at a fine scale can be effectively used for these conservation measures. We suggest that seed sampling strategies should include collection from seed trees over 100 m apart to avoid collecting samples from relatives. Moreover, according to the recommendations of Hoban and Schlarbaum (2014), the range-wide population structure should also be considered to capture more diversity/ alleles, especially rare alleles. The authors suggest that when a population structure exists, dispersed sampling across the range is optimal to capture more alleles than constrained sampling in one region because of spatially restricted alleles (Hoban and Schlarbaum 2014;Hoban et al. 2018). Including the rangeedge populations in the seed sampling strategy is important because these populations often maintain unique alleles, which promote divergence in comparison to the core of the range, where populations share common alleles (Hampe and Petit 2005). According to Sękiewicz et al. (2018), a subtle population genetic structure was detected for C. atlantica despite the restricted species distribution range (see Fig. 1A in Sękiewicz et al. 2018). Currently, our investigation did not reveal a geographically dependent genetic structure due to the detection of intensive gene admixture among populations (Online Resource 5); however, these results may be due to the low number of populations studied. Considering the previous results (Sękiewicz et al. 2018), populations from the Aghbar Forest, Alouss and Tazalt stands (near Alouss) likely represent unique genetic assemblages, where private alleles were noted; therefore, they should be considered in seed collections.
Currently, C. atlantica is struggling with threats caused by climate change, habitat loss and strong human pressure. The predictions for the species are not very optimistic, as was indicated by SDM modelling. Without a broad-scale conservation strategy focusing on in situ and ex situ conservation, there is a serious concern about the survival of the species in the western High Atlas.