Alain Billionnet
email: alain.billionnet@ensiie.fr.

Quantifying extinction probabilities of endangered species for phylogenetic conservation prioritization may not be as sensitive as might be feared

Keywords: Biodiversity conservation, phylogenetic diversity, protected species, extinction risk, optimization, greedy algorithm

Quantifying extinction probabilities of endangered species for phylogenetic conservation prioritization may not be as sensitive as might be feared

Introduction

An increasing number of species are threatened with extinction [START_REF] Stork | Re-assessing current extinction rates[END_REF]. The causes for this growing threat include habitat loss, degradation, and fragmentation, agriculture, human over-population, deforestation, poaching and hunting, invasive species, and climate change (Collen et al. 2013). The loss of biodiversity that would result from a massive extinction would have disastrous consequences [START_REF] Cardinale | Biodiversity loss and its impact on humanity[END_REF]). The international community agrees that the protection of threatened and endangered species is a priority issue but the resources available to implement these protections are limited (Convention on Biological Diversity 2011). It is therefore crucial to carefully select the species on which the protection efforts will be concentrated in order to preserve biodiversity as much as possible. There are many ways to assess biodiversity (see, for example, [START_REF] Moreno | Measuring biodiversity in the Anthropocene: a simple guide to helpful methods[END_REF]) and also the very extensive book by [START_REF] Maclaurin | What Is Biodiversity? The University of[END_REF] on this subject). In this paper, the measurement of biodiversity is based on the notion of phylogenetic diversity. We consider that the survival of a protected species is assured -which is not the case for an unprotected species -and we consider the classic conservation problem which consists in determining, among a set of species T, the "best" subset of species to protect within a limited budget. More precisely, we consider the problem of determining, among a set of species T, k species to be protected so as to maximize the resulting expected phylogenetic diversity of T, k being a fixed value. The phylogenetic diversity is a widely used criterion for expressing the amount of biodiversity (evolutionary history) of a set of species captured by a phylogeny. The reader can refer to [START_REF] Maclaurin | What Is Biodiversity? The University of[END_REF]Sterelny 2008, Chap.7 andPellens &Grandcolas 2016) for a very comprehensive presentation of the role of phylogenetic diversity, in large sense, in the field of biological conservation. The expected phylogenetic diversity requires the precise knowledge of the extinction probabilities of the species concerned and these probabilities are generally not easy to quantify [START_REF] Mace | Quantification of extinction risk: IUCN's system for classifying threatened species[END_REF][START_REF] Mooers | Converting endangered species categories to probabilities of extinction for phylogenetic conservation prioritization[END_REF].

In this note, we use the measure of phylogenetic diversity created by Faith (Faith 1992a(Faith , 1992b) and we study the influence, on the solutions of the considered conservation problem, of the values attributed to the extinction probabilities.

The expected phylogenetic diversity

A rooted phylogenetic tree can be considered as a directed tree, ) , , , (

 T A N H 
, where

N is the set of nodes, A, the set of arcsor directed edges or branches -,T, the set of leaves representing the species, and  , the vector of the length of each branch. The branches of A are directed from the root of H to its leaves. The nodes of the phylogenetic tree that are neither root nor leaves are called internal nodes. The internal nodes have one predecessor and at least two successors. For all A a  , we denote by a L the set of species for which there is a path from the terminal node of a to these species. As mentioned in Section 1, we use the measure of phylogenetic diversity created by Faith: The phylogenetic diversity of the set T of species associated with the tree H is equal, by definition, to  A a a

 . Intuitively, it represents the total amount of evolutionary history embodied in a set of species, since the time of the most recent common ancestor of the set [START_REF] Chao | Phylogenetic diversity measures based on Hill numbers[END_REF]. In this paper, we refer to PD to indicate this measure. The reader is referred to [START_REF]Biodiversity and evolutionary history: useful extensions of the PD phylogenetic diversity assessment framework[END_REF][START_REF] Faith | The PD phylogenetic diversity framework: linking evolutionary history to feature diversity for biodiversity conservation[END_REF] and also [START_REF] Collen | Conservation prioritization in the context of uncertainty[END_REF][START_REF] Isaac | Mammals on the EDGE : Conservation priorities based on threat and phylogeny[END_REF]) for a comprehensive discussion about the use of PD in measuring biodiversity. [START_REF] Mimouni | Phylogenetic diversity and its conservation in the presence of phylogenetic uncertainty: a case study of cladoceran communities in urban waterbodies[END_REF] point out that PD is an established index with a large enough body of literature regarding its computational and combinatorial aspects. For example, [START_REF] Hartmann | Phylogenetic diversity: From combinatorics to ecology[END_REF] review algorithmic, mathematical, and stochastic results concerning PD. As mentioned in the introduction, the book of [START_REF] Maclaurin | What Is Biodiversity? The University of[END_REF] and the book edited by Pellens and Grandcolas (2016) deal in detail with the importance of phylogenetic diversity, and especially Faith's phylogenetic diversity, for biodiversity conservation.

Suppose that an extinction probability, i p , is associated with each species, i, of T. The probability that the information associated with an arc a of the tree is retained is equal to the probability that at least one of the species "reachable" from ai.e. belonging to a Lsurvives. This way, we can define (from the Faith's PD) the expected phylogenetic diversity (ePD) of T. It is equal to the sum on all arcs of T of the probability that the information associated with this arc is conserved, multiplied by the length of this arc (Witting & Loeshchke 1995). More precisely, the ePD of

T is equal to    A a a a P ) 1 (  where    a L i i a p P
. We suppose here that the extinction probabilities of the species are independent. The ePD is recognized as a relevant criterion for evaluating the biodiversity associated with a set of species. It accounts for complementarity among species [START_REF] Faith | Phylogenetic diversity, functional trait diversity and extinction: avoiding tipping points and worst-case losses[END_REF]. As an example, consider the hypothetical ultrametric phylogenetic tree of Figure 1 and assume that the 7 species considered all have the same extinction probability. Let us examine the impact of the following 2 events (of same probability): the extinction of the 3 species 1, 6 and 7, the 4 other species remaining alive, and the extinction of the 3 species 2, 3 and 4, the 4 other species remaining alive. The first event would cause the loss of 22 units while the second would cause the loss of only 13 units.

The use of ePD requires the knowledge of the extinction probabilities, i p , for each species i, and these probabilities are generally not easy to quantify [START_REF] Mooers | Converting endangered species categories to probabilities of extinction for phylogenetic conservation prioritization[END_REF]. We consider that the studied species are divided into several groups more or less threatened and that an extinction probability is associated with each group. It is possible, for example, to consider the 5 IUCN Red List categories: Least Concern, Near Threatened, Vulnerable, Endangered and Critically Endangered. Note that ranking species in these categories is a difficult task [START_REF] Mace | Quantification of extinction risk: IUCN's system for classifying threatened species[END_REF]. The extinction probability of the species i, i p , is then equal to the extinction probability assigned to the group of threat to which it belongs. [START_REF] Redding | Incorporating phylogenetic measures into conservation prioritization[END_REF] and [START_REF] Isaac | Mammals on the EDGE : Conservation priorities based on threat and phylogeny[END_REF] have proposed ways to transform the Red List categories to extinction probabilities. [START_REF] Mooers | Converting endangered species categories to probabilities of extinction for phylogenetic conservation prioritization[END_REF] have shown that the choice of transformation of Red List categories to different extinction probabilities can affect the resulting species ranking for the metrics EDGE [START_REF] Isaac | Mammals on the EDGE : Conservation priorities based on threat and phylogeny[END_REF]) and HEDGE [START_REF] Steel | Hedging our bets: The expected contribution of species to future phylogenetic diversity[END_REF][START_REF] Mooers | Converting endangered species categories to probabilities of extinction for phylogenetic conservation prioritization[END_REF][START_REF] Forest | Phylogeny, extinction and conservation: embracing uncertainties in a time of Phil[END_REF][START_REF] Nunes | The price of conserving avian phylogenetic diversity: a global prioritization approach[END_REF]. In this study, we examine the impact of the extinction probability value attributed to each group (and thus to each species) on the optimal protection strategy associated with the classical conservation problem outlined above and precisely defined below.

Determining the best set of species to protect

The problem that we consider here can be stated as follows: given a phylogenetic tree ) , , , (

 T A N H 
and extinction probabilities, p i , for each species, i, of T, determine a set of k species of T whose protection maximizes the ePD of T. Remember that we assume that the protection of a species ensures its survival. In other words, its probability of extinction decreases from p i to 0. This basic problem is a particular case of the Noah's Ark problem which was first defined by [START_REF] Weitzman | The Noah's ark problem[END_REF]. It is shown in (Hartman & Steel 2006 ) and [START_REF] Pardi | Algorithms on phylogenetic trees[END_REF]) that the considered conservation problem can be solved easily by a greedy algorithm. The algorithm begins with an empty set, S , and sequentially adds to S the species, i, which maximizes the ePD of the set } {i S  until the cardinal of S is k.

Link with the HEDGE metric: The well-known HEDGE metric defines species priority ranks for a fixed tree and fixed extinction probabilities. The HEDGE score of the species i is equal to the increase of the ePD of T obtained by protecting this species with respect to the ePD of T obtained in protecting no species (Figure 2). One way of approaching the considered problem -determine the species of T to be protected so as to maximize the resulting ePD of T-consists in calculating the HEDGE score of each species and then retaining the k species with the highest scores. It is to be noted, however, that the HEDGE score does not take into account the protection actions actually taken. Indeed, if the species that has the highest HEDGE score is effectively protected (its extinction probability becomes equal to 0) the HEDGE scores of some other species can decrease significantly. To overcome this drawback, a natural extension of the HEDGE method, called I-HEDGE, has been proposed by [START_REF] Jensen | I-HEDGE: determining the optimum complementary sets of taxa for conservation using evolutionary isolation[END_REF]. This method consists in gradually calculating the HEDGE score of each species taking into account at each stage the fact that a new species becomes protected. The new species that becomes protected is the one with the highest score, this one being calculated taking into account the species already chosen to be protected. This method, which therefore integrates the notion of phylogenetic complementarity, corresponds exactly to the greedy algorithm proposed by [START_REF] Pardi | Algorithms on phylogenetic trees[END_REF] and Hartman & Steel (2006), and mentioned above.

Figure 2 illustrates the computation of the EDGE and I-HEDGE scores.

Main result

We present in this section experimental results which show that the choice of the values attributed to the extinction probabilities associated with the 5 groups of species considered does not have an immense influence on the solution of the problem under consideration although this choice may significantly affect the set of species to be protected as a priority.

Consider the two scenarios defined by the probabilities the optimal set of species to be protected in scenario  . As we have seen above, these sets can be determined by the greedy algorithm. For each subset S of T and for each scenario  , we denote by ) , ( S T ePD  the ePD of T generated by the protection of the species of S and calculated with the probabilities of the scenario  .

We are interested first of all in the case where one decides to protect the species of 1 * S .

We then compute the resulting ePD of T with the two sets of extinction probabilities, that is to say The first of these gaps provides the relative error made (on the maximum ePD) in considering that the real scenario is the scenario 2 when in reality it is the scenario 1. The second gap provides the relative error made in considering that the real scenario is the scenario 1 when in reality it is the scenario 2.

We first calculated gap 1 and gap 2 on 10,000 randomly generated instances corresponding to non ultrametric trees and also on 10,000 randomly generated instances corresponding to ultrametric trees. The description of all these instances and the corresponding detailed results are presented in Section 5.

For the non ultrametric trees, the largest of the gaps obtained (out of 20,000) is equal to 2.63% (Table 5). Therefore, we can say that, at least under these experimental conditions, the possible error committedon the maximum ePD valuesby the choice of transformation of categories to different extinction probabilities is relatively limited. However, it should be noted that the set of k species, optimal in scenario 1, may be quite different from the set of k species optimal in scenario 2. We measured their dissimilarity by the cardinal of their symmetric difference divided by the cardinal of their union. This criterion therefore varies from 0 to 1. For example, the dissimilarity between the two optimal sets corresponding to the largest gap is equal to 0.70 which corresponds to two very different sets (Table 5). Note that [START_REF] Mooers | Converting endangered species categories to probabilities of extinction for phylogenetic conservation prioritization[END_REF] had already observed that the choice of transformation of Red List categories to different extinction probabilities can affect significantly the resulting species ranking for the two well-known metrics EDGE and HEDGE.

For the ultrametric trees, the largest of the gaps obtained (out of 20,000) is slightly higher than the one obtained for non-ultrametric trees since it is equal to 4.13% (Table 5).

Here too, the dissimilarity between the two optimal sets corresponding to this gap is important. It is equal to 0.59 (Table 5).

We then calculated gap 1 and gap 2 on 10,000 instances generated from the phylogenetic tree associated with the 52 Madagascar lemurs (species and subspecies) and studied by [START_REF] Pardi | Algorithms on phylogenetic trees[END_REF]. In this case, the phylogenetic tree considered (structure, length of the branches and category to which each of the 52 lemurs belongs) is definitively fixed. It includes 97 nodes in total. Each of the 10,000 instances is generated by associating to the 52 lemurs 2 different scenarios of extinction probabilities. The description of these instances and the corresponding detailed results are presented in Section 5. As in experiments on randomly generated phylogenetic trees we found that 1) the possible error committedon the maximum ePD valuesby the choice of transformation of categories to different extinction probabilities is relatively limited (0.89%) and 2) the two optimal sets corresponding to the largest error are very different (dissimilarity equal to 0.42).

Description of the instances and detailed results

Randomly generated non ultrametric phylogenetic trees

The 10,000 instances considered were generated in the following way:

Construction of the hypothetical phylogenetic trees: First, the number of non-leaf nodes is uniformly and randomly chosen between 50 and 1,000. Then, for each of these nodes, the number of children is uniformly and randomly generated between 2 and max d , and max d is itself drawn at random uniformly between 2 and 4.

Lengths of the branches:

The length of each branch is uniformly and randomly generated between 1 and max  , and max 

is itself drawn at random uniformly between 5 and 20.

Assignment of extinction probabilities:

The set of species is divided into five categories and the category to which each species belongs is drawn at random in two different ways. In the first case, this group is drawn at random uniformly in the set {1,2,3,4,5}. The numbers of species in each group are then comparable. In the second case, the group is drawn at random to approximately obtain the following distribution of the species in the five groups, as in [START_REF] Mooers | Converting endangered species categories to probabilities of extinction for phylogenetic conservation prioritization[END_REF]): 2%, 4%, 9%, 9%, 76%. Considering the case 1 or the case 2 is chosen at random in an equiprobable manner. To define the extinction probabilities of an unprotected species we were inspired by the description of the 5 IUCN Red List categories (IUCN 2012):

For each tree considered, two sets of extinction probabilities (corresponding to two hypothetical scenarios) are drawn at random uniformly in the intervals given in Table 1. are drawn at random uniformly in the intervals given in Table 2. These intervals represent a variation around the values of the extinction probabilities proposed by [START_REF] Pardi | Algorithms on phylogenetic trees[END_REF]. Maximum number of species that can be protected: As in experiments on randomly generated phylogenetic trees, the maximum number of species that can be protected is equal to the greatest integral number less than or equal to  multiplied by the number of species ( 52), and

 is drawn at random uniformly in the set {0.1, 0.2, 0.3, 0.4, 0.5}.

For the three types of phylogenetic trees considered (non-ultrametric, ultrametric, Madagascar lemur), statistics on the 10,000 generated instances, characteristics of the instance corresponding to the largest gap and solutions of the conservation problem considered for this instance and for the two associated scenarios are presented in Tables 345.

With regard to the phylogenetic tree associated with Madagascar lemurs we also measured the impact that possible errors on the lengths of the branches could have on the solution of the considered conservation problem: determine k species to protect in order to maximize the resulting ePD. For this we have definitively assigned to each species the survival probability given by [START_REF] Pardi | Algorithms on phylogenetic trees[END_REF] and we have randomly generated two sets of values for the lengths of the branches, thus obtaining two scenarios. The value of each branch is uniformly and randomly generated in the interval [0.75 a , 1.25 a ] where  a is the value of the branch given by Pardi. We therefore consider that a significant error ( % 25 

) can be made on the lengths of the branches. We have thus built 10,000 different instances. For each instance, we have determined the 2 optimal sets S *1 and S *2 corresponding to the two scenarios and we have calculated the two relative gaps, gap 1 and gap 2 , defined in Section 4.

Here, ) , ( S T ePD  is the ePD of T generated by the protection of the species of S and calculated with the branch lengths of the scenario  ( =1 or 2). The maximum of these 2 gaps over the 10,000 instances generated, that is to say the maximum error that can be made by choosing the wrong scenario, is equal to 1.68%. This error is larger than that obtained by considering uncertainties on the survival probabilities. Moreover the two optimal sets of species corresponding to the largest gap are very different since their dissimilarity is equal to 0.64. If we admit that an error of % 50  can affect the length of the branches, the maximum relative gap goes to 5.28% and the dissimilarity to 0.57.

Conclusions and perspectives

To determine which of a given set of endangered species should be protected as priorities, the criterion of the expected phylogenetic diversity (ePD) is recognized as a relevant criterion. However, its use requires explicit knowledge of the extinction probabilities of the species considered when they are protected and also when they are not protected. We assume here that the extinction probability of the protected species is zero. We were interested in the influence of the conversion of the extinction risk of a species to an extinction probability of this species. The extinction risk of a species is given by the fact that it belongs to a certain category (such as those of the IUCN Red List). We showed experimentally, on randomly generated ultrametric and non-ultrametric phylogenetic trees as well as on a real tree, that the choice of the extinction probabilities values is not as important as might be feared. More precisely, if we consider 2 scenarios, corresponding to two different conversions, the ePD generated by the protection of a set of k species, optimal in the scenario 1 and calculated with the probabilities of this scenario is not very distant from the ePD generated by the protection of a set of k species, optimal in the scenario 2 and calculated with the probabilities of the scenario 1. In other words, the error committed by retaining the scenario 2 when the real scenario would be the scenario 1 is not very important, at least in our experimental conditions. The terms "scenario 1" and "scenario 2" can be exchanged in the preceding two sentences. It should be noted, however, that the set of k species, optimal in the scenario 1, may be quite different from the set of k species optimal in the scenario 2.

We considered in this study that the survival probabilities of the protected species were equal to 1. It would be interesting to study other values for these probabilities and to measure the impact of the values retained (if the species is protected and if it is not protected) on the optimal ePD. It would also be interesting to study the impact of the uncertainties that inevitably exist on the structure of the phylogenetic trees considered and the values attributed to their branches. In section 5, we have given some preliminary results that show that the consequences of uncertainties on branch lengths may not be negligible.

  A hypothetical ultrametric phylogenetic tree associated with 7 species; it has 10 branches. The length of each branch is indicated next to the branch.

Figure 2 .

 2 Figure2. A hypothetical phylogenetic tree with 8 species. The lengths of the branches are indicated next to each branch. The extinction probabilities of each unprotected species (p i ) are shown under each species. The HEDGE score for each species (HEDGE i ) is shown in brackets under the corresponding species. The I-HEDGE score for each species (I-HEDGE i ) is shown in square brackets under the corresponding species. According to these last values, the 3 species whose protection maximizes the ePD of T are 1, 6 and 8, and the corresponding ePD is equal to 55.6 (30.8416+10.7520+8.0064+6.0000).

  probability of the species of the group g in the scenario .

Table 1 .

 1 Randomly generated trees: possible values of the extinction probability of each category.

	Category Interval in which the extinction
		probability is chosen at random
	1	[0.50,1.00]
	2	[0.20,0.50]
	3	[0.10,0.20]
	4	[0.05,0.10]
	5	[0.00,0.05]

Table 2 .

 2 Madagascar lemurs: possible values of the extinction probability of each category.

	Category Interval in which the extinction
		probability is chosen at random
	1	[0.90,1.00]
	2	[0.65,0.85]
	3	[0.40,0.60]
	4	[0.15,0.35]
	5	[0.00,0.10]

Table 3 .

 3 Statistics on the 10,000 randomly generated instances.

	Non-ultrametric trees	Ultrametric trees	Madagascar lemurs tree

Table 4 .

 4 Characteristics of the instances corresponding to the largest relative gap.

		Non-	Ultrametric	Madagascar
		ultrametric trees	trees	lemurs tree
	Number of non-leaf nodes	727	147	-
	Number of species	1420	290	-
	d max	4	4	-
	max 	7	9	-
	Length of the longest path	72	-	-
	from the root to a species
	Length of the shortest path	16	-	-
	from the root to a species
	Number of species	297, 293, 267,	55, 64, 50,	-
	in each group	274, 289	58, 63
	Scenario 1: Extinction	0.86, 0.24, 0.14,	0.51, 0.47, 0.12,	0.99, 0.76, 0.48,
	probabilities by category:	0.07, 0.01	0.06, 0.01	0.32, 0.05
	Scenario 2: Extinction	0.50, 0.46, 0.13,	0.99, 0.22, 0.18,	0.95, 0.67, 0.59,
	probabilities by category:	0.07, 0.03	0.07, 0.03	0.16, 0.01
		0.1	0.2	0.3
		(142 species).	(58 species)	(15 species).
	Scenario 1: ePD obtained	11003.39	8262.40	94.07
	by protecting no species
	Scenario 2: ePD obtained	11237.61	7748.15	98.70
	by protecting no species
	ePD obtained by protecting	13485.00	10520.00	143.18
	all the species

Table 5 .

 5 Results associated with the solution of the problem in the two scenarios for the instances described in Table4.

		Non-ultrametric	Ultrametric	Madagascar
		trees	trees	lemurs tree
	Scenario 1: Maximum ePD	12176.37	9407.09	131.21
	Scenario 1: Distribution by category	category 1: 142 category 1: 31	category 1: 5
	of the species to be protected		category 2: 27	category 2: 6
				category 3: 2
				category 4: 2
	ePD obtained in scenario 1	11856.61	9265.58	130.87
	by protecting the species of the
	set optimal in scenario 2:
	gap 1 (%)	2.63	1.50	0.26
	Scenario 2: Maximum ePD	12001.33	9605.16	132.97
	Scenario 2: Distribution by category	category 1: 66	category 1: 55	category 1: 4
	of the species to be protected	category 2: 76	category 2: 3	category 2: 6
				category 3: 5
	ePD obtained in scenario 2	11914.76	9208.83	131.79
	by protecting the species of the
	set optimal in scenario 1
	gap 2 (%)	0.72	4.13	0.89
	Dissimilarity of the two optimal sets	0.70	0.59	0.42

Maximum number of species that can be protected: Finally, the maximum number of species that can be protected is equal to the greatest integral number less than or equal to  multiplied by the number of species, and  is drawn at random uniformly in the set {0.1, 0.2, 0.3, 0.4, 0.5}.

We also conducted other experiments slightly different from the previous ones by modifying the way extinction probabilities are attributed. In these experiments, the extinction probabilities of the species may vary within a category. More precisely, the extinction probability of each species is drawn at random uniformly in the range corresponding to its category. For example, if species 27 belongs to category 2, its extinction probability is drawn at random uniformly in the interval [0.20, 0.50]. We again generated 10,000 instances and we obtained results comparable to those we had previously obtained: The mean of the 20,000 relative gaps is equal to 0.19%, the corresponding standard deviation is equal to 0.13%, and the largest relative gap is equal to 1.21%. The instance corresponding to this largest gap concerns 67 species divided into five categories respectively formed of 16, 12, 10, 15 and 14 species. In this instance, the maximum number of species that can be protected is equal to 6 and the dissimilarity of the two optimal sets obtained is equal to 0.50.

Randomly generated ultrametric phylogenetic trees

The 10,000 instances considered were generated in exactly the same way as for nonultrametric trees, and then the lengths of the branches linked to a species in each resulting tree are modified to obtain an ultrametric tree. Let L be the length of the longest path connecting the root to a species and let l i be the length of the path connecting the root to the species i. For each species i, the length of the arc connected to this species is modified by adding to it the value L-l i .

Phylogenetic tree associated with Madagascar lemurs [START_REF] Pardi | Algorithms on phylogenetic trees[END_REF] The 10,000 instances considered were generated in the following way:

Assignment of extinction probabilities: The set of species is divided into five categories. The category to which each species belongs is given by [START_REF] Pardi | Algorithms on phylogenetic trees[END_REF]. For each instance considered, two sets of extinction probabilities (corresponding to two hypothetical scenarios)


Marine Pouget

Sami Youssef

Jérémy Migliore

Marianick Juin

Frédéric Médail

Alex Baumel

Phylogeography sheds light on the central-marginal hypothesis in a Mediterranean narrow endemic plant

Keywords: Genetic diversity, geographical range, ecological niche, pleistocene, Provence, conservation, migration, population distinctiveness

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

A fundamental issue in ecology and conservation biology is how evolutionary processes shape variation across species' whole geographical range [START_REF] Lesica | When are peripheral populations valuable for conservation ?[END_REF][START_REF] Brown | On the relationship between abundance and distribution of species[END_REF][START_REF] Gaston | The structure and dynamics of geographic range[END_REF][START_REF] Guo | Incorporating latitudinal and central-marginal trends in assessing genetic variation across species ranges[END_REF]. According to the "central-marginal" hypothesis [START_REF] Eckert | Genetic variation across species' geographical ranges: the central-marginal hypothesis and beyond[END_REF], genetic diversity decreases and genetic differentiation increases toward the periphery of the range. The ecological and evolutionary features of populations living at range limits of species [START_REF] Sexton | Evolution and ecology of species range limits[END_REF] are of particular importance in the context of global change [START_REF] Sexton | Gene flow increases fitness at the warm edge of a species' range[END_REF]. Low habitat quality at the range margins may lead to higher rates of extinction, while increased isolation reduces the rate of recolonization as well as of gene flow.

However, the decline in genetic diversity at range limits is not a ubiquitous trend, neither in multispecific comparisons [START_REF] Eckert | Genetic variation across species' geographical ranges: the central-marginal hypothesis and beyond[END_REF][START_REF] Eckstein | Genetic structure among and within peripheral and central populations of three endangered floodplain violets[END_REF][START_REF] Hardie | Evolutionary ecology at the extremes of species' ranges[END_REF][START_REF] Guo | Incorporating latitudinal and central-marginal trends in assessing genetic variation across species ranges[END_REF], nor within intra-specific case studies [START_REF] Hoban | Range-wide distribution of genetic diversity in the North American tree Juglans cinerea: a product of range shifts, not ecological marginality or recent population decline[END_REF]. As emphasized by several authors [START_REF] Garner | Genetic diversity across a vertebrate species' range: a test of the central-peripheral hypothesis[END_REF][START_REF] Eckert | Genetic variation across species' geographical ranges: the central-marginal hypothesis and beyond[END_REF][START_REF] Moeller | Population genetics and the evolution of geographic range limits in an annual plant[END_REF], looking for the causes of a central-marginal structure of genetic diversity is particularly challenging because of the diverse alternative models and the effects of current demographic and ecological factors combined with persistent historical influences. Exploring historical influences, e.g. whether marginal populations are younger than core populations [START_REF] Eckert | Genetic variation across species' geographical ranges: the central-marginal hypothesis and beyond[END_REF], has been seen as relevant in the Northern Hemisphere, where Pleistocene climatic oscillations are still perceptible in the layout of the current patterns of widespread species distribution, through persistence in refugial areas and post-glacial migration roads [START_REF] Taberlet | Comparative phylogeography and postglacial colonization[END_REF][START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF][START_REF] Petit | Glacial refugia: hotspots but not melting pot of genetic diversity[END_REF][START_REF] Hewitt | Genetic consequences of climatic oscillations in the Quaternary[END_REF][START_REF] Hampe | Conserving biodiversity under climate change: the rear edge matters[END_REF][START_REF] Soltis | Comparative phylogeography of unglaciated eastern North America[END_REF]. Thus, range shifts and distance from areas ensuring the persistence of populationsnotably refugia [START_REF] Stewart | Refugia revisited: individualistic responses of species in space and time[END_REF] -are often invoked to explain the coarse geographical changes in genetic variation [START_REF] Garner | Genetic diversity across a vertebrate species' range: a test of the central-peripheral hypothesis[END_REF][START_REF] Hoban | Range-wide distribution of genetic diversity in the North American tree Juglans cinerea: a product of range shifts, not ecological marginality or recent population decline[END_REF]and references in Hardie and Hutchings, 2010).

Yet neither the historical influences shaping genetic variation, nor the central-marginal hypothesis, have been invoked in connection with Mediterranean endemic plants. Our aim here was to explore these questions, for three reasons. Firstly, as underlined recently by Nieto [START_REF] Feliner | Southern European glacial refugia: A tale of tales[END_REF], the post-glacial colonization paradigm pattern is insufficient as a basis for exploring Mediterranean evolutionary patterns and processes during the Pleistocene.

Secondly, in the highly structured landscapes of the Mediterranean, narrow endemic plants are so frequent [START_REF] Médail | Ecological characteristics and rarity of endemic plants from southeast France and Corsica: implications for biodiversity conservation[END_REF][START_REF] Casazza | Ecological characteristics and rarity of endemic plants of the Italian Maritime Alps[END_REF] that they are one of the main features of Mediterranean biodiversity [START_REF] Thompson | Ecological differentiation of Mediterranean endemic plants[END_REF]. Thirdly, because their distribution is structured in small areas, narrow endemic plants could make it easier to test the central-marginal hypothesis at the scale of the whole range; moreover, they constitute key candidates to study persistence and speciation [START_REF] Stebbins | Endemism and speciation in the california flora[END_REF][START_REF] Chown | Speciation and rarity: separating cause from consequence[END_REF].

Recent studies of some narrow endemic Mediterranean plants (e.g. [START_REF] Medrano | Geographical structuring of genetic diversity across the whole distribution range of Narcissus longispathus, a habitat-specialist, Mediterranean narrow endemic[END_REF][START_REF] Molins | Molecular phylogeography of Thymus herba-barona ( Lamiaceae ): Insight into the evolutionary history of the flora of the western Mediterranean islands[END_REF][START_REF] Mayol | Patterns of genetic variability and habitat occupancy in Crepis triasii (Asteraceae) at different spatial scales: insights on evolutionary processes leading to diversification in continental islands[END_REF][START_REF] Nicoletti | Spatial Genetic Structure of Campanula sabatia, a Threatened Narrow Endemic Species of the Mediterranean Basin[END_REF] support the hypothesis of their long-term persistence during the Pleistocene and their fluctuating genetic diversity.

Nevertheless, the link between historical influences, geographical range and genetic diversity in Mediterranean narrow endemic plants remains relatively unelucidated in the context of the central-marginal hypothesis. To explore this link, we focus on the Caryophyllaceae Arenaria provincialis Chater and Hallid. (Chater and Halliday, 1964), which inhabits the calcareous outcrops of Southern Provence, in south-eastern France.

Recent studies on geographical distribution and ecological niches [START_REF] Véla | Exemple d'une nouvelle évaluation du statut de menace suivant les critères de l'UICN version 3.1. : le cas de l'endémique provençale Arenaria provincialis Chater & Halliday (= Gouffeia arenarioides DC., Caryophyllaceae)[END_REF][START_REF] Youssef | Factors Underlying the Narrow Distribution of the Mediterranean Annual Plant Arenaria provincialis (Caryophyllaceae)[END_REF] provide the basis for examining central and marginal populations of A.

provincialis throughout its entire range, taking into account both abundance and ecological aspects. According to the central-marginal hypothesis, populations of A. provincialis living near the centre of the distribution are expected to exhibit more genetic diversity than populations living further afield (Fig. 1). Moreover, at range limits, populations of A. provincialis experience contrasting ecological conditions [START_REF] Youssef | Factors Underlying the Narrow Distribution of the Mediterranean Annual Plant Arenaria provincialis (Caryophyllaceae)[END_REF] that could raise conservation issues if, concurrently, their genetic diversity decreased and their isolation increased with respect to more central populations.

Here, occurrence data are used to show that abundance is not uniformly distributed across the range of A. provincialis, which confirms the importance of exploring the central-marginal hypothesis. First, we test the consistency of our hypotheses on a central-marginal structure of genetic diversity using multilocus markers, through analysis of amplified fragment length polymorphism (AFLP) markers. Second, we infer the historical influences throughout the distribution from the phylogeography, based on nucleotidic variation of four chloroplastic loci analysed in a Bayesian framework. Finally, we use ecological data to link marginality in the ecological gradient with range position and phylogeography, thus revealing the conservation significance of our results.

MATERIALS AND METHODS

Study system and occurrence data

In south-eastern France, the winter annual plant Arenaria provincialis [START_REF] Chater | Arenaria L[END_REF] is distributed within an area of less than 145 km² along an altitudinal gradient and in different rocky habitats [START_REF] Véla | Exemple d'une nouvelle évaluation du statut de menace suivant les critères de l'UICN version 3.1. : le cas de l'endémique provençale Arenaria provincialis Chater & Halliday (= Gouffeia arenarioides DC., Caryophyllaceae)[END_REF][START_REF] Baumel | Ecological magnitude and fine scale dynamics of the Mediterranean narrow endemic therophyte, Arenaria provincialis (Caryophyllaceae)[END_REF]. Its southern limits occur in the driest context of steep scree slopes at low altitudes and short distances from the coastline, whereas its northern limits lie in colder conditions near the summits of the outcrops.

Variation in abundance of Arenaria provincialis throughout its entire range was estimated from the variation in density of occurrence on samples collected during several field campaigns. Data were filtered to ensure that occurrences were at least 100 metres apart, to avoid redundancy and overestimation of abundance due to oversampling. 591 occurrences were used to model abundance variation by applying an isotropic Gaussian kernel function (density.ppp function; "spatstat" package; R Development Core Team, 2011).

Occurrence data were also used to set discrete, non-overlapping geographical units to organize measures of genetic diversity and for phylogeography analysis (see below).

Preliminary analyses distinguished 5 geographical units (A-E; Fig. 1a; Table 1). These geographical units were delineated via non-hierarchical clustering analysis (kmeans function, "stats" package; R Development Core Team, 2011) from coordinates of A. provincialis occurrences. The barycentre of the distribution is situated in the northern part of area C.

DNA extraction and molecular methods

Forty populations of Arenaria provincialis were sampled throughout its entire geographical range (Table 1) between April 2009 and May 2010. Leaves were dried with silica gel, and then stored at -20°C. Total DNA extraction followed the CTAB protocol of [START_REF] Doyle | A rapid DNA isolation method for small quantities of fresh tissues[END_REF].

The polymorphism of the chloroplast genome was examined on two or three randomly sampled individuals per population (n = 94) for four plastid regions: internal part of matK gene, trnK intron plus 5' part of matK, tnrL-trnF and trnT-trnL intergenic spacer (supplementary material Table S1). Sequence characterization of the four chloroplastic sequences and the combined matrix are reported in Table 2. According to the most recent phylogeny published on this genus [START_REF] Youssef | Factors Underlying the Narrow Distribution of the Mediterranean Annual Plant Arenaria provincialis (Caryophyllaceae)[END_REF], A. provincialis and A. cinerea are sister species. The same cpDNA regions were sequenced for three samples of A. cinerea sampled at three different locations, as well as one sample of A. ciliata and one sample of A. hispida (Table 1; Genbank accession number in supplementary material Table S2).

Polymerase Chain Reactions (PCR) were performed in a 50 µl volume reaction that contained 2 µl of each primer (5 mM), 1.25 units of Taq polymerase (Q-Biogen), 5 µl of incubation buffer, 6 µl of MgCl 2 (10 mM), 4 µl of dNTP (25 mM) and 5 µl of DNA (50 ng.µl -1 ). Amplification was performed using PTC-200 thermocycler (MJ Research, Inc.; cycling conditions in supplementary material Table S3). Direct sequencing (ABI 3730 xl, Applied Biosystems) was carried out by LGC genomics services (Germany). After quality checking and alignment, the cpDNA data set contained 94 sequences.

Variation in genetic diversity throughout the range was investigated by multilocus genotyping based on the Amplified Fragment Length Polymorphism technique (AFLP) developed by [START_REF] Vos | AFLP: a new technique for DNA fingerprinting[END_REF]. Two primer combinations for the selective PCR (fluorescent dyes in brackets, Eurofins MWG Operon, Ebersberg, Germany) were used: EcoRI (6-FAM) -AAGG ⁄ MseI CCAG, EcoRI (6-FAM)-AACG ⁄ MseI-CAAG. Selective amplification was performed in 20µL volumes with 5 pmol of each primer, 0.65 mM of MgCl2, 0.5 mM of dNTPs, 1 unit of Taq DNA polymerase (Q-Biogen) and 5 µL of 100x diluted preamplification. The selective amplification thermocycle profile was 94°C for 2 min, 10 cycles of 94°C for 30 s, 65°C for 30 s (step -0.7°C per cycle), 72°C for 1 min, followed by 20 cycles at 94°C for 30 s, 56°C for 30 s, 72°C for 1 min and 72°C for 5 min. PCR products were separated and quantified on an ABI 3730xl DNA analyzer (Applied Biosystems, Foster City, CA). The reliability of AFLP markers was checked by repeating the complete analysis from DNA amplification to AFLP screening on 10 samples for each pair of primers.

More than forty polymorphic markers were generated, the minimum threshold for the resolution of a genetic structure with AFLPs [START_REF] Gaudeul | Genetic diversity and differentiation in Eryngium alpinum L. (Apiaceae): comparison of AFLP and microsatellite markers[END_REF]. Amplified labelled fragments were revealed on an ABI 3730 DNA Analyzer (PE Applied Biosystems) and fingerprint profiles were scored with Genemapper 3.7 (PE Applied Biosystems). The results of the scoring were exported as a presence/absence matrix and used for further analysis. Two individuals per population were studied for all populations. However, AFLP failed for some samples, and the final AFLP data set comprised 73 genotypes.

Variation in AFLP across the geographical range

Genetic diversity was estimated by the percentage of polymorphic loci, the number of private markers and the Shannon entropy index, which were computed from AFLP data for each of the 5 geographical units (A-E) described above (Table 3), using GenAlEx 6.41 (Genetic Analysis in Excel; [START_REF] Peakall | GenAlEx 6: genetic analysis in Excel. Population genetic software for teaching and research[END_REF]. Genetic differentiation was estimated by the mean of genotype divergence within each geographical unit. Genotype divergences were based on jaccard index (dist.binary function; "ade4" package; R Development Core Team,

2011).

Phylogeography of A. provincialis

A network of relationships among cpDNA haplotypes was established from a NeighborNet analysis [START_REF] Bryant | Neighbor-Net: an agglomerative method for the construction of phylogenetic networks[END_REF] using SplitsTree4 [START_REF] Huson | Application of phylogenetic networks in evolutionary studies[END_REF]. Arenaria cinerea, A. hispida, and A. ciliata were included in the analysis so as to identify the branch representing the oldest divergence in A. provincialis. To compute pairwise distances between haplotypes, a distance using a GTR + gamma substitution model was chosen according to jModeltest [START_REF] Posada | jModelTest: phylogenetic model averaging[END_REF] with empirical estimation of nucleotide frequencies (Fig. 2).

The phylogeography of A. provincialis was based on the Bayesian MCMC methods of BEAST 1.7 [START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF], using all sequences obtained for A. provincialis (four cpDNA genes and n=94 samples). All five geographical units (A-E) were used to set a geographic state for each sample. The ancestral geographic states were inferred at each node using the Bayesian Stochastic Search Variable Selection (BSSVS; [START_REF] Lemey | Bayesian phylogeography finds its roots[END_REF] (Fig. 3). Without a priori knowledge about the history of A. provincialis, we decided not to specify prior rate for the ancestral geographic states; thus, geographic state transitions (dispersal) are estimated solely from the genealogical processes. This method was recently used to build discrete phylogeographies in intra-specific studies (e.g. [START_REF] Winkler | Tales of the unexpected: phylogeography of the arctic-alpine model plant Saxifraga oppositifolia (Saxifragaceae) revisited[END_REF]. We used Bayes factor (BF) analysis to identify well-supported geographic state transitions with strong posterior support (i.e. BF values >3, [START_REF] Lemey | Bayesian phylogeography finds its roots[END_REF].

The maximum clade credibility (MCC) tree produced by BEAST 1.7 was used to map the ancestral distributions reconstructed for each node via the BSSVS method. We ran the MCMC analysis for 10 million steps, sampling states every 1000 generations. The coalescent model for molecular data was set with a GTR + gamma substitution model, an empirical estimation of nucleotide frequencies, a coalescent-constant size tree prior, a random starting tree, a strict molecular clock and default program prior distribution on other model parameters. The coalescent model for the discrete geographical data was set with a symmetric substitution model and the BSSVS was then enabled. Convergence and acceptable mixing of the sampled parameters were checked using TRACER (part of BEAST) (supplementary material Fig. S1); all ESS values were higher than 200. TREEANNOTATOR (part of BEAST) was used to choose the MCC tree after a 10% burn-in.

The most probable location of each node was labelled on the MCC tree with the corresponding geographic state. For the main phylogroups and the root nodes, the distribution of posterior probabilities of geographic states was also shown to account for uncertainties in ancestral reconstructions (Fig. 3). The history of migration events within the distribution of A. provincialis was synthetized on Figure 4 according to ancestral geographic state reconstruction. [START_REF] Youssef | Factors Underlying the Narrow Distribution of the Mediterranean Annual Plant Arenaria provincialis (Caryophyllaceae)[END_REF] found that the flora of open rocky habitats in Southern Provence is organized along two ecological gradients, according to an Outlying Mean Index analysis (OMI; [START_REF] Dolédec | Niche separation in community analysis: a new method[END_REF] based upon 623 floristic relevés and 14 variables : altitude, slope, exposition, mean annual temperature and precipitation (BIOCLIM variables; [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF], cover percentage of high and low shrubs, high and low trees, herbaceous plants and total plant cover, as well as the proportion of rocks, pebbles and gravel. The first gradient is mainly climatic, related to altitude and distance from the coastline, from thermophilous coastline habitats to colder and wetter habitats located at high altitude in the study area. The second gradient is in vegetation cover, from scrublands to scree slopes, measured within a 100 m² circular area. This vegetation structure is mainly determined by substrate granulometry and slope steepness.

Correlation between phylogeography and ecological gradients

In our analysis, each location, corresponding to a leaf of the MCC tree of cpDNA sequences, was reported with respect to its contribution to the two ecological gradients described above using the dotchart.phylog function ("adephylo" and "ade4" packages; R Development Core Team, 2011). The contributions were centred; therefore the mean value on each gradient corresponded to the barycentre of the ecological distribution of A. provincialis. Finally, the phylogenetic auto-correlation on the two gradients (N=94) was tested using the Abouheif test with the abouheif.moran function ("adephylo" package, R Development Core Team, 2011), based on the Moran statistic and a randomization procedure (1000 permutations) [START_REF] Pavoine | Testing for phylogenetic signal in phenotypic traits: new matrices of phylogenetic proximities[END_REF].

RESULTS

Abundance variation across the geographical range

The density graph obtained via the Kernel method shows, in different intensities of grey, the smoothing of occurrence density inside a 500x500 metre plot (Fig. 1b), with a maximum of 17 occurrences. Density is higher and more continuous in the south-western part of the distribution and lower, as well as more fragmented, in the north-eastern part. The core of abundance is not in line with the geographic barycentre, but located in the south-western part (mainly B area). Moreover, abundance decreases on the northern and eastern margins but not on the southern and western margins.

Variation in AFLP across the geographical range

When AFLP data were collated according to the 5 geographical units (Table 3), the C area was distinguished as having the highest number of private markers (13) and the highest Shannon entropy index (0.33). The B area was also characterized by high diversity, while the lowest diversity was observed in the A area for all indexes (Table 3). Medium diversity was observed in both areas D and E. Genetic differentiation, estimated by mean pairwise AFLP genotype divergence within each geographical unit, was lowest in the C area.

Phylogeography of A. provincialis

When the four cpDNA sequences for A. provincialis (Table 2) were combined, the alignment reached 2,151 bp in length. Twenty eight haplotypes (H1 to H28; Table 1) were identified for the 94 samples of A. provincialis from 59 variable sites (51 parsimony informative positions), and 5 indels. The trnL-trnF DNA region presents the highest nucleotide diversity (0.011), in spite of having the fewest number of nucleotides (267 bp) and haplotypes (only 9).

Conversely, the longest region, trnK-matK (830 bp), is characterized by the highest number of haplotypes ( 16) and the lowest nucleotide diversity (0.003).

NeighborNet analysis revealed two main phylogroups within A. provincialis. The outgroup taxa (A. cinerea, A. cilliata and A. hispida) are connected to the branch separating these two main phylogroups (Fig. 2). The geographical distribution of cpDNA haplotypes revealed a strong phylogeographical signal: the first phylogroup (H1 to H15) was mainly situated in the northern part of the distribution range (groups D and E), while the second phylogroup (H15 to H28) was restricted to the south (areas A, B and C). The two phylogroups are separated by the Huveaune river valley, but a few haplotypes of the first phylogroup were also found in the southern part of the distribution range (group C, haplotypes H5, H11 and H15).

The MCC tree provided by BEAST and based on the 94 samples (Fig. 3) is congruent with the result of the NeighborNet analysis. The root of the MCC tree is located on the branch separating the two main phylogroups. The BSSVS method used to reconstruct the ancestral areas of Arenaria provincialis identified area C as ancestral (highest posterior probability for the root, Fig. 3). After the initial divergence of the two phylogroups, the history of A.

provincialis was characterized by migrations to E and B areas. The first phylogroup originated from area E and migrated eastwards to area D. The second phylogroup originated from area B, migrated westwards to area A and then migrated back to area C. There were at least three dispersal events back towards former areas (e.g. from D to E (H4, H5), from E to C (H5), or from B to C; Fig. 3 and Fig. 4), but no migration event from the southern areas to the north occurred after the divergence of the phylogroups. BF analysis revealed that the geographic state transitions with strongest posterior support were between C and B (BF=24) and between C and E (BF= 22).

Correlation between ecological and phylogenetic structures

The Abouheif tests exhibit significant phylogenetic autocorrelation on the two ecological gradients provided by OMI analysis for all haplotypes (p value = 0.001) as well as for the two phylogroups separately (p value = 0.001). Descriptive results of the Abouheif tests are given as supplementary material (Table S4). Fig. 5 shows that samples closely related in phylogeny tend to share the same position in the ecological gradients.

Across the first gradient, samples from sites with above average climatic conditions (first OMI axis) are mainly in the first phylogroup, while those below average are mainly in the second phylogroup. In the former areas of both phylogroups, i.e. C and E (H16, H17, H18, H1, H2, and H3), climatic conditions differed but were not extreme. The most marginal values, corresponding to either the coldest and wettest or the warmest and driest climatic conditions, were for locations at the periphery, where the most recent haplotypes were found (D versus A areas) and for which migrations from former areas (B, C and E) were inferred (see above results).

Across the second gradient, that of vegetation cover, the distribution of samples was also phylogenetically structured. Samples belonging to the second phylogroup are mainly distributed in sites characterized by low vegetation cover, high slope values and high gravel cover values (corresponding to screes), whereas samples belonging to the first phylogroup are mainly distributed in sites with low vegetation and low gravel cover values and high rock cover values (representing lapiaz). Again, the most marginal values were observed in areas A, B and D, while values for the C and E areas were at the centre of the gradient.

DISCUSSION

A hotspot of biodiversity, the Mediterranean Basin has been described as a mosaic of refugial areas harboring numerous endemic taxa [START_REF] Médail | Glacial refugia influence plant diversity patterns in the Mediterranean Basin[END_REF]. But few studies have investigated precisely the phylogeography of narrow endemic plants, especially outside the major regional Mediterranean biodiversity hotspots [START_REF] Médail | Hot-spots analysis for conservation of plant biodiversity in the Mediterranean basin[END_REF]. Furthermore, annual plants are generally not well represented in phylogeographical studies. Therefore our study focusing on a winter annual plant of southern Provence, a region that remains relatively unexplored in molecular studies, represents a contribution towards more comprehensive knowledge of plant evolution in the Mediterranean region.

Unexpected chloroplast DNA variation

In a recent study of its ecological niche and phylogenetic position [START_REF] Youssef | Factors Underlying the Narrow Distribution of the Mediterranean Annual Plant Arenaria provincialis (Caryophyllaceae)[END_REF], Arenaria provincialis was shown to be the only one of the six Arenaria species of Provence able to colonize scree slopes. This suggests a probable relationship between its history and the emergence of lapiaz and scree slopes derived from the intense frost erosion processes occurring during the Pleistocene glacial periods.

Here, NeighborNet analysis based on the combined sequences of four cpDNA loci detected a well-resolved genealogy with two phylogroups in A. Provincialis, themselves connected to the outgroups A. cinerea, A. hispida and A. ciliata (Fig. 2). After the divergence of the two main phylogroups of A. provincialis, the cpDNA haplotypic diversity accumulated, leading to the current high numbers of haplotypes ( 28) observed within a small area (less than 145 km²).

In comparison, a study based on Arenaria humifusa, endemic to the Arctic Arc, showed only two substitutions and one indel for 1,578 bp in two cpDNA loci, revealing three haplotypes distributed over more than 5,000 km² [START_REF] Westergaard | Glacial survival may matter after all: nunatak signatures in the rare European populations of two west-arctic species[END_REF].

Our study, like other Mediterranean case studies [e.g. in insular endemics, [START_REF] Molins | Molecular phylogeography of Thymus herba-barona ( Lamiaceae ): Insight into the evolutionary history of the flora of the western Mediterranean islands[END_REF][START_REF] Mayol | Patterns of genetic variability and habitat occupancy in Crepis triasii (Asteraceae) at different spatial scales: insights on evolutionary processes leading to diversification in continental islands[END_REF], confirms that strong genetic structuring is plausible in any part of the distribution range of Mediterranean endemics, and therefore that the size of a distribution range is not correlated with levels and structure of genetic diversity. The well-known climatic impact of the last glacial maximum (21,000 years bp) on species distributions [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF] has also been shown for vegetation in the western Mediterranean area [START_REF] Cheddadi | Putative glacial refugia of Cedrus atlantica deduced from Quaternary pollen records and modern genetic diversity[END_REF][START_REF] Médail | Glacial refugia influence plant diversity patterns in the Mediterranean Basin[END_REF]Fady and Connors, 2010). However, the high cpDNA diversity and polymorphism of A. provincialis we observed in such a restricted area support the idea that the preservation and accumulation of genetic diversity during the Ice Ages could be a common pattern in the Mediterranean, as underlined recently by Nieto [START_REF] Feliner | Southern European glacial refugia: A tale of tales[END_REF]. Our study yields new evidence that extreme or azonal habitats, i.e. rocky terrain, cliffs and steep slopes, may have promoted the persistence of endemic plants and enhanced specific traits [START_REF] Thompson | Ecological differentiation of Mediterranean endemic plants[END_REF]. Arenaria provincialis, with its phylogeographic structure, therefore affords an excellent opportunity to relate long-term persistence and patterns of genetic diversity to the central-marginal hypothesis [START_REF] Eckert | Genetic variation across species' geographical ranges: the central-marginal hypothesis and beyond[END_REF][START_REF] Sexton | Evolution and ecology of species range limits[END_REF].

Spatial distribution and the central-marginal hypothesis

The kernel model of density (Fig. 1b) clearly shows that abundance of Arenaria provincialis is not homogeneous across the geographical range. Area B has the highest abundance. In the northern part of its distribution, A. provincialis is markedly less abundant and more fragmented. A bias in field inventories is unlikely to explain this pattern, since inventories of A. provincialis have been made for several years and its northern and eastern range limits have been intensively investigated. The north-eastern parts of the distribution correspond to isolated patches of open rocky habitats most frequently embedded in mattoral or forest matrixes. A decrease in abundance is observed towards the northern and eastern range limits.

Genetic diversity showed a different pattern from that of contemporary abundance: measured by the polymorphism and Shannon entropy of AFLP markers (Table 3), it increased towards the centre (areas B and C), whereas AFLP genotype divergence, used here to describe genetic differentiation, was significantly lower in area C (Table 3). Molecular markers thus show that the populations of A. provincialis with the highest diversity and that are the least differentiated are near the central part of the distribution (mostly C and secondary B areas), whereas marginal populations with the lowest genetic diversity are in areas A, D, and E, near the western, eastern and northern parts of the distribution respectively. The structure of genetic diversity seems roughly consistent with the hypotheses formulated here (see introduction), and we cannot rule out the contemporary effects of isolation and small effective population size, notably at the northern and eastern limits. However, we make the point below that this pattern of contemporary genetic diversity is also well correlated with Pleistocene range dynamics and expansion.

Past migrations from the centre to the periphery

The structured phylogeography of A. provincialis and the reconstruction of ancestral areas, based on the BSSVS method of [START_REF] Lemey | Bayesian phylogeography finds its roots[END_REF], enabled us to reconstruct a step-by-step history of the migrations responsible for the current range of A. provincialis. The ancestral area of A. provincialis was thus found to be situated near the centre of the distribution (C area), with evidence that after migrations to the west (B area) and north (E area), migrations occurred toward the eastern and western periphery (A and D areas, Fig. 3 and Fig. 4).

Moreover, phylogeography and the reconstruction of ancestral states revealed that backward migrations also occurred to C and E areas (Fig. 3 and Fig. 4). BF analysis revealed strong posterior support for dispersal between C and B and between E and C areas, thus confirming central position C as both an ancestral and a crossroads area. In this phylogeographic context, it should be noted that the frequency of AFLP private markers (Table 3) supported a longer persistence of A. provincialis in area C, and secondarily in area B. Consequently, area C, which is central in the distribution, appears to have been crucial in the history of A.

provincialis.

The history of A. provincialis therefore reveals that the central-marginal pattern is strongly associated with a history of expansion from former areas (Fig. 4). Such correlations between genetic diversity and distance from refugia have been found elsewhere (e.g. [START_REF] Garner | Genetic diversity across a vertebrate species' range: a test of the central-peripheral hypothesis[END_REF][START_REF] Hoban | Range-wide distribution of genetic diversity in the North American tree Juglans cinerea: a product of range shifts, not ecological marginality or recent population decline[END_REF], but over a coarser scale and in the context of glaciations. To the best of our knowledge, this correlation has never before been observed for a Mediterranean narrow endemic plant. In the light of the extensive cpDNA nucleotidic variation and marked phylogeographic structure revealed here (28 haplotypes in two phylogroups), it is likely that the rocky habitats situated in the central parts (C, and thereafter B and E areas) of the distribution range have sheltered A. provincialis for a long time (i.e. probably since the early Pleistocene). Originating from these central-ancestral areas, the populations situated today at the range limits (A and D areas) are younger and have more limited AFLP diversity (Table 3).

An extensive survey of the literature yielded no conclusive evidence from recent geological history of the focus area that might have helped to calibrate precisely in time the haplotype genealogy of A. provincialis. Regarding past environmental change, paleoecological data are also very scarce for rocky outcrops of the focus area, and ancestral distribution modelling would be difficult to realize for a narrow endemic plant living in azonal habitats. Of the plants common to the rocky habitats of A. provincialis, the pioneer tree Juniperus phoenicea has been hypothesised to have undergone expansion during arid and cold periods of the Pleistocene, such as the Last Glacial Maximum [START_REF] Boratynski | High level of genetic differentiation of Juniperus phoenicea (Cupressaceae) in the Mediterranean region: geographic implications[END_REF]. Similarly, it has been proposed that Cheirolophus intybaceus reached southern France from the Iberian Peninsula during the mid-late Pleistocene [START_REF] Garnatje | Balearic insular isolation and large continental spread framed the phylogeography of the western Mediterranean Cheirolophus intybaceus s.l. (Asteraceae)[END_REF]. The sole relevant hypothesis, therefore, is that a decrease in vegetation cover due to climate cooling and drying during the major glaciation phases beginning around 0.9 Mya (De [START_REF] Beaulieu | Vegetation dynamics in northwestern Mediterranean regions: Instability of the Mediterranean bioclimate[END_REF] may have favored the expansion of pioneer plants like A. provincialis.

Conservation significance of the phylogeographic pattern

The habitat of A. provincialis is structured around two main ecological gradients mainly determined by altitude and vegetation cover [START_REF] Baumel | Ecological magnitude and fine scale dynamics of the Mediterranean narrow endemic therophyte, Arenaria provincialis (Caryophyllaceae)[END_REF][START_REF] Youssef | Factors Underlying the Narrow Distribution of the Mediterranean Annual Plant Arenaria provincialis (Caryophyllaceae)[END_REF]Fig. 5A). As supported by the significant phylogenetic autocorrelation on these two ecological gradients (Abouheif test), the populations having the most divergent haplotypes are now living at the extremes of its two ecological gradients (Fig. 5C). During the history of its expansion, A. provincialis reached the top of its elevation gradient in southern Provence (i.e. a range of 1,000 metres), where its populations are currently experiencing contrasting climatic conditions (Axis 1, Fig. 5A and Fig. 5C). For example, A. provincialis is the sole Arenaria species living both near the littoral and near the summit of the Sainte Baume mountain range (almost 1,000m above sea level), in contrast to its congeners A. grandiflora and A. aggregata which are present only in the second location. Yet populations located near the littoral, and belonging mainly to the second phylogroup, are living in the most open and steep habitats of the species, such as huge south-facing seaside screes (Axis 2, Fig. 5 A,C). This means that today's populations of A. provincialis are organized across an ecological range explained both by the history of the species and by the geographical structure of its genetic diversity.

In conclusion, population distinctiveness based on molecular variation and ecological data is well supported here: populations from the two phylogroups are not genetically, nor ecologically, exchangeable [START_REF] Crandall | Considering evolutionary processes in conservation biology[END_REF]. Further study will be needed in order to evaluate the probable adaptive polymorphism of A. provincialis suggested by the strong correlation between phylogeography and ecological gradients [START_REF] Moritz | Strategies to protect biological diversity and the evolutionary processes that sustain it[END_REF]. Our results point to key issues in terms of conservation. Firstly, the two phylogroups of A. provincialis have different ecologies, and exhibit different patterns of fragmentation and abundance (Fig. 1, Fig. 5). The populations situated near the south-western and north-eastern range limits contribute significantly to this fragmented structure. Secondly, the central populations have persisted the longest and been the source of expansions, but they are not currently characterized by the highest abundance.

Exploring the abundance, genetic diversity and phylogeography of a Mediterranean narrow endemic plant, this study highlights the fact that current range size and abundance patterns are not sufficient to predict the organization of genetic diversity. However, our study reveals the crucial role played in A. provincialis' history by the areas centrally situated in its distribution. Haplotypes are clustered by geographical groups (from A to E). For population codes see Table 1. Chloroplast DNA haplotypes detected (with the number of individuals carrying a given haplotype in brackets) for each of the 40 populations of Arenaria provincialis sampled relative to geographical group memberships, codes, GPS coordinates (decimal degrees) and altitude (geographical location of haplotype in Fig. 1).

FIGURES AND TABLES

Table 2. Characterization and diversity of the four chloroplast DNA sequences separately and for the combined matrix (Total) examined for Arenaria provincialis.

Fig. 1 .

 1 Fig. 1. (a) Geographical distribution of the 28 haplotypes detected for the 94 individuals of Arenaria provincialis combining four cpDNA regions (internal part of matK gene, trnK intron plus 5' part of matK, tnrL-trnF and trnT-trnL intergenic spacer for a total of 2,150 bp).

  Abundance of Arenaria provincialis estimated by an isotropic gaussian kernel function. 591 occurrences were used, at least 100 metres apart (see Materials and Methods for more details). The maximum abundance was 17 occurrences within a 500x500 metre plot.

Fig. 2 .

 2 Fig.2. NeighborNet diagram generated using SplitsTree4 connecting the 28 haplotypes detected for the 94 individuals of Arenaria provincialis combining the polymorphism of four cpDNA locus (internal part of matK gene, trnK intron and 5' part of matK, tnrL-trnF and trnT-trnL intergenic spacer for a total of 2,150 bp). Each circle represents a haplotype, with size proportional to number of individuals (see table1).

Fig. 3 .

 3 Fig. 3. Maximum clade credibility tree based on sequence variation at four cpDNA loci of Arenaria provincialis and reconstruction of ancestral areas according to BSSVS method. Nodes are named according to the highest posterior probability of location state. Node histograms for the three principal nodes show posterior probability distributions of location state.

Fig. 5 .

 5 Fig. 5. (a) Projection of the first two axes of the OMI analysis following Youssef et al. (2011) realized on 14 variables (described in the text) and 254 stations.

Table 1 .

 1 Crosses represent the 156 occurrences of Arenaria provincialis used inYoussef et al.

Table 3 .

 3 Summary statistics of AFLP data for the geographical groups of populations of Arenaria provincialis. Number of investigated individuals (N), number of AFLP fragments (Frag tot ) per population, number of fragments that only occur in one population (private), percentage of fragments exhibiting intrapopulational polymorphism (%poly), Shannon diversity index (Shan), and genetic differentiation (Diff) are listed.

Table 1 .

 1 Table 2.

		fragment	haplotypes	polymorphic	haplotype	nucleotide
		lenght	number	sites	diversity (h)	diversity (π)
	MatK	711	14	19	0.90		0.00787
	trnK-matK	830	16	12	0.88		0.00297
	trnL-trnF	267	9	12	0.85		0.01099
	trnT-trnL	359	13	8	0.77		0.00613
	total	2150	28	51	0.94		0.00626
	Table 3.
		A	B	C	D		E
	N	11	24	12	16		10
	Frag tot	133	212	215	186		160
	Private	1	9	13	2		4
	% poly	50,2	81,67	80,48	72,11		61,35
	Shan	0,20±0,02	0,28±0,01	0,33±0,02	0,26±0,01	0,25±0,02
	Diff	0,697±0,1	0,703±0,1	0,644±0,08 0,716±0,09	0,701±0,1

Acknowledgements: We thank the Conservatoire Botanique National Méditerranéen de Porquerolles for providing unpublished data about Arenaria taxa locations and allowing us to collect Arenaria material. We are grateful to Errol Véla for his help on the field and botanical assessment. Our thanks to Marjorie Sweetko for English language editing. This research was funded by the CNRS, the French National Office of Forests (ONF) and the General council of Bouches du Rhône district (CG13). Marine Pouget was supported by the Provence-Alpes-Côte d'Azur region (PACA).

Table S1. Primers and references for amplification of chloroplast loci.

Table S2. Genbank accession number of the 28 haplotypes of Arenaria provincialis and Outgroup.

Table S3. Cycling condition of the PCR.

Table S3.