Anthropogenic deforestation and climate dryness as drivers of demographic decline and genetic erosion in the southernmost European fir forests

A better understanding of long-term effects of climate and historical anthropogenic changes is needed to define effective conservation measures of endangered forest inhabiting managed landscapes. Diversification and distribution of Mediterranean firs are attributed to the global climate change during the Miocene and Quaternary as well as to the effects of human-driven deforestation. We evaluated the impact of climate change and anthropogenic activities in shaping the genetic diversity and structure of Abies pinsapo Boiss. (Pinaceae), a relict fir endemic from SW Spain. We genotyped a total of 440 individuals from 44 populations by using two different molecular markers (cpSSRs and nSSRs). Overall, low genetic structure was found; however, incipient differentiation appeared within mountain ranges. Analyses suggest that the effects of isolation by distance and lithological or topographical diversity were not enough to structure the populations of the different mountain ranges. The combined role of genetic drift in the small populations and the anthropogenic action associated with forest management has shaped the current genetic pattern of this fir species in the study area. Demographic inference analyses pointed to a very recent synchronic divergence (eleventh–sixteenth century) of the ancestral A. pinsapo population into its current scattered distribution range. Although population bottlenecks were supported by several analyses, the conservation of this endangered species seems not to be limited by lacking genetic diversity, while threats of current climate change and habitat loss must be regarded.


Introduction
The current genetic diversity and structure of populations across their distribution range have been strongly influenced by climate changes over the course of the Quaternary ice age (Hewitt 2000). During interglacial periods, distribution ranges of cold-adapted species, such as boreal and alpine trees, shrank massively, sometimes causing vicariant distributions where small remnant populations became genetically isolated and severely bottlenecked (Hewitt 2000). These species expanded their refugium ranges during glacial periods reducing the genomic variability at the leading edge due to founder effects, although some genetic admixture between refugium populations increased diversity in newly colonized areas (Petit et al. 2003). Therefore, repeated cycles of range expansion and contraction result in a heterogeneous mosaic of genome diversity across extant populations that determine Communicated by Oliver Gailing.
Juan Luis García-Castaño and Salvador Talavera: Posthumous author. the evolutionary potential of relict tree species (Hampe and Petit 2005;Hampe and Jump 2011). Identifying those areas of reduced and increased diversity and evaluating the impact of historical demographic processes that underlie them are of utmost importance if we aim to guarantee their persistence in an increasingly warming world.
Beside climate changes, human-driven deforestation also impairs evolutionary processes that shape the genetic diversity of populations by extirpating trees (and their genes) while slowing down gene flow as fragmentation increases (Di Marco and Santini 2015;Frankham et al. 2017). Overall, deforestation exacerbates the genetic erosion of tree populations already depauperated by several cycles of range expansion and contraction (Dale 1997). However, we lack reliable, empirical estimates on the relative impact of climate driven quaternary range changes vs. human-driven deforestation in shaping the observed patterns of genetic variation in relict alpine trees highly threatened by climate warming. This knowledge gap can be addressed by applying chloroplast and nuclear genetic markers to depict spatial genetic patterns and to infer their historical demographic trends. Chloroplast markers (cpSSRs) are paternally inherited in most conifers, haploid, and show slow-paced mutation rate which, overall, make them little variable, while nuclear markers (nSSRs) are biparentally inherited, diploid, and show a fast-paced mutation rate, which make them highly variable (Avise 2004).
Here we evaluate the impact of climate change and anthropic activities in underlying the genetic diversity and structure of the Tertiary relict fir Abies pinsapo Boiss. (Pinaceae). Colonization and diversification patterns of firs in the Mediterranean are shaped by global climate change during the Miocene, while Quaternary climatic cycles would prompt secondary contact zones and posterior isolation of species (Balao et al. 2020). Further, human-driven deforestation threatens Mediterranean firs with a geographic scattered distribution, many of them endemic to few mountain ranges. Currently, this species occupies three mountain ranges at the southwestern region of the Iberian Peninsula: Sierra de Grazalema (ca. 400 ha), Sierra de las Nieves (ca. 1000 ha), and Sierra Bermeja (ca. 50 ha). Scattered stands and isolated trees, assumed to represent the remains of former larger populations, are rather common throughout the current A. pinsapo range. Modeling studies, supported by pollen fossil records, evidence that the distribution range of A. pinsapo has massively shrunk during the Holocene (Alba-Sánchez et al. 2010 due to a combination of climate warming and an intense anthropic activity that includes severe deforestation driven by a flourished naval industry during the sixteenth century, an indiscriminate logging during the nineteenth century onwards for timber, and uncontrolled livestock that hampered natural regeneration during the twentieth century (Becerra Parra 2006). Previous genetic studies on a few populations suggested evidence for bottlenecks and genetic drift, as well as low among-population genetic differentiation (Sánchez-Robles et al. 2014b;Cobo-Simón et al. 2020). Based on an exhaustive sampling across 44 populations that span the complete distribution range and by combining cpSSR and nSSR markers, we aim to: (1) quantify the spatial distribution of the genetic diversity. Based on known distribution range shifts and an intense deforestation documented in the area Alba-Sánchez et al. 2019;Cobo-Simón et al. 2020), we expect low levels of genetic diversity and increased values of genetic structure particularly when dealing with small and geographically isolated populations, presumably subject to severe bottlenecks in the past; (2) test whether geographic features explain observed genetic patterns. To that end we apply landscape genetic models to test for isolation by distance (i.e., an increased genetic distance as the geographic distance increases caused by reduced levels of gene flow, Wright 1943). Extant populations are distributed across mountain ranges and, therefore, we expect that, in addition to geographic distance, orographic features and effective pollen dispersal distance would limit gene flow (Sánchez-Robles et al. 2014a). To support this hypothesis, we further identify genetic clusters that typically arise when populations become genetically isolated and we depict the location of geographic barriers to gene flow (Balkenhol et al. 2016);and (3) seek for, quantify and date past bottlenecks that might have contributed to shape current genetic patterns. On this account, we use cpSSRs and nSSRs to date demographic bottlenecks inferred by applying Approximate Bayesian Computation models (Cornuet et al. 2018). Bottlenecks caused by climate changes are expected to happen as the southern distribution range limit retreated after the last glaciation ca. 10 kya. Finally, we discuss the evolutionary consequences of past bottlenecks for relict and cold adapted trees in the light of increasing climate warming.

cpSSR and nSSR molecular analyses
The genotyping error rate for cpSSRs was estimated by genotyping a replicate taken at random from each population (44 samples in total; 10% of the total). Repetitions were carried out starting from the same DNA extraction. Allele combinations across the four cpSSR loci were treated as a haplotype, considering different haplotypes as a unique size variant combination of each cpSSR locus. The haplotype variation present in the different populations, as well as in the three mountain ranges as a whole, was explored by calculating the number of haplotypes (Nh) and the number of private haplotypes (Nh priv ) using the software Arlequin v3.5 (Excoffier et al. 2005) and SPAGeDi v1.4 (Hardy and Vekemans 2002).
Regarding the nuclear markers, several tests were carried out to detect the different types of errors present in the nSSR data: (i) Micro-Checker v2.2.3 (van Oosterhout et al. 2004) was used to estimate the rate of null alleles (alleles that are present but that cannot be genotyped); (ii) Pedant v1.0 (Johnson and Haydon 2007) was used to estimate the rate of dropout alleles (e 1 ; homozygous individuals genotyped as heterozygous due to random failure to amplify an allele) as well as false alleles (e 2 ; poorly genotyped alleles) using 10,000 iterations; and (iii) Fstat v2.9.3.2 (Goudet 2001) was used to test linkage disequilibrium between pairs of loci applying the Bonferroni correction for multiple  comparisons. Additionally, the different mountain ranges, as well as the different populations within them, were genetically characterized by estimating the mean number of alleles and effective alleles (N e ) corrected for the sample size (Nielsen et al. 2003), the genetic diversity based on expected heterozygosity (H e ; Nei 1978), and the inbreeding index (F IS ). Such diversity estimates were obtained with SPAGeDi and Fstat. ANOVAs were used to analyze mountain range differences, assessed as a fixed effect, for N e , H e , and F IS (the latter ones, previously arcsine-square root transformed as their values, after rescaling for F IS , are 0-1 limited,); Tukey's HSD post-hoc tests (α = 0.05) were applied to test for pairwise differences. These analyses were performed with the JMP v6 software (SAS Institute Inc.). Furthermore, Spearman correlation tests were applied for both markers to detect the possible relationship between altitude and genetic diversity (based on the expected heterozygosity: H e ). In a similar way, we proceeded to detect the relationship between H e and the degree of isolation of a population (measured as the average distance to the other populations within the same mountain range), which was not correlated with the altitude (Spearman's ρ = 0.026, P = 0.865, n = 44). These analyses were also performed with the JMP v6 software (SAS Institute Inc.).

Genetic structure
To evaluate different hypotheses about the partition of the genetic variance, analyses of molecular variance were performed (AMOVAs; Arlequin v3.5). Groupings were carried out to represent the different mountain ranges and the different possible pairwise combinations between them. Additionally, based on the presence of peridotite, a grouping represented the different lithology outcrop areas. Variance components were tested for significance using an exact nonparametric test with 10,000 permutations.
For both markers, the global genetic structure was analyzed by Bayesian assignment analysis using the Structure v2.3.3 (Pritchard et al. 2000) software. Correlated allelic frequencies were considered (since due to migration or existence of common ancestors it was expected that there would be similar allele frequencies among populations) and 10 independent repetitions for each K were employed (K = 1-9). The number of iterations was 1 × 10 6 , with a burn-in period of 150,000. To determine the optimal number of clusters (K), we checked the Structure results with the web-based Structure Harvester (Earl and vonHoldt 2012), which implements the Evanno et al. (2005) method using the adhoc statistics ΔK.
Furthermore, a phylogeographic test was carried out in order to study the genetic differentiation under different evolutionary models (Pons and Petit 1996). For this test, both markers were used independently (cpSSRs and nSSRs). The test compared observed randomly permuted N ST (pN ST = G ST ) values. G ST estimates are based on allele frequencies and follow an infinite allele mutational model (IAM) for microsatellites, whereas N ST takes into account the distance between different haplotypes, following a stepwise mutational model (SMM). Significantly higher N ST than G ST values indicate the presence of phylogeographic structure, i.e., different alleles are more related at a within-than at an among-population level (Lynch and Crease 1990;Pons and Petit 1996). The test was performed using SPAGeDi through 20,000 permutations. N ST was computed using a matrix of genetic distances based on the sum of the squared differences (R ST ; Hardy et al. 2003). In order to detect the existence of phylogeographic signals at different scales, the test was applied independently for the three mountain ranges and for pairs of mountains (i.e., pairwise approach).
In order to explore the relationships between the different chlorotypes found, a haplotype network (Minimum Spanning Network) was constructed in ARLEQUIN using a matrix of F ST distances based on the number of the different alleles found.

Isolation by distance and barriers to gene flow
The possible existence of isolation by distance (IBD) was analyzed by a Mantel test using the R package ecodist v2.0.7 (Goslee and Urban 2007). The Nei's genetic distance matrix (Nei 1972) was compared with the geographic distance matrix. The test was carried out using Spearman's nonparametric correlation with 10,000 permutations.
Likewise, genetic barriers (i.e., geographic areas with pronounced genetic discontinuity among populations) were explored with the Barrier v2.2 (Manni et al. 2004) software. The boundaries among the mountain ranges associated with the highest rates of genetic differentiation were detected based on a geographic distance matrix, and using the Delaunay triangulation and the Voronoi polygons. To evaluate barrier robustness, the population pairwise genetic distance matrix based on Jost's D index (Jost 2008) was bootstrapped 1,000 times with the R package boot v1.2-43 (Canty and Ripley 2019) and, subsequently, run in Barrier.
Lastly, to explore and characterize the complex topology resulting from ancient and recent genetic interactions, we used a model based on Graph Theory (Dyer and Nason 2004). The genetic covariation between the different populations was analyzed for each marker by using the R package popgraph v1.2 (Dyer 2013). For each marker, the relationships between populations (nodes) were taken into account at a significance level of 0.050. Subsequently, we assessed the congruence among the connectivity patterns (Structural Congruence), as well as among the pairwise distances (Distance Congruence) of the graphs obtained for each marker (Dyer et al. 2010).

Demographic analysis
Finally, to trace the demographic history of A. pinsapo, the approximate Bayesian computation (ABC) analysis procedure (Beaumont et al. 2002) implemented in DIYABC v2.1 (Cornuet et al. 2014) was performed based on both cpSSR and nSSR datasets. The different cpSSR and nSSR loci were treated separately and the default values of the priors were used for all parameters. The summary statistics used were mean gene diversity, F ST , and the classification index. An estimate of 20 years was used for the generation time in order to calculate the number of generations (Arista and Talavera 1994;Cobo-Simón et al. 2020). Four evolutionary scenarios were tested with DIYABC, based on the results obtained in the different analyses. Briefly, we tested a synchronic and asynchronic divergence of the three mountain ranges (scenarios 1 and 2) and we also tested the presence of bottlenecks just for Sª Bermeja (scenario 3) and for the three ranges (scenario 4). For more details, see Supplementary Material, Appendix B and Table S4.
Additionally, Migraine software ) was used to detect changes in population size (i.e., bottlenecks or expansion) for the three mountain ranges using maximum likelihood estimation (MLE). For each range and genetic marker (i.e., cpSSRs and nSSRs), we ran the software under the OnePopVarSize demographic model (a single population model with a single continuous past variation in population size) in order to estimate the bounds of the demographic parameters: current effective population size (2Nµ current), scaled duration of population size change (Dg/2N), ancestral effective population size (2Nµ ancestral), and the current population/ancestral population ratio (N-ratio). A generalized stepwise mutation model (GSM) was used for nSSR, whereas a strict stepwise mutation model was set for cpSSR. Preliminary parameter bounds were estimated with 4 iterations, 500 point estimates and 2000 runs per point. Final inferences were performed setting 8 iterations with 500 point estimates and 2000 runs per point. We evaluated significant changes in population size checking if 95% confidence intervals (CIs) excluded N-ratio = 1.0.

Genetic diversity and spatial distribution of the Spanish fir
The overall total error rate for the cpSSR markers was 1.14%. For the nSSR ones, we found null alleles for most of the loci used when all populations were considered as a whole (Table S2) We recorded 30 different cpSSR haplotypes and 86 nSSR alleles (Table S2). The number of effective alleles (N e ) showed no significant difference among mountain ranges for cpSSRs (F 2,41 = 0.85, P = 0.433) but did show a significant one for nSSRs (F 2,41 = 3.31, P = 0.047), between Sª Bermeja (N e = 2.11) and Sª de las Nieves (N e = 2.70). Genetic diversity showed no significant difference among mountain ranges both for cpSSRs (F 2,41 = 1.65, P = 0.205) and nSSRs (F 2,41 = 3.15, P = 0.053). Furthermore, we found significant differences for the inbreeding coefficient (F IS ; F 2,41 = 4.66, P = 0.015) between Sª Bermeja (F IS = − 0.16) and the other two mountain ranges (Sª de Grazalema: F IS = 0.01; Sª de las Nieves: F IS = − 0.01). For more details, see Table 2.
The global analysis of molecular variance (AMOVA ;  Table S3) showed that most of the genetic diversity was found within populations (cpSSRs: 80.4%; nSSRs: 90.1%). Mountain ranges played a significant role in structuring the Table 2 Global and within mountain ranges genetic diversity obtained for cpSSRs and nSSRs n: number of individuals sampled; N: mean number of haplotypes; N e : number of effective alleles (Nielsen et al. 2003); Nh: number of haplotypes; Nh priv : number of private haplotypes; H e : genetic diversity corrected for the size of the samples (Nei 1978 Although the clustering analysis (Figs. S1 and S2) showed that sampled individuals could be grouped according to K = 4 (cpSSRs) and K = 2 (nSSRs) following the ΔK approach, these clusters did not fit any well-defined geographic pattern, i.e., clusters contained individuals sampled across all populations and mountain ranges. For nSSRs, the likelihood increased as K increased, whereas the mean likelihood reached a maximum in K = 4-5 for the cpSSR analysis.
The combination of alleles from the different cpSSR loci resulted in a total of 30 different haplotypes. Fig. S3 shows the relationship among the different haplotypes found, not appreciating defined patterns of relationship among them across the different mountain ranges.
The Monmonier's maximum-difference, conducted on the two markers, showed a genetic border that separated Sª Bermeja from the other two mountain ranges, Sª Grazalema and Sª de las Nieves (Fig. 2a).
The phylogeographic test resolved that the most probable mutation model that occurs for both markers, considering all populations together, was the infinite allele model (IAM). However, the significance for both markers varied considerably, being this model more appropriate for the cpSSR markers (N ST = 0.143 vs. pN ST = 0.147; P = 0.493) than for the nSSR ones (N ST = 0.116 vs. pN ST = 0.061; P = 0.051). When the test was conducted for pairs of mountain ranges, the presence of genetic structure for the cpSSR markers was only observed when Sª de las Nieves and Sª Bermeja were compared (N ST = 0.197 vs. pN ST = 0.092; P = 0.036), indicating that, for these markers, Sª Bermeja is more differentiated from Sª de las Nieves than from Sª de Grazalema.

Environmental factors that shape the distribution of the genetic diversity across the landscape
We observed a positive and significant correlation between geographic distance and genetic distance between populations, suggesting an overall isolation by distance pattern (cpSSRs: Spearman's ρ = 0.286; nSSRs: Spearman's ρ = 0.267; both P < 0.001). On the one hand, altitude and genetic diversity (H e ) were not correlated (cpSSRs: Spearman's ρ = 0.036, P = 0.815; nSSRs: Spearman's ρ = 0.096, P = 0.534). On the other hand, whereas cpSSRs revealed that the level of geographic isolation of a population and its genetic diversity (H e ) were negatively correlated The genetic covariance networks obtained at significance level of 0.050 for both markers showed congruence in their structure (Structural Congruence: CDF = 0.024) and distances (Distance Congruence: r = 0.076; P = 0.020). The combined network of both markers (Fig. 2b) showed a greater connectivity (number of connections between populations) between Sª Grazalema and Sª de las Nieves than between these mountain ranges and Sª Bermeja. Four population connections were observed between Sª Bermeja and Sª de las Nieves, which were weaker than the only connection found between Sª Bermeja and Sª Grazalema (which still remained at the P = 0.015 significance level). Only the connections between populations present within each mountain range, and between populations of Sª de las Nieves and Sª Grazalema, remained when the significance level increased to a higher value (P = 0.010), remaining Sª Bermeja isolated from the other two mountain ranges.
The results of the maximum likelihood estimation under the single population model with a single continuous past variation in population size were mostly congruent with ABC analyses. Significant population contractions between the ancestral and the current A. pinsapo populations were detected for the three mountain ranges for cpSSR genetic markers. Additionally, significant population bottlenecks were detected for Sª Bermeja and Sª de Grazalema for nSSR genetic markers, whereas N-ratio was not significantly different from 1.0 in Sª de las Nieves (see Table S7).

Contrasted genetic diversity patterns among markers
The global pattern observed in our study supports a high population genetic diversity for plastidial and nuclear markers, which contrasted with the long-term isolation and the relict character of this species (Hampe and Petit 2005;Hampe and Jump 2011). Similar genetic diversity values have been previously shown in A. pinsapo (Terrab et al. 2007;Dering et al. 2014;Sánchez-Robles et al. 2014b;Cobo-Simón et al. 2020), as well as in other relict Mediterranean firs, such as those from north Africa (A. marocana, A. tazaotana, and A. numidica;Terrab et al. 2007;Sánchez-Robles et al. 2014b) and the eastern Mediterranean fir group (Awad et al. 2014;Hrivnák et al. 2017). These studies found relatively high within-population genetic diversity, as a presumable consequence of recurrent admixture events among isolated populations (Aleksić and Geburek 2014). This is contrary to what has been observed in conifers with wide distribution range, where high genetic diversity relies on their vast population sizes (Vendramin et al. 1999;Bouillé and Bousquet 2005).
On the one hand, genetic diversity for the plastidial markers was similar to that obtained by Sánchez-Robles et al. (2014b) and Cobo-Simón et al. (2020). The larger sampling in our study did not result in an increase in cpSSR genetic diversity, suggesting that populations show similar diversity values along the complete range. Regarding the nSSR markers, the genetic diversity obtained in our study revealed that, overall, the pinsapo fir exhibits a high genetic diversity (overall mean H e = 0.523; overall mean N e = 2.63), which was largely congruent with previous results (mean H e = 0.596; mean N e = 2.83), using the same markers (Cobo-Simón et al. 2020).
Remarkably, there was incongruence in terms of genetic diversity between cpSSRs and nSSRs, which showed different patterns for the different mountain ranges. This could be due to the fact that genetic diversity can be considered from two different methodological perspectives for haploid and diploid genomes: genetic identity and heterozygosity. Considering that plastidial markers are more sensitive to genetic drift than nuclear ones (Birky et al. 1983;Birky 1995), populations with a small effective size would be much more affected by stochastic effects than large populations. Therefore, the former ones could differentiate more easily by the fixation of different alleles (Wright 1982;Slatkin 1985). Moreover, being other factors equal (as mutation rates), this would increase the global genetic diversity when all populations are considered as a whole (Wahlund 1928). In this way, Sª Bermeja, with the smallest size in comparison to the rest of mountain ranges, could be the most affected one by drift processes. This is supported by the greatest genetic differentiation among its populations. Moreover, its F IS values, the lowest ones in comparison to those from Sª de Grazalema and Sª de las Nieves, might also be explained by drift, but by selection based on inbreeding depression as well (Table 2). By contrast, due to their larger sizes, Sª de Grazalema and Sª de las Nieves mountain ranges would be less affected; congruently, populations within these ranges were genetically more similar (lower F ST values).
It is noteworthy that a potential relationship between genetic diversity and altitude, which might be based on lithology, slope, land-use history, or climate change refuges (Feurdean et al. 2013), was not found. In addition, a relationship of isolation by distance was detected for the plastidial  markers, resulting in lower genetic diversity values (this effect was not detected for the nuclear ones).
Regarding the nSSRs, it is remarkable the high presence of null alleles found for most loci, as well as the presence of endogamous effects when all populations were considered as a whole or grouped into the mountain range they belonged to (with the exception of Sª Bermeja). However, these effects disappeared when populations were considered independently. This lack of congruence could be due to the Wahlund effect (Wahlund 1928), where an excess of observed homozygotes (greater expected than observed heterozygosity) results from merging of populational units rather than to the presence of null alleles or inbreeding effects themselves, indicating that the mountain ranges as a whole could be composed by small nuclei isolated from each other, which is consistent with the isolation effect found.

Low genetic structure but incipient differentiation
The results obtained by the AMOVA, the Bayesian analyses and the phylogeographic tests indicated a low genetic structure when all populations were considered globally, which is congruent with a previous study using AFLP markers (Sánchez-Robles et al. 2014b). However, there are some signs of genetic structure at a mountain range level. Overall, the cpSSR markers showed greater genetic differentiation (F ST = 0.196) than the nSSR markers (F ST = 0.099). Additionally, the restricted pollen dispersal of the genus Abies in general (Mazier et al. 2008) and A. pinsapo in particular (average pollen dispersal distance of 113-227 m; Sánchez-Robles et al. 2014a) would increase the drift effect due to the limited gene flow, also for cpSSRs.
The results obtained by the analyses of molecular variance (AMOVAs), population covariance graphs and the Monmonier algorithm were congruent with the highest differentiation found in Sª Bermeja, which had the smallest population size of the three mountain ranges and, therefore, where drift could act faster (Wright 1982;Nora et al. 2011). The F ST value found for both markers were two times higher for Sª Bermeja than for Sª de las Nieves and Sª de Grazalema populations, suggesting that Sª Bermeja populations, despite their proximity (see Fig. 1), were more differentiated than those from Sª de las Nieves and Sª de Grazalema, both with a similar genetic differentiation value. Moreover, the phylogeographic tests for the cpSSR markers indicated that Sª Bermeja was genetically closer to Sª de Grazalema than to Sª de las Nieves. This relationship was also indicated by the genetic covariance graphs (Fig. 2b), where population connections between these two mountain ranges showed greater significance values. However, we found this phylogeographic signal between both mountain ranges only for one marker. This, together with the high number of connections, although less significant, present between Sª Bermeja and Sª de las Nieves populations in the genetic covariance graphs, suggest that this differentiation is very weak, and that it may indicate a practically contemporary separation of Sª Bermeja from the other two mountain ranges.

Recent divergence and demographic decline of the pinsapo fir
In this regard, the results obtained by the ABC analyses point to 'synchronic divergence with general bottleneck' scenario as the most probable scenario for the cpSSR markers. For nSSRs, the 'synchronic divergence' among the three mountain ranges was the most probable one. However, according to the posterior probability obtained (95% CI), there was no significant difference between both scenarios. However, population bottlenecks were also suggested by the Migraine analysis of nSSR and cpSSR but just for Sª de Grazalema and Sª Bermeja mountain ranges. Sª de las Nieves would have maintained its effective population size in spite of the population decline. The ABC-and Migraine-based estimates for divergence times and effective population sizes should be considered a crude approximation since these algorithms do not take into account continuous migration or pollen-mediated gene flow, which counteracts differentiation, underestimating the assessed divergence and bottleneck time values (van Loo et al. 2015;Hrivnák et al. 2017).
Evidences of genetic isolation and genetic drift found for geographically isolated populations within each mountain range indicated a restricted gene flow. Furthermore, the fact that no significant effect of inbreeding was found when populations were considered individually, even for the most isolated and less diverse populations, supported that the current distribution of A. pinsapo could have originated from a relatively recent fragmentation (given the longevity of the species), probably due to unfavorable environmental conditions, recurrent fires and/or human activity (Becerra Parra 2006;Soto García 2006;Linares et al. 2009;Esteban et al. 2010), still maintaining levels of heterozygosity that are high enough to allow these symptoms to manifest.
A. pinsapo would have suffered expansion and contractions in its distribution range due to climatic oscillation along the Holocene (Alba-Sánchez et al. 2010). The Abies macro-and microfossil records found in localities in the south of the Iberian Peninsula, where it is not present today, leave evidence of the presence of the Spanish fir after the Holocene Thermal Maximum (c. 6000 yr BP), even in the last two millennia (Pantaléon-Cano et al. 2003;Olmedo-Cobo et al. 2017;Pardo-Martínez et al. 2021). Moreover, a palaeoecological study for the recent Holocene indicates an initial phase of great stability during the Islamic period, followed by significant degradation after Christian reconquest (ca. 1487-1530 cal AD). The Modern Era (1530-1800 cal AD) marked the beginning of forest management with cycles of conservation and deforestation related to the use of the Spanish fir for the construction of road network and iron foundry, among other uses (Guzmán Álvarez et al. 2013). After that, there was a progressive reduction in its representation until current date (Alba-Sánchez et al. 2018).
Finally, all western Mediterranean firs are narrow endemics and many populations have undergone considerable reductions in historical times due to an intensive management that has stretched well into the twentieth century. For A. pinsapo, many studies and recent historical photographs (Barbey 1931;Ceballos and Vicioso 1933;Gil Albarracín 2002;Garrido Domínguez 2006;Becerra Parra 2006;among others) suggest that the pinsapo fir landscapes have suffered depletion since the middle of the eighteenth century and, especially, during the nineteenth century, as a consequence of higher demand for coal, firewood and beams for construction. Additionally, the overgrazing derived from the numerous goat and sheep herds would have avoided further regeneration, added to the action of fires in reducing the area of the species. Nevertheless, the recent (last half of the twentieth century) expansion of several fragmented populations of A. pinsapo from scattered remaining stands, following the implementation of conservation measures (Guzmán Álvarez et al. 2013), could have diluted some populational differentiations and bottlenecks (e.g., Sª de Grazalema).

Implications for conservation
Our results showed that A. pinsapo populations have apparently been able to maintain relatively high levels of plastidial DNA diversity. Moreover, their genetic variation is little structured and the evolutionary "uniqueness" of single populations (in terms of endemic haplotypes; see Petit et al. 2003;Hampe and Petit 2005) is, therefore, relatively low.
The observed molecular patterns have several implications for the conservation and management of the species. The 44 A. pinsapo populations, taken as a whole, represent two delimited major evolutionarily units (ESUs): Sierra de Grazalema-Sierra de las Nieves vs. Sierra Bermeja (Moritz 1994). Most stands could serve as sources of plant material for restoration programs, although populations from Horno de la Miera, Sª del Pinar, and Sª de Zafalgar (from Sª de Grazalema), Los Sauces, Cubero, Caucón, Lajar, Froncaile, Yedra, Navas, Pilones, Cañada del Cuerno, Cañadas bajo, Cañada de Enmedio, and Cañada de las Ánimas (from Sª de las Nieves), and stands from the southern area of Sª Bermeja seem to be the most suitable ones in terms of haplotype richness (for more details, see Fig. S4). Genetic factors such as inbreeding depression may be of relatively little importance for population performance. Hence, conservation measures should mainly dealing to diminish the probability of stochastic events such as fires and the prevention of further habitat loss due to climate change and/or human impact.
The endangered situation of this species under a climate change scenario might be enhanced, regarding conservation and adaptive management insights. Recent climate change has been related to increasing drought and A. pinsapo forest decline. Recurrent mortality events, mainly at low-elevation sites within the elevation range of A. pinsapo have been noted, whilst the high-elevation sites yielded a recent growth enhancement and increasing forest cover. Given the limited migration potential of this climate relict, adaptive management programs should consider assisted migration, upward in elevation or searching for potential sites by species distribution modelling for the conservation of this species. Furthermore, despite it is not plausible the existence contrasting lineages inhabiting every single population within the mountain ranges, the existence of gene frequency shifts associated with drought-induced mortality suggests some selective pressures that might enhance the differentiation (see for instance, Kuparinen et al. 2010;Cobo-Simón et al. 2021). Finally, collections of plant material for ex-situ activities or restoration programs could help by maximizing the number of populations rather than the number of individuals within populations.

Conclusions
This study contributes to a better understanding of the effects of historical climatic and recent anthropogenic changes on forest species inhabiting managed landscapes. Analyses suggest the existence of significant bottlenecks. This, together with genetic drift and isolation by distance, underlies the current genetic diversity distribution of the pinsapo fir across mountain ranges.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.