# Signature of mid-Pleistocene lineages in the European silver fir (Abies alba Mill.) at its geographic distribution margin

Caroline Scotti-Saintagne, Thomas Boivin, Marie Suez, Brigitte Musch, Ivan Scotti, Bruno Fady

► To cite this version:

Caroline Scotti-Saintagne, Thomas Boivin, Marie Suez, Brigitte Musch, Ivan Scotti, et al. Signature of mid-Pleistocene lineages in the European silver fir (Abies alba Mill.) at its geographic distribution margin. Ecology and Evolution, 2021, 11 (16), pp.10984-10999. 10.1002/ece3.7886. hal-03318278

HAL Id: hal-03318278

https://hal.inrae.fr/hal-03318278

Submitted on 9 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Signature of mid-Pleistocene lineages in the European silver fir (Abies alba Mill.) at its geographic distribution margin

Caroline Scotti-Saintagne<fnref n="1" /> | Thomas Boivin<fnref n="1" /> | Marie Suez<fnref n="1" /> | Brigitte Musch<fnref n="2" /> |
Ivan Scotti<fnref n="1" /> | Bruno Fady<fnref n="1" /> ID

## Abstract

In a conservation and sustainable management perspective, we identify the eco-
logical, climatic, and demographic factors responsible for the genetic diversity pat-
terns of the European silver fir (Abies alba Mill.) at its southwestern range margin
(Pyrenees Mountains, France, Europe). We sampled 45 populations throughout the
French Pyrenees and eight neighboring reference populations in the Massif Central,
Alps, and Corsica. We genotyped 1,620 individuals at three chloroplast and ten nu-
clear microsatellite loci. We analyzed within- and among-population genetic diver-
sity using phylogeographic reconstructions, tests of isolation-by-distance, Bayesian
population structure inference, modeling of demographic scenarios, and regression
analyses of genetic variables with current and past environmental variables. Genetic
diversity decreased from east to west suggesting isolation-by-distance from the Alps
to the Pyrenees and from the Eastern to the Western Pyrenees. We identified two
Pyrenean lineages that diverged from a third Alpine-Corsica-Massif Central lineage
0.8 to 1.1 M years ago and subsequently formed a secondary contact zone in the
Central Pyrenees. Population sizes underwent contrasted changes, with a contrac-
tion in the west and an expansion in the east. Glacial climate affected the genetic
composition of the populations, with the western genetic cluster only observed
in locations corresponding to the coldest past climate and highest elevations. The
eastern cluster was observed over a larger range of temperatures and elevations. All
demographic events shaping the current spatial structure of genetic diversity took
place during the Mid-Pleistocene Transition, long before the onset of the Holocene.
The Western Pyrenees lineage may require additional conservation efforts, whereas
the eastern lineage is well protected in in situ gene conservation units. Due to past
climate oscillations and the likely emergence of independent refugia, east-west ori-
ented mountain ranges may be important reservoir of genetic diversity in a context
of past and ongoing climate change in Europe.

<footnote n="1">¹INRAE, Ecologie des Forêts
Méditerranéennes (URFM), Avignon, France</footnote>

<footnote n="2">²ONF, UMR 0588 BioForA, Orléans, France</footnote>

<footnote>**Correspondence**

Bruno Fady, Ecologie des Forêts
Méditerranéennes (URFM), INRAE, Avignon
Cedex 9, France.
Email: bruno.fady@inrae.fr</footnote>

<footnote>**Funding information**

Eranet Biodiversa, Grant/Award Number:
ANR 12 EBID 000 3; Parc naturel régional
des Pyrénées Ariégeoises; European
Commission, Grant/Award Number: 676876</footnote>
KEYWORDS

admixture, conservation, demography, French Pyrenees, genetic diversity, keystone species, Phylogeography, Quaternary climate, range-edge

# 1 | INTRODUCTION

As the world's biodiversity continues to decline (Secretariat of the Convention on Biological Diversity, 2020), the Post-2020 Biodiversity Framework of the Convention on Biological Diversity (CBD) and the EU Biodiversity Strategy for 2030 have called for transformative changes that will halt biodiversity decline and maintain ecosystem, species, and genetic diversity for the benefit of people (Hoban et al., 2020).

Among these transformative changes is the need to better consider and protect genetic diversity. It is a key component of resilience and adaptability in its own right and it can be harnessed for the benefit of people using nature-based innovative management suitable at times of climate uncertainty (Fady et al., 2016; Hoban et al., 2020). Knowledge of how genetic diversity is structured is important for delineating conservation strategies and identifying areas where conservation efforts should be focused (Moritz & Faith, 1998). Populations with the highest genetic diversity are likely to display increased resilience and those harboring different ancestral alleles in past refuge zones are likely adapted to different environments, a crucial information for prioritizing population conservation (Fady et al., 2016; Fernandez-Fournier et al., 2021; Hampe & Petit, 2005).

A significant number of phylogeographic studies have demonstrated the role of Holocene and the Last Glacial Maximum in the distribution of species, populations, and their genetic diversity in the world's temperate regions, emphasizing climate change and dispersal as the main drivers of observed and modeled patterns. Current species geographic distribution in the Northern Hemisphere is often the result of a rapid south-to-north colonization during the warm Holocene period that followed the Last Glacial Maximum (Hewitt, 1999). The consequence for European plants is a strong latitudinal gradient of genetic diversity, with often higher values in the south than in the north (Petit et al., 2005), although reverse clines have been observed, for example, for Artic-Alpine plants which also demonstrate lower diversity margin effects (Hirao et al., 2017). In the Mediterranean basin, where major glacial refugia were located, genetic diversity shows an additional, longitudinal gradient, particularly in trees and insects, with higher diversity in the east than in the west (Fady & Conord, 2010), interpreted as resulting from both east-west Holocene migrations and differential local drift during the Last Glacial Maximum (see e.g., the widespread forest tree *Pinus nigra*, Scotti-Saintagne et al., 2019).

Additionally, less widespread evidence indicates that earlier Pleistocene climate and geological events have also shaped current within-species genetic diversity (Riddle, 2016). While Holocene warming has undoubtedly affected genetic diversity patterns in many species, earlier events of colonization and successive vicariance at different times (Kropf et al., 2006; Schmitt et al., 2006) may also

have had a strong impact on current genetic diversity (Avise, 2000). This may be particularly true for long-lived, sessile organisms such as forest trees. In Europe, forest tree taxa have repeatedly shifted their distributions and disappeared locally throughout the Pleistocene (Magri et al., 2017), resulting in earlier and deeper genetic differentiation patterns than those determined by the last glacial cycle, such as with *Cupressus sempervirens* (Bagnoli et al., 2009) or *Quercus cerris* (Bagnoli et al., 2016) for example.

Past climate changes have had a significant impact on the genetic structure of forest tree populations, and current human-induced, strong climatic changes can be expected to modify it again. Focusing on populations at the warm margin of their range is of particular relevance as they represent limits to species distributions, which habitats may be strongly altered by climatic change (Crawford, 2008). Silver fir (*Abies alba* Mill.) is a target European conifer species for conservation programs. It is an economically important species as well as a keystone species which maintains high biodiversity forest ecosystems (Mauri et al., 2016). While not threatened *per se*, widespread mortality has been observed recently, particularly at habitat and geographic margins (Cailleret et al., 2014; Camarero et al., 2011).

The distribution range of silver fir is wide but patchy, from southern Poland to northern Spain (Wolf, 2003). The genetic diversity of *Abies alba* is high, indicating that despite geographic patchiness, gene flow is high, but generally declining from east to west (Fady & Conord, 2010; Mosca et al., 2012). At a broad scale, genetic diversity is structured into several lineages resulting from past climate changes. The Italian Peninsula populations were found to have diverged from Balkan populations during the middle Pleistocene (Belletti et al., 2017; Piotti et al., 2017), while the Eastern and Western Alps lineages derived from Balkan and Italian refugia are believed to have merged during the Holocene (Cheddadi et al., 2014; Liepelt et al., 2009). At their south-westernmost Alpine margin, silver fir populations belong to two different genetic groups, the Southern Alps cluster which also includes the Massif Central populations and the Pyrenees cluster (Liepelt et al., 2009; Sancho-Knapik et al., 2014). Evidence of local adaptation to drought and frost has been reported in both the Alps and the Pyrenees (Csilléry et al., 2020; Matías et al., 2016; Roschanski et al., 2016).

Compared to the Alps, where the spatial distribution of genetic diversity has been intensively studied, little information is available for the Pyrenees. The geographic distribution of silver fir there is fragmented, as a result of both past climate changes and human impact during the Antiquity (grazing or logging for industrial activity) and until the industrial revolution (Ejarque et al., 2009; Galop & Jalut, 1994; Pérez-Sanz et al., 2013). This fragmentation has led to the isolation of old-growth (sensu FAO, 2015) native forests, many of them in marginal habitats where biodiversity conservation issues are most pressing in Europe (Sabatini et al., 2018).
The Pyrenees are the south-westernmost large alpine-type ecosystem in Europe making it a hotspot of diversity and endemism (Nagy et al., 2003; Ozenda, 1985; Villar & Dendaletche, 1994) as well as a glacial refugium for many mountain species (Duriez et al., 2007; Horn et al., 2009; Taberlet et al., 1998). The Pyrenees' longitudinal structure, stretching 450 km from west to east between the Atlantic Ocean and the Mediterranean Sea, confers them unique ecological characteristics. Along the Pyrenees, the landscape changes from oceanic and lower elevation in the west to Mediterranean and variable elevation and geomorphology toward the east (Ninot et al., 2007). These longitudinally contrasted environments have led to a west-east genetic divergence attributed to secondary refugia in many animal species, such as the mountain ringlet butterfly *Erebia epiphoron* (Schmitt et al., 2006), the ground-dwelling spider *Harpactocrates ravastellus* (Bidegaray-Batista et al., 2016), the pine processionary moth *Thaumetopoea pityocampa* (Rousselet et al., 2010), and plant species such as snapdragons *Antirrhinum* sp. (Liberal et al., 2014), *Rhododendron ferrugineum* (Charrier et al., 2014), the clover *Trifolium alpinum* (Lauga et al., 2009), the European oaks (Petit et al., 2002), the European beech *Fagus sylvatica* (Magri et al., 2006), and silver fir (Matías et al., 2016; Sancho-Knapik et al., 2014). In these animal and plant phylogeography studies, identified lineages are always shared between France and Spain at any given longitude, indicating that the east-west ridge of high peaks in the Central Pyrenees did not act as a barrier to migration.

Here, we analyze the phylogeography of silver fir with a Northern Pyrenean focus, through detailed landscape-scale sampling and using neighboring silver fir populations (located in the higher latitude French mountain ranges of the Massif Central, the southern Alps, and Corsica) as external references. We addressed the following questions:

1. What is the evolutionary origin of the Pyrenean lineages? How and when did they diverge from each other and from Alpine lineages?

2. What are the locations of the geographical boundaries between genetic lineages along the Pyrenees, and is there admixture between them?

3. Which demographic events related to past climate oscillations have left the most significant imprint on genetic diversity and can this information be used for furthering *in-situ* gene conservation efforts in France?

For this, we sampled the Northern Pyrenees only, in France, thus assuming that, as for all current mid-latitude Northern Hemisphere conifers, *Abies alba* ancestors expanded southward across Europe from circumboreal high latitudes and diversified during the general cooling of the late Eocene and Neogene period, starting some 25 million years ago and culminating during the Pleistocene (Barnosky, 1987; Rundel, 2019; Xiang et al., 2015). We expect that our results will provide a phylogeographic picture valid for the entire Pyrenees (both France and Spain) as range-wide evidence from

other species does not demonstrate the existence of separate lineages on each side of the central high peak ridge of the Pyrenees.

## 2 | METHODS

### 2.1 | Population sampling and DNA extraction

Sampling consisted of an average of 30 trees per population from a total of 52 populations (Table S1.1 in Appendix S1). Forty-four populations were sampled spanning a 400 km east/west axis in the French Pyrenees between the Mediterranean and the Atlantic coasts. The Western Pyrenees populations (PYR_W) were sampled in the area between the Bay of Biscay and the Pic d'Anie, where elevation does not exceed 2,000 meters above sea level. The Central Pyrenees populations (PYR_C) were sampled between the Pic d'Anie and the Puymorens pass, in the highest part of the Pyrenees with peaks higher than 3,000 m. The Eastern Pyrenees populations (PYR_E) were sampled in the area east of the Puymorens pass down to the Mediterranean coast; their highest elevation (Puigmal peak, 2,910 m) is intermediate between that of the Central and Atlantic Pyrenees. To compare the genetic diversity observed in the Pyrenees with surrounding mountain ranges, we also sampled four populations in the Massif Central: three in the Western Alps and one in Corsica (collectively defined as "Alpine populations" as they belong to the same lineage, Liepelt et al., 2009).

We sampled wild populations large enough (1.5–3 ha) to allow sampling trees at least 20 meters apart from each other. Sampling consisted of three cambium disks per tree, frozen in liquid nitrogen in portable devices during collection. DNA extraction followed the protocol of Demesure et al. (1995).

### 2.2 | Microsatellite genotyping

Among the 52 sampled populations, 50 were analyzed at both chloroplast (cpSSR) and nuclear (nuSSR) microsatellites and two from Western Alps at only one type of marker (VES at nuSSR and LURE at cpSSR).

Three cpSSRs, Pt30141 (Liepelt et al., 2001), Pt15169, and Pt71936 (Vendramin & Ziegenhagen, 1997), were amplified by polymerase chain reaction (PCR) using the protocol of Vendramin and Ziegenhagen (1997). cpSSR pt15169 displays an 18 base pair insertion-deletion (indel), whereas the other two display mononucleotide indel patterns. Genotypes were scored using an ABI PRISM 3100 automatic sequencer (Applied Biosystems) at ONF in Orléans (France), and allele size was assigned using GeneMapper v3.7. In addition, we used ten nuclear microsatellites (nuSSR): six from Cremer et al. (2006) (SF1, SFb4, SFb5, SF50, SF78, and SF333) and four from Hansen et al. (2005) developed for *Abies nordmanniana* and transferable to *A. alba* (NFF3, NFH15, NFF7, NFH3). One nuSSR (SF1) displays a trinucleotide pattern, whereas all the others are dinucleotides. The ten nuSSR were genotyped in two multiplexes
using the Multiplex PCR Kit (QIAGEN). The final PCR volume was optimized to 10 µl, and the PCR mix for both multiplexes was 5 µl of the QIAGEN Multiplex PCR Master Mix, 2 µl of DNA (10 ng/µl), and 0.2 µM primers. For both multiplexes, the PCR thermal profile was the following: an initial step at 95°C for 15 min, followed by 35 cycles at 95°C for 30 s, 57°C for 90 s, and 72°C for 90 s, with a final 30-min extension step at 60°C. PCR products were scored at INRAE in Avignon (France) using an AB 3730 XL automatic sequencer (Applied Biosystems) with LIZ-600 as internal size standard.

## 2.3 | Statistical analysis

Statistical computing of data was made using R version 3.6.3 (R Core Team, 2020) on RStudio (RStudio Team, 2019).

### 2.3.1 | Genetic diversity

Estimation of the null allele frequencies at each locus was performed in each population using the maximum likelihood approach and the EM algorithm of Dempster et al. (1977) available in the Genepop v4.6.9 R package (Rousset, 2008). The proportion of alleles that are unique (in only one individual in the population) $A_u$ was computed for each population using the StrataG v2.4.905 R package (Archer et al., 2017). The number of alleles $N_A$, the number of effective alleles $N_{A_e}$ (Nielsen et al., 2003), the genetic diversity $H_e$ corrected for sample size (Nei, 1978), and the inbreeding coefficient $F_i$ (for nuSSR) were estimated using Spagedi1.4b (Hardy & Vekemans, 2002). Allelic richness and private allelic richness were estimated using rarefaction methods, as implemented in ADZE (Szpiech et al., 2008), to take into account heterogeneity in sample sizes. The number of rarefied gametes was 50 for nuSSRs and 22 for cpSSRs. The distribution of genetic diversity among the sampled populations was explored using a principal correspondence analysis (PCA) performed on the allelic frequencies at the individual level using the R package Adegenet (Jombart, 2008).

### 2.3.2 | Isolation by distance and phylogeographic signals

Pairwise differentiation between populations and between mountain ranges was computed using GENEPOP 4.7.5 (Rousset, 2008). A phylogenetic tree using the distance matrix of $F_{ST}$ between populations was built with Poptree (Takezaki et al., 2010) using the neighbor-joining method (Saitou & Nei, 1987). We used the genotypic data of four Mediterranean firs, A. pinsapo (Boiss.), A. cilicica ((Antoine & Kotschy) Carrière), A. cephalonica (Loud.), and A. bormuelleriana (Mattf.) as outgroup, obtained from the study of Awad et al. (2014) produced on the same sequencer and using common Western Alps samples as reference for correctly assigning allele length. Tree topology was tested by performing 10,000 bootstrap

tests (Felsenstein, 1985). Linear regressions of pairwise population statistics (both $F_{ST}$ and $R_{ST}$, Slatkin (1995)) on geographical distances were performed to test for the presence of a phylogeographic signal, checking whether $R_{ST}$ remained higher than $F_{ST}$ after permuting alleles sizes over allelic states (Hardy et al., 2003) using SPAGeDi 1.4 (Hardy & Vekemans, 2002). Geographic distances between populations were computed using geographic Distance Matrix Generator version 1.2.3 (Ersts, American Museum of Natural History, Center for Biodiversity and Conservation, http://biodiversityinformatics.amnh.org/open_source/gdmg).

### 2.3.3 | Effective migration surface

Spatial population structure in the Pyrenees was explored by estimating effective migration surfaces (EEMS; Petkova et al., 2015). The EEMS method approximates a demographic model which evolves under equilibrium in time (isolation-by-distance) based on the stepping-stone model (Kimura & Weiss, 1964). EEMS allows to test for deviations from isolation-by-distance and to test whether the relationship between genetic diversity and geography differs among habitats due to a combination of environmental and historical factors (gene flow barriers and/or admixture zones). We modeled the forty-four Pyrenean populations on a dense regular grid, with migration between neighboring demes, and each population belonging to a different deme. Preliminary runs were made to obtain acceptance proportions of the MCMC between 20%–50%. Fine-tuning was performed by modifying the proposed variances as follows: mSeedsProposalS2 = 0.05, qSeedsProposalS2 = 0.01, mEffctProposalS2 = 2.5, qEffctProposalS2 = 0.002, and mrateMuProposalS2 = 0.05. Since the choice of the grid could influence results, we followed the authors' instructions by averaging the estimates over two different grids (200 and 300 demes). For each, we performed three analyses, each with a different random seed. All EEMS analyses were run for 10 million iterations, sampling every 10,000 iterations with a burn-in of 500,000. The six runs were combined and graphs were made using the accompanying R package rEEMSplots (http://www.github.com/dipetkov/eems).

### 2.3.4 | Blind approach to explore population genetic structure

We automatized runs of Structure 2.3.4 (Pritchard et al., 2000) using StrAuto (Chhatre & Emerson, 2017) to infer population structure and to assign individuals to populations. Due to possible gene flow between populations, we modeled (i) ancestry assuming that individuals have a mixed ancestry and (ii) allele frequencies assuming that frequencies in the different populations are likely to be correlated. Considering the number of genetic groups already observed in silver fir (Liepelt et al., 2009; Matías et al., 2016), the number of clusters expected in this study was three (Alps, Western Pyrenees, and Eastern Pyrenees); however, we tested the possibility of having
an under-structuring and therefore estimated the probability of ob-
serving from *K* = 1 to *K* = 8 clusters. For each value of *K*, we per-
formed ten independent iterations with a burn-in of 50,000 Markov
chain Monte Carlo (MCMC) steps, followed by a run of 500,000
steps. We used Structure Harvester 0.6.94 (Earl & vonHoldt, 2012)
and Clumpak (Kopelman et al., 2015) to compile and visualize results.
To infer the most likely number of genetic clusters, we selected the
value of *K* which gave the highest peak for Δ*K* (Evanno et al., 2005)
and checked that the curve of the likelihood for this value of *K*
reached a plateau as suggested by Pritchard et al. (2000).

2.3.5 | Demographic history

To infer the demographic history of silver fir, we compared scenarios using the coalescent/approximate Bayesian computation (ABC) framework implemented in DIYABC 2.1.0 (Cornuet et al., 2010) and both nuSSR and cpSSR data. To remove biases caused by gene flow between populations (which DIYABC cannot model), we only modeled the population having the highest ancestry coefficient in each cluster identified by the STRUCTURE analysis: Borc for the Western Pyrenees (Q = 0.91); Lesp for the Eastern Pyrenees (Q = 0.88); and Vent for the Southern Alps (Q = 0.93). An alternative sampling scheme would have been to pool a subset of individuals across

populations belonging to the same genetic cluster, but this could
have impacted the detection and the quantification of population
size changes (Chikhi et al., 2010). In addition, we included a popu-
lation (Stla, Q = 0.44) representing the contact zone between the
genetic lineages in the Pyrenees. To confirm that the results were
not biased by our choice of populations, we duplicated the analy-
sis by choosing another set of 4 populations: Issa (Q = 0.92), Bele
(Q = 0.87), Sigu (Q = 0.60), and Punt (Q = 0.86) for the Western
Pyrenees, Eastern Pyrenees, admixed Pyrenees and Southern Alps,
respectively (Figure 1).

In scenario S1 (Simple split), a single event of divergence took
place at time T₁ from a common ancestor to give rise to the Alps,
the Western Pyrenees, and the Eastern Pyrenees lineages, followed
by a later secondary contact at time T₂ between the two lineages
from the Pyrenees to give rise to the admixed lineage (Figure 2).
In scenario S2 (Hierarchical split), two events of divergence were
considered. The Alps and Pyrenees lineages diverged first at T₁ fol-
lowed by a divergence at T₂ between the Western and the Eastern
Pyrenees and then at T₃ by a secondary contact between lineages
from Pyrenees. In scenario S3 (Alps–Eastern Pyrenees Split followed
by two admixtures), at T₁ a first event of divergence gave rise to
both the Alps and the Eastern Pyrenees lineages, then at T₂ the two
lineages came into contact to give rise to the Western Pyrenees lin-
age. At T₃, a secondary contact took place between the Pyrenean

<figure><img src="image_3.png" /><figcaption>FIGURE 1 Location of the 52 silver fir populations sampled at the western edge of the geographic distribution of *Abies alba* (Mill.) in Europe and their membership to the three genetic clusters identified. The blue-green area represents the geographic distribution of silver fir (see http://www.euforgen.org/species/abies-alba/)</figcaption></figure>

populations belonging to the same genetic cluster, but this could
have impacted the detection and the quantification of population
size changes (Chikhi et al., 2010). In addition, we included a popu-
lation (Stla, Q = 0.44) representing the contact zone between the
genetic lineages in the Pyrenees. To confirm that the results were
not biased by our choice of populations, we duplicated the analy-
sis by choosing another set of 4 populations: Issa (Q = 0.92), Bele
(Q = 0.87), Sigu (Q = 0.60), and Punt (Q = 0.86) for the Western
Pyrenees, Eastern Pyrenees, admixed Pyrenees and Southern Alps,
respectively (Figure 1).

In scenario S1 (Simple split), a single event of divergence took
place at time T₁ from a common ancestor to give rise to the Alps,
the Western Pyrenees, and the Eastern Pyrenees lineages, followed
by a later secondary contact at time T₂ between the two lineages
from the Pyrenees to give rise to the admixed lineage (Figure 2).
In scenario S2 (Hierarchical split), two events of divergence were
considered. The Alps and Pyrenees lineages diverged first at T₁ fol-
lowed by a divergence at T₂ between the Western and the Eastern
Pyrenees and then at T₃ by a secondary contact between lineages
from Pyrenees. In scenario S3 (Alps–Eastern Pyrenees Split followed
by two admixtures), at T₁ a first event of divergence gave rise to
both the Alps and the Eastern Pyrenees lineages, then at T₂ the two
lineages came into contact to give rise to the Western Pyrenees lin-
age. At T₃, a secondary contact took place between the Pyrenean
<figure><img src="image_4.png" /></figure>

**FIGURE 2** Demographic scenarios tested to infer past divergence events between four populations of silver fir. Four demographic scenarios are tested: Scenario S1 (Simple split), Scenario S2 (Hierarchical split), Scenario S3 (Alps-East Pyrenees Split followed by two admixtures), and Scenario S4 (Alps-West Pyrenees Split followed by two admixtures). In the four scenarios, fluctuation of the population sizes was allowed along final branches. Additional details are given in the main text

lineages to give rise to the admixed population of the Pyrenees. Finally, scenario S4 (Alps-West Pyrenees Split followed by two admixtures) is as in S3, but the Alps and Western Pyrenees lineages diverged first and came into contact later to give rise to the Eastern Pyrenees population. In the four scenarios, fluctuation of population sizes was allowed along final branches (at time $T_a$, $T_b$, $T_c$, $T_d$ for the Alps, Eastern Pyrenees, Western Pyrenees, and admixed Pyrenees populations, respectively).

We defined the prior distributions of the demographic parameters based on results obtained in the literature. All details of our approach are explained in Table S3.1. Both nuSSRs and cpSSRs were simulated using the generalized stepwise mutation model (Estoup et al., 2002; Zhivotovsky et al., 1997). The mean mutation rate of nuSSRs ($\mu Nu_{SSR}$) was drawn from a log-uniform distribution constrained between $10^{-4}$ and $10^{-3}$, while for cpSSR, $\mu Cp_{SSR}$ was sampled between $10^{-8}$ and $10^{-3}$. Priors for times of demographic events were drawn from a normal distribution of mean equal to 25,000 and standard deviation equal to 12,500. Priors for demographic sizes were drawn from a normal distribution of mean equal to 5,000 and standard deviation of 2,500. We applied the procedure "evaluate scenario-prior combination" (Cornuet et al., 2010), which is based on both a principal component analysis and a test of rank, to choose a subset of the available summary statistics, that did not under- or overestimate genetic distance in simulated data sets, compared to the empirical one. Eight summary statistics were chosen (four single-population summary statistics, and four two-population summary statistics, Table S3.2 in Appendix S3).

For each scenario, one million simulations were performed. The most-likely scenario was selected according to the results of a

polytomous weighted logistic regression (Beaumont, 2010). Due to the high number of summary statistics, we replaced them with components of a linear discriminant analysis. The most-likely scenarios were considered successfully recovered if the maximum posterior probability of the scenario was outside the 95% confidence intervals estimated for all other candidate scenarios (Cornuet et al., 2008). Confidence in scenario choice was evaluated by calculating type I error (classification error) and type II error from 1,000 pseudo-observed data sets selected within the 1% simulations closest to observed statistics for each scenario (Cornuet et al., 2010). The posterior distribution of demographic parameters was finally estimated for the best scenario using a linear regression on the 1% closest simulations and applying a *logit* transformation to parameter values. Since the effective population size $N_e$ depends on an unknown mutation rate value, we compared demographic fluctuations across time between genetic clusters by estimating the posterior distribution of the ratio between the current and past population diversity parameters using the ratio = $\theta_0/\theta_1$, with $\theta_0 = N_0\mu0$ (present) and $\theta_1 = N_1\mu1$ (past) (Barthe et al. (2017). Thus, by assuming that the mutation rate is constant through time ($\mu0 = \mu1$), the estimate of the demographic size ratio is no longer dependent on the (unknown) value of the mutation rate.

Bias and precision of parameter estimations were computed (mean relative bias (MRB), relative root mean square (RMSE), 50% and 95% coverage and Factor 2). The goodness-of-fit of the combination of the best model and parameter posterior distributions was appreciated using the model-checking procedure in which the observed values of four additional summary statistics (Table S3.2 in Appendix S3) are compared to the values obtained from 1,000
<figure><img src="image_5.png" /><figcaption>**FIGURE 4** Principal component analysis (PCA) based on individual allelic frequencies at eight nuSSR and one cpSSR haplotype in silver fir populations sampled in six mountain ranges. Colors represent mountain ranges: Massif Central (purple), Western Alps (dark blue), Corsica (light blue), Western Pyrenees (orange), Central Pyrenees (peach), and Eastern Pyrenees (red). The code names of the populations are indicated in boxes at the center of gravity of their allelic distribution. The ellipses summarize the dispersion of the point cloud of each population with the center of the ellipse being the center of gravity of the population</figcaption></figure>

<figure><img src="image_6.png" /></figure>

replicated data generated under the scenario which gave the best probability. We did not use the same summary statistics to avoid risk of overestimating the quality of the fit.

### 2.3.6 | Correlating between genetic diversity and environmental factors

Our goal here was neither to look for signatures of selection and adaptation, nor to reveal causal relationships between neutral genetic diversity and environmental variables. Instead, it was to test whether environmental characteristics in the Northern Pyrenees (mainly past and present climatic variables) could be used as a proxy

to model the organization of the genetic diversity in the Northern Pyrenees.

Environmental characteristics were obtained for each Pyrenean population. Slope, altitude, latitude, longitude, and nineteen bioclimatic variables from contemporary periods (1970–2000) and older periods (Mid-Holocene approx. 6,000 years before present (BP), Last Glacial Maximum approx. 22,000 years BP, Last interglacial approx. between 120,000 and 140,000 years BP) were extracted from the WorldClim database (Hijmans et al., 2005). The descriptions of the environmental variables are reported in Appendix S4, Table S4.1 and the values for each population in Table S4.2.

For each genetic variable (Appendix S1, Table S1.2, Table S1.3), Pearson correlations with environmental variables were computed
<figure><img src="image_7.png" /><figcaption>FIGURE 5 Structure analysis of 51 silver fir populations sampled in six mountain ranges in France. Each vertical bar represents a silver fir tree, organized by the longitude of the sampled populations in Western Pyrenees, Central Pyrenees, Eastern Pyrenees, Massif Central, Western Alps, and Corsica. Colors represent group membership. The proportion of each color in a bar represents the probability of assignment of the individual to each of the groups. Details of the populations are given in Appendix S1 Table S1.1</figcaption></figure>

by applying a local false discovery rate (lfdr, Fdrtool R package, Strimmer, 2008) equal to 0.05. The presence of multicollinearity was tested and confirmed by calculating the variance inflator factor (VIF, car R Package, Fox & Weisberg, 2018) and using a threshold of VIF > 10 to consider multicollinearity. Since all tested variables caused multicollinearity, we then applied a forward stepwise regression which prevents multicollinearity problems in generalized linear models. The Bayesian information criterion (BIC) was used to find the most parsimonious regression models able to explain the genetic variables. For each genetic variable to explain, the best explicative environmental variables were finally included in a linear model to estimate the adjusted R-square of the regression model and to estimate the regression coefficient and the *p* value of each retained environmental variables.

Environmental variables and data, genetic diversity parameters and many method descriptions and parameters are available in the Appendices (S1 to S5). All appendices are deposited open access at the following address: https://doi.org/10.15454/X6ZJU0.

3 | RESULTS

3.1 | Genetic diversity

The SF333 and SF50 nuSSR markers had a frequency of null alleles greater than 0.05 in more than half of the populations analyzed and were discarded. The combination of the three cpSSRs generated 72 haplotypes (eight variants for pt71936, 15 for pt30141, and three for ssr15169). The haplotype frequency distribution was L-shaped, with a large number of low-frequency haplotypes (60 with frequency lower than 0.01) and one common haplotype (frequency = 0.48). Population and mountain-range genetic diversity levels are reported in Appendix S1, Table S1.3. Overall, genetic diversity was higher in Alpine populations (H<sub>e</sub> at nuSSR between 0.72 and 0.74) than in the Pyrenees (0.55–0.65), while genetic diversity in the Pyrenees was highest in the eastern part of the range (Figure 3, Figure S1.1 in Appendix S1).

Figure 4 displays the first two axes of an individual genotype-
based PCA. The first two axes of the PCA explain 16.4% of the total
variance (9.1% and 7.3% for axes 1 and 2, respectively). The centers
of gravity of the Western Alps, the Massif Central, and Corsica popu-
lations lie in the upper right quadrant of the PCA plot, while the cen-
ters of gravity of the Pyrenean populations are scattered throughout
the remaining quadrants. Population and mountain-range dot scat-
ters overlap extensively, suggesting weak genetic differentiation.

Pairwise genetic differentiation (F<sub>ST</sub>, R<sub>ST</sub>) between mountain ranges was significant at both nuSSR and cpSSR and globally stronger between the Alps and the Pyrenees than among zones within the Pyrenees (Table S1.4 in Appendix S1). The F<sub>ST</sub> matrix-based phylogenetic tree had a deep split separating Alpine and Pyrenean populations (Figure S1.2 in Appendix S1). The Pyrenean branch had a ladder shape, suggesting isolation-by-distance from east to west in the Pyrenees, as confirmed by the significant positive slope of the linear regressions of pairwise statistics (both F<sub>ST</sub> and R<sub>ST</sub>) on geographical distances. However, there was no significant difference between the slopes of the linear regressions of pairwise statistics based on allele identity (F<sub>ST</sub>) and allele size (R<sub>ST</sub>) meaning that populations could have been isolated for a relatively short period of time (in comparison to the mutation rate) or that gene flow between populations was strong. The EEMS confirmed the progressive colonization of the Pyrenees. However, in some parts, the presence of deviations from the step-by-step model was observed as either a stronger restriction of gene flow or an excess of gene flow (Appendix S1 - Figure S1.3).

3.2 | Genetic cluster inference

The most likely number of genetic clusters among all sampled pop-
ulations was *K* = 3 (deltaK = 247) (Figure S2.1a in Appendix S2).
Two geographically distinct genetic clusters were detected in the
Pyrenees (Western and Eastern Pyrenees clusters) while one cluster
grouped the Alpine populations (Figure 5). Western Alps populations
were more genetically homogeneous than populations from Massif
Central or from Corsica, which had respectively 15% and 10% of
individuals admixed with the Pyrenees lineage (vs. 0% in Western Alps). Western and Eastern Pyrenees samples belonged to two distinct ancestral parental populations, whereas Central Pyrenees populations (Juze, Bout, Stla, Beth) were split between the two clusters with a high degree of admixture. The two-cluster genetic structure of the Pyrenees was confirmed when the Alpine populations were removed from the analysis (Figures S2.1b and S2.2 in Appendix S2).

## 3.3 | Demographic history

Evaluations of the simulated data sets to test for the scenarios of divergence are provided in Appendix S3, Figure S3.1. The differences between the probabilities associated with each scenario were small, even for with scenario S2 the most probable (Figure S3.2, in Appendix S3), but only slightly more supported than the others. Redoing the analysis using a different set of four populations confirmed our first results (Figure S3.5, in Appendix S3). Further support to scenario S2 came from the phylogenetic tree obtained using $F_{ST}$ values (Figure S1.2 in Appendix S1). However, the demographic inference lacked both accuracy and power, with a high posteriori error rate ($p = .69$) for all four scenarios. For scenario S2, type I error (rejection of S2 when true) and the type II error (acceptation S2 when not true) were respectively 0.56 and 0.65. Model adequacy of the best scenario S2 is nonetheless acceptable since among thirty tests, only four observed summary statistics (the mean genetic diversity between pairs of Pyrenees clusters) deviated significantly from its simulated distribution ($p$-value < .05; Table S3.7 in Appendix S3). Parameter estimates were also reliable with a relative median of the absolute errors ranging from 0.185 to 0.416 (Table S3.8 in Appendix S3).

Considering that the four scenarios gave similar estimates for temporal events (T1, T2, and T3), population size, admixture, and mutation rates (see Appendix S3, Tables S3.3 to S3.6), we conservatively decided to estimate model parameters as the mean of the medians of all scenarios. The 95% coverage of the posterior probability corresponds to the lowest and the highest limit obtained for the four scenarios. The Pyrenean and Alpine lineages diverged from an ancestral population 27,875 generations ago (95% confidence interval [CI] = 47,900–11,600). The Pyrenean lineage then separated 24,100 generations ago (95% CI = 39,300–12,900) into the Eastern and Western Pyrenean lineages. These came into contact 21,775 generations ago (95% CI = 37,900–11,400) to establish an admixture zone. Considering that silver fir is sexually mature at 30 to 40 years (Mauri et al., 2016), the mean of the medians for the divergence between the Alpine and the Pyrenean lineages was 836,250–1,115,000 years ago. The mean of the medians for the divergence between the two Pyrenean lineages was 723,000–964,000 years ago and the mean of the medians for their merging into a secondary contact zone located in the Central Pyrenees 653,250–871,000 years ago. For each of the four populations, we detected a synchronic population resizing 461,550–615,400 years ago. Thus, all major events of lineage divergence took place during the Mid-Pleistocene Transition (0.7–1.25 million years ago), long before the Last Glacial Maximum (LGM) and

the onset of the Holocene. All those estimates were confirmed by the second approximate Bayesian computing analysis using four different populations (Figure S3.5, Table S3.9 to S3.11).

In the Eastern Pyrenees, the mode of posterior distribution of the ratio of current to past size ($N_{current}/N_{past}$) was greater than 1 over 95% of its distribution for both nuSSR and cpSSR (mode of $N_{current}/N_{past}$ = 182 and 235 for nuSSR and cpSSR, respectively) indicating a population expansion event (Figure 6). In contrast, in the Western Pyrenees, the $N_{current}/N_{past}$ ratio was smaller than 1 over 55% and 95% of its distribution for nuSSR and cpSSR respectively (modes of $N_{current}/N_{past}$ = 0.69 and 0.05 for nuSSR and cpSSR, respectively) indicating a population contraction event. In the Western Alps and the admixture zone of the Pyrenees (Figure S3.4 in Appendix S3), the results were less consistent between the two types of markers, and a demographic equilibrium ($N_{current}/N_{past}$ = 1) cannot be excluded from the 95% posterior distribution envelope.

## 3.4 | Correlations between environmental variables and genetic parameters in the Pyrenees

The Pearson correlations gave 206 significant correlation tests with a local false discovery rate lower than 0.05 (Table S4.3 in Appendix S4). However, only 47 remained significant after the forward stepwise variable selection applied using a generalized linear model (Table S4.4 in Appendix S4). There were 12 significant correlations with precipitation seasonality (BIO15), nine with longitude, eight with mean temperature during the wettest quarter (BIO8), seven with mean temperature during the driest quarter (BIO9), four with temperature seasonality (BIO4), two with precipitations during the driest month (BIO14), and one significant correlation each with isothermality (BIO3), minimum temperature during the coldest month (BIO6), precipitations during the driest quarter (BIO17), precipitations during the warmest quarter (BIO18), and latitude. There was no significant difference (Pearson's chi-squared test) between the number of significant bioclimatic explanatory variables related to temperature and those related to precipitation (21 and 16, respectively) and the absence of difference remained when nuSSR and cpSSR were considered separately. Considering the period, Pearson's chi-squared tests indicated that cpSSR variables were more strongly associated with middle Holocene (6,000 years ago) bioclimatic variables, while the nuSSR variables were more strongly associated with last interglacial (LIG, 120,000 years ago) bioclimatic variables (Figure S4.4 in Appendix S4).

Among the significant regressions, high values of the membership coefficient for the western genetic cluster ($Q_w$) were negatively associated with both longitude and mean temperatures during the LIG wettest quarter (Figure 7). Populations belonging to the western cluster ($Q_w$ > 0.75) were associated with lower temperature (on average 3.94°C, standard deviation $SD = 2.13$°C), whereas the populations belonging to the eastern cluster ($Q_w$ < 0.25) were associated with a larger range of temperatures (on average 6.3°C, $SD = 4.02$°C).
<figure><img src="image_8.png" /></figure>

FIGURE 6 Comparison of the mode of the ratio between current and ancestral population sizes ($N_{current}/N_{past}$) in the western and eastern genetic clusters in French silver fir populations. Left columns: nuSSR; right column: cpSSR. The prior distribution is shown in gray dashed line, while the red dotted line is the posterior distribution. When the posterior distribution (in red) overlaps the prior distribution (mean of the prior = $log_{10}(N_{current}/N_{past}) = 0$, that is, the mode of the prior is $10^0 = 1$, that is, $N_{current} = N_{past}$), we cannot exclude the hypothesis of a demographic equilibrium. This is the case in the Central Pyrenees (Pyr_C) at nuSSR and cpSSR. On the other hand, we can reject the equilibrium hypothesis when the mode of the posterior distribution is greater than 1 over 95% of its distribution. This is the case in the Eastern Pyrenees at both nuSSR and cpSSR where the ratio of current and ancestral population sizes is greater than zero, indicating a strong population expansion about 15,000 generations ago (see Table S3.4 for the timing of the event).

<figure><img src="image_9.png" /></figure>

FIGURE 7 Plot of the regression of membership coefficient for the western genetic cluster $Q_w$ in French silver fir populations according to the best explicative environmental variables. (a) 2D plot, $Q_w$ according to the longitude (model 1). (b) 2D plot, $Q_w$ according to the wettest quarter during LIG. (c) 2D plot, residual of the model 1 according to the mean temperature at the wettest quarter during LIG. Red dots are for values of $Q_w > 0.75$ and correspond to the western genetic cluster, dark orange dots for $Q_w < 0.25$ and correspond to the eastern genetic cluster, pale orange dots for $0.25 < Q_w < 0.75$ and correspond to the admixed cluster

4 | DISCUSSION

4.1 | Past genetic divergence and demography of silver fir in the Pyrenees

The ABC-based estimates of timing of divergence and effective population sizes should be considered as approximations since they are closely related to the mutation rate, unknown in silver fir. Here, we used published mutation rates from other plants as priors. Also, the main limitation of DIYABC is the assumed absence of migration

among populations after they have diverged, which results in a probable under estimation of divergence time. In fact, our STRUCTURE analysis shows that admixture is still present in some easternmost and westernmost populations, outside of the admixture zone, indicating recent gene flow, either long distance or due to neighboring exotic plantation. Further, to convert generation time estimates into calendar time estimates, we assumed that sexual maturity is reached between 30 and 40 years in silver fir (Mauri et al., 2016), a reasonable approximation but subject to environmental variations. Finally, the use of only eight nuSSR and three cpSSR loci could explain in
part why the very different scenarios of divergence we tested were not more strongly separated. While being aware of these limitations, the estimation of the demographic parameters remained consistent across the two repeated analyses and dated the origin of the two Pyrenean lineages to the middle Pleistocene, 0.7 M years ago, when Quaternary glacial cycle duration shifted from lasting an average of 41,000 years to an average of 100,000 years (Pisias & Moore, 1981).

Pyrenean lineages may have diverged from a single Alpine ancestor but also independently from two Alpine ancestral lines, as previously suggested for the evergreen shrub *Rhododendron ferrugineum* (Charrier et al., 2014). Silver fir was present between Eastern France (Jura) and the Pyrenees during the Pleistocene since pollen was detected in central France in the southern Loire river basin (Gauthier et al., 2005; Visset et al., 2003 and see Figure S5.1 in Appendix S5), but the lack of data from southwestern France makes it impossible to draw further conclusions.

Our results support the multiple refugia scenario, defined as "Refugia within Refugia" by Gómez and Lunt (2006) to explain biodiversity in the Iberian Peninsula. Even though we could not determine their exact location, the signal of a progressive colonization from east to west presupposes two glacial refugia: one for the eastern lineage, in the Eastern Pyrenees where the genetic diversity is highest, and a second more hypothetical refugium in the Central Pyrenees or in the pre-Pyrenees chain as suggested for the perennial plant *Borderea pyrenaica* Miégev. (Segarra-Moragues et al., 2007). The hypothesis of a possible refuge in the Western Pyrenees must be taken with caution, since it is not supported by fossil data. In addition, migrations during the Quaternary period led to secondary contact zones where the Eastern and Western Pyrenean lineages merged into admixed populations in the Central Pyrenees.

## 4.2 | Distribution and environmental drivers of genetic diversity of silver fir in the Pyrenees

Our study shows that genetic diversity of silver fir decreases from east to west, Alpine populations being more diverse than populations in the Northern Pyrenees. This confirms the observations made on coarser sampling grids (Fady & Conord, 2010; Sancho-Knapik et al., 2014; Vendramin et al., 1999). A decrease in genetic diversity between the Alps and the Pyrenees is also observed in several plants (Charrier et al., 2014; Despres et al., 2002) and animals (Caizergues et al., 2003). It is explained by harsher (colder and drier) climatic conditions on resident populations during the LGM in the west than in the east and/or long-distance colonization during the Holocene or previous interglacials from favorable refugial habitats (Conord et al., 2012). This continental trend remains true at the Northern Pyrenees range scale for *Abies alba* and is probably explained by isolation-by-distance from east to west, as suggested by autocorrelation analyses, EEMS, and the historical demography signal (expansion in eastern and contraction in western populations). However, we also detected a contact zone allowing gene flow between genetic clusters in the Northern Central Pyrenees which could explain the

absence of phylogeographic signal (absence of difference between $F_{ST}$ and $R_{ST}$) and the deviation from a stepping-stone model in EEMS.

Silver fir displays a very clear genetic structure along the range of the Pyrenees as already observed in other plant and animal species (see references in the introduction), demonstrating the roles of east-west land barriers in models of simple processes of differentiation, colonization, and mixing (Wallis et al., 2016). The geographical position of the Pyrenees between the Atlantic Ocean and the Mediterranean Sea, as well as the contrasted topography along the range (lower altitude in the west), resulted in the current contrasted regional climates but also in regional contrasts in glaciation patterns during the Quaternary periods (Calvet, 2004). This can explain why longitude is the environmental variable that explains most of the genetic diversity and structure in silver fir populations, dividing the Northern Pyrenees into two large blocks: the Western and the Eastern Pyrenees.

Additional factors can explain patterns of genetic diversity and in particular the distribution of ancestry coefficients. The relationships observed between the degree of membership to genetic clusters and the temperature during the wettest quarter at interglacial periods suggest a significant role of the past climate in the distribution of genetic diversity. Temperature during the wettest quarter is currently four degrees higher in the Eastern than in the Western Pyrenees (8°C vs. 4°C, Table S4.2), and it was five degrees higher during the last interglacial period, although mean temperature was much colder then, and probably more selective for trees (6.3°C vs. 1.2°C, Table S4.2). The strength of such environmental contrasts could have contributed to the maintenance of a strong differentiation between Pyrenean lineages. Whether long-lasting differentiation in different environments has led to differential selection remains to be tested, but indications that regional environmental gradients lead to diversifying selection can be found in *Abies alba* studies from both the Spanish Pyrenees (Matías et al., 2016) and the Alps (Brousseau et al., 2016; Csilléry et al., 2020; Roschanski et al., 2016; Vitasse et al., 2009).

## 4.3 | Conservation perspectives

The conservation of forests and forest habitats is part of the global strategy to conserve biodiversity, under the legal responsibility of each nation following the framework of the CBD (https://www.cbd.int/). Since 2016, recommendation WCC-2016-Rec-104 of the International Union for Conservation of Nature (IUCN) "Integrating autochthonous forest genetic diversity into protected area conservation objectives" recommends that geographic entities that have *in situ* genetic diversity conservation as their main purpose should be recognized as protected areas (IUCN, 2016). Prioritization is essential to minimize biodiversity loss (Brooks et al., 2006), and in Europe, it is done under the guidance of EUFORGEN, the European Forest Genetic Resources Programme (http://www.euforgen.org/). Identifying which forests deserve the status of gene conservation unit (GCU) warrants the identification of well-defined evolutionary
lineages and the selection of forests where management is compat-
ible with the occurrence of genetic processes governed by natural
selection (Koskela et al., 2013; Lefèvre et al., 2020).

Major efforts are made in Europe for the conservation of the
genetic resources of silver fir, with a total of 263 GCUs (http://www.
euforgen.org/species/abies-alba/), located in seven environmen-
tal zones (http://portal.eufgis.org), following the classification of
Metzger et al. (2013) modified by de Vries et al. (2015). Seven GCUs
currently represent the genetic diversity of the French Pyrenees, all
analyzed in this study.

Four GCUs (the forests of Hares, Fanges, Arques, and Canigou) from the Eastern Pyrenees belong to the eastern lineage with a membership coefficient of 0.77 to 0.88. One GCU (the forest of Saint Lary) from the Central Pyrenees is admixed between the two lineages. The remaining two GCUs (the forests of Hautacam and Barousse), also from the Central Pyrenees, belong to the western lineage with a membership coefficient of 0.77 to 0.80. While the Eastern Pyrenees (in the geographical sense) and the eastern lineage are properly covered, with GCUs from forests at both low and high elevation, the Western Pyrenees (in the geographical sense) are not covered and the western lineage is poorly covered, with low elevation populations and the western margin of the distribution lacking. Geographical edges of distribution ranges demonstrate original local adaptations (Hampe & Petit, 2005; Parisod & Joost, 2010). Despite its low genetic diversity, the western margin of the geographical distribution of silver fir in the Pyrenees may harbor local adaptations that need further evaluation and protection.

Our study shows an absence of strict overlap between genetic
clustering and geographic (topography and climate) structure in
the Pyrenees and suggests that geography alone cannot be used
as proxy for sampling for conservation. As genetic differentiation
can occur over short distance along elevation gradients in silver fir
(Csilléry et al., 2020; Roschanski et al., 2016), we advise that GCUs
should be as large as possible (typically larger than 100 ha) to make
gene flow possible across different elevations. In this perspective,
including populations of admixed origin as GCUs (such as the Central
Pyrenees Saint Lary forest) can significantly contribute to safe-
guarding evolutionary processes important for management under
a changing climate, both for local resilience and as a resource for
plantation forests (Hamilton & Miller, 2016; vonVoldt et al., 2017).

A remote and rather inaccessible area until recently, the
Pyrenees represent one of the last refuges in Europe for some of
its most spectacular and endangered mammals and birds, including
the brown bear (Ursus arctos arctos L.), the Pyrenean ibex (Capra
pyrenaica pyrenaica Schinz), and the bearded vulture (Gypaetus bar-
batus L.). Many strategies are developed to protect these species in
the Pyrenees such as the strict protection of old-growth silver fir,
beech, and pine forests (Sabatini et al., 2018) under the Biodiversity
Strategy of the European Union, the Natura 2000 network, national
forest strategies, and regional initiatives. Creating and managing
GCUs for silver fir and other forest tree species, maintaining their
resilience, will further help conserve these and other mountain plant
and animal species present in the Pyrenees.

5 | CONCLUSION

Although much emphasis has been placed on the Holocene period, this study suggests that mid-Pleistocene events have also durably shaped the current spatial structure of genetic diversity in natural forest tree populations. East-west oriented mountain ranges, perpendicular to the main climate change directions of the Pleistocene, are suitable model systems to study the phylogeography and the demographic fate of mid- to high-elevation mountain species, possibly emerging from different refugia. Geographical variables range-wide remain insufficient and potentially misleading to capture patterns of genetic diversity for sampling with conservation purposes. Genetic studies are needed at the scale of a whole mountain range to allow the identification of all local particularities (e.g., in response to past environmental changes) for building and/or refining conservation strategies in mountain forest tree species.

ACKNOWLEDGMENTS

We thank Pierre Gonin and Grégory Sajdak (CNPF Toulouse, France) for coordinating the collection of material in the Pyrenees and Massif Central, and Norbert Turion and Olivier Gilg (INRAE UEFM Avignon, France), for collecting the material in the Southern Alps and in Corsica. We also thank Bénédicte Le Guerroué (ONF Orléans, France) for DNA extraction and cpSSR genotyping and Anne Roig and Matthieu Lingrand (INRAE URFM Avignon, France) for nuSSR genotyping. This study received funding from the Parc naturel régional des Pyrénées Ariégeoises and the Eranet Biodiversa project TipTree (Project ANR 12 EBID 000 3). This publication is also part of a project that has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 676876 (GenTree). We also acknowledge the support of the French National Commission on Forest Genetic Resources (CRGF, https://agriculture.gouv.fr/la-politique-nationale-de-conservation-des-ressources-genetiques-forestieres).

CONFLICT OF INTEREST

The authors have no competing interests to declare.

AUTHOR CONTRIBUTION

**Caroline Scotti-Saintagne:** Data curation (lead); Formal analysis (lead); Methodology (lead); Writing-original draft (equal); Writing-review & editing (equal). **Thomas Boivin:** Formal analysis (equal); Methodology (equal); Writing-review & editing (equal). **Marie Suez:** Formal analysis (supporting); Methodology (supporting). **Brigitte Musch:** Conceptualization (equal); Funding acquisition (equal); Writing-review & editing (equal). **Ivan Scotti:** Methodology (equal); Writing-review & editing (equal). **Bruno Fady:** Conceptualization (lead); Funding acquisition (lead); Writing-original draft (equal); Writing-review & editing (equal).

DATA AVAILABILITY STATEMENT

Genotypic data are archived and publicly accessible open access at:
https://doi.org/10.15454/HVPAON.
ORCID

Bruno Fady

<figure><img src="image_10.png" /></figure>

https://orcid.org/0000-0003-2379-7617

REFERENCES

Archer, F. I., Adams, P. E., & Schneiders, B. B. (2017). Stratag: An r package for manipulating, summarizing and analysing population genetic data. *Molecular Ecology Resources*, 17, 5-11.

Avise, J. C. (2000). *Phylogeography: The history and formation of species*. Harvard University Press.

Awad, L., Fady, B., Khater, C., Roig, A., & Cheddadi, R. (2014). Genetic structure and diversity of the endangered fir tree of Lebanon (*Abies cilicica* Carr.): Implications for conservation. *PLoS One*, 9(2), e90086. https://doi.org/10.1371/journal.pone.0090086

Bagnoli, F., Tsuda, Y., Fineschi, S., Bruschi, P., Magri, D., Zhelev, P., Paule, L., Simeone, M. C., González-Martínez, S. C., & Vendramin, G. G. (2016). Combining molecular and fossil data to infer demographic history of *Quercus cerris*: Insights on European eastern glacial refugia. *Journal of Biogeography*, 43, 679-690.

Bagnoli, F., Vendramin, G. G., Buonamici, A., Doulis, A. G., Gonzalez-Martinez, S. C., La Porta, N., Magri, D., Raddi, P., Sebastiani, F., Fineschi, S., & Fineschi, S. (2009). Is Cupressus sempervirens native in Italy? An answer from genetic and palaeobotanical data. *Molecular Ecology*, 18(10), 2276-2286.

Barnosky, C. W. (1987). Response of vegetation to climatic changes of different duration in the late Neogene. *Trends in Ecology & Evolution*, 2(8), 247-250. https://doi.org/10.1016/0169-5347(87)90007-3

Barthe, S., Binelli, G., Hérault, B., Scotti-Saintagne, C., Sabatier, D., & Scotti, I. (2017). Tropical rainforests that persisted: Inferences from the Quaternary demographic history of eight tree species in the Guiana shield. *Molecular Ecology*, 26, 1161-1174. https://doi.org/10.1111/mec.13949

Beaumont, M. A. (2010). Approximate Bayesian computation in evolution and ecology. *Annual Review of Ecology, Evolution and Systematics*, 41, 379-405. https://doi.org/10.1146/annurev-ecols ys-102209-144621

Belletti, P., Ferrazzini, D., Ducci, F., De Rogatis, A., Mucciarelli, M., & Mucciarelli, M. (2017). Genetic diversity of Italian populations of *Abies alba*. *Dendrobiology*, 77, 147-159. https://doi.org/10.12657/ denbio.077012

Bidegaray-Batista, L., Sánchez-gracia, A., Santulli, G., Maiorano, L., Guisan, A., Vogler, A. P., & Arnedo, M. A. (2016). Imprints of multiple glacial refugia in the Pyrenees revealed by phylogeography and palaodistribution modelling of an endemic spider. *Molecular Ecology*, 25, 2046-2064. https://doi.org/10.1111/mec.13585

Brooks, T. M., Mittermeier, R. A., da Fonseca, G. A. B., Gerlach, J., Hoffmann, M., Lamoreux, J. F., Mittermeier, C. G., Pilgrim, J. D., & Rodrigues, A. S. L. (2006). Global biodiversity conservation priorities. *Science*, 313, 58-61. https://doi.org/10.1126/science.1127609

Brousseau, L., Postolache, D., Lascoux, M., Drouzas, A. D., Källman, T., Leonarduzzi, C., Liepelt, S., Piotti, A., Popescu, F., Roschanski, A. M., Zhelev, P., Fady, B., & Vendramin, G. G. (2016). Local adaptation in European firs assessed through extensive sampling across altitudinal gradients in southern Europe. *PLoS One*, 11(7), 1-26. https://doi.org/10.1371/journal.pone.0158216

Cailleret, M., Nourtier, M., Amm, A., Durand-Gillmann, M., & Davi, H. (2014). Drought-induced decline and mortality of silver fir differ among three sites in Southern France. *Annals of Forest Science*, 71(6), 643-657. https://doi.org/10.1007/s13595-013-0265-0

Caizergues, A., Bernard Laurent, A., Brenot, J. F., Ellison, L., & Rasplus, J. Y. (2003). Population genetic structure of rock ptarmigan *Lagopus mutus* in Northern and Western Europe. *Molecular Ecology*, 12, 2267-2274.

Calvet, M. (2004). The Quaternary glaciation of the Pyrenees. *Developments in Quaternary Sciences*, 2(1), 119-128.

Camarero, J. J., Bigler, C., Linares, J. C., & Gil-Pelegrín, E. (2011). Synergistic effects of past historical logging and drought on the decline of Pyrenean silver fir forests. *Forest Ecology and Management*, 262, 759-769. https://doi.org/10.1016/j.foreco.2011.05.009

Charrier, O., Dupont, P., Pornon, A., & Escaravage, N. (2014). Microsatellite marker analysis reveals the complex phylogeographic history of *Rhododendron ferrugineum* (Ericaceae) in the Pyrenees. *PLoS One*, 9(3), e92976.

Cheddadi, R., Birks, H. J. B., Tarroso, P., Liepelt, S., Gömöry, D., Dullinger, S., Meier, E. S., Hülber, K., Maiorano, L., & Laborde, H. (2014). Revisiting tree-migration rates: *Abies alba* (Mill.), a case study. *Vegetation History and Archaeobotany*, 23, 113-122. https://doi.org/10.1007/s00334-013-0404-4

Chhatre, V. E., & Emerson, K. J. (2017). StrAuto: Automation and parallelization of STRUCTURE analysis. *BMC Bioinformatics*, 18(1), 192. https://doi.org/10.1186/s12859-017-1593-0

Chikhi, L., Sousa, V. C., Luisi, P., Goossens, B., & Beaumont, M. A. (2010). The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes. *Genetics*, 186(3), 983-995. https://doi.org/10.1534/genetics.110.118661

Conord, C., Gurevitch, J., & Fady, B. (2012). Large-scale longitudinal gradients of genetic diversity: A meta-analysis across six phyla in the Mediterranean basin. *Ecology and Evolution*, 2(10), 2600-2614. https://doi.org/10.1002/ece3.350

Cornuet, J.-M., Ravigné, V., & Estoup, A. (2010). Inference on population history and model checking using DNA sequence and microsatellite data with the software DIYABC (v1.0). *BMC Bioinformatics*, 11, 401.

Cornuet, J.-M., Santos, F., Beaumont, M. A., Robert, C. P., Marin, J.-M., Balding, D. J., Guillemaud, T., & Estoup, A. (2008). Inferring population history with DIY ABC: A user-friendly approach to approximate Bayesian computation. *Bioinformatics*, 24(23), 2713-2719. https://doi.org/10.1093/bioinformatics/btn514

Crawford, R. M. M. (2008). *Plants at the margin. Ecological limits and climate change* (494 pp.). Cambridge University Press.

Cremer, E., Liepelt, S., Sebastiani, F., Buonamici, A., Michalczyk, I. M., Ziegenhagen, B., & Vendramin, G. G. (2006). Identification and characterization of nuclear microsatellite loci in *Abies alba* Mill. *Molecular Ecology Notes*, 6, 374-376. https://doi.org/10.1111/j.1471-8286.2005.01238.x

Csilléry, K., Buchmann, N., & Fady, B. (2020). Adaptation to drought is coupled with slow growth, but independent from phenology in maternal silver fir (*Abies alba* Mill.) populations. *Evolutionary Applications*, 13, 2357-2376.

de Vries, S. M. G., Alan, M., Bozzano, M., Burianek, V., Collin, E., Cottrell, J., Ivankovic, M., Kelleher, C. T., Koskela, J., Rotach, P., Pietto, L., & Yrjänä, L. (2015). *Pan-European strategy for genetic conservation of forest trees and establishment of a core network of dynamic conservation units*. Rome, Italy: European Forest Genetic Resources Programme (EUFORGEN), Bioversity International.

Demesure, B., Sodzi, N., & Petit, R. J. (1995). A set of universal primers for amplification of polymorphic non-coding regions of mitochondrial and chloroplast DNA in plants. *Molecular Ecology*, 4, 129-134. https://doi.org/10.1111/j.1365-294X.1995.tb00201.x

Dempster, A., Laird, N., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via EM Algorithm. *Journal of the Royal Statistical Society Series B*, 39, 1-38.

Despres, L., Loriot, S., & Gaudeul, M. (2002). Geographic pattern of genetic variation in the European globeflower *Trollius europaeus* L. (Ranunculaceae) inferred from amplified fragment length polymorphism markers. *Molecular Ecology*, 11, 2337-2347. https://doi.org/10.1046/j.1365-294X.2002.01618.x

Duriez, O., Sachet, J. M., Menonie, E., Pidancier, N., Miquel, C., & Taberlet, P. (2007). Phylogeography of the Capercaillie in Eurasia: What is the conservation status in the Pyrenees and Cantabrian
Mounts? *Conservation Genetics*, 8, 513–526. https://doi.org/10.1007/s10592-006-9165-2

Earl, D. A., & vonHoldt, B. M. (2012). Structure harvester: A website and program for visualizing STRUCTURE output and implementing the Evanno method. *Conservation Genetics Resources*, 4, 359–361. https://doi.org/10.1007/s12686-011-9548-7

Ejarque, A., Julià, R., Riera, S., Palet, J. M., Orengo, H. A., Miras, Y., & Gascón, C. (2009). Tracing the history of highland human management in the eastern Pre-Pyrenees: An interdisciplinary palaeoenvironmental study at the Pradell fen, Spain. *The Holocene*, 9(8), 1241–1255. https://doi.org/10.1177/0959683609345084

Estoup, A., Jarne, P., & Cornuet, J.-M. (2002). Homoplasy and mutation model at microsatellite loci and their consequences for population genetics analysis. *Molecular Ecology*, 11, 1591–1604. https://doi.org/10.1046/j.1365-294X.2002.01576.x

Evanno, G., Regnaut, S., & Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. *Molecular Ecology*, 14(8), 2611–2620. https://doi.org/10.1111/j.1365-294X.2005.02553.x

Fady, B., & Conord, C. (2010). Macroecological patterns of species and genetic diversity in vascular plants of the Mediterranean basin. *Diversity and Distributions*, 16(1), 53–64. https://doi.org/10.1111/j.1472-4642.2009.00621.x

Fady, B., Cottrell, J., Ackzell, L., Alia, R., Muys, B., Prada, A., & González-Martínez, S. C. (2016). Forests and global change: What can genetics contribute to the major forest management and policy challenges of the twenty-first century? *Regional Environmental Change*, 16(4), 927–939. https://doi.org/10.1007/s10113-015-0843-9

FAO. (2015). Global Forest Resources Assessment. Terms and definitions. In Food and agriculture organization of the united nations (Ed.), *Forest resources Assessment Working Paper*, Vol. 180 (pp. 36). FAO.

Felsenstein, J. (1985). Confidence limits on phylogenies: An approach using the bootstrap. *Evolution*, 39, 783–791.

Fernandez-Fournier, P., Lewthwaite, J. M. M., & Mooers, A. Ø. (2021). Do we need to identify adaptive genetic variation when prioritizing populations for conservation? *Conservation Genetics*, 22, 205–216. https://doi.org/10.1007/s10592-020-01327-w

Fox, J., & Weisberg, S. (2018). *An R companion to applied regression*, 3rd ed. Sage.

Galop, D., & Jalut, G. (1994). Differential human impact and vegetation history in two adjacent Pyrenean valleys in the Ariège basin, southern France, from 3000 BP to the present. *Vegetation History and Archaeobotany*, 3(4), 225–244.

Gauthier, G., Vannière, B., & Martineau, R. (2005). Histoire des feux et pratiques agraires du Néolithique à l'âge du Fer en région Centre: Implications territoriales démographiques et environnementales. *Gallia Préhistoire*, 47, 167–186.

Gómez, A., & Lunt, D. H. (2006). Refugia within refugia: Patterns of phylogeographic concordance in the Iberian Peninsula. In S. Weiss, & N. Ferrand (Eds.), *Phylogeography in Southern European Refugia: Evolutionary Perspectives on the Origins and Conservation of European Biodiversity* (pp. 155–182). Springer.

Hamilton, J. A., & Miller, J. M. (2016). Adaptive introgression as a resource for management and genetic conservation in a changing climate. *Conservation Biology*, 30(1), 33–41. https://doi.org/10.1111/cobi.12574

Hampe, A., & Petit, R. J. (2005). Conserving biodiversity under climate change: The rear edge matters. *Ecology Letters*, 8(5), 461–467. https://doi.org/10.1111/j.1461-0248.2005.00739.x

Hansen, O. K., Vendramin, G. G., Sebastiani, F., & Edwards, K. J. (2005). Development of microsatellite markers in *Abies nordmanniana* (Stev.) Spach and cross-species amplification in the *Abies* genus. *Molecular Ecology Notes*, 5(4), 784–787.

Hardy, O. J., Charbonnel, N., Fréville, H., & Heuertz, M. (2003). Microsatellite allele sizes: A simple test to assess their significance on genetic differentiation. *Genetics*, 163(4), 1467–1482. https://doi.org/10.1093/genetics/163.4.1467

Hardy, O. J., & Vekemans, X. (2002). spagedi: A versatile computer program to analyse spatial genetic structure at the individual or population levels. *Molecular Ecology Notes*, 2(4), 618–620. https://doi.org/10.1046/j.1471-8286.2002.00305.x

Hewitt, G. M. (1999). Post-glacial re-colonization of European biota. *Biological Journal of the Linnean Society*, 68, 87–112. https://doi.org/10.1111/j.1095-8312.1999.tb01160.x

Hijmans, R. J., Cameron, S., Parra, J. L., Jones, P. G., & Jarvis, A. (2005). Very high-resolution interpolated climate surfaces for global land areas. *International Journal of Climatology*, 25, 1965–1978. https://doi.org/10.1002/joc.1276

Hirao, A. S., Watanabe, M., Tsuyuzaki, S., Shimono, A., Li, X., Masuzawa, T., & Wada, N. (2017). Genetic diversity within populations of an arctic-alpine species declines with decreasing latitude across the Northern Hemisphere. *Journal of Biogeography*, 44, 2740–2751. https://doi.org/10.1111/jbi.13085

Hoban, S., Bruford, M., D'Urban Jackson, J., Lopes-Fernandes, M., Heuertz, M., Hohenlohe, P. A., Paz-Vinas, I., Sjögren-Gulve, P., Segelbacher, G., Vernesi, C., Aitken, S., Bertola, L. D., Bloomer, P., Breed, M., Rodríguez-Correa, H., Funk, W. C., Grueber, C. E., Hunter, M. E., Jaffe, R., ... Laikre, L. (2020). Genetic diversity targets and indicators in the CBD post-2020 Global Biodiversity Framework must be improved. *Biological Conservation*, 248, 108654. https://doi.org/10.1016/j.biocon.2020.108654

Horn, A., Stauffer, C., Lieutier, F., & Kerdelhué, C. (2009). Complex postglacial history of the temperate bark beetle *Tomicus piniperda* L. (Coleoptera, Scolytinae). *Heredity*, 103, 238–247. https://doi.org/10.1038/hdy.2009.48

IUCN. (2016). *IUCN resolutions, recommendations and other decisions* (p. 106). IUCN.

Jombart, T. (2008). Adegenet: A R package for the multivariate analysis of genetic markers. *Bioinformatics*, 24(11), 1403–1405. https://doi.org/10.1093/bioinformatics/btn129

Kimura, M., & Weiss, G. H. (1964). The stepping stone model of population structure and the decrease of genetic correlation with distance. *Genetics*, 49(4), 561–576. https://doi.org/10.1093/genetics/49.4.561

Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N., & Mayrose, I. (2015). Clumpak: A program for identifying clustering modes and packaging population structure inferences across K. *Molecular Ecology Resources*, 15, 1179–1191.

Koskela, J., Lefèvre, F., Schueler, S., Kraigher, H., Olrik, D. C., Hubert, J., Longauer, R., Bozzano, M., Yrjänä, L., Alizoti, P., Rotach, P., Vietto, L., Bordács, S., Myking, T., Eysteinsson, T., Souvannavong, O., Fady, B., De Cuyper, B., Heinze, B., ... Ditlevsen, B. (2013). Translating conservation genetics into management: Pan-European minimum requirements for dynamic conservation units of forest tree genetic diversity. *Biological Conservation*, 157, 39–49. https://doi.org/10.1016/j.biocon.2012.07.023

Kropf, M., Comes, H. P., & Kadereit, J. W. (2006). Long-distance dispersal vs vicariance: the origin and genetic diversity of alpine plants in the Spanish Sierra Nevada. *New Phytologist*, 172, 169–184. https://doi.org/10.1111/j.1469-8137.2006.01795.x

Lauga, B., Malaval, S., Largier, G., & Regnault-Roger, C. (2009). Two lineages of *Trifolium alpinum* (Fabaceae) in the Pyrenees: Evidence from random amplified polymorphic DNA (RAPD) markers. *Acta Botanica Gallica*, 156(3), 317–330.

Lefèvre, F., Alia, R., Bakkebø Fjellstad, K., Graudal, L., Oggioni, S. D., Rusanen, M., Vendramin, G. G., & Bozzano, M. (2020). Dynamic conservation and utilization of forest tree genetic resources: Indicators for in situ and ex situ genetic conservation and
*forest reproductive material.* European Forest Genetic Resources Programme (EUFORGEN), European Forest Institute.

Liberal, I. M., Burrus, M., Suchet, C., Thébaud, C., & Vargas, P. (2014). The evolutionary history of *Antirrhinumin* the Pyrenees inferred from phylogeographic analyses. *BMC Evolutionary Biology*, 14, 146. https://doi.org/10.1186/1471-2148-14-146

Liepelt, S., Cheddadi, R., de Beaulieu, J.-L., Fady, B., Gömöry, D., Hussendörfer, E., Konnert, M., Litt, T., Longauer, R., Terhörne-Berson, R., & Ziegenhagen, B. (2009). Postglacial range expansion and its genetic imprints in *Abies alba* (Mill.) - A synthesis from palaeobotanic and genetic data. *Review of Palaeobotany and Palynology*, 1-2(153), 139-149. https://doi.org/10.1016/j.revpalbo.2008.07.007

Liepelt, S., Kuhlenkamp, V., Anzidei, M., Vendramin, G. G., & Ziegenhagen, B. (2001). Pitfalls in determining size homoplasy of microsatellite loci. *Molecular Ecology Notes*, 1, 332-335. https://doi.org/10.1046/j.1471-8278.2001.00085.x

Magri, D., Di Rita, F., Aranbarri, J., Fletcher, W., & González-Sampériz, P. (2017). Quaternary disappearance of tree taxa from Southern Europe: Timing and trends. *Quaternary Science Reviews*, 163, 23-55. https://doi.org/10.1016/j.quascirev.2017.02.014

Magri, D., Vendramin, G. G., Comps, B., Dupanloup, I., Geburek, T., Gömöry, D., Latałowa, M., Litt, T., Paule, L., Roure, J. M., Tantau, I., Van Der Knaap, W. O., Petit, R. J., & De Beaulieu, J.-L. (2006). A new scenario for the Quaternary history of European beech populations: Palaeobotanical evidence and genetic consequences. *New Phytologist*, 171, 199-221. https://doi.org/10.1111/j.1469-8137.2006.01740.x

Matías, L., Gonzalez-Díaz, P., Quero, J. L., Camarero, J. J., Lloret, F., & Jump, A. S. (2016). Role of geographical provenance in the response of silver fir seedlings to experimental warming and drought. *Tree Physiology*, 36(10), 1236-1246. https://doi.org/10.1093/treephys/tpw049

Mauri, A., de Rigo, D., & Caudullo, G. (2016). *Abies alba* in Europe: Distribution, habitat, usage and threats. In J. San Miguel Ayanz, D. de Rigo, G. Caudullo, T. Houston Durrant, & A. Mauri (Eds.), *European Atlas of Forest Tree Species*. Publ. Off. EU, Luxembourg, pp. E01493b+.

Metzger, M. J., Bunce, R. G. H., Jongman, R. H. G., Sayre, R., Trabucco, A., & Zomer, R. (2013). A high-resolution bioclimate map of the world: A unifying framework for global biodiversity research and monitoring. *Global Ecology and Biogeography*, 22(5), 630-638. https://doi.org/10.1111/geb.12022

Moritz, C., & Faith, D. P. (1998). Comparative phylogeography and the identification of genetically divergent areas for conservation. *Molecular Ecology*, 7(4), 419-429. https://doi.org/10.1046/j.1365-294x.1998.00317.x

Mosca, E., Eckert, A. J., Liechty, J. D., Wegrzyn, J. L., La Porta, N., Vendramin, G. G., & Neale, D. B. (2012). Contrasting patterns of nucleotide diversity for four conifers of Alpine European forests. *Evolutionary Applications*, 5(7), 762-775. https://doi.org/10.1111/j.1752-4571.2012.00256.x

Nagy, L., Grabherr, G., Körner, C., & Thompson, D. B. A. (2003). In L. Nagy, G. Grabherr, C. Körner, & D. B. A. Thompson (Eds.), *Alpine biodiversity in Europe*, Vol. 167. Springer.

Nei, M. (1978). Estimation of average heterozygosity and genetic distance from a small number of individuals. *Genetics*, 89(3), 583-590. https://doi.org/10.1093/genetics/89.3.583

Nielsen, R., Tarpy, D. R., & Reeve, H. K. (2003). Estimating effective pattern number in social insects and the effective number of alleles in a population. *Molecular Ecology*, 12(11), 3157-3164. https://doi.org/10.1046/j.1365-294X.2003.01994.x

Ninot, J. M., Carillo, M., Font, X., Carreras, J., Ferré, A., Masalles, R. M., Soriano, I., & Vigo, J. (2007). Altitude zonation in the Pyrenees. A Geobotanic Interpretation. *Phytocoenologia*, 37(3-4), 371-398.

Ozenda, P. (1985). *La végétation de la chaîne alpine dans l'espace montagnard européen*. Masson.

Parisod, C., & Joost, S. (2010). Divergent selection in trailing- versus leading-edge populations of *Biscutella laevigata*. *Annals of Botany*, 105(4), 655-660. https://doi.org/10.1093/aob/mcq014

Pérez-Sanz, A., González-Sampériz, P., Moreno, A., Valero-Garcés, B., Gil-Romera, G., Rieradevall, M., Tarrats, P., Lasheras-Álvarez, L., Morellón, M., Belmonte, A., Sancho, C., Sevilla-Callejo, M., & Navas, A. (2013). Holocene climate variability, vegetation dynamics and fire regime in the central Pyrenees: The Basa de la Mora sequence (NE Spain). *Quaternary Science Reviews*, 73, 149-169. https://doi.org/10.1016/j.quascirev.2013.05.010

Petit, R. J., Csaikl, U. M., Bordács, S., Burg, K., Coart, E., Cottrell, J., van Dam, B., Deans, J. D., Dumolin-Lapègue, S., Fineschi, S., Finkeldey, R., Gillies, A., Glaz, I., Goicoechea, P. G., Jensen, J. S., König, A. O., Lowe, A. J., Madsen, S. F., Mátyás, G., ... Kremer, A. (2002). Chloroplast DNA variation in European white oaks: Phylogeography and patterns of diversity based on data from over 2600 populations. *Forest Ecology and Management*, 156(1-3), 5-26. https://doi.org/10.1016/S0378-1127(01)00645-4

Petit, R. J., Hampe, A., & Cheddadi, R. (2005). Climate changes and tree phylogeography in the Mediterranean. *Taxon*, 54(4), 877-885. https://doi.org/10.2307/25065474

Petkova, D., Novembre, J., & Stephens, M. (2015). Visualizing spatial population structure with estimated effective migration surfaces. *Nature Genetics*, 48, 94-100. https://doi.org/10.1038/ng.3464

Piotti, A., Leonarduzzi, C., Postolache, D., Bagnoli, F., Spanu, I., Brousseau, L., Urbinati, C., Leonardi, S., & Vendramin, G. G. (2017). Unexpected scenarios from Mediterranean refugial areas: Disentangling complex demographic dynamics along the Apennine distribution of silver fir. *Journal of Biogeography*, 44(7), 1547-1558. https://doi.org/10.1111/jbi.13011

Pisias, N. G., & Moore, T. C. (1981). The evolution of Pleistocene climate: A time series approach. *Earth and Planetary Science Letters*, 52(2), 450-458.

Pritchard, J. K., Stephens, M., & Donnelly, P. (2000). Inference of population structure using multilocus genotype data. *Genetics*, 155, 945-959. https://doi.org/10.1093/genetics/155.2.945

R Core Team (2020). *R: A language and environment for statistical computing*. R Foundation for Statistical Computing.

Riddle, B. R. (2016). Comparative phylogeography clarifies the complexity and problems of continental distribution that drove AR Wallace to favor islands. *Proceedings of the National Academy of Sciences*, 113(29), 7970-7977. https://doi.org/10.1073/pnas.1601072113

Roschanski, A. M., Csilléry, K., Liepelt, S., Oddou-Muratorio, S., Ziegenhagen, B., Huard, F., Ullrich, K. K., Postolache, D., Vendramin, G. G., & Fady, B. (2016). Evidence of divergent selection for drought and cold tolerance at landscape and local scales in *Abies alba* Mill. in the French Mediterranean Alps. *Molecular Ecology*, 25, 776-794. https://doi.org/10.1111/mec.13516

Rousselet, J., Zhao, R., Argal, D., Simonato, M., Battisti, A., Roques, A., & Kerdelhué, C. (2010). The role of topography in structuring the demographic history of the pine processionary moth, *Thaumetopoea pityocampa* (Lepidoptera: Notodontidae). *Journal of Biogeography*, 37, 1478-1490.

Rousset, F. (2008). *Genepop'007: A complete reimplementation of the Genepop software for Windows and Linux*. *Molecular Ecology Resources*, 8, 103-106. https://doi.org/10.1111/j.1471-8286.2007.01931.x

RStudio Team. (2019). *RStudio: Integrated development for R*. RStudio Inc.

Rundel, P. W. (2019). A neogene heritage: Conifer distributions and endemism in mediterranean-climate ecosystems. *Frontiers in Ecology and Evolution*, 7, article 364. https://doi.org/10.3389/fevo.2019.00364

Sabatini, F. M., Burrascano, S., Keeton, W. S., Levers, C., Lindner, M., Pötzschner, F., Verkerk, P. J., Bauhus, J., Buchwald, E., Chaskovsky, O., Debaive, N., Horváth, F., Garbarino, M., Grigoriadis, N., Lombardi, F., Marques Duarte, I., Meyer, P., Midteng, R., Mikac, S., ... Kuemmerle,
T. (2018). Where are Europe's last primary forests? *Diversity and Distributions*, 24, 1426-1439. https://doi.org/10.1111/ddi.12778

Saitou, N., & Nei, M. (1987). The neighbor-joining method: A new method for reconstructing phylogenetic trees. *Molecular Biology and Evolution*, 4, 406-425.

Sancho-Knapik, D., Peguero-Pina, J., Cremer, E., Camarero, J., Fernández-Cancio, Á., Ibarra, N., Konnert, M., & Gil-Pelegrín, E. (2014). Genetic and environmental characterization of *Abies alba* Mill. Populations at its western rear edge. *Pirineos. Revista De Ecología De Montaña*, 169, e007.

Schmitt, T., Hewitt, G. M., & Muller, P. (2006). Disjunct distributions during glacial and interglacial periods in mountain butterflies: *Erebia epiphron* as an example. *Journal of Evolutionary Biology*, 19, 108-113. https://doi.org/10.1111/j.1420-9101.2005.00980.x

Scotti-Saintagne, C., Giovannelli, G., Scotti, I., Roig, A., Spanu, I., Vendramin, G. G., Guibal, F., & Fady, B. (2019). Recent, Late Pleistocene fragmentation shaped the phylogeographic structure of the European black pine (*Pinus nigra Arnold*). *Tree Genetics & Genomes*, 15, 76. https://doi.org/10.1007/s11295-019-1381-2

Secretariat of the Convention on Biological Diversity. (2020). Global Biodiversity Outlook 5. Montreal.

Segarra-Moragues, J. G., Palop-Esteban, M., González-Candelas, F., & Catalán, P. (2007). *Nunatak* survival vs. *tabula rasa* in the Central Pyrenees: A study on the endemic plant species *Borderea pyrenaica* (Dioscoreaceae). *Journal of Biogeography*, 34, 1893-1906.

Slatkin, M. (1995). A measure of population subdivision based on microsatellite allele frequencies. *Genetics*, 139, 1463.

Strimmer, K. (2008). fdrtool: A versatile R package for estimating local and tail area- based false discovery rates. *Bioinformatics*, 24, 1461-1462. https://doi.org/10.1093/bioinformatics/btn209

Szpiech, Z. A., Jakobsson, M., & Rosenberg, N. A. (2008). adze: A rarefaction approach for counting alleles private to combinations of populations. *Bioinformatics*, 24, 2498-2504. https://doi.org/10.1093/bioinformatics/btn478

Taberlet, P., Fumagalli, L., Wust-Saucy, A.-G., & Cosson, J. F. (1998). Comparative phylogeography and postglacial colonization routes in Europe. *Molecular Ecology*, 7(4), 453-464. https://doi.org/10.1046/j.1365-294x.1998.00289.x

Takezaki, N., Nei, M., & Tamura, K. (2010). POPTREE2: Software for constructing population trees from allele frequency data and computing other population statistics with Windows interface. *Molecular Biology and Evolution*, 27(4), 747-752. https://doi.org/10.1093/molbev/vmp312

Vendramin, G. G., Degen, B., Petit, R. J., Anzidei, M., Madaghiele, A., & Ziegenhagen, B. (1999). High level of variation at *Abies alba* chloroplast microsatellite loci in Europe. *Molecular Ecology*, 8(7), 1117-1126. https://doi.org/10.1046/j.1365-294x.1999.00666.x

Vendramin, G. G., & Ziegenhagen, B. (1997). Characterization and inheritance of polymorphic plastid microsatellites in *Abies*. *Genome*, 40, 857-864.

Villar, L., & Dendaletche, C. (1994). Pyrenees. France, Spain and Andorra. In S. D. Davis, V. H. Heywood, & A. C. Hamilton (Eds.), *Centres of plants diversity, a guide and strategy for their conservation*, Vol. 1 (pp. 61-64). Information Press.

Visset, L., Cyprien, A. L., Carcaud, N., & Ouguerram, A. (2003). De la fin du Mésolithique au Néolithique, l'évolution du paysage végétal dans le bassin de la Loire océanique et moyenne. In *Approche archéologique de l'environnement et de l'aménagement du territoire lignien*. Notes: Acte de Colloque. Fédération Archéologique du Loiret, 117-123.

Vitasse, Y., Delzon, S., Bresson, C. C., Michalet, R., & Kremer, A. (2009). Altitudinal differentiation in growth and phenology among populations of temperate-zone tree species growing in a common garden. *Canadian Journal of Forest Research*, 39, 1259-1269. https://doi.org/10.1139/X09-054

vonVoldt, B. M., Brzeski, K. E., Wilcove, D. S., & Rutledge, L. Y. (2017). Redefining the role of admixture and genomics in species conservation. *Conservation Letters*, 11, 1-6.

Wallis, G. P., Waters, J. M., Upton, P., & Craw, D. (2016). Transverse alpine speciation driven by glaciation. *Trends in Ecology & Evolution*, 31(12), 916-926. https://doi.org/10.1016/j.tree.2016.08.009

Wolf, H. (2003). *Euforgen Technical guidelines for genetic conservation and use of Silver fir* (*Abies alba*). Bioversity International.

Xiang, Q. P., Wei, R., Shao, Y. Z., Yang, Z. Y., Wang, X. Q., & Zhang, X. C. (2015). Phylogenetic relationships, possible ancient hybridization, and biogeographic history of *Abies* (Pinaceae) based on data from nuclear, plastid, and mitochondrial genomes. *Molecular Phylogenetics and Evolution*, 82, 1-14. https://doi.org/10.1016/j.ympev.2014.10.008

Zhivotovsky, L. A., Feldman, M. W., & Grisheckin, S. A. (1997). Biased mutations and microsatellite variation. *Molecular Biology and Evolution*, 14, 926-933. https://doi.org/10.1093/oxfordjournals.molbev.a025835

## SUPPORTING INFORMATION

Additional supporting information may be found online in the Supporting Information section.

**How to cite this article:** Scotti-Saintagne, C., Boivin, T., Suez, M., Musch, B., Scotti, I., & Fady, B. (2021). Signature of mid-Pleistocene lineages in the European silver fir (*Abies alba* Mill.) at its geographic distribution margin. *Ecology and Evolution*, 00, 1-16. https://doi.org/10.1002/ece3.7886