New conservation viewpoints when plants are viewed at one level higher. Integration of phylogeographic structure, niche modeling and genetic diversity in conservation planning of W Mediterranean larkspurs

Protection and management of closely related endangered species and subspecies at a very narrow regional scale is the origin of multiple dysfunctional conservation decisions. These include artificially increased IUCN risk assessment categories and derived consequences: poor effectiveness in allocating public and private funds or repeat of unnecessary actions/ facilities. Data provided by the revisited study of a group of W Mediterranean larkspurs (Delphinium ser. Fissa), including new data on demography, niche modeling, genetic diversity and phylogeography, contributed to a new and wider analysis of causes of threat. Although current IUCN Red List regulations did not allow for assessments at levels higher than a specific rank, scientific information suggests that in some cases this could be a better approach for sound scientifically-based biodiversity conservation planning and action. © 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Biodiversity Conservation, as a discipline of crisis, by definition (Soul e, 1985), usually works with uncertainty and limited knowledge to solve increasingly complex conservation problems. Unfortunately, the availability of better scientific knowledge about endangered species is scarce, with gaps existing in critical information, such as incomplete sampling, only partly known species distribution ranges and very few data on biological attributes (demography, genetic variability, breeding system and reproductive features, and mutualisms, among others). In fact, there is still much biodiversity to be discovered (Forest et al., 2015;Deli c et al., 2017;Treurnicht et al., 2017).
Overall, in a well-balanced framework of biodiversity management, conservation needs can push for evolutionary, phylogenetic and taxonomic updates, at least in speciose groups (Forest et al., 2015). Numerous daily conservation decisions concerning regional administrations (under the classical pressure pair: local populations' problems/locally needed solutions) show that effective, long-term and sustainable conservation actions could only be adopted following sound scientific knowledge. This cannot usually be acquired as rapidly as needed to counteract stochastic decreases of population viability and, as a consequence, not all conservation efforts are correctly designed to obtain successful results (L opez-Pujol et al., 2013).
Nevertheless, a well combined and intercommunicative network of biodiversity researchers and managers may provide very useful win-win advances on both sides (e.g. Center for Biodiversity and Conservation Science, University of Queensland, Australia, at https://cbcs.centre.uq.edu.au/). Additionally, some controversial conservation policy decisions that are not based on sound scientific data could be avoided, and researchers could focus significant efforts on satisfying the technical needs of managers without giving up higher knowledge challenges. Thus, conservation-oriented studies can provide the perfect approach to specific or general biodiversity scientific problems.
Delphinium L. ser. Fissa Pawl. (Delphinieae, Ranunculaceae) is a group of hemicryptophyte larkspurs of narrowly dissected leaves represented in the W Mediterranean area by four closely related taxa, including two species; D. bolosii C. Blanch e & Molero and D. mansanetianum Pitarch, Peris & Sanchis, and two subspecies of the widely distributed species D. fissum Waldst. & Kit., D. fissum subsp. sordidum (Cuatrec.) Amich, E. Rico & J. S anchez and D. fissum subsp. fissum. All entities, except the last, are endemic to the Iberian Peninsula (Blanch e and Molero, 1983) ( Fig. 1A and A1). Common systematic practice in this group still follows a traditional combination of morphological and cytogenetic features (Blanch e, 1991;Warnock, 1995;Ramírez-Rodríguez et al., 2016). Fine molecular taxonomic delimitation at species or subspecies level in perennial larkspurs is still not available, although genus/tribe levels have been investigated by Jabbour and Renner (2012a). The four taxa have small geographic ranges, some of them being recently described so their real distribution may be poorly known, and all have been generally perceived as an endangered biodiversity component of the W Mediterranean flora (protected by law at several levels; Table A1).
However, some disagreements on i) inter-and intraspecific taxonomic delimitation (see synonyms list in Table A2); ii) levels of threat (see distinct results of IUCN risk assessment depending on regional approach or scientific team in Table A1) or iii) the effective conservation actions undertaken (see differential law consideration depending on taxa or countries in Table  A1), all suggest that a ground level, local view of this group may not help to fully understand natural relationships or present vs. past gene exchanges among populations, levels and structure of genetic diversity (and taxonomy in accordance with biology), nor to precisely identify levels of threat or define conservation strategies based on the best available scientific knowledge.
Delphinium ser. Fissa may be considered, therefore, as a suitable case study to explore if, from a slightly higher level, a better resulting overview is obtained (i.e. by considering a multidisciplinary and higher taxa approach to evolution and conservation issues). In our case, there is agreement that the series taxonomic rank is the most natural one in the genus Delphinium. The morphological treatments of the classical monographs (Pawlowsky, 1963;Munz, 1967;Malyutin, 1987;Warnock, 1995) are also generally supported by the clade resulting structure after the molecular phylogeny study (Jabbour and Renner, 2012a), where series resolution is the best compartmental interpretation of diversity patterns within the genus (each taxonomic group will better fit the systematic assemblage as an appropriate unit for the challenge of an upperstep approach).
This complex of four W Mediterranean Delphinium taxa has received increasing attention during the last 20 years. New scientific data are available (both from field exploration and lab research as well as from gray literature derived from conservation actions) that offer a significant increase of information about known populations as well as new genetic data and analyses. However, some gaps in sound biological knowledge still exist. Thus, we can derive the following goals: 1) to update evolutionary knowledge and phylogeography; 2) to infer potential distribution areas of taxa and to test for niche divergence/ conservatism; and 3) to propose new conservation guidelines under a global view of the whole series Fissa to check if a slightly higher approach can help in providing useful tips to better tackle conservation challenges.

Study species and sampling strategy
The four Delphinium taxa studied herein belong to a monophyletic clade (Jabbour and Renner, 2012a). The scope of ser. Fissa is described in Blanch e and Molero (1983) and the full nomenclatural details and synonyms for the studied taxa are shown in Table A2. Fifteen populations were surveyed in this study (Table 1, Fig. A1). Representative voucher specimens are kept in the herbarium BCN (Universitat de Barcelona). For D. bolosii, all known populations were sampled, including the recently discovered ones (ALO, CAM and MUR). Three populations of D. fissum subsp. fissum, from the W part of its range (the Maritime Alps) were studied, whereas a representative selection of the whole range of D. fissum subsp. sordidum (five populations) was sampled. Finally, the single known population of D. mansanetianum was also studied. All 15 populations were scored for both allozyme and cpDNA analyses (sample sizes are indicated in Table 1). For allozymes, fresh leaf samples were collected in the field and stored at 4 C until protein extraction. Leaf fragments for DNA survey were collected and placed in plastic bags filled with silica gel and kept desiccated at À80 C in an ultra-freezer until DNA extraction. Chromosome counts for many populations were mainly compiled from Bosch et al. (2016). In addition, for cytogenetic research, seed material was collected from ALO and CAM populations of D. bolosii and flower buds from the D. mansanetianum population (MOS). All field collections obtained the necessary permissions from conservation authorities and were accurately performed using nondestructive methods.
Loci were numbered consecutively and alleles at each locus were labelled alphabetically, beginning with the most anodal form in both cases. Allozyme frequencies at each locus were calculated for each population. To estimate the levels of genetic diversity, the following statistics were computed: P 95 , the percentage of polymorphic loci when the most common allele had a frequency of <0.95; A, the mean number of alleles per locus; PA, the number of private alleles; AR, the mean allelic richness using a 'rarefaction method' that compensates for uneven population sample sizes (El Mousadik and Petit, 1996); H o , the observed heterozygosity; and H e , the unbiased expected heterozygosity. All parameters were calculated using GenAlEx v. 6.5 (Peakall and Smouse, 2006) with the exception of AR, which was computed with FSTAT v. 2.9.3 (Goudet, 1995). To test whether populations were under Hardy-Weinberg equilibrium (HWE), Wright's (Wright, 1965) F IS was estimated following the method of Weir and Cockerham (1984) with GENETIX v. 4.05 (Belkhir et al., 1996). Statistical significance of F IS values for each locus per population was tested by permutation tests (10,000 randomizations), also using GENETIX.
Genetic structure was assessed through five different methods. First, the Bayesian algorithm implemented in Structure v. 2.3.4 (Pritchard et al., 2000) was used under an admixture ancestry model with correlated allele frequencies. The burn-in period and Markov chain Monte Carlo (MCMC) were set to 50,000 and 500,000 iterations, respectively, and K was run from 1 to 15 (with 10 replicates). The most likely value of K was determined by the DK statistic of Evanno et al. (2005), with the aid of STRUCTURE HARVESTER v. 0.6.94 (Earl and vonHoldt, 2012). For the most likely K, CLUMPP v. 1.1.2 (Jakobsson and Rosenberg, 2007) was used to combine the results of the 10 replicates of the best K. To plot the output result produced by CLUMPP, we used the program DISTRUCT v. 1.1 (Rosenberg, 2004). Second, an analysis of molecular variance (AMOVA) was run with the aid of GenAlEx v. 6.1, in order to estimate the influence of individuals within populations, populations within species, and species on the observed genetic variation. Third, pairwise genetic distance between populations was calculated using Nei et al. (1983) genetic distance (D A ); distance matrices were converted into UPGMA (unweighted pair-group method using arithmetic averages) and NJ (neighbor-joining) trees (after 1000 bootstrap replicates), employing the program Populations v. 1.2.30 (Langella, 1999), and visualized with TreeView v. 1.6 (Page, 1996). Fourth, a principal coordinate analysis (PCoA) was conducted using GenAlEx based on genotypic distances. Lastly, Barrier v. 2.2 software (Manni et al., 2004), based on Monmonier's (Monmonier, 1973) algorithm, was run in order to identify genetic barriers between populations, with significance tested by means of 1000 bootstrap matrices of D A . Geographical coordinates of populations were used to obtain a Voronoï tessellation where barriers were delineated. To gain insight into the patterns of recent (i.e. within the last few generations) gene flow between populations, migration rates were estimated using the program BayesAss v. 1.3 (Wilson and Map in Fig. A1; N ¼ number of studied individuals/population. Rannala, 2003). We ran 3 Â 10 6 MCMC iterations, with a burn-in of 999,999 iterations and a sampling frequency of 2000 by setting delta at 0.15 (the default value).

cpDNA
Two intergenic spacer regions from plastidial DNA (atpB-rbcL and tnrV-ndhC) were amplified using the primers described in Heinze (2007). Polymerase chain reaction (PCR) amplifications were made in a total of 50 mL including 1X PCR buffer, 4 mM MgCl 2 , 0.04 mM dNTP, 0.16 mM of each primer, 0.5 U of Taq Polymerase and 50e100 ng of DNA. The temperature profile for both regions started with an initial denaturation step of 3 min at 92 C. Then 30 cycles of 30 s at 94 C of DNA denaturation, 30 s at 55 C for primer annealing and 30 s at 72 C of elongation (for the atpB-rbcL spacer), or 30 s at 94 C, 30 s at 50 C and 45 s at 72 C (for the tnrV-ndhC spacer), with a final step of 3 min at 72 C. Amplified products were sent to Macrogen Inc. (Korea) for sequencing. A statistical parsimony haplotype network was constructed using TCS v. 1.21 (Clement et al., 2000). We also evaluated the genetic structure using the Bayesian clustering method implemented in BAPS v. 6.0. (Corander et al., 2008). Ten replicates from K ¼ 2 to K ¼ 15 were run and the most likely K value was chosen according to the highest log marginal likelihood values. The evaluation of genetic variability through different species was determined using AMOVA and was carried out using Arlequin v. 3.5.1.3 software (Excoffier and Lischer, 2010). SAMOVA was constructed using SAMOVA v. 1.0 software (Dupanloup et al., 2002) running from K ¼ 2 to K ¼ 15. As for allozymes, we used the Monmonier's maximum difference algorithm implemented in Barrier software to identify the five strongest possible barriers between populations. The algorithm was applied on D S Nei genetic distances (Nei, 1972) obtained from GenAlEx.

Chromosome counts
Germinated seeds in Petri dishes were treated for 2 h at 4 C with a colchicine solution (0.05%) and washed with distilled water for 5 min. The fixation was carried out with the Farmer reagent (absolute ethanol and glacial acetic acid [3:1 v/v]), softening with acid hydrolysis with 1N hydrochloric acid (at 60 C for 10 min), and staining by immersion in 2% acetic orcein at room temperature for 24 h. Fixed flower buds were directly stained with aceticecarmine 1%. Observations were made with a Zeiss Axioplan microscope, and karyotypes were analyzed from a minimum of five metaphase plates of good quality. See detail of procedure and calculated parameters in Bosch (1999).

Ecological niche modeling
In order to evaluate the present potential distribution of Delphinium ser. Fissa, we performed Ecological Niche Modeling (ENM) using the maximum entropy algorithm implemented in MaxEnt v. 3.3 (Phillips et al., 2006), which requires presenceonly data. This algorithm has been proved to have a high robust performance and accuracy, in comparison to other ENM methodologies (Pearson et al., 2007;Wilting et al., 2010). Delphinium mansanetianum was deliberately excluded from the ENM because it occurs in a single location. The georeferenced occurrence records used are illustrated in Fig. 1A and detailed lists are kept at the University of Barcelona (due to confidentiality commitments with conservation authorities), but will be available from the authors upon request and agreement. To avoid spatial autocorrelation bias, presence points were filtered keeping only one within each pixel of 30 arc-sec (ca. 1 km). The filtering was done automatically by MaxEnt software in the case of ENM, and for niche comparison analysis by the R script used (see details below). As model variables, a set of 19 bioclimatic layers was downloaded from the WorldClim website (http://www.worldclim.org/, Hijmans et al., 2005) under current conditions (ca. 1960À1990). The information for the 19 variables was extracted from 1000 points randomly generated in the study area ( Fig. 1A) with the aid of ArcGIS v. 10.2 (ESRI, 2013). After a Pearson pairwise correlation analysis, we retained a set of nine relatively uncorrelated (r ! j0.85j) variables to perform the models: mean diurnal range (bio2), isothermality (bio3), temperature seasonality (bio4), minimum temperature of the coldest month (bio6), mean temperature of the wettest quarter (bio8), mean temperature of the driest quarter (bio9), precipitation seasonality (bio15), precipitation of the driest quarter (bio17), and precipitation of the coldest quarter (bio19). Variables were selected based on their relative contribution to the models (percentage contribution, permutation importance and jackknife of regularized gaining train), making a consensus selection to include the most influential variables for the three taxa. Given the low number of occurrences for D. fissum subsp. sordidum (13 after filtering) and D. bolosii (seven after filtering), we tested whether these species were suitable for modeling following the jackknife method developed by Pearson et al. (2007). Both taxa proved adequate for further analysis. To get the definitive ENM models, MaxEnt was run 50 times using the 0 bootstrap 0 method for D. fissum subsp. sordidum and D. bolosii, and 0 subsample 0 for D. fissum subsp. fissum (with 25% of the occurrence data set chosen randomly to test the model). The parameters selected to run the models were set as recommended (Phillips et al., 2006): auto features, regularization multiplier ¼ 1, maximum number of iterations ¼ 500, convergence threshold ¼ 0.00001 and maximum number of background points ¼ 10,000. The AUC values were used as a metric to assess models' performance. The ENM predictions were visualized in ArcGIS and exported as a binary presence/absence distribution according to the maximum sensitivity plus specificity logistic threshold, which is very robust with all types of data (Liu et al., 2016).

Niche comparisons
To aid in species delimitation and to test whether the closely related, allopatrically-distributed Delphinium taxa are occupying different climatic niches (evaluating the same nine bioclimatic variables as for the ENM), we used the approach based on environmental space through a PCA-env developed by Broennimann et al. (2012). As proposed by Silva et al. (2016), to construct the PCA-env the background climatic areas for each taxon were extracted from a minimum convex polygon of 0.3 degrees of buffer size. The observed occurrences were smoothed with a kernel density function and then projected to the environmental space of 500 Â 500 grid-cell resolution, where each cell represents unique climatic conditions. To explore whether the niches of the different Delphinium taxa studied have been conserved or have diverged, we used the data extracted from the PCA-env to calculate the niche overlap metric of Schoener's D (Schoener, 1970;Warren et al., 2008) (from 0 representing no overlap to 1 complete overlap), and to perform two statistical tests (Warren et al., 2008;Broennimann et al., 2012); the niche equivalency (only considering the occurrences) and niche similarity (considering both the occurrences and background data). The PCA-env analysis was conducted using the original R code of Broennimann et al. (2012) with posterior  (Phillips et al., 2006). modifications of Silva et al. (2016), with slight self-adjustments. Finally, the values of climatic variables from presence records were extracted using ArcGIS to compare each variable separately.

Allozymes
Fifteen consistent loci were scored: Aat-1, Aco-1, Aco-2, Acp-3, Adh, Dia-2, Idh, Mdh-1, Mdh-2, Me, 6Pgd-1, 6Pgd-3, Pgi-1, Pgm-1 and Pgm-2. Of these, nine were polymorphic across populations according to the 95% allele-frequency criterion, and 11, considering those harboring two or more alleles regardless of their frequencies (Table 2). Considering all populations as a whole, a total of 33 alleles were observed (Table 2). Of these, 28 were detected in D. bolosii, while 25, 20 and 18 alleles were found in D. fissum subsp. sordidum, D. fissum subsp. fissum and D. mansanetianum populations, respectively. The population with the largest number of alleles was RBO of D. bolosii, with 22, and the poorest populations were ESC and MAG of D. fissum subsp. sordidum, both with 16 alleles. Although MUR (D. bolosii) had one allele less than the conspecific RBO, the former was the most singular population as it harbored three private alleles, all diagnostic for D. bolosii. In fact, D. bolosii showed six diagnostic alleles, followed by four in D. fissum subsp. fissum (the other two taxa had no exclusive alleles; Tables 2 and 3). Overall, low levels of genetic diversity were found for all studied taxa (Table 3) with D. bolosii and D. fissum subsp. sordidum showing the highest and the lowest levels of variability, respectively. Indeed, the most genetically rich populations were those of D. bolosii, regarding polymorphic loci (P 95 ¼ 33.33 in ULL2), allelic richness (A ¼ 1.467 in RBO; AR ¼ 1.400 in MUR), and expected heterozygosity (H e ¼ 0.117 in ULL2) (Table 3). Mean fixation indices (F) for almost all populations (with the single exception of BOU) were not significantly different from zero, showing a general concordance with the Hardy-Weinberg expectations; in accordance with this, only 10 out of 46 F-values at individual loci were significantly different from zero (data not shown).
According to Evanno's approach, K ¼ 3 and K ¼ 5 were the most likely numbers of genetic clusters for Delphinium ser. Fissa in the W Mediterranean (Fig. A2). In both cases, the only taxon with its own cluster was D. fissum subsp. fissum, whereas populations of D. bolosii showed the same genetic make-up: two pairs of populations belonging to two different clusters, with Table 2 Allele frequency for 11 polymorphic loci of populations of Delphinium ser. Fissa in W Mediterranean. the remaining two populations (ALO, CAM) showing a high degree of admixture (Fig. 2). Populations of D. fissum subsp. sordidum plus D. mansanetianum were clustered following a geographical pattern: W Iberian vs. E Iberian populations (Fig. 2). AMOVA results confirmed the low correspondence between species and genetic clusters: the among-taxa component (i.e. when the populations were grouped into species) was only one third (Table 4). The segregation of both D. bolosii and D. fissum subsp. sordidum in several clearly-differentiated genetic subgroups, as well as the distinctiveness of D. fissum subsp. fissum, were also supported by both the UPGMA/NJ trees (Fig. 3) and the PCoA analysis (Fig. 4). The first genetic barrier detected, using Monmonier's algorithm, mostly separated D. fissum subsp. fissum from D. bolosii (with 86% bootstrap support; Fig. 5), with subsequent barriers separating D. bolosii populations as well as D. bolosii itself from D. fissum subsp. sordidum. Also, although weak, the separation between E and W Iberian populations was evidenced (Fig. A3). The BayesAss analysis indicated a near absence of recent gene flow between populations, except for a few population pairs (Table 5).

cpDNA data
The aligned sequences of the two cpDNA regions represented 1622 bp (atpB-rbcL, 762 bp; trnV-ndhC 860 bp) and resulted in three haplotypes ( Fig. 6A and B). Haplotype 2 was shared between all four studied taxa, haplotype 1 was exclusive to D. fissum subsp. sordidum and haplotype 3 was restricted to a single population of D. fissum subsp. fissum (BOU; Fig. 6A). No Table 3 Main parameters of genetic diversity at the population level for the four studied taxa of Delphinium ser. Fissa, based on 15 allozyme loci.  The run with highest lnPr (X j K) out of a total of 10 is shown.
intrapopulation variation was observed. The parsimony network (Fig. 6C) showed that the maximum number of mutational steps to occur between different haplotypes can be inferred to be six. The Bayesian analysis of cpDNA population structure revealed two genetically distinct groups (Fig. 7). One of these groups was characterized by populations carrying haplotype 1 and the second group included haplotypes 2 and 3. AMOVA analyses   attributed 62.29% of the genetic variation to the among-taxa component (Table 4). From the spatial genetic analysis with SAMOVA, the maximal value of F ct was determined as 1 (K ¼ 3). Subdivision groups were: 1) ESC, HER, MAG, VIL; 2) ALO, CAM, RBO, ULL1/ULL2, MUR, SAD, COR, SAV, MOS; and 3) BOU. Regarding the Barrier analysis (Fig. 7), the first genetic barrier we can infer was between three populations of the W range of D. fissum subsp. sordidum (ESC, HER, MAG) and the COR population of the E range of this taxon. The second barrier revealed was located between populations of D. fissum subsp. fissum, separating SAV from BOU. The third barrier separated the westernmost population of D. fissum subsp. fissum (BOU) from MUR (the easternmost population of D. bolosii). Finally, fourth and fifth barriers isolated the ALO population from contiguous ones of D. bolosii (data not shown). The Bayesian analysis of cpDNA population structure revealed two genetically distinct groups (Fig. 7).

Chromosome counts
Two additional counts, from ALO and CAM populations of D. bolosii not previously studied, gave 2n ¼ 16 chromosomes (Fig. 8), showing normal balanced eudiploid cytotypes. Karyotypes presented typically bimodal patterns with two submetacentric and six subtelocentric chromosome pairs (Table A4). The chromosome number of 2n ¼ 16 is reported here for the first time in D. mansanetianum (MOS).

Ecological niche modeling
Niche models generated by MaxEnt analyses had AUC >0.93 for the three Delphinium taxa (Table A3), indicating that the models performed well at predicting current species distribution. Jackknife tests revealed that the highest informative variables for the models differed greatly between taxa (Table A3). The predicted areas of suitable habitat for the three taxa (Fig. 1B) were broadly congruent with their respective occurrence points. Nevertheless, the predicted suitable habitats of D. fissum subsp. sordidum and D. bolosii covered a somewhat greater area than their actual ranges known to date; for example, suitable habitats of D. bolosii appeared in the central-east part of the Iberian Peninsula, over 200 km away from its nearest known occurrences (Fig. 1).

Niche comparisons
In the PCA-env analysis, the two principal components (PC) together explained 75.69% of the total variation (59.24% from PC1 and 16.45% from PC2; Fig. 9). The variables most closely correlated with PC1 and also with PC2, were precipitation variables: precipitation of the driest quarter (bio17), and precipitation seasonality (bio15), respectively (data not shown). The individual PCA-env plots for each taxon (Fig. 9B-D ) revealed that D. fissum subsp. fissum presents the highest climatic niche breadth in relation to the occurrence density, and also the widest climatic background areas (Fig. 9B). All three studied taxa showed narrow niche occupancies in comparison with their background climates ( Fig. 9C and D). A clear niche differentiation was shown by the multiple niche PCA-env plot displaying 50% of occurrence density (Fig. 9E), where a non-overlap pattern was observed. However, when 100% of occurrence density was plotted (Fig. 9F), marginal niches were partly overlapped, especially between D. bolosii and D. fissum subsp. fissum.
The observed overlap values of D index were low in all niche comparisons (D ¼ 0.000e0.009; Table 6). According to the niche equivalency and similarity tests, there was no significant niche conservatism pattern for any of the taxa comparisons. Conversely, niche divergence between D. fissum subsp. fissum and D. fissum subsp. sordidum was observed in the equivalency test (Table 6). For D. bolosii, there was no clear pattern of niche divergence or conservatism, with respect to the other taxa, probably due to the methodological limitations of using fewer than 10 occurrences.

Genetic diversity and gene flow
Levels of genetic diversity are relatively low (from moderate to extremely low) when compared with endemic and narrowly-distributed plants (Hamrick and Godt, 1990), but also with endemic and endangered species in the Mediterranean Basin ). However, differential diversities can be found, both within species [H e of largest populations (ULL 1/2) of D. bolosii is ca. 0.110e0.120, i.e. almost twice that of the smallest ones (CAM, ALO), with ca. 0.060e0.070; Table 3], and between species (the mean H e value for D. bolosii is 7e8 times larger than that of D. fissum subsp. sordidum and D. mansanetianum; Table 3). In terms of allelic richness, these findings can be attributed partly to population size effects (N e and percentage of vegetative propagation) and partly to allele loss due to extreme/peripheral area effects (by accumulation of bottlenecks and isolation from the core genetic pool of the species) (Frankham et al., 2010) on a hypothetical westward/ southward migration route (further details can be found in the section "Differentiation and speciation"). At this point it should be noted that the less genetically variable species (D. fissum subsp. sordidum) is situated in the westernmost extreme of the series Fissa range and, among its populations, the most isolated and peripheral one is MAG, which is also the most depauperated population both in terms of allele richness and heterozygosity ( Fig. 1 and Table 3). Notably, the only species with two cytotypes (D. bolosii) is the most genetically diverse, which, again, highlights the importance of this species in the diversi-  Table 5), indicating that, at least within the last few generations, some gene flow between populations has taken place. The populations most involved in recent interchange of genes are those of D. bolosii (especially RBO, which acts as a source of migrants). These findings do not agree with the generic assessments made by Ramírez-Rodríguez and Amich (2017), assuming that D. fissum populations >10 km apart are genetically isolated. On the contrary, recent discoveries of new plant stands in remote locations not previously known (see below) suggest some recent mobility of groups of individuals, including intermediate spots connecting previously known populations (e.g. spots along river gorges in Mura creek, Montsant and Segre rivers for D. bolosii ;Bosch, 2009;Turmo, 2017). Genetic similarities, even between populations from different taxa (Figs. 2 and  3), in addition, suggest ancestral connections among populations. STRUCTURE (Fig. 2) and Barrier (Fig. 7) analyses help in identifying population clustering, drawing a consistent pattern of diversity and connectivity. The first barrier detected with allozymes might be interpreted as a sharp genetic break (bootstrap support ¼ 86%) between the E cluster corresponding to D. fissum subsp. fissum which, in turn, and according to cpDNA, is genetically subdivided by the Rhone River valley, a well-known biogeographic barrier (M edail and Diadema, 2009), and the whole W Iberian cluster, where a more diverse group of populations (¼ D. bolosii) acts as a diversity "reservoir" (higher genetic diversity, more private alleles, higher number of total alleles). Within this suggested "reservoir", the population MUR seems to be the main link in the chain (L opez-Pujol et al., 2014), given its high genetic variability but also singularity (in the UPGMA, it appears as sister to the remaining populations of the Iberian cluster; Fig. 3). Leaving aside D. bolosii, the rest of the Iberian populations are structured in two genetic subclusters (W and E populations), following a geographical pattern, with only a minor difference between allozymes and cpDNA (MAG population is included within either the W or the E group).  Delphinium mansanetianum is genetically most closely related to the COR population of D. fissum subsp. sordidum, its geographically closest population (they are separated by ca. 140 km, and both are included within the same mountain range, the Iberian System), although at present the species is not exchanging genes, either with this population nor the others (Table  5). Ancestral gene flow from/to D. mansanetianum can, however, be suggested, on the basis of its genetic affinities, including with the S populations of D. bolosii (located not far away d ca. 160 km).
If contemporary gene flow among the populations of our study system is actually taking place, two conservation implications have to be noted: 1) levels of threat in connected groups must be lowered according to IUCN (2017) guidelines, as demographic or genetic rescue is possible and subpopulation isolation is not complete; and 2) reintroduction, reinforcements and translocations must carefully take into account compartmentalization of gene diversity in order not to break isolation when it exists, to avoid outbreeding depression and the loss of population genetic identity (Frankham et al., 2010).

Chromosome number
The previously known chromosome number for D. bolosii was 2n ¼ 18 (Blanch e, 1991;Bosch, 1999), but recently, a population with 2n ¼ 16 (MUR) (Blanch e et al., 2014;Bosch et al., 2019) was identified. The two new 2n ¼ 16 eudiploid counts reported here (ALO, CAM) add evidence that such species previously considered exclusively dysploid (2n þ 2) (Blanch e et al.,  The niche overlap metric (Schoener's D) represents the overlap level between climatic niches compared (D ¼ 0 no overlap; D ¼ 1 complete overlap). The niche equivalency (eq) and niche similarity (sim) tests result significant (P < 0.05) when niche overlap is less equivalent/similar than randomly expected (niche divergence; D), or more equivalent/similar than randomly expected (niche conservatism; C). NA ¼ No Data. 1997) must be better interpreted as sharing two cytotypes (2n ¼ 16 or 18, with, respectively, two or one submetacentric long chromosome pairs). This new view suggests some evolutionary divergence at the cytogenetic level in this species, in contrast to the reported constant 2n ¼ 16 or 32, for the remaining European taxa of ser. Fissa (Table A4 and Bosch et al., 2016), even for D. mansanetianum, are reported here as 2n ¼ 16 for the first time, opening a new line of research into karyotype evolution and the significance of such changes, both at intra-and interpopulation levels. These two cytotypes in D. bolosii could also be interpreted as a reservoir of diversity and an important link in the patterns of diversification of the series Fissa in the W Mediterranean (Blanch e et al., 2014). This change in chromosome number (2n ¼ 16 vs. 18) is considered as the basis of diversification of the genus Delphinium (Jabbour and Renner, 2012a).

Niche modeling and niche comparison
The obtained models seem to finely accord with real distribution patterns of the three taxa examined, and their predictive character has been (blind) verified by new population discoveries published in recent times; such results were known by our team only after completion of modeling a few months ago, and lie precisely over the expected areas: D. bolosii in Segre River, Noguera, Catalonia (Turmo, 2017); D. fissum subsp. sordidum in Mines of Santo Adrião, Tr as-os-Montes, Portugal  and D. fissum subsp. fissum in Gorges de Caranç a, Pyren ees Orientales, France (Andrieu et al., 2010). Indeed, niche modeling has been proved to be a helpful tool in plant conservation planning, both in understanding distribution of rare plant species in the Iberian Peninsula (Draper, 2010), as previously performed for D. fissum subsp. sordidum (Rus et al., 2018) and for locating new populations of several rare endangered species in the mountains of NE Iberian Peninsula (Buira, 2016).
In Delphinium ser. Fissa, in addition, the identification of the main influential factors explaining distributions through the PCA-env supports the delimitation of subspecies/species through habitat specialization where morphology or genetic signals are weak. Although niches cannot be regarded as totally different among the studied taxa of ser. Fissa, there are signals of niche divergence, especially regarding extreme climatic factors, such as the precipitation of the driest quarter (bio17) or the precipitation seasonality (bio15), the highest contributing variables to the PC1 (59.24%). For example, occurrences of D. bolosii are characterized by much larger precipitations during the driest quarter compared to those of D. fissum subsp. sordidum (mean ± SD ¼ 120.7 ± 13.4 mm vs. 63.8 ± 20.0 mm) but, at the same time rather lower compared to D. fissum subsp. fissum (162.6 ± 45.2 mm; Fig. A4). Regarding the precipitation seasonality, populations of D. fissum subsp. sordidum show the highest values of coefficient of variation (39.5 ± 7.0; D. fissum subsp. fissum ¼ 21.4 ± 6.6; D. bolosii ¼ 24.9 ± 1.2; Fig. A4). This pattern may highlight an ecological differentiation among Delphinium taxa, the Iberian populations of D. bolosii and D. fissum subsp. sordidum being better adapted to drier climates with a lower water supply than populations of D. fissum subsp. fissum in the Alps. Interestingly, we detected that the three taxa shared a relatively high proportion of climatic backgrounds (i.e. analog climates). However, the three entities seem to be occupying different climatic spaces, thus actively selecting distinct environments. This fact could be attributed to a true niche divergence scenario (Atwater et al., 2018), as a consequence of the action of a probable ecological speciation process (Peterson, 2011).
The complex niche evolution patterns detected in Delphinium ser. Fissa, as reported for other Iberian groups like the genus Aquilegia (Jaime et al., 2015), may be adaptive responses to the quick and deep climate changes in the Iberian Peninsula during the Pleistocene (ca. 2.59e0.01 My). The diversification of Delphinium ser. Fissa seems very recent, probably starting after the onset of the late Neogene climate cooling, as indicated by molecular clock dating analysis that reported speciation within the series Fissa around mid-Pleistocene, characterized by abrupt climatic changes (Jabbour and Renner, 2012a).

Differentiation and speciation
The Delphinium ser. Fissa group is assumed to be of E Mediterranean/Irano Turanian origin and arrived in the W Mediterranean following post-Messinian crisis paths, after 5.33 My (Blanch e and Molero, 1983;Blanch e, 1991;Ramírez-Rodríguez et al., 2016). A careful review of populations' precise locations, following field exploration, reveals a particular concentration of populations in soil bags among rocky outcrops following river valleys and narrow gorges, which can indicate migration routes and/or preserved stands, as a result of refugia outputs of such habitats in the Mediterranean (M edail and Diadema, 2009). Barrier analysis presented here confirms the role of mountain massifs and river valleys as steps and nodes of an EeW migration, reinforcing the hypothesis of an expansion from the Alps to the Iberian Peninsula (Orellana et al., 2007b;Ramírez-Rodríguez et al., 2016).
Differentiation was manifested through multiple layers of diversity: geographic and ecological differentiation (temperature but, especially, rainfall; see above), plus a layer of chromosome diversity (at population level) and allele diversity (see also above), and a further layer of flower morphology and pollination biology diversity (Guerrant, 1982;Blanch e, 1991;Bosch, 1999;Jabbour and Renner, 2012b), by a combination of both allopatric and sympatric speciation but with scarce genetic differentiation consequences (low allozyme and cpDNA diversity). Thus ecological differentiation should have played an important role, as early Delphinium biosystematics work suggested (Blanch e, 1991;Pawlowsky, 1963), in speciation, whereas genetic data have been unable to neatly resolve the ser. Fissa species group. As is shown by Jabbour and Renner (2012a), diversification of the terminal clade branches (i.e. last node driving to current ser. Fissa species) is poorly resolved, and has been dated (by molecular clock techniques) at 1.39 My. Supplementary evidences of low genetic differentiation between species come from the easy production of hybrids under greenhouse conditions, well known by larkspur gardeners and reportedly from Delphinium ser. Fissa and ser. Pentagyna (Bosch, 1999;Ramírez-Rodríguez and Amich, 2017) after breaking allopatric isolation by human transport. Such pattern of mainly ecological and relatively recent speciation is well known from other insect/hummingbird pollinated Ranunculaceae such as Aquilegia (Hodges et al., 2004;Kramer and Hodges, 2010). In this genus, Li et al. (2014) suggest that natural selection has facilitated the formation of distinct genetic variation patterns and habitat adaptation has been driving the ecologically based evolution of reproductive isolation.

Conservation
Data reported here could contribute to a better understanding of the conservation issues of this set of related species and subspecies that is more aligned with scientific evidence, and provide some guidelines for designing appropriate conservation action plans.
Delphinium ser. Fissa in the W Mediterranean is a group of four taxa, which (excepting D. fissum subsp. fissum, the only one non-endemic) are considered as threatened (from "CR" d critically endangered to "EN" d endangered, following 3.1 IUCN categories) at country/region level and are protected accordingly (Table A1). Similar status is also applied to representatives of the series in the E Mediterranean region, such as Delphinium fissum subsp. caseyi (B.L. Burtt) C. Blanch e & Molero, a Cypriot endemic evaluated as CR by Kadis et al. (2006). In all cases, danger is mainly reported as a combination of rarity (small population sizes, low subpopulation number, reduced Area of Occupancy [AOO] or Extent of Occurrence [EOO]) linked to habitat alteration or destruction risks. No significant biotic threats are reported under good habitat conservation, even considering the high dependency of these species on fragile mutualistic assemblages, such as insect pollination (Bosch, 1999). Although competition with blackberry bushes' flower resources (pollen, attractiveness) was observed in D. bolosii, the main threats are due to anthropogenic causes (changes in agricultural management) (Orellana et al., 2008). Even strong interanual variations in some adaptive characteristics (such as germination rate in Delphinium fissum subsp. sordidum, see Herranz et al., 2010) could be interpreted as intrinsic to the species and not recorded as threatening. Low to moderate levels of genetic diversity and heterozygosity, as reported here, did not imply a loss of fitness, at least in the surveyed cases like D. bolosii populations (Orellana et al., 2007a), where well-balanced and structured genetic compartmentalization and multilocus genotypes distribution were found at population level (Bosch, 2009: Fig. 11;L opez-Pujol et al., 2014: Fig. 6). Except for very extreme cases, such as D. mansanetianum, showing a uniquely isolated population with only 26 reproductive individuals (Bosch et al., 2005), extremely low allozyme diversity (Table 3), and below any minimum viable population (MVP) standard (Jamieson and Allendorf, 2012), good adaptation to moderate levels of genetic diversity could be assumed for the overall group, with working connections among certain populations, as gene flow data show. Allelic richness in combination with levels of overall genetic diversity and threat due to demographic limitations or habitat destruction could help to identify priority units for conservation, such as the Relevant Genetic Units for Conservation (RGUC) applied to other European endangered taxa belonging to Androcymbium (Caujap e-Castells and Pedrola-Monfort, 2004), Boleum (P erez- Collazos et al., 2008) or Astragalus (Peñas et al., 2016).
Demographic monitoring essentially reported: 1) long-term stability of population censuses; 2) significant fluctuations of annual rates of reproductive/vegetative individuals following available resources and climate conditions (Orellana et al., 2008;Bosch et al., 2009;Rocha et al., 2011;Blanch e et al., 2014;Ramírez-Rodríguez et al., 2016Bosch et al., 2019, and data therein and unpublished results); and 3) a combined vegetative/sexual strategy of perpetuation. Total numbers of emerged individuals can also fluctuate, sometimes even dramatically, but a significant part of any population consists of a subterraneous dormant rhizome bank contributing to further demographic recoveries. These fluctuations must not be considered as reduction in terms of Red List Assessments (IUCN, 2017) and no consistent patterns of decline have been reported for any of the four considered larkspur taxa in the last 15 years. Resulting from recent field work, the number of known individuals and subpopulations is continuously increasing (except for D. mansanetianum), including the rediscovery, in locations where they had been apparently extinct for more than a century, of both D. bolosii in Sant Llorenç Massif (Blanch e and Bosch, 2014) and D. fissum subsp. fissum in the Pyrenees (Andrieu et al., 2010;Lewin, 2015), and the publication of new records for D. fissum subsp. sordidum for the W Iberian Peninsula (Ramírez-Rodríguez et al., 2016. We attributed these findings to the focus on these particular species derived from their appearance in Red Books, thus attracting the attention of plant explorers, as well as from regular population monitoring by government agencies in charge of conservation. The D. fissum group could thus be regarded as being composed of units under incomplete genetic isolation. Despite its old origin, this situation could be similar to that of neoendemics or groups of very recent origin and still poorly differentiated. IUCN criteria for risk evaluation (IUCN, 2012(IUCN, , 2017 may not be sufficiently refined for this pattern of diversity but it is a controversial question. At the same time, it seems desirable to preserve phylogenetic diversity (see key discussion in Vane- Wright et al., 1991;Faith et al., 2004;Forest et al., 2015) also at species level, and in active evolution processes, as the guarantee for future self-sustainable survival, that could be named "terminal clade conservation" (if appropriate, then, supraspecific taxa could also be considered as conservation units, as infraspecific are). But active management of genetically poorly isolated populations (with potential hybridization capacity) should also be considered a risky practice, to include prevention of growing closely related species/populations in the same botanical garden or in the field as a result of introductions (e.g. Ward, 2014, for Delphinium pavonaceum).
The above text is a general depiction of the main conservation issues for the ser. Fissa in the W Mediterranean, taken as a whole and considering the relevant available scientific data. However, this approach is not accepted by the IUCN Red List Authority which clearly states that "Assessments of higher taxa (i.e. above the species level) may NOT be included on the IUCN Red List" (IUCN, 2017); a well-established principle (certainly recommended for practical reasons) but one that needs to be open to wider and more plastic interpretations when scientific considerations can be argued (i.e. when gene flow among subpopulations and taxa can be inferred, as is the case). However, as the species level is an intrinsic characteristic of the IUCN Red List criteria, the IUCN cannot be criticized for this. On the contrary, following criteria similar to or modified from those of the IUCN could lead to a proposal for some new, better criteria to evaluate supraspecific groups. Further comparative integration of similar approaches in other taxonomic groups will require additional attention in order to produce such a new evaluation system. As only species or infraspecific taxa can currently be assessed by IUCN Red List tools, the two following survival factors can be remarked on taking account of the available data reported both in this paper and in the literature, thus: a) Intrinsic factors. In summary, a noticeable increase of AOO and EOO for each species (and for the whole group) is clear, which should produce an obligatory improvement in the conservation status for A, B, C or D criteria (IUCN, 2017: Table 2.1), imposing category changes due to improved knowledge (i.e. a New information transfer reason; not due to a genuine recovery of conservation condition; see IUCN, 2017: 12). Such category changes are also required because the total population for each species is over the threshold for CR or EN categories, and considering the rescue effect of geographically close locations to be taken into account in regional assessments (Gardenfors et al., 2001;Blanch e and Bosch, 2014). Certain cases, such as D. mansanetianum, however, must be maintained in the CR category (D criteria), due to the extremely low number of individuals (<25 adult individuals ; Bosch et al., 2005; and unpublished data). b) Extrinsic factors. Reported documented threats to the four taxa all come from anthropic disturbances determining the loss of individuals and decrease both of habitat extent and habitat quality. These include the overcollection of plants as ornamentals, frequentation, outdoor sports effects such as climbing, running, walking or bullfight exhibitions, overpasture and flock management, changes in forest/agricultural practices, pesticides, opening or widening of roads, fires and floods (Kadis et al., 2006;Blanch e et al., 2007;Andrieu et al., 2010;Rocha et al., 2011;Li et al., 2014;Ramírez-Rodríguez et al., 2017; and according to many unpublished results from field experience). Reports of these episodes and incidents are increasing, and stochastic events should also be expected, including climate change expected effects, e.g. temperature increases inducing seeds to additional dormancy (Herranz et al., 2010). Thus, a perception of risk must not be overlooked and a careful analysis should be made of the results accumulated after some years of application of the recovery plan for the population COR of D. fissum susbp. sordidum in Castilla-La Mancha (Spain), which began in 2002 (DOCM, 2002;Herranz et al., 2010).
Considering these two main drivers, the necessary threat category changes should be adopted in the near future and the derived practical conservation measures appropriately undertaken for each taxon without forgetting the IUCN recommendation that interactions among sub-units should be carefully considered when planning conservation actions (IUCN, 2017: 7). New assessments from updated field information (demography, current threats, etc.) will be an excellent opportunity to explore new conservation cooperative approaches between local, regional, national and international administrations.
In fact, recent contributions to tackling unclear conservation status designation among distinct territories have been developed by complementary approaches (probabilistic, demographic) (Nantel et al., 2018, and references therein) and we encourage the use of an integrative approach to determining the threat status of species and groups of species. However, an increasing amount of literature underlines the shifting of ecological niche, population density, reproductive performance, and significant differences in many other population-related features, between the center and the margin(s) of species' ranges. This has been recently explored under a theoretical approach (Pironon et al., 2017), and with real data for 11 Mediterranean plants (Papuga et al., 2018), among others. Often, the populations at the margins and those at the center of a species' range are included within different administrative units. This factor should be taken into consideration to ensure that significant conservation units' maintenance is compatible with a single threat assessment and a unique conservation strategy.
The results presented herein could lead in the future to a supranational recovery program for Delphinium ser. Fissa in Western Europe. National or regional borders should not necessarily be a limitation on proposing technical, supranational rules. For instance, the Council of Europe, through the Bern Convention, financed the drafting of a European recovery strategy for Cypripedium calceolus (Terschuren, 1999). Although countries are not obliged to comply with a supranational plan, it can play an important role as a guideline for the development of national or regional recovery plans for the targeted species.