Barbu Cm

Allienne Jf

Dr N Bech

Pyrenean ptagmigans decine under climatic and human influences through the Holocene

Keywords: demographic history, climate change, changes in land use, Holocene, habitat fragmentation, Lagopus muta

In Europe, the Quaternary is characterized by climatic fluctuations known to have led to many cycles of contraction and expansion of species geographical ranges. In addition, during the Holocene, historical changes in human occupation such as colonization or abandonment of traditional land uses can also affect habitats. These climatically or anthropically induced geographic range changes are expected to produce considerable effective population size change, measurable in terms of genetic diversity and organization. The rock ptarmigan (Lagopus muta) is a small-bodied grouse occurring throughout Northern hemispheric arctic and alpine tundra. This species is not considered threatened at a continental scale, but the populations in the Pyrenees are of concern because of their small population size, geographical isolation and low genetic diversity. Here, we used 11 microsatellites to investigate genetic variations and differentiations and infer the overall demographic history of Pyrenean rock ptarmigan populations. The low genetic variability found in these populations has been previously thought to be the result of a bottleneck that occurred following the last glacial maximum (i.e., 10 000 years ago) or more recently (i.e., during the last 200 years). Our results clearly indicate a major bottleneck affecting the populations in the last tenth of the Holocene. We discuss how this decline can be explained by a combination of unfavorable and successive events that increased the degree of habitat fragmentation.

INTRODUCTION

In Europe, the Quaternary is characterized by climatic fluctuations known to have led to many cycles of contraction and expansion of species geographical ranges [START_REF] Hewitt | Some genetic consequences of ice ages, and their role in divergence and speciation[END_REF][START_REF] Taberlet | Comparative phylogeography and postglacial colonization routes in Europe[END_REF]. These range changes, induced by climatic cycles, are expected to produce considerable reorganization in effective population size and so in genetic structure patterns [START_REF] Hewitt | Some genetic consequences of ice ages, and their role in divergence and speciation[END_REF]. The last climatic transition occurred at the end of the last glacial maximum between the Pleistocene and the Holocene 10000 years ago. Following this transition period, species had to take colonization routes toward refuges where they found habitats corresponding to their ecological requirements [START_REF] Taberlet | Comparative phylogeography and postglacial colonization routes in Europe[END_REF]. Species adapted to temperate climate were in general confined to refuges during glacial periods, whereas the cold-adapted species were in refugia during interglacial periods [START_REF] Stewart | Refugia revisited: individualistic responses of species in space and time[END_REF].

There is growing evidence that the Earth's climate has been warming since the beginning of the 1900s (IPCC, 2001). This recent climate warming is expected to have a strong impact on the distribution of biomes, forcing species to shift their range. Indeed, several studies have shown that at temperate latitudes, cool climate species have already started to shift their ranges northward or have entered mountaintop refugia [START_REF] Willis | PALEOECOLOGY: the refugial debate[END_REF][START_REF] Walther | Trends in the upward shift of alpine plants[END_REF][START_REF] Parmesan | Ecological and evolutionary responses to recent climate change[END_REF]. In mountainous area, such a warming-induced upward shift has been observed, for example, in the beech forest (Fagus sylvatica) of the Montseny mountains (northeast Spain), which has shifted upward by 70m, at the highest altitudes in the last 55 years (Pen˜uelas and Boada, 2003). This uppermost limit of the forest, referred to as the tree line, is often a clear vegetation boundary separating the landscape into forested and treeless zones (Del [START_REF] Barrio | Thermal seasonality of the high mountain belts of the Pyrenees[END_REF][START_REF] Korner | A world-wide study of high altitude treeline temperatures[END_REF][START_REF] Shoo | Detecting climate change induced range shifts: where and how should we be looking?[END_REF]. Vegetation shifts to higher altitudes have a direct impact on the distribution and diversity of highmountain biota, and the isolation and contraction of these habitats are hypothesized to cause population decline and genetic bottlenecks.

In addition to the effects of climate warming on the distribution of biomes, several studies have shown that changes in human activities strongly influence the processes of vegetation succession [START_REF] Grace | Impacts of climate change on the tree line[END_REF][START_REF] Carcaillet | Holocene tree-limit and distribution of Abies alba in the inner French Alps: anthropogenic or climatic changes?[END_REF][START_REF] Chauchard | Patterns of land-use abandonment control treerecruitment and forest dynamics in Mediterranean mountains[END_REF][START_REF] Jetz | Projected impacts of climate and land-use change on the global diversity of birds[END_REF]Batllori and Gutie´rrez, 2008;[START_REF] Pompe | Climate and land use change impacts on plant distributions in Germany[END_REF][START_REF] Popy | A weak upward elevational shift in the distributions of breeding birds in the Italian Alps[END_REF]Ame´ztegui et al., 2010). Among these, changes in land use, including intense exploitation and land abandonment, are major factors affecting the maintenance and distribution of species. For example, land overuse and deforestation have been shown to cause population decline and genetic bottlenecks in various ecosystems, from alpine to tropical environments [START_REF] Regnaut | Population genetics of capercaillie (Tetrao urogallus) in the Jura and the Pyrenees: a non-invasive approach to avian conservation genetics[END_REF][START_REF] Goosens | Genetic signature of anthropogenic population collapse in orang-utangs[END_REF].

The rock ptarmigan (Lagopus muta) is one of the cool climate species that found refuge at northern latitudes or in mountaintop habitats following Quaternary temperature increases [START_REF] Caizergues | Population genetic structure of rock ptarmigan Lagopus mutus in Northern and Western Europe[END_REF]. Indeed, the rock ptarmigan is a small-bodied grouse with a large circumpolar distribution that also occurs in southern mountains of Europe. In this respect, ptarmigans are considered among the most effective sentinel species for evaluating responses to environmental changes. In southern Europe, the species inhabits alpine habitats (2000-2900m above sea level (a.s.l.)) dominated by ericaceous shrubs mixed with other dwarf alpine plants and screes. The rising of trees and shrub lines in alpine areas associated with the warming of temperatures would likely result in less alpine habitat in smaller patches. In addition, most of the rock ptarmigan habitats in European mountains are affected by livestock grazing, ski resort infrastructures and more generally mountain sports. In the French Pyrenees, small and threatened populations of L. muta pyrenaica survive in mountaintops of a restricted range occurring at the southern margin of the repartition area [START_REF] Novoa | Effect of weather on the reproductive rate of rock Ptarmigan Lagopus muta in the eastern Pyrenees[END_REF]. These Pyrenean populations are characterized by two genetic clusters separated by a valley [START_REF] Bech | Population genetic structure of rock ptarmigan in the 'sky islands' of French Pyrenees: implications for conservation[END_REF]. Moreover, these populations have a very low genetic diversity relative to populations occurring in the Alps and in Scandinavia [START_REF] Caizergues | Population genetic structure of rock ptarmigan Lagopus mutus in Northern and Western Europe[END_REF][START_REF] Bech | Population genetic structure of rock ptarmigan in the 'sky islands' of French Pyrenees: implications for conservation[END_REF]. According to [START_REF] Caizergues | Population genetic structure of rock ptarmigan Lagopus mutus in Northern and Western Europe[END_REF], this low genetic diversity could result from a bottleneck arisen during the Pyrenees colonization after the last glacial maximum or from a recent reduction in the effective population size.

Based on four supplementary microsatellite markers (total=11 markers) and more individuals (i.e., 45% supplementary; total =170), this study aims to reinvestigate the genetic structure and reconstruct the demographic history of Pyrenean rock ptarmigan populations. This study attempts to test if Pyrenean populations experienced a demographic bottleneck and estimate its date and magnitude. Combining these results with the compilation of historical data on human colonization and trends in tree line within Pyrenees, we discuss the relative impact on ptarmigan populations of the end of the last glacial maximum -10000 years before present (YBP) -versus more recent events such as global warming and changes in agropastoral activities (<200 YBP).

MATERIALS AND METHODS Population genetics of rock ptarmigans

Sampling and genotyping. Rock ptarmigan samples (n= 170) have been collected over the last 20 years in various regions of the central and eastern Pyrenees. Some of these samples (n=117) have been analyzed previously, evidencing a marked genetic structure of the population [START_REF] Bech | Population genetic structure of rock ptarmigan in the 'sky islands' of French Pyrenees: implications for conservation[END_REF]. Liver and wings were collected during the hunting season (from late September through late October) in Luchon (n=15), Melles (n=13), Ustou (n= 19), Aston (n=45) and Campcardos (n= 16) (Figure 1). In the eastern part of the Pyrenees, where the species is not currently hunted, we collected molted feathers or muscle samples from dead birds found during census in July and August or from birds captured for radio tracking (from July to early October) in Caranc¸a-Puigmal (n=35) and Canigou (n= 27). Based on the geographical distribution of the 170 samples and the moderate dispersal abilities of the rock ptarmigan -radio-tracking data showed a mean natal dispersal distance o6km for young birds (Novoa, unpublished data), we grouped individuals into seven populations, each representing a particular geographical location. DNA was extracted and genotyped as detailed by [START_REF] Bech | Population genetic structure of rock ptarmigan in the 'sky islands' of French Pyrenees: implications for conservation[END_REF] but using the following 11 microsatellite markers developed for other avian species: LLSD7, LLSD4 and LLSD10 [START_REF] Piertney | Isolation and characterization of hypervariable microsatellites in the red grouse Lagopus lagopus scoticus[END_REF]; LLST7 [START_REF] Piertney | Local genetic structure in red grouse (Lagopus lagopus scoticus): evidence from microsatellite DNA markers[END_REF]; TTT2 [START_REF] Caizergues | Isolation and characterization of microsatellite loci in black grouse (Tetrao tetrix)[END_REF]; BG6, BG15 and BG19 [START_REF] Piertney | Polymorphic microsatellite DNA markers in black grouse (Tetrao tetrix)[END_REF]; Aru1F25 and Aru1G4 [START_REF] Ferrero | Sixteen new polymorphic microsatellite markers isolated for red-legged partridge (Alectoris rufa) and related species[END_REF]; and MNT45 [START_REF] Reed | Microsatellite loci for genetic mapping in the turkey (Meleagris gallopavo)[END_REF].

Genetic analysis. We used the Micro-Checker software [START_REF] Oosterhout | Micro-checker: software for identifying and correcting genotyping errors in microsatellite data[END_REF] to identify the presence of null alleles or scoring errors resulting from stuttering. Deviation from Hardy-Weinberg expectancies was analyzed using a Markov chain-based test (1000 permutations) as implemented in GenePop version 3.4 [START_REF] Raymond | Genepop (version 1.2): population genetics software for exact tests and ecumenicism[END_REF].

Linkage disequilibria were analyzed using the global tests with Bonferroni correction to adjust the level of significance for multiple testing in FSTAT version 2.9.3.2 [START_REF] Goudet | FSTAT, a program to estimate and test gene diversities and fixation indices (version 2.9.3)[END_REF]. We computed allelic richness, expected heterozygosity and Fis values for all loci and populations and finally the pairwise Fst values according to [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] using the same software. However, for highly variable multilocus data sets, Fst of 1 is often statistically impossible to obtain [START_REF] Hedrick | A standardized genetic differentiation measure[END_REF]. Therefore, we calculated the upper most Fst value for our microsatellite data set, using both FSTAT and RECODE DATA version 0.1 [START_REF] Meirmans | Using the AMOVA framework to estimate a standardised genetic differentiation measure[END_REF]. This allowed to calculate G'st values, a standardized measure of genetic differentiation calculated by dividing obtained Fst values by the maximum possible Fst.

We also examined the population structure using the Bayesian approach implemented in STRUCTURE version 2.2 software [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF]. We chose the admixture model with correlated allele frequencies among populations. Each simulation (i.e., K= 1-7) was replicated 20 times as recommended by [START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF], with a 5 X 10 5 burn-in period followed by 10 6 steps. The number of independent populations (K) was determined with Structure Harvester version 0.6.1 [START_REF] Earl | STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method[END_REF].

Past demographic analysis. To assess changes in the effective population size, we used the method of [START_REF] Storz | Testing for genetic evidence of population expansion and contraction: an empirical analysis of microsatellite DNA variation using a hierarchical Bayesian model[END_REF] implemented in the MSVAR 1.3 software. This likelihood-based Bayesian method is known for producing robust reconstructions of the demographic history of populations [START_REF] Girod | Inferring population decline and expansion from microsatellite data: a simulation-based evaluation of the MSVAR method[END_REF]. MSVAR analyzes the information from the full allelic distribution using a coalescence-based framework under the stepwise mutation model. This model assumes a stable population of size N1 (ancestral size) began to decrease or increase exponentially some time in history (ta generations ago) to the current population size N0 (current size). Given the allelic distribution and prior information, this method allows estimating N0, N1 and the time T since the population started changing in size. Prior distributions for N0, N1 and T were assumed to be lognormal. The means and standard deviations of these prior distributions were themselves drawn from prior (or hyperprior) distributions. For all simulations, we used a generation time of 3.5 years, the mean age of reproducing females in this population [START_REF] Novoa | Demographic traits of rock ptarmigan in the french Alps and PyreneesSandercock[END_REF]. For each genetic cluster identified from STRUCTURE, as previously described, we performed four independent runs (replicates) with different random seeds and uninformative hyperpriors (Supplementary Information S1). To prevent autocorrelation in the Markov chains, they were run for 5 X 10 9 steps and thinned by 5 X 10 4 to get a total of 100000 draws from the posterior distribution. The first 10% of the chains were discarded to avoid bias in parameter estimates resulting from the starting conditions. We used the statistic of [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF] (G&R diagnostic) to test for the convergence of the four independent chains for each cluster. Smooth posterior distributions were estimated from the raw MSVAR data using the local regression as implemented in the Locfit R package [START_REF] Loader | Local Regression and Likelihood[END_REF]. Finally, we performed a Bayes factor analysis to test, within the Bayesian framework, the relative probability every 50 years over the last 20000 years, thus including the whole Holocene and the end of the Pleistocene [START_REF] Kass | Bayes factors[END_REF]. Bayes factors 43 indicate positive evidence [START_REF] Kass | Bayes factors[END_REF].

Environmental data

Agropastoral activities. The historic changes of agropastoral activities were documented based on archival data compiled from the following socioeconomic and local studies: [START_REF] Susplugas | Le sol et la vegetation dans le Haut-Vallespir (Pyrenees-Orientales) Montpellier[END_REF][START_REF] Becat | Atlas de Catalunya Nord[END_REF][START_REF] Chevalier | La vie humaine dans les Pyrenees ariegeoises[END_REF][START_REF] Adph | Hommes et troupeaux des Pyrenees: Images des Hautes-Pyrenees[END_REF][START_REF] Bosch De La Trinxeria | Records d'un excursionista Biblioteca Selecta[END_REF][START_REF] Rendu | La Montagne d'Enveig. Une estive pyre´ne´enne dans la longue dure´e[END_REF]. These documents allowed reconstructing the regional history of livestock in various parts of the Pyrenees from 1750 to recent times.

Temperatures. Mean annual temperatures, determined locally, were used to estimate temperature changes according to both altitude (adiabatic thermal gradient) and time (temperature increase because of climate change). These temperature data, provided by Meteo-France, were obtained from the Pic du Midi meteorological station (the oldest station in the Pyrenees: 31T; 4757712° N, 266865° E; 2877m a.s.l.) and from the Tarbes meteorological station, located 35km north at lower elevation (31T; 4790610° N, 262272° E; 310m a.s.l.). Temperature increase was estimated at the Pic du Midi station by calculating the difference between the mean annual temperature for the 1980-2009 and 1882-1911 periods (the period of 30 years is considered long enough to smooth out year-to-year variations). To assess if changes in temperature have been consistent in the two stations, we performed a Pearson's correlation test between temperature data from stations of Tarbes and Pic du Midi using GraphPad InStat software version 3.10.

RESULTS

Genetic analysis

Several markers had null alleles, but all were randomly distributed across loci and populations, indicating a good quality of the data. Across all markers and populations, we found no significant Fis values (P-value threshold after Bonferroni correction, P=0.0036) and no evidence for linkage disequilibrium (P-value threshold after Bonferroni correction, P=0.0024), suggesting no departure from the Hardy-Weinberg equilibrium. The microsatellite marker set reveals a polymorphism characterized by a heterozygosity ranging from 0 to 0.736 and an allelic richness ranging from 1 to 5.992. The details of raw microsatellite data are available in Supplementary Information S2.

Measures of genetic differentiation detected a significant genetic structure within the seven populations -mean Fst value=0.045 and upper most Fst value=0.5728. We show in Table 1 both measures of genetic differentiation: Fst and G 0 st. This genetic structure was corroborated by Bayesian assignment of individuals to K population clusters according to the allelic distribution, revealing a highest posterior probability for two clusters (K=2; Evanno's delta K value=250) located on either side of the Segre valley (Figure 1 and Supplementary Information S3). Individuals from Canigou and Caranca-Puigmal were mainly assigned to the first cluster (mean±s.d.: 0.907±0.001 and 0.812±0.001, respectively), whereas individuals in Campcardos, Aston, Ustou, Melles and Luchon were mainly assigned to the second cluster (mean±s.d.: 0.622±0.002, 0.745±0.001, 0.813±0.001, 0.873±0.002 and 0.756±0.001, respectively). These results are consistent with previous conclusions by [START_REF] Bech | Population genetic structure of rock ptarmigan in the 'sky islands' of French Pyrenees: implications for conservation[END_REF].

Population genetic decline

Even using different prior and hyperprior sets, we observed a good convergence of the four independent chains in each cluster (G&R diagnostic=1.05 and 1.06 for clusters 1 and 2, respectively). Posteriors indicated Pyrenean ptarmigans underwent major population decline, with an effective ancestral population size (N1) about two orders of magnitude larger than the current population (N0) for both clusters (Figure 2). Precisely, the posterior median for the effective population size plunged in the first cluster from 16491 to 100 individuals and in the second cluster from 13715 to 178 (details shown in Supplementary Information S4). The posterior distribution for the estimated time since the population started to decrease (T) exhibited a median value of 786 and 525 years for clusters 1 and 2, respectively (Figure 3 and Supplementary Information S4). The Bayes factor provided a positive evidence (Bayes factors >3; [START_REF] Kass | Bayes factors[END_REF] for a period ranging from 200 to 2000 YBP, highlighting a period outside of both initial hypotheses (i.e., approximately 10000 YBP or o200 YBP) (Figure 4).

History of agropastoral activities

Historical statistics on agropastoral activities in the Pyrenees are rare and scattered, but all socioeconomic studies of the area concerning the eighteenth and nineteenth centuries suggest a high level of mountain land use at the beginning of the nineteenth century and a decline in regional livestock from 1850 onward (Figure 5) [START_REF] Susplugas | Le sol et la vegetation dans le Haut-Vallespir (Pyrenees-Orientales) Montpellier[END_REF][START_REF] Becat | Atlas de Catalunya Nord[END_REF][START_REF] Chevalier | La vie humaine dans les Pyrenees ariegeoises[END_REF][START_REF] Adph | Hommes et troupeaux des Pyrenees: Images des Hautes-Pyrenees[END_REF][START_REF] Bosch De La Trinxeria | Records d'un excursionista Biblioteca Selecta[END_REF][START_REF] Rendu | La Montagne d'Enveig. Une estive pyre´ne´enne dans la longue dure´e[END_REF].

The numbers of sheep reported for the nineteenth century suggest that Pyrenean upland areas were subject to major grazing pressure for several decades. In the three administrative regions (Hautes-Pyrenees, Ariege and Pyrenees-Orientales) encompassing most of our sampling sites, administrative documents of the nineteenth century report 232000 sheep in Ariege in 1810 [START_REF] Chevalier | La vie humaine dans les Pyrenees ariegeoises[END_REF], 435364 in the Hautes-Pyre´ne´es in 1866 [START_REF] Adph | Hommes et troupeaux des Pyrenees: Images des Hautes-Pyrenees[END_REF] and 300 000 in the Pyrenees-Orientales in 1850 [START_REF] Becat | Atlas de Catalunya Nord[END_REF], suggesting that approximately 1 000 000 sheep were grazing in the mountains of these regions, an area of 4166 km 2 in the middle of the nineteenth century. In 2008, 213400 sheep have been recorded in the same area. On a more anecdotal note, Bosch de la Trinxeria (1887), in his records of travel in the eastern Pyrenees, reported a herd of 16000 sheep grazing an alpine pasture of approximately 800Ha from mid-June to mid-September of 1880. This represents a density of 2000 sheep per 100Ha, a grazing pressure at least four times greater than those currently reported in studies addressing the effects of grazing at extreme densities on birds [START_REF] Fuller | Changes in sheep numbers in Britain: implications for bird populations[END_REF][START_REF] Loe | Positive short-term effects of sheep grazing on the alpine avifauna[END_REF].

Temperature change

At the Pic du Midi meteorological station, the mean annual temperature (±s.e.) between 1882 and 1911 was 1.50±0.12°C, whereas between 1980 and 2009 it was 0.50±0.16°C. The mean temperature then increased by approximately 1.0°C in 100 years, although larger variations have been observed on the scale of decade (Figure 5). Results from Pearson's correlation test indicated that temperatures recorded in the Tarbes and the Pic du Midi stations have increased at a similar manner since 1940 (N=56, r 2 =0.63, P<0.0001), suggesting a coherent trend of the altitudinal gradient over the years.

DISCUSSION

In agreement with a previous study [START_REF] Bech | Population genetic structure of rock ptarmigan in the 'sky islands' of French Pyrenees: implications for conservation[END_REF], we found a significant genetic structure for the ptarmigans across the Pyrenees. Based on a larger sample of individuals and loci, these results corroborated the presence of two genetic clusters located on either side of the Segre valley. Indeed, given the low dispersal abilities of rock ptarmigans, the 18-km width of the Segre valley represents a current barrier to the gene flow for rock ptarmigans [START_REF] Bech | Population genetic structure of rock ptarmigan in the 'sky islands' of French Pyrenees: implications for conservation[END_REF]. These Pyrenean ptarmigan populations were thought to have undergone a severe population decline probably during the Pyrenees colonization after the last glacial maximum and/or more recently [START_REF] Caizergues | Population genetic structure of rock ptarmigan Lagopus mutus in Northern and Western Europe[END_REF]. Our results clearly indicated that Pyrenean ptarmigans experienced a major decrease during the Holocene, most likely between 200 and 2000 years ago. This scenario takes place outside of both previously proposed hypotheses (i.e., approximately 10 000 YBP or <200 YBP). Our inferences are consistent over two different genetic clusters and with different prior and hyperprior sets. Moreover, the low level of population genetic differentiation within clusters suggests that artifactual bottlenecks are unlikely and that the population decline cannot result from a genetic substructure [START_REF] Chikhi | The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes[END_REF]. We cannot reject the occurrence of an additional older bottleneck as the model assumes the existence of a single bottleneck. Older decline events may have occurred, but our genetic analysis highlights a major population decrease during the last tenth of Holocene. In the following section, we put in perspective the genetic decline we observe with major events affecting the Pyrenean landscape during the Holocene.

9300-7600 YBP: the early Holocene

Although it is difficult to address the history of the Pyrenean landscape during the Holocene, the palynological and anthracological data clearly showed that the tree line probably reached its highest levels (i.e., >2460m a.s.l.) during this period [START_REF] Cunill | Holocene treeline changes on the south slope of the Pyrenees: a pedoanthracological analysis[END_REF][START_REF] Cunill | Holocene highaltitude vegetation dynamics in the Pyrenees: a pedoanthracology contribution to an interdisciplinary approach[END_REF]. This seems to be a generalized phenomenon in Southern Europe because the tree line also reached its highest level in the Alps at the same period [START_REF] Ali | Pinus cembra L. (arolla pine), a common tree in the inner French Alps since the early Holocene and above the present tree line: a synthesis based on charcoal data from soils and travertines[END_REF][START_REF] Talon | Reconstruction of Holocene high-altitude vegetation cover in the French southern Alps: evidence from soil charcoal[END_REF]. The high elevation of the tree line during this period likely resulted in a loss and fragmentation of habitat for the rock ptarmigan. Considering that the current tree line was at 2100m a.s.l. (Del [START_REF] Barrio | Thermal seasonality of the high mountain belts of the Pyrenees[END_REF][START_REF] Korner | A world-wide study of high altitude treeline temperatures[END_REF], we estimate that the amount of suitable rock ptarmigan habitat in the Pyrenees during this period was 80% less than the current area. As a consequence, we expect that the species distribution range reached then its highest degree of fragmentation during the Holocene. Even if the climatic cycles of this period did not seem the triggering factor of the bottleneck, we suggest that the associated high level of fragmentation would have negatively affected the rock ptarmigan populations.

5100-2200 YBP: the beginning of the Bronze Age

During this period, the pedoanthracology coupled to the palynology showed successive phases of fires that led to the largest Pinus decline and a clear subsequent increase in herbaceous plants and shrubs. The increase of both herbaceous plants and shrubs, such as Genista, is evidence of recurrent burns to maintain pastures [START_REF] Cunill | Holocene highaltitude vegetation dynamics in the Pyrenees: a pedoanthracology contribution to an interdisciplinary approach[END_REF]. This corresponds to the intensification of anthropic activities in the Pyrenean mountains at the end of the Neolithic and beginning of the Bronze Age [START_REF] Galop | La foret, l'homme et le troupeau dans les Pyrenees: 6000 ans d'histoire de l'environnment entre Garonne et Mediterranee[END_REF][START_REF] Miras | Prehistoric impact on landscape and vegetation at high altitudes: an integrated palaeoecological and archaeological approach in the eastern Pyrenees (Perafita valley, Andorra)[END_REF][START_REF] Cunill | Holocene treeline changes on the south slope of the Pyrenees: a pedoanthracological analysis[END_REF]. The Bronze Age was followed, in the Pyrenees, by the Roman period, which marked the beginning of human establishment in the landscape. Even if the varying degrees of Romanization likely led to a great regional disparity between Pyrenean valleys [START_REF] Cunill | Holocene highaltitude vegetation dynamics in the Pyrenees: a pedoanthracology contribution to an interdisciplinary approach[END_REF], human activities, not reaching intensities later seen, probably limited the habitat fragmentation further favoring the reestablishment of local Pyrenean ptarmigan populations and facilitating an increase of their genetic diversity.

1400-500 YBP: the medieval age

At the beginning of the medieval period, the Pyrenees were marked by a generalized intensification of human activities that profoundly affected the Pyrenean landscape. Among these activities were the Grand Transhumances organized by feudal and church powers. This is attested by the sedimentary macrocharcoal records that show periodic burns undoubtedly intended for the creation and preservation of pasturing spaces [START_REF] Marugan | Historia del Pallars del origens als nostres dies[END_REF]Curnill et al., 2013). Even if moderate sheep densities can favor the expansion of several bird species [START_REF] Loe | Positive short-term effects of sheep grazing on the alpine avifauna[END_REF], intense livestock grazing is generally considered to be detrimental to ground-nesting birds [START_REF] Fuller | Changes in sheep numbers in Britain: implications for bird populations[END_REF] and may affect the reproduction of grouse [START_REF] Myrberget | The effect of sheep grazing on a willow grouse breeding habitat[END_REF][START_REF] Baines | The implications of grazing and predator management on the habitats and breeding success of black grouse Tetrao tetrix[END_REF][START_REF] Jouglet | Impact du paturage ovin estival sur l'habitat et les effectifs du tetras-lyre (Tetrao tetrix) dans les Hautes-Alpes[END_REF][START_REF] Storch | Genetic correlates of spatial population structure in central European capercaillie Tetrao urogalus and black grouse T. tetrix: a project in progress[END_REF]. Despite the absence of studies specifically linking declines in ptarmigan numbers to grazing pressure, anecdotal observations indicate that intensive sheep grazing negatively affects the nesting attempts of rock ptarmigans [START_REF] Miquet | Sur des destructions de nids de Lagopede alpin Lagopus mutus dues aux ongules domestiques[END_REF]. On this basis, and in agreement with our genetic dating estimates, we suggest this significant intensification of human establishment in the Pyrenean Mountain has marked the beginning of a decline of the overall reproductive performance of the rock ptarmigan populations and resulted in a marked reduction of the effective population size during this period. As the population size reduction is the main reason for loss of genetic diversity [START_REF] Frankham | Introduction to Conservation Genetics[END_REF], we suggest this population reduction, coupled to a habitat still fragmented, has limited the gene flow and created genetically isolated small populations.

500-150 YBP: a continued pressure

According to the palynology, this period is marked by the highest recorded levels of Poaceae (Cerealia) and by the minimal values of the Pinus pollen. This suggests a great human impact on the Pyrenean landscape [START_REF] Cunill | Holocene highaltitude vegetation dynamics in the Pyrenees: a pedoanthracology contribution to an interdisciplinary approach[END_REF]. According to archived data, this human pressure showed, at the beginning of the 1800s (i.e., 300 YBP), an expansion of agropastoral activities that resulted in herds of several thousands of sheep, which did not favor the rock ptarmigan populations. In our historical compilation, we show that at the beginning of the 1850s, in the Pyrenees, the agropastoral activities began to switch from a situation of overgrazing to land abandonment ;

The land abandonment, since the 1850s, may have had a positive effect on grouse reproduction, at least temporarily. However, it has been shown in Mediterranean mountains that abandonment of traditional land uses probably triggered the decrease of seminatural open fields in altitude, a suitable habitat for the rock ptarmigan [START_REF] Dirnbock | A regional impact assessment of climate and land-use change on alpine vegetation[END_REF][START_REF] Penuelas | A global change-induced biome shift in the Montseny mountains (NE Spain)[END_REF][START_REF] Chauchard | Patterns of land-use abandonment control treerecruitment and forest dynamics in Mediterranean mountains[END_REF]. The sudden switch from massive overgrazing to almost total abandonment then likely did not allow a rebound of the rock ptarmigan population. In addition, according to the temperature data we present, the beginning of the 1900s in the Pyrenees corresponds to the beginning of a warming period that likely favored the expansion of subalpine forests (Batllori and Gutie´rrez, 2008), of which the uppermost limit was no longer limited by grazing [START_REF] Chauchard | Patterns of land-use abandonment control treerecruitment and forest dynamics in Mediterranean mountains[END_REF]. Indeed, in the contemporary period, warming is considered to be one of the main drivers of the general shift of alpine plants to higher altitudes [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF][START_REF] Walther | Trends in the upward shift of alpine plants[END_REF][START_REF] Parmesan | Ecological and evolutionary responses to recent climate change[END_REF][START_REF] Sekercioglu | Climate change, elevational range shifts and bird extinctions[END_REF][START_REF] Lenoir | A significant upward shift in plant species optimum elevation during the 20th century[END_REF]. In the Pyrenees, this temperature increase (i.e., 0.10°C per decade) seems to be the result of cloud cover or the greenhouse effect (Bu¨cher and Dessens, 1991) and is consistent with the increase of temperature recorded globally and at high altitudes [START_REF] Vuille | Mean annual temperature trends and their vertical structure in the tropical Andes[END_REF][START_REF] Pepin | A global comparison of surface and free-air temperatures at high elevation[END_REF]. Thus, at the end of the nineteenth century, the concomitant occurrence of agropastoral abandonment and warming likely favored the densification and altitudinal progression of forests at the expense of seminatural open fields in altitude (i.e., rock ptarmigan habitat), with land uses finally reaching very low levels during the twentieth century [START_REF] Garcia-Ruiz | Land-use changes and sustainable development in mountain areas: a case study in the Spanish Pyrenees[END_REF][START_REF] Davasse | Forets charbonniers et paysans dans les Pyrenees de l'est du moyen age à nos jours[END_REF]. Recently, several observations of such an upward shift (i.e., 200m upward) have been recorded during the same period for various plant species in the Alps and Pyrenees [START_REF] Walther | Trends in the upward shift of alpine plants[END_REF][START_REF] Lenoir | A significant upward shift in plant species optimum elevation during the 20th century[END_REF]. Based on the increase of recorded temperatures in the Pyrenees, we suggest that alpine habitats continued to shift upward and fragment. As a consequence, alpine species, such as the rock ptarmigan, took refuge in sky islands [START_REF] Haila | A conceptual genealogy of fragmentation research: from island biogeography to landscape ecology[END_REF], which continue to contract at the beginning of the twenty-first century.

The decline of Pyrenean ptarmigans likely results from a combination of unfavorable events occurring during the Holocene. Accessing to the demographic history of threatened species and populations from genetic data is often complex because genetic signal could have been influenced by different evolutionary events [START_REF] Salmona | Signature of a pre-human population decline in the critically endangered reunion island endemic forest bird Coracina newtoni[END_REF][START_REF] Girod | Inferring population decline and expansion from microsatellite data: a simulation-based evaluation of the MSVAR method[END_REF]. Nevertheless, our time estimates from genetic analysis, together with historical data, suggest the increasing anthropic pressure intensification, from the medieval age through the nineteenth century, has been a major factor in rock ptarmigan decline. Since 1850, the reversal from heavy land use to almost total abandonment, together with a rise in temperatures, suggests a continued or even increased pressure on alpine populations and rock ptarmigans in particular. In addition, because the method employed in our genetic analysis more easily identifies ancient population collapses rather than very recent ones [START_REF] Salmona | Signature of a pre-human population decline in the critically endangered reunion island endemic forest bird Coracina newtoni[END_REF][START_REF] Girod | Inferring population decline and expansion from microsatellite data: a simulation-based evaluation of the MSVAR method[END_REF] and because the modal values of the posterior distribution indicate very recent peaks, a strong and recent decline cannot be excluded.

The rock ptarmigan populations, in the Alps, do not present genetic impoverishment, suggesting that this population did not experience a genetic bottleneck even though the events mentioned above also likely occurred in the Alps. We suggest that the bottleneck, detected in this study, is mainly due to the fragmentation of the Pyrenean rock ptarmigan habitat. Fragmentation, more than habitat size reduction, characterizes the division of the two genetic clusters we studied. The higher elevation of the Alps prevented high levels of fragmentation of remaining suitable alpine rock ptarmigan habitats; consequently, populations in the Alps did not experience population declines and associated loss of genetic diversity.

According to recent studies [START_REF] Revermann | Habitat at the mountain tops: how long can Rock Ptarmigan (Lagopus muta helvetica) survive rapid climate change in the Swiss Alps? A multi-scale approach[END_REF] and the predictions of Meteo-France [START_REF] Nogues-Bravo | Exposure of global mountain systems to climate warming during the 21st Century[END_REF][START_REF] Nogues-Bravo | Climate change in Mediterranean mountains during the 21st century[END_REF], we can expect, in the near future, a further upward shift of the Pyrenean tree line and consequently an upward shift of the suitable habitat of the rock ptarmigan. However, the future tree line position is challenging to predict in the Pyrenees as increasing temperature may also induce more frequent summer droughts, negatively affecting the tree line dynamics by reducing growth and seedling viability (Batllori and Gutie´rrez, 2008). According to our results on the impact of anthropogenic disturbance, the future trend of seminatural open fields in altitude is even more difficult to predict because of the uncertain evolution of the temporal and spatial intensity of anthropogenic disturbance during the next century. The protection of seminatural open fields in altitude (i.e., highest life zones) is crucial for preserving natural systems and biodiversity. But, among altitudinal biodiversity, some species considered are early indicators of how ecosystems respond to climate change and land use [START_REF] Koumoundouros | Current genetic isolation and fragmentation contrasts with historical connectivity in an alpine lizard (Cyclodomorphus praealtus) threatened by climate change[END_REF]. As prediction of habitat trend is complex, observation of these high-altitude sentinel species should inform our conservation policies. For these conservation policies, our results bring a hope: the strong anthropogenic impact we demonstrate on alpine species habitat can be used to preserve the seminatural open fields of altitude and their associated biodiversity FIGURE 1 a: Study area map representing geographical location of sampled populations and their membership in genetic clusters highlighted using STRUCTURE software. The black and grey colors, in the pie charts, indicate the average membership coefficient of each population to genetic cluster 1 and genetic cluster 2, respectively. b: Cross section linking the major mountain peaks east of the Pyrenees (white dashed curve on the site map a). This cross section passes through mountain peaks and crosses the Se`gre valley in its narrowest part. The black area represents the current rock ptarmigan habitat. We used a Digital Elevation Model provided by the NASA shuttle radar topography mission (SRTM 90 m Digital Elevation; file: strm_37_03) implemented in MapInfo v. 7.5 software to generate this figure. We used the universal transverse Mercator coordinate system (UTM, Zone 31 North, WGS84).

FIGURE 2

Current and ancestral population sizes. Posterior distribution for current (N0) and ancestral (N1) effective population size. The grey dashed lines represent the prior distribution for N0 and N1 (for more details see Supplementary Information S1). The solid lines correspond to the different runs for the cluster 1 (black line) and cluster 2 (bold black line). This figure reveals that the current population sizes are smaller than the ancestral population sizes for all studied populations. This indicates that bottlenecks occurred in each population.

FIGURE 3

Time since bottleneck (T). Posterior distribution for time since population decrease (a: log scale; b: scale in years before present). The solid lines correspond to cluster 1 (black line) and cluster 2 (black bold line). The grey dashed lines in a represent the different prior distributions for T. This figure reveals that bottlenecks have occurred concomitantly in both clusters.

FIGURE 4

Bayes factor for different periods. A recent scenario assuming the changes in climatic and land-uses corresponds to the period from current to 2000 years before present. The scenario assuming the Pyrenees colonization after the last glacial maximum corresponds to the period around 10,000 years before present. We pooled the four independent runs for each cluster in order to compute BF values of T. These BF values were computed every 50 years over the last 20,000 years. Black line: BF distribution for the cluster 1 and black bold line: BF distribution for the cluster 2. We considered BF values > 3 (i.e., above the horizontal dasheddotted) as significant positive evidence [START_REF] Kass | Bayes factors[END_REF].

FIGURE 5

Historic livestock (dashed lines) and temperature (solid line) trends in the Pyrenees. The solid curve corresponds to the mobile mean (with a period of 10 years) of the temperature, provided by Meteo France and obtained from the meteorological station at Pic du Midi (31T; 4757712 N, 266865 E; 2877 m a.s.l.). The dashed curves with grey triangles and crosses represent respectively the historic volution of sheep herds in the central Pyrenees (A.D.P.H., 1982) and in the eastern Pyrenees [START_REF] Becat | Atlas de Catalunya Nord[END_REF][START_REF] Chevalier | La vie humaine dans les Pyrenees ariegeoises[END_REF].

TABLE 1

Pairwise Fst values for each population comparison (below diagonal) and their significance level (above diagonal) P-value threshold is adjusted with the Bonferroni correction, P=0.002. The Fst values that reveal a significant genetic differentiation are in bold. G'st values, a standardized measure of genetic differentiation, are in parentheses.

ACKNOWLEDGEMENTS

We are very grateful to M Conesa, V Parmain, JP Me´taille´, D Galop, C Gilbert and S Beltran for constructive discussion and three anonymous referees for helpful comments. We thank S Somot and A Braun (Meteo-France) for providing predicted temperature values. We also thank V Bech, Karthik Sethuraman and the English Manager Science Editing for editing the English language in the paper. This study was supported by the Office National de la Chasse et de la Faune Sauvage, the Ministere de l'Enseignement Superieur et de la Recherche and the Centre National de la Recherche Scientifique (CNRS).

DATA ARCHIVING

Data were deposited in the Dryad repository: doi:10.5061/ dryad.55425.

CONFLICT OF INTEREST

The authors declare no conflict of interest.