Population genetic structure and management perspectives for Armeria belgenciencis, a narrow endemic plant from Provence (France)

1Institut Méditerranéen de Biodiversité et d’Ecologie marine et continentale (IMBE), Aix Marseille Univ, Avignon Université, CNRS, IRD, Faculté des Sciences et Techniques St-Jérôme, Marseille, France 2Conservatoire botanique national méditerranéen de Porquerolles (CBNMed), 34 avenue Gambetta, 83400 Hyères, France 3Conservatoire d’Espaces Naturels (CEN), L’Astragale, 888 chemin des Costettes, 83340 Le Cannet-des-Maures, France *Corresponding author: alex.baumel@imbe.fr REGULAR PAPER


INTRODUCTION
Conserving threatened species requires to properly design conservation units. Different strategies can be applied to conserve rare endemic plants (e.g., Van Rossum et al. 2017), but all need to be based on a robust knowledge and delineation of the entities targeted for in or ex situ conservation (Ryder 1986;Médail & Baumel 2018;Levin 2019). Delineating conservation units is never straightforward, and it often requires the judicious use of available data. However, data are often too restricted to support the exclusiveness or distinctiveness of population sets (Crandall et al. 2000). In the case of incipient species, insufficient isolation and divergence may have precluded the formation of morphologically discrete and reciprocally monophyletic entities for the molecular loci considered. In such cases, combining systematics with population genetic approach is crucial to consider the balance between recognizing ongoing diversification (Moritz 2002;Daïnou et al. 2014) and the risk of over-splitting taxa into populations that may compromise both management and rescue, and ultimately their long-term persistence (Frankham et al. 2012).
The genus Armeria (Plumbaginaceae), known as "thrift", exemplifies this scenario. Integrated by about eighty small perennial polycarpic mostly Mediterranean species, its centre of diversity is in the Iberian Peninsula (Nieto Feliner 1990). Molecular systematic studies in Armeria showed the recurrence of hybridization due to the weakness of reproductive barriers and frequent range fluctuation (Fuertes Aguilar et al. 1999;Tauleigne-Gomes & Lefèbvre 2005;Piñeiro et al. 2011;Nieto Feliner et al. 2019). Therefore, the delineation of taxa within this genus is based on a correlation criterion between ecological, geographical and morphological isolation which often results in taxa consisting of few populations (Nieto Feliner 1990). Due to frequent introgressive hybridization (Nieto Feliner et al. 2019) and sometimes also the lack of information on DNA barcode markers, population genetic approaches based on multilocus markers are highly relevant to set conservation units since they are expected to provide a finer reconstruction of the evolutionary history of Armeria taxa.
We report here the results of a study on the genetic diversity structure of Armeria belgenciensis Donadille ex Kerguélen and geographically closely related populations of Armeria arenaria (Pers.) Schult. Armeria belgenciensis is a narrow endemic represented by only one population, which recently experienced a strong demographic decline caused by repeated disturbances in its habitat and was the object of a demographic rescue. Within two kilometres of A. belgenciensis the presence of Armeria arenaria subsp. peirescii Baumel, Auda & Médail represented also by one small population, poses the question of whether the two populations are genetically isolated and should be considered as separated taxa for conservation. According to the literature (Tison et al. 2014), A. arenaria subsp. peirescii flowers one to two months earlier than A. belgenciensis. However, their spatial proximity, strong morphological similarity and sharing of the same ITS ribosomal DNA genotype suggested hybridization between them (Baumel et al. 2009).
We applied a population genetic approach to evaluate the isolation between A. belgenciensis and A. arenaria subsp. pereiscii. Considering that Armeria taxa can exchange genetic material when they meet, our approach is based on the assumptions that genetic differentiation among taxa is expected to be higher than within taxa and their isolation has to be supported by their genetic clustering. Our sampling was therefore designed to estimate (1) genetic admixture between A. belgenciensis and A. arenaria subsp. pereiscii and (2) genetic differentiation among them and close neighbour populations belonging to A. arenaria (A. arenaria subsp. bupleuroides Godron & Gren. and A. arenaria subsp. pradetensis Médail, Baumel & Auda). We also monitored flowering phenology of A. belgenciensis and A. arenaria subsp. peirescii to test for their hypothetical pre-zygotic isolation.
Our genotyping method is based on Amplified Fragment Length Polymorphism (AFLP) because AFLP markers have a higher genome coverage than single locus DNA barcode markers while having a lower cost than genotyping by sequencing methods. Usefulness of AFLP markers for phylogeography or systematics of recently diverged taxa has been proven in several case studies (Despres et al. 2002;Duminil & di Michele 2009;Hardion et al. 2012;De Riek et al. 2013;Medrano et al. 2014;Gutiérrez-Larruscain et al. 2019).
We discuss the implications of our results and propose practical recommendations for the conservation of these two threatened taxa.

Studied taxa
Armeria belgenciensis is known from a single location in southern Provence (Solliès-Toucas, Var, France). This population experienced a strong demographic bottleneck from about 1000 individuals in the 1960s (Donadille 1969) to less than 150 individuals in 2004 and 34 in 2007 (Conservatoire botanique national méditerranéen de Porquerolles (CBNMed) unpublished data). A population rescue implemented by CBNMed and Conservatoire d'Espaces Naturels Provence-Alpes-Côte d'Azur (CEN-PACA) since 2009 has helped to restore the census size to nearly 458 individuals in 2019. Within 2 km of the A. belgenciensis population, one small population of A. arenaria was discovered in 2002 and described in 2009 as A. arenaria subsp. peirescii (Baumel et al. 2009). In 2019, the census size of A. arenaria subsp. peirescii was 158 individuals. Armeria belgenciensis is located on a dolomitic sandy plateau with an important tree cover mainly consisting of Pinus pinaster with P. halepensis and Quercus ilex. The P. pinaster trees are at most 50-60 years old and they regenerate. The area belongs to seven landowners and there is no planned forest management. In 2009 the area with A. belgenciensis has been placed under legal protection ("Arrêté préfectoral de protection de biotope") and is regularly monitored by CBNMed. Small clearings were cut to remove some pines around the main stand of A. belgenciensis and their needles were raked. Armeria belgenciensis is located in clearings and forest edges, with other psammophilous species such as Arenaria modesta, Helichrysum stoechas, Iberis ciliata, Iberis linifolia, Silene otites or Eu-phorbia seguieriana. In contrast, A. arenaria subsp. pereiscii grows in an open steep rocky habitat within a vegetation dominated mainly by Quercus coccifera, Rosmarinus officinalis and Thymus vulgaris.

Sampling
Leaves from two hundred individuals were collected from six populations ( fig. 1). According to our first objective described in introduction, the main sampling effort concerned A. belgenciensis and A. arenaria subsp. pereiscii and thus 50 individuals were sampled for each of those two populations (BEL and PER) allowing to comprehensively sample their spatial distribution for detecting admixed individuals (supplementary file 1). According to our second objective the sampling was completed with 50 individuals from two nearby populations of A. arenaria subsp. bupleuroides (situated at "Aiguilles de Valbelle", named VAL and "Méounes-les-Montrieux" named MEO), 25 individuals from a population of A. arenaria subsp. pradetensis situated near the Mediterranean coastline (situated at "Le Pradet" named PRA) and 25 individuals of A. arenaria subsp. bupleuroides (situated at "Gigaro" named GIG) from an isolated population, 52 km away from A. belgenciensis. The sampling coordinates are included in the data sets deposited in Dryad Digital Repository.

DNA extraction, AFLP and genetic diversity analysis
Total DNA was extracted using NucleoSpin Plant II Kit (Macherey & Nagel, Germany) from c. 15 mg of silica dried leaves. DNA concentrations were measured using a photometer (Biophotometer, Eppendorf, Germany) and homogenized at 10 ng/µl. The AFLP protocol was conducted on 100 ng of DNA according to Vos et al. (1995) with modifications (see supplementary file 2). Electrophoresis was done using an ABI 3730xl DNA Analyzer.
The marker bands produced by the four primer combinations were entered automatically using the RawGeno application (Arrigo et al. 2009) with stringent criteria to select markers (minimum size peaks greater than 100 UFR, minimum width of 1.2 bp and maximum of 1.75 bp). Despite the stringency of this automatic procedure, it generates an error rate that may be greater than a manual selection. To reduce this error rate, reproducibility of the markers was tested on 12 samples. Markers not reproducible in 100% of the 12 replicates, i.e., having an error rate equal to or greater than 8.33%, were removed. In addition, markers with a frequency lower than 5% or higher than 95% were also removed. Individuals with very few markers, for which the AFLP procedure was not optimal, were removed as well. Pl. Ecol. Evol. 153 (2), 2020 The genetic diversity structure was investigated using four different methods. First, STRUCTURE (Pritchard et al. 2000) was used to investigate clustering models. From one to 10 genetic groups (K = 1-10) were tested allowing admixture at the individual level, with correlated allele frequencies and without prior information on geographical origin. Ten replicate runs per K were performed, each having a burn-in period of 200 000 generations and a chain length of 1 000 000 generations. The results of the 10 replicates were averaged with CLUMPACK (Kopelman et al. 2015) for each K model. STRUCTURE HARVESTER (Earl & vonHoldt 2012) was used to obtain likelihood values across the multiple values of K as well as to apply the delta K criterion to select the K model representing the optimal partition of the genetic data set (supplementary file 3).
Second, we used the model-free Discriminant Analysis of Principal Components (DAPC) implemented in the adegenet R package (Jombart et al. 2010). The DAPC uses a k-means algorithm, principal components and discriminant analysis to identify the best clustering of the data set. For DAPC the choice of the K clustering solution is based on the Bayesian information criterion measure of goodness of fit (R package adegenet, find.cluster and dapc functions; Jombart et al. 2010). Because the challenges for estimating the optimal K clusters with STRUCTURE (Janes et al. 2017) also apply for DAPC, all the results of STRUCTURE and DAPC from K = 1 to K = 6 are shown in supplementary files 4-6. The retained STRUCTURE and DAPC clusters were depicted using the compoplot function of the adagenet R package. The main pattern of genetic structure revealed by STRUCTURE and DAPC was used to perform a principal coordinate analysis (PCoA) based on inter-individual distances computed with the Dice coefficient (dist.binary function, ade4 R package; Dray & Dufour 2007).
Third, we used the Ward hierarchical clustering method (R package, hclust functions with the correction of Murtagh & Legendre 2014). This last method uses the inter-individual distances computed with the Dice coefficient to generate a dendrogram showing groups that minimize within-group dispersion. Finally, an AMOVA was performed to estimate overall and pairwise differentiations among populations with GeneAlex 6.5 (Peakall & Smouse 2012).
Within population genetic diversity was analysed to estimate the unbiased Nei genetic diversity (HE), the Shanon index (SH) and the standardized index of association (RBARD) with the poppr R package (Kamvar et al. 2014).

Flowering phenology monitoring
Occasional field observations reported a summer flowering phenology for A. belgenciensis and a spring phenology for A. arenaria subsp. peirescii (Baumel et al. 2009;Tison et al. 2014

Genetic diversity structure
After the filtering steps of AFLP raw data, 328 markers were selected for a total of 187 individuals. To identify the main pattern of genetic structure we used the models of STRUC-TURE assigning the 187 thrift genotypes to genetic clusters. According to the delta K criterion (supplementary file 3), the optimal solutions were the K = 2 and K = 4 models ( fig. 2; but see all in supplementary file 4). For the K = 2 model, BEL and PER populations constituted the first group and the remainder the second group. For the K = 4 model, BEL and PER populations are split to form the two first groups whereas VAL, MEO and PRA constituted the third group and GIG the fourth. The BIC measure of goodness of fit (supplementary file 5) identified six as the optimal number of groups for the DAPC analysis. These six groups corresponded to the six sampled populations albeit with high admixture between VAL and MEO. The STRUCTURE patterns are congruent with the results of DAPC for two and six groups, respectively ( fig. 3). The DAPC revealed a higher level of admixture between BEL and PER than STRUCTURE. The Ward algorithm generated a dendrogram with a hierarchical pattern of clustering that is similar to the STRUCTURE and DAPC results ( fig. 2). The PCoA of genotypic dissimilarities revealed a pattern of individual grouping that is congruent with STRUCTURE and DAPC admixture plots ( fig. 3). The mean PhiST was 0.2 and the minimum pairwise PhiST was 0.13 between VAL and MEO, followed by the PhiST between BEL and PER, whereas the two highest values, 0.4 and 0.42, were obtained between BEL and GIG (table 1).
Overall patterns of genetic structure indicate a close relationship between the BEL and PER populations, which form a distinct genetic group from other populations. Moderate differentiation (PhiST = 0.14), as well as multivariate methods DAPC and PCoA, supported the existence of genetic exchanges between BEL and PER although STRUCTURE identified a lesser extent of admixture ( fig. 2). According to the multivariate methods, the gene flow is symmetrical with the same proportion of introgressed individuals in both populations whereas STRUCTURE supports an asymmetrical introgression of BEL into PER.

Genetic diversity
Within-population genetic diversity exhibited similar values in the six populations (table 2). HE ranged between 0.2 and 0.26 and SH between 3.14 and 3.91. Private markers to sin-  gle populations were relatively few. The normalized index of association (RBARD) was significant and positive but small in all populations.

Flowering phenology
BEL showed a higher number of capitula and flowers per blooming individual. The mean number of capitula per plant was 2.5 (SE 0.13) for BEL versus 2 (SE 0.09) for PER. The mean number of flowers by plant was 8.5 (SE 0.65) for BEL versus 4.9 (SE 0.34) for PER. BEL has also larger confidence intervals for these variables (fig. 4). The flowering period of PER started earlier than BEL and had a shorter duration ( fig. 4). Over the three years of flowering censuses, we observed that PER started blooming in May whereas BEL started in June. The blooming peak was in the beginning of June for PER and in the beginning of July for BEL. In 2017 the frequency of visits was not sufficient to observe flowering overlap between the two populations. However, in 2018 and 2019 an overlap was clearly observed between the end of the PER blooming period and the beginning of the BEL blooming period. This overlap, occurring at the end of June/beginning of July, lasted 8 days in 2019 (from June 24 to July 1) and 15 days in 2018 (from June 8 to June 22).

DISCUSSION
Since the Mediterranean Basin hotspot includes many rare and vulnerable narrow endemic plant species, the number of available molecular systematic studies is not sufficient to understand the complexity of the evolution of many groups and thus to properly design conservation units (Médail & Baumel 2018). Our multilocus molecular approach was able to distinguish a hierarchical genetic differentiation structure with taxonomic and conservation implications.

Genetic structure and taxonomic implications
The main result of this study based on 328 AFLP markers is that different methods of clustering revealed a close relationship between Armeria belgenciensis (BEL) and A. arenaria subsp. peirescii (PER) as suggested by a previous morphometric study (Baumel et al. 2009). The fact that GIG, the geographically most distant population, was grouped with VAL, MEO and PRA in all analyses stresses the distinctiveness of BEL and PER, even though the latter was described as a subspecies of A. arenaria. BEL and PER must therefore be conserved independently of A. arenaria and we propose to gather these two populations under A. belgenciensis. This proposal has the advantage of establishing a taxonomical scheme based on a robust analysis of genetic differentiation and it also implies the legal protection of a second population of A. belgenciensis. VAL, MEO and PRA also form an original genetic group despite the distance that separates these populations. Within this group, the PRA population (A. arenaria subsp. pradentensis) is moderately differentiated according to PhiST and forms a genetic group according to the clustering methods. Its genetic isolation is not complete ( fig. 2) but it makes sense to keep it as a separate conservation unit, and possibly maintaining its subspecies status based on three arguments: its geographical isolation, its specific rocky habitat located on the waste of an old copper mine and its different morphology (Baumel et al. 2009). Previous works emphasized the evolution of heavy-metal tolerance in Armeria (e.g., Vekemans & Lefèbvre 1997).
Our results also raise the question of the status of the most remote A. arenaria population (GIG), which appears to be an isolated and independent genetic group in all analyses, and occurs on a distinct semi-halophilous habitat, near the coastline. On the basis of these results we are facing a systematic dilemma. On the one hand, maintaining A. arenaria subsp. pradetensis would require describing the GIG populations as an independent subspecies. On the other hand, if we followed the main split in the Ward dendrogram, or the STRUCTURE solution for K = 2, only two taxa should be recognized: A. belgenciensis and A. arenaria subsp. bupleuroides. Deciding between these two treatments will re- quire an ad hoc study, based on a comprehensive sampling of A. arenaria over its entire range.

Genetic diversity of Armeria belgenciensis
The only known population of Armeria belgenciensis (BEL) has recently experienced a very strong decline (c. 95% of the census size), which has urged a reinforcement procedure to save the population. It was expected that the genetic diversity of this population would be small relative to other Armeria populations. For example, intra-population diversity in Armeria caespitosa (Ortega) Boiss., which is endemic to mountainous areas in central Spain encompassing small populations with often less than 500 individuals, was found to be generally low with mean values of 0.18 for HE and 0.24 for SH (García-Fernández et al. 2013). For Armeria pungens (Link) Hoffmanns. & Link, a coastal species located in Spain, Portugal and Corsica-Sardinia, the diversity was also quite low with mean values of 0.11 for HE and 0.12 for SH (Piñeiro et al. 2007). These values were even lower for Armeria maderensis Lowe (HE = 0.06, SH = 0.11), a rare endemic of the mountains of Madeira (Piñeiro et al. 2011 Baumel 2018). This result indicates that some factors may have limited the impact of the strong recent bottleneck. We suggest the role of the soil seed bank, the rescue performed by the CBNMed and the gene flow from the PER population (see above) as hypothetical mitigating processes.
The RBARD values were significant and positive but small in all populations, rejecting the null hypothesis of no linkage among markers and indicating little deviation from panmictic mixing (Kamvar et al. 2014). Because clonality or selfing are very unlikely in these plants (Baker 1966), this positive and significant values for RBARD are probably caused by a small Wahlund effect. An interesting point is that PER has a lower demographic size, lower genetic diversity (SH), and higher RBARD index. Although the PER population seems stable since its discovery (CBNMed unpublished records), our analyses support a low effective demographic size compared to the greater genetic diversity and higher rate of genetic mixing that characterizes the BEL population.

Implications for the conservation of Armeria belgenciensis
To our knowledge no forest fires occurred in the area during the last 50 years, meaning that the woody vegetation (forest and matorral) could have acted as a physical barrier between BEL and PER for several generations. Incongruences among the clustering methods limit our discussion about the extent and the direction of gene flow. Progress in this direction will require more data and estimating parameters of an isolation by migration model. However, the BEL and PER populations are not strongly differentiated and are probably connected by low gene flow. Their flowering overlap ( fig. 4) is narrow but sufficient to allow cross pollination in early summer by the numerous generalist pollinators regularly observed visiting Armeria (García-Camacho & Escudero 2009;Nieto Feliner et al. 2019). Moreover, a hiking trail passing close to both populations may facilitate ectozoochorous dispersal of the fruits (enclosed in their calyx) between the populations (Tauleigne-Gomes & Lefèbvre 2005). Gene flow, even low, may contribute to their persistence in the long term because the small size of these two populations (458 individuals for A. belgenciensis and 158 individuals for A. arenaria subsp. peirescii) are well under the effective number of 1000 individuals recommended to maintain evolutionary potential (Frankham et al. 2014). A growing body of studies highlights the pivotal role of biological connectivity for the genetic rescue of small fragmented populations as well as to minimize the extinction risk under future environmental change (Frankham 2015;Ralls et al. 2018). Our recommendation is, therefore, to keep these two populations separated as two different management units for seed sourcing or for any further population recovery plan, but to not act against their connectivity, for example by installing fences or removing trails passing nearby. Increasing connectivity between the BEL and PER populations is not pertinent right now but could be considered in the future to prevent the harmful effects of small demographic sizes, i.e., inbreeding, allele effect or limited adaptive potential (Hamilton & Miller 2016