Little hope for the polyploid endemic Pyrenean Larkspur (Delphinium montanum): Evidences from population genomics and Ecological Niche Modeling

Abstract Species endemic to restricted geographical ranges represent a particular conservation issue, be it for their heritage interest. In a context of global change, this is particularly the case for plants which belong to high‐mountain ecosystems and, because of their ecological requirements, are doomed to survive or disappear on their “sky islands”. The Pyrenean Larkspur (Delphinium montanum, Ranunculaceae) is endemic to the Eastern part of the Pyrenees (France and Spain). It is now only observable at a dozen of localities and some populations show signs of decline, such as a recurrent lack of flowering. Implementing population genomics approach (e.g., RAD‐seq like) is particularly useful to understand genomic patterns of diversity and differentiation in order to provide recommendations in term of conservation. However, it remains challenging for species such as D. montanum that are autotetraploid with a large genome size (1C‐value >10 pg) as most methods currently available were developed for diploid species. A Bayesian framework able to call genotypes with uncertainty allowed us to assess genetic diversity and population structure in this system. Our results show evidence for inbreeding (mean G IS = 0.361) within all the populations and substantial population structure (mean G ST = 0.403) at the metapopulation level. In addition to a lack of connectivity between populations, spatial projections of Ecological Niche Modeling (ENM) analyses under different climatic scenarios predict a dramatic decrease of suitable habitat for D. montanum in the future. Based on these results, we discuss the relevance and feasibility of different conservation measures.


| INTRODUC TI ON
Extinction risk associated to climate change is expected not only to increase but also to accelerate for every degree rise in global temperatures (Urban, 2015). In this context, the number of threatened species by reduced population sizes, habitat degradation, and habitat fragmentation is increasing. Small and fragmented populations are particularly prone to be affected by drift and inbreeding depression and generally display decreased levels of genetic diversity (e.g., Bergl et al., 2008;Dixo et al., 2009). The increase of anthropogenic and environmental pressures, both being intimately bound, represents an unprecedented threat for vulnerable species. In particular, populations from species of high-mountain ecosystems might not be able to move upward rapidly enough to keep within their thermal tolerances and might be doomed on their "sky islands" (see McCormack et al., 2009, e.g., De Gabriel Hernando et al., 2021, Kidane et al., 2019. Thus, species endemic to specific mountain ranges represent a major conservation issue, if only for their heritage interest. It is critical to characterize the distribution of genetic diversity within-and between-their populations to evaluate the inbreeding risk and study population structure to better understand how isolated they really are, in order to undertake effective conservation actions (Calevo et al., 2021).
The threatened Pyrenean larkspur Delphinium montanum DC, 1815 perfectly illustrates this situation. This plant (Family Ranunculaceae), endemic to the Eastern part of the Pyrenees (France and Spain), is now only observable at a dozen of localities (two already got extinct over the 20th century) and its global population size was estimated to range between 7800 and 10,300 individuals in total (Aymerich et al., 2020). It occurs between alpine and subalpine levels, between 1600 and 2400 m (Lopéz-Pujol et al., 2007).Delphinium montanum is vulnerable as a result of increasing levels of multiple threats. Some of its populations are located nearby touristic sites (i.e., mountain resorts or hiking trails) and are negatively impacted by trampling (Simon et al., 2001). Flowers and fruits predation by herbivores, particularly Pyrenean chamois/isards (Rupicapra pyrenaica pyrenaica) was also shown to be responsible for the seed bank's decrease in some populations (Aymerich, 2003;Simon et al., 2001). Importantly, differences in population dynamics are observable from one population to another. Some populations are stable, whereas others show decrease in population size (4-Pedraforca) and/or absence of flowering individuals since 2011 (2-Nohèdes). Pioneering work of Simon et al. (2001) and Lopéz-Pujol et al. (2007) also suggested relatively low effective size as well as, evidence of inbreeding within-populations and significant population structure. However, these two studies were based on a limited number of genetic markers (i.e., 7 and 14 allozyme loci, respectively) and populations studied (i.e., 2 and 7, respectively).
Implementing a population genomic approach on D. montanum could be particularly fruitful to inform conservation strategies by providing a substantial increase in the resolution of the genetic marker dataset and genome complexity reduction protocols such as Reduction site-associated DNA sequencing (or RAD-seq, Baird et al., 2008) allow to get such information on tens to hundreds of individuals. However, the features of the genome of the species make Delphinium montanum particularly challenging from a technical and methodological point of view. D. montanum is an autotetraploid species (Lopéz-Pujol et al., 2007;Simon et al., 2001). Our recent flow cytometry estimates of 10.32 pg/1C for this species (i.e., about twice the 1C-value reported for several other diploid Delphinium species) corroborates previous knowledge and confirms the large genome size of D. montanum (Bertrand et al. unpublished data). A large genome size can be considered has a limitation as itself, be it because it requires higher sequencing efforts and costs to carry out conservation studies. The difficulty to infer genotype and allele frequencies on polyploid species with RAD-seq techniques as well as the restricted set of methods available to analyze such data (see Dufresne et al., 2014;Meirmans et al., 2018;Meirmans & van Tienderen, 2013) is also likely to explain the relatively low number of population genomics studies dealing with RAD-seq like techniques on nonmodel polyploid organisms (Ahmad et al., 2021;Clevenger et al., 2015;van de Peer et al., 2017).
Out the handful of studies that have addressed population genomics questions on polyploid organisms without reference genome, some have considered the genotypes as diploid-like data, an assumption that allows to use classical variant calling tools with default settings (e.g., Brandrud et al., 2017, see also Brandrud et al., 2019;Brandrud et al., 2020;Záveská et al., 2019). Consensus allelic information may also be coded in the form of ambiguous sites at heterozygous positions (e.g., Wagner et al., 2021). At last, a couple of methods are able to incorporate ploidy-aware calling algorithms and use genotype likelihoods to infer tetraploid genotypes (e.g., Ahmad et al., 2021;Brandrud et al., 2020;Karbstein et al., 2020;Záveská et al., 2019). In studies for which both kind of methods were used, the latter one generally outputted a higher number of variable sites but signals were found to be overall congruent (Brandrud et al., 2020;Záveská et al., 2019).
The recent studies of Záveská et al. (2019), Brandrud et al. (2020) and Ahmad et al. (2021) used the RAD-seq loci from the catalog generated in the popular pipeline Stacks (Catchen et al., 2011(Catchen et al., , 2013 to generate a synthetic reference onto which the raw reads are mapped back before to infer tetraploid genotypes with the approach implemented in EBG (Empirical Bayes Genotyping in Polyploids, Blischak et al., 2018, see also Blischak et al., 2016). As a putative alternative, Clark et al. (2019) proposed a package called PolyRAD to call genotypes with uncertainty from polyploid data in an R environment.
PolyRAD can deal with multiple levels of ploidy and is able to use for example population structure, linkage disequilibrium, and/or selffertilization as priors to estimate genotype probabilities. To our best knowledge, this promising and convenient method has been rarely used to implement a population genomics approach on natural populations of polyploid nonmodel organisms.
In this study, we first aimed at reevaluating the conservation status of Delphinium montanum following such population genomics approach. We used a protocol called normalized Genotyping By Sequencing (nGBS) that provided us with thousands of genetic markers (SNPs) from which we inferred patterns of diversity and differentiation with an unprecedented resolution in order to quantify inbreeding and level of connectivity between populations. We then used an Ecological Niche Modeling (ENM) approach to better understand the bioclimatic features that constraint the current geographic distribution of D. montanum. Based on several Global Climate Models (CGM), we also predicted the putative evolution of the geographic distribution of the species in the future, at different time periods (2011-2040, 2041-2070, and 2071-2100) and under different, more or less pessimistic but still realistic scenarios of evolution of global warming (or RCPs for Representative Concentration Pathways).

| Geographic distribution and population sampling
The Pyrenean Larskpur (Delphinium montanum) geographic distribution ranges from Serra del Cadí-Pedraforca (Figure 1,localities 4,5,8,and 9) up to the Puigmal massif (localities 1, 3, and 6) via Tosa d'Alp massif (locality 7). There, the distance between the most distant sites is about 60 km, and only one locality is located outside this area: Nohèdes (locality 2), in the Madres massif. This latter locality is separated from the Puigmal massif populations by the pit of Cerdagne-Conflent, which could have formed a historical barrier to gene flow. In terms of population size, Cadí-Pedraforca is the most important with up to 70% of the global population size (5450-6700 individuals). The sites of the Puigmal massif represent 25%-30% of the global population (2200-2700 individuals). Localities of Madres massif (2) and Tosa d'Alp (7) are numerically low (<300, from which none are actually flowering and 35 individuals, respectively).
We sampled a total of 106 individuals from 9 localities representative of the geographic distribution of Delphinium montanum (in the Eastern Pyrenees: 3 sites in France and 6 in Spain) between July and October 2020 (Table 1, Figure 1, see also Appendix S1). At each locality, we collected on individual distant from at least 1 meter from each other, about 1-2 cm 2 of leaf tissues that were stored in 90%

| Molecular procedures
Genomic DNA extraction and genotyping were subcontracted to LGC Genomics GmbH (Berlin, Germany). Genotyping-By-Sequencing (GBS) was performed following a specific protocol called nGBS for "normalised GBS". This double digest Restriction site-Associated DNA seq (or RAD-seq) like protocol relies on two restriction enzymes (PstI and ApeKI, in our case) to reduce D. montanum genome complexity and includes a normalization step that aims at avoiding repetitive regions. The resulting 106 individually barcoded libraries were sequenced in paired-end mode (2 × 250 bp) on an Illumina NovaSeq 6000 with an expectation of a minimum number of 1.5 million read pairs per sample.

| Genomic data processing
We used Stacks v.2.41 (Catchen et al., 2011(Catchen et al., , 2013 to build loci from Illumina reads, de novo (i.e., without aligning reads to a reference genome). Stacks consists of a wrapper of several scripts that are usually run sequentially including: process_radtags to demultiplex and clean reads, denovo_map.pl to build loci within individuals, create a catalog and match all samples against it and populations to further filter the SNPs obtained at the population level and compute basic population genetic statistics. We first optimized several key parameters of the pipeline: -m (the minimum number of identical raw reads required to form a putative allele), -M (the number of mismatches allowed between alleles to form a locus), and -n (the number of mismatches allowed between loci during construction of the catalog) on a subset of 12 individuals representative of the whole data set (i.e., geographic origin and coverage). As recommended by several authors, we varied M and n (fixing M = n) while keeping m = 3 (see Paris et al., 2017;Rochette & Catchen, 2017). The combination -m 3, -M 3, and -n 3 was found to be the most suitable to maximize the number of SNPs, assembled and polymorphic loci in our case and was used to run Stacks on the whole data set (see Appendix S2).
As D. montanum is a tetraploid species, we also specified the option "-X ustacks:--max_locus_stacks 5" (in the ustacks script's) to allow the number of alleles per site to reach the number of 5 instead or 3 by default. However, as the populations script is unable to deal with polyploid genotypes, we followed a different approach to call variants and complete genomic data processing.
We then used polyRAD (Clark et al., 2019), a method implemented in an R-package to call genotypes with uncertainty from sequencing data in polyploids. The Bayesian genotype caller implemented in polyRAD is able to import read depth from the catalog and matches files from Stacks through a specific function called poly-  Population genetic statistics were computed based on 5095 SNPs. Statistically significant values are indicated in bold (p < .05). a H E : unbiased expected heterozygosity measured following Nei (1987) in Genodive.
binomial distribution) and adjusted its parameter to 6 (from a range of 2 to 14). This information was used to call genotypes while taking into account population genetic structure. The most likely genotypes of the remaining 5 095 loci were then exported in several formats such as genind objects or Structure (.str).

| Population genomics analyses
To get an overview of the overall genetic diversity and differentiation among individuals and populations, we first performed a Principal Component Analysis (PCA) based on a matrix of 106 individuals (as rows) and 5095 SNPs (as columns) coded as a genind object with the R-package "adegenet" (Jombart, 2008). We then used Genodive v.3.04 (Meirmans, 2020) to compute expected and  Rousset, 1997) and its statistical significance was tested with a Mantel test between matrices of pairwise G ST and geographical distance.
To further investigate population structure and characterize putative migration/admixture event, we used sNMF (Frichot et al., 2014) as implemented in the R-package "LEA" (Frichot & François, 2015) to estimate individual ancestry coefficients based on sparse non-negative matrix factorization algorithms. sNMF is particularly suitable for our genotype dataset as it can deal with polyploid data and is fast enough to be applied on hundreds of individuals and thousands of markers with reasonable computation time. The number of genetic clusters was varied from K = 1 to 10, and analyses were run with 10 replicates at each value of K.

| Ecological niche modeling
In order to assess the future trends of spatial dynamics of D.  (Fick & Hijmans, 2017) in a way that they consist of downscaled Global Climate Models (GCM) output temperature and precipitation estimates at a maximal resolution of 30 arc seconds (approximately 1 km). In this study, we used the most recent version of CHELSA as the methodology and bias correction it implements have been shown to outperforms WORLDCLIM in mountainous environments (Bobrowski et al., 2021;Karger et al., 2017). Current climate data correspond to time aver- in considering a total of 9 out the 106 occurrences initially available.
We then examined raw values for the 19 bioclimatic variables extracted at each occurrence point and noticed that the variable bio19 ("precipitation of coldest quarter") displayed obviously aberrant (i.e., negative) values. This variable bio19 was thus discarded from subsequent analyses. From the 18 remaining variables (bio1 to bio18), we selected the less correlated ones (Person correlation coefficient <0.75) thanks to the R-package "caret" (Kuhn, 2019) and kept five variables: bio5 ("mean daily maximum air temperature of the warmest month"), bio7 ("annual range of air temperature"), bio14 ("precipitation amount of the driest month") and bio15 ("precipitation seasonality") for further analyses.
The calibration area for the models was created as a buffer of 1.5° around the minimum convex polygon encompassing all occurrences. We used the maximum entropy method (implemented in MaxEnt ver. 3.4.1, Phillips et al., 2006Phillips et al., , 2017 to calibrate model and evaluated model performance thanks to the package "ENMeval" (Muscarella et al., 2014) as implemented in ENMwizard. We evaluated models using a geographic partition scheme of type "block" and to another (this holds true for the number of total alleles A, too, see Table 1). The number of private alleles Ap is consistent with the relative level of geographic isolation of the populations.
The degree of genetic differentiation was found to be relatively high overall (

| Suitable habitats for Delphinium montanum and prediction of its spatial dynamics
ENMeval analyses identified RM = 2 and a combination of Quadratic and Hinge feature (QH) classes as the best-performing parameters F I G U R E 2 (a) Principal Component Analysis (PCA) displaying the two first axes (PC1 and PC2) representing 14.85% and 9.21% of the total genetic variance. PCA was computed based on the 106 individuals of D. montanum genotyped at 5 095 SNPs whose colors represent sampling localities. (b) Values of the cross-entropy criterion for a number of clusters ranging from K = 1 to 10 (10 sNMF runs each). The optimal number of K was found to be 6. (c) Barplot of ancestry coefficients obtained from sNMF for 106 individuals for K = 2 and K = 6, based on 5095 SNPs for calibrating the final ENM following the "LowAIC" and M = 2 and a combination of Product and Hinge feature (PH) classes following the "AUC" optimality criteria. The first model had the best corrected Akaike Information Criterion (AICc) score (152.25) and an AICc weight of 0.33, whereas the second one had an AICc score of 174.63 and an AICc weight of 4.53 × 10 −6 . These two models had mean omission rates of 0.375 and 0.25 for the 10th percentile and the lowest presence training thresholds, respectively, as well as a mean test AUC of 0.988. Values of these diagnostic metrics are provided for all candidate models in Appendix S4. Both models (RM = 2/PH and RM = 2/QH) considered bio15 ("precipitation seasonality") as being the most important explanatory variable with a contribution of 88.6% and 90.0%, respectively, followed by bio7 ("Temperature annual range (bio5-bio6))" with 6.2% and 7.1% and bio5 "Maximum temperature of the warmest month" with 5.2% and 3.0%; the contribution of bio14 "Precipitation of the driest month," being 0% with both models. Current suitable area for D. montanum was found to be 2670 and 2365 km 2 following "LowAIC" and "AUC" selection model criteria, respectively (the values were identical whatever the threshold method used). All future consensus spatial projections predict a dramatic decrease of the suitable area (Table 2 and Figure 4). For RCPs 7.0 and 8.5, the suitable area would have completely disappeared for the time period 2071-2100 (but not for RCP 2.6). We also noticed some slightly different outcomes depending on the selection model criteria with RCP 8.5 scenarios giving sometimes higher values than RCP 7.0, and even RCP 2.6 (e.g., time period 2011-2040, "Low AIC" criterion).

| Conservation genomics on a polyploid species with large genome size
The

| A set of populations isolated on sky islands
Our results confirm that populations of D. montanum can be seen  (2011-2040, 2041-2070, and 2071-2100) and various Representative Concentration Pathways (RCP 2.6, 7.0, and 8.5). by a matrix of unsuitable habitats and that gene flow has little to no chance to be efficiently maintained within this metapopulation. All populations considered in this study were found to display deficit in heterozygotes (G IS > 0.315, p < .05). These results confirm previous findings of Simon et al. (2001) and López-Pujol et al. (2007) and show that the autotetraploid status of the species does not impede relatively low levels of heterozygosity (H O < 0.17).

F I G U R E 4 Geographic areas holding suitable climatic conditions for
Polyploidy is expected to be associated with increased genetic variation, a characteristic that provides more diversity to adapt (especially to extreme environments such as mountainous ones) and allows buffering effect against deleterious mutations (see van de Peer et al., 2021 and references therein). From a conservation point of view, polyploidy does not seem to be associated with increased risk of extinction even though plant with big genome size may be more prone to be under threat (Pandit et al., 2011;Vinogradov, 2003). As two other closely related congeneric spe-

| Climate change as the main threat for Delphinium montanum
The two previous studies of Simon et al. (2001) and López-Pujol et al. (2007) have listed several threats to explain the decline of D. montanum, at least in some particular populations. According to Simon et al. (2001), natural disturbance of habitats due to rock falls and avalanches, trampling because of an overall increase in human visits and animal predation (e.g., caterpillar, Pyrenean chamois/isards and perhaps red and roe deers) as well as competition for insect pollinators with simultaneously blooming species (e.g., Aconitum spp.) have thus been considered as more or less serious problems for D. montanum.
The most recent study of López-Pujol et al. (2007)  further supporting that this phenomenon could be nonspecific to D. montanum and perhaps caused by bioclimatic or at least environmental modifications. There, output from Regional Climate Models (RCMs) predict an increase in mean annual temperature that could reach 3.1 to 4.5°C (+1.5°C already observed since the mean of the period 1961-1990) associated with a decrease in precipitation of up to 25% and perhaps more compared to present-day conditions (Lespinas et al., 2014). According to the same study, precipitations would mostly decrease in summer (−40%) and in spring (−20%) but would remain similar the rest of the year. This in accordance with our ENM predictions suggesting that bio15 "Precipitation seasonality," the annual coefficient of variation of rainfall may be by far, the most important variable in explaining the bioclimatic niche of D. montanum.

| RECOMMENDATIONS FOR CONS ERVATION
Ex-situ conservation would be a way to preserve the genetic diversity of D. montanum, but also a bridge for future in-situ conservation We can identify three main methods of in-situ conservation that could be envisaged for D. montanum (see Godefroid et al., 2011;Maschinski & Albretcht, 2017;Ren et al., 2014): reinforcement, translocation and managed relocation (or assisted migration). The first one is reinforcement, which consists in adding ex-situ grown individuals in their natural population of origin to increase the total number of individuals, increase the spatial occupancy of the species and thereby reduce extinction risk. On the one hand, the increase in population size provided by reinforcement could contribute to coun- For the very peculiar situation of 2-Nohèdes's locality, which seems to be the population the most at risk of extinction at the moment, there is still the question of the lack of flowering. To test whether this problem may be associated with the absence of a specific environmental stimulus, we may propose to either (i) experimentally modify temperature/water supply, in situ, (ii) translocate some individuals to localities where D. montanum seem healthy (e.g., Cadí range) and see if flowering can be restored and/or (iii) experimentally modify temperature/water supply in the lab to better understand the conditions required for flowering and figure out the limiting factor in the wild. Alternatively, it may be due to the fact that that this locality is formed by immature individuals and that the local population itself may be recovering from an important mortality event. At the moment, we have no further arguments to support this scenario.

| CON CLUS ION
Overall, our results confirm that the remaining populations of D. montanum show evidence of decreased genetic diversity and evidence of inbreeding at a genomic scale. The degree of genetic structure observable is consistent with a set of strongly isolated populations no longer able to maintain gene flow between their sky islands. In this context, in-situ conservation methods are likely to have an effect on the status of the species. Reinforcement would help counteracting decrease in local population sizes, and translocation would help increasing genetic diversity, buffering inbreeding depression and perhaps, bring the variation required for D. montanum to adapt. This being said, ex-situ crossings are first recommended to verify the absence of outbreeding depression as local genetic pools were found to differ substantially. As our results support that environmental change may represent the main risk of extinction for D. montanum all these measures may be however considered as of limited interest on the long term, if the ecological niche of the species disappears. They also thank Olivier François, Lindsey Clark, Patrick Meirmans, and Neander Marcel Heming for technical advice and insightful discussions.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.