diversity Article Reusing Old and Producing New Data Is Useful for Species Delimitation in the Taxonomically Controversial Iberian Endemic Pair Petrocoptis montsicciana/ P. pardoi (Caryophyllaceae) Neus Nualart 1 , Sonia Herrando-Moraira 1 , Eduardo Cires 2,3 , Moisès Guardiola 4 , Emilio Laguna 5 , David Pérez-Prieto 1 , Llorenç Sáez 4 and Jordi López-Pujol 1, * 1 Botanic Institute of Barcelona (IBB, CSIC-Ajuntament de Barcelona), Passeig del Migdia, s/n, 08038 Barcelona, Spain; nnualart@ibb.csic.es (N.N.); sonia.herrando@ibb.csic.es (S.H.-M.); perezprieto@hotmail.com (D.P.-P.) 2 Departamento de Biología de Organismos y Sistemas, Universidad de Oviedo, Área de Botánica. C/Catedrático Rodrigo Uría s/n, 33071 Oviedo, Spain; cireseduardo@uniovi.es 3 Instituto de Recursos Naturales y Ordenación del Territorio (INDUROT), Campus de Mieres, C/Gonzalo Gutiérrez Quirós s/n, 33600 Mieres, Spain 4 Systematics and Evolution of Vascular Plants (UAB)—Associated Unit to CSIC, Departament de Biologia Animal, Biologia Vegetal i Ecologia, Facultat de Biociències, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain; moises.guardiola@uab.cat (M.G.); Llorens.Saez@uab.cat (L.S.)  5  Generalitat Valenciana, Servicio de Vida Silvestre—CIEF (Centro para la Investigación y Experimentación Forestal), Avda. Comarques del País Valencià, 114, 46930 València, Spain; laguna_emi@gva.es Citation: Nualart, N.; * Correspondence: jlopez@ibb.csic.es Herrando-Moraira, S.; Cires, E.; Guardiola, M.; Laguna, E.; Abstract: Petrocoptis montsicciana and P. pardoi are two Iberian endemic taxa of Caryophyllaceae family Pérez-Prieto, D.; Sáez, L.; López-Pujol, with an unclear taxonomic delimitation, being variously treated as independent species, subspecies J. Reusing Old and Producing New or even synonyms. In the present study, allozyme raw data obtained in the early 2000s have been Data Is Useful for Species reused with improved tools to survey genetic structure, and complemented with modeling and niche Delimitation in the Taxonomically Controversial Iberian Endemic Pair comparative analyses to shed light on species delimitation. Genetic structure was investigated using Petrocoptis montsicciana/P. pardoi four approaches: Bayesian clustering, Monmonier’s algorithm, Principal Coordinate Analysis (PCoA), (Caryophyllaceae). Diversity 2021, 13, and Analysis of Molecular Variance (AMOVA). Ecological niche differences have been assessed 205. https://doi.org/10.3390/ through Ecological Niche Modeling (ENM) using MaxEnt, and Principal Component Analysis using d13050205 both occurrence records and background climate (PCA-env). Genetic analysis confirms the distinction between both taxa, and the scenario of a progenitor–derivative (P–D) is suggested. In agreement Academic Editor: Julio Peñas de Giles with genetic data, niche analysis shows clear differences between their climate regarding species occurrences and background spaces. Climate divergence could be explained, at least partially, by the Received: 15 April 2021 abundance of rocks where species live although differences at the microclimate instead of the regional Accepted: 10 May 2021 climate should be explored in future research. Given the genetic distinction between P. montsicciana Published: 12 May 2021 and P. pardoi, both taxa should be regarded as separate ‘Management Units’ (MUs). Publisher’s Note: MDPI stays neutral Keywords: chasmophytic species; conservation; endemic species; genetic diversity; niche analysis with regard to jurisdictional claims in published maps and institutional affil- iations. 1. Introduction Petrocoptis A. Braun ex Endl. (Caryophyllaceae) is one of the flagship genera of Copyright: © 2021 by the authors. higher plants of the Iberian Peninsula, given to a series of reasons that include: (1) its Licensee MDPI, Basel, Switzerland. alleged endemic status; (2) the conservation concerns of several of its members (some are This article is an open access article narrow endemics and highly threatened); and (3) the beauty of the plants and the habitats distributed under the terms and where they grow (Figure S1). Petrocoptis is one of the 27 endemic genera in the Iberian conditions of the Creative Commons flora, with nine species [1], although it should be included within Silene L. according to Attribution (CC BY) license (https:// some authors [2]. However, recent molecular studies (based both on nuclear and plastid creativecommons.org/licenses/by/ sequences) have confirmed its independence from Silene [3], having probably diverged 4.0/). approximately in Late Miocene, ca. 10.6 Mya [4]. Most of Petrocoptis taxa have restricted Diversity 2021, 13, 205. https://doi.org/10.3390/d13050205 https://www.mdpi.com/journal/diversity Diversity 2021, 13, 205 2 of 14 ranges with small, isolated populations, which are the main reasons for their inclusion in red lists and red data books as well as protection catalogues ([4] and references therein). In addition to stochasticity associated with the small size of populations, rock-climbing also represents a very important threat for this genus [5]. Commonly known as clavel de roca (rocky carnation in Spanish), Petrocoptis representatives are perennial herbs with red, pink-purple or white pentamerous flowers that grow on rocky cracks and fissures in calcareous vertical rock walls, overhangs, narrow gorges and caves (Figure S1). The delimitation of species within Petrocoptis is unclear, with the number of species ranging from four to 12 species depending on the taxonomical treatment (see Table 1 in [4]). The two easternmost-distributed taxa (Petrocoptis montsicciana O. Bolòs & Rivas Mart. and P. pardoi Pau) are a good example of such taxonomic inconsistency, having been treated as independent species (e.g., [6,7]), as two subspecies of P. crassifolia Rouy [8] or even as synonyms, both under P. pardoi [9] or under Silene pardoi (Pau) Mayol & Rosselló subsp. pardoi [2]. Petrocoptis montsicciana and P. pardoi are, indeed, morphologically hardly distinguishable (Figure S1), with very small differences in the shoots and in the hairs of the strophiole according to Flora iberica [7], or with slightly different size of the seeds according to [6,8]. However, the most recent treatment of the genus, which is based only on molecular markers, suggest that—pending additional studies to fully resolve infrageneric relationships—both P. montsicciana and P. pardoi might deserve the rank of species [4]. These two species were also the object of a full genetic diversity study using allozymes [10], which suggested that P. montsicciana and P. pardoi were genetically close but different species, likely constituting a progenitor–derivative species pair. In cases where species delimitation is problematic, such as in the P. montsicciana/P. par- doi pair, a multidisciplinary approach (also known as “integrative taxonomy” [11,12]) is highly advisable, with the combination of morphological, molecular and ecological data having provided notable success in several case studies (e.g., [13,14]). Herein for the first time in this genus we are combining niche (including both modeling and niche comparative analyses) with genetic analyses, to get additional insights into the relationships between P. montsicciana/P. pardoi. The genetic analyses will consist of re-analyzing the raw allozyme data of [10] because tools to survey genetic structure have improved significantly during the last two decades, mostly thanks to the appearance of programs that employ Bayesian clustering algorithms and edge detection techniques [15]. We have chosen a representative of each of these two families of methods, as well as other additional approximations not employed in [10], to estimate genetic divergence at the population level between the two species. Reusability of data is a growing policy priority in science, often linked to the poli- cies of data sharing and open data [16,17], and it is particularly helpful in the context of the current COVID-19 pandemics, because field sampling and carrying new experiments could be extremely difficult [18] even if field investigation remains crucial for the conservation of threatened species, especially for those with a small distribution range [19,20]. 2. Materials and Methods 2.1. Genetic Analysis The spatial genetic structure for the pair P. montsicciana/P. pardoi has been assessed using the same raw genetic data as that of [10]. These raw data are the genotype matrix of a total of 223 individuals belonging to seven populations (Figure 1), four of P. montsicciana (N = 140) and three of P. pardoi (N = 83), and using 16 allozyme loci (Aat, Aco-1, Aco- 2, Adh, Dia-2, Dia-3, Mdh-1, Mdh-4, Me, 6Pgd-1, 6Pgd-2, Pgi-2, Pgm-1, Pgm-2, Prx-1 and Prx-2). Population genetic structure was investigated using four different approaches. First, the Bayesian clustering method implemented in the software Structure v.2.3.4 [21] was run assuming the admixture model with correlated allele frequencies. Based on exploratory runs, we restricted the number of assayed clusters (K) to range K = 1–8 and each K was estimated 20 times with a length of burn-in period of 105 iterations and 106 Markov chain Monte Carlo (MCMC) replications. The optimal number of clusters was determined using the ∆K statistical approach by [22], with the aid of Structure Harvester Bayesian clustering method implemented in the software Structure v.2.3.4 [21] was run assuming the admixture model with correlated allele frequencies. Based on exploratory runs, we restricted the number of assayed clusters (K) to range K = 1–8 and each K was Diversity 2021, 13, 205 3 of 14 estimated 20 times with a length of burn-in period of 105 iterations and 106 Markov chain Monte Carlo (MCMC) replications. The optimal number of clusters was determined using the ΔK statistical approach by [22], with the aid of Structure Harvester v.0.6.94 [23]. Pro- grams v.0.6.94Clumpp v.1.1.2 [24] [23]. Programs and Distruct Clumpp v.1.1.2v.1.1 [24][25] andwere used for Distruct the[25] v.1.1 post-processing were used for of the the results obtained from post-processing of theStructure. Second,from results obtained the location Structure.of potential Second, genetic barriers the location between of potential populations genetic barrierswasbetween exploredpopulations through the Monmonier’s was explored throughalgorithm implemented algorithm the Monmonier’s in Barrier v.2.2 [26]. Theinsignificance implemented Barrier v.2.2of[26]. barriers was tested of The significance forbarriers nSSR data by bootstrapping was tested for nSSR data 1000 by Nei’s DA genetic bootstrapping distances 1000 Nei’s D[27] A matrices genetic that distanceswere [27]previously matrices obtained that were with Microsatellite previously obtained Analyzer (MSA) v.4.05 with Microsatellite [28]. Third, Analyzer (MSA)GenAlEx v.6.1Third, v.4.05 [28]. [29] was used to GenAlEx perform v.6.1 a Principal [29] was used to Coordinate AnalysisCoordinate perform a Principal (PCoA) based on codominant Analysis (PCoA) based genotypic distances on codominant both at distances genotypic the indi- vidual both at and populationand the individual levels. Fourth,levels. population and Fourth, last, an andAnalysis last, anofAnalysis Molecular Variance of Molecular (AMOVA) using theusing Variance (AMOVA) program ARLEQUIN the program v.3.11 [30]v.3.11 ARLEQUIN was performed. AMOVA AMOVA [30] was performed. was car- was under ried carriedtwounder two hypotheses: hypotheses: (1) all populations (1) all populations belong tobelong a singletospecies; a singleand species; and (2) taking (2) taking into account into theaccount existencetheofexistence of two two different different species species and and therefore therefore three establishing establishing hierar- three hierarchical chical levels among levels of analysis: of analysis: amonggroups, taxonomic taxonomic amonggroups, among populations populations within species within and speciespopulations. within and within populations. Location Figure 1.1.Location of seven of the the seven studied studied populations populations of the of the pair pair Petrocoptis Petrocoptis montsiccianamontsicciana (red (red points)/P. pardoi (bluepardoi points)/P. points) for the (blue genetic points) foranalysis plusanalysis the genetic the rest of their plus theoccurrences rest of theirused in the niche occurrences usedanal- in ysis. (A) P. the niche pardoi:PP-AGV, analysis. Aiguaviva (A) P. pardoi: PP-AGV,rocks, Valencian Aiguaviva Community, rocks, 40◦ 47W, 40°47′ N, 00°09′ Valencian Community, 0 N, N00 ◦ 090PP- = 22; W, BAR, N = 22;Barranc PP-BAR, de Barranc la Mare dede la Déu rocks, Mare Valencian de Déu Community, rocks, Valencian 40°44′ N, 40 Community, 00°11′ W, 00 ◦ 440 N, N◦=1130; PP-CBD, 0 W, N = 30; Cantal Badat rocks, Valencian Community, 40°46′ N, 00°10′ W, N = 31. (B) P. montsicciana: PP-CBD, Cantal Badat rocks, Valencian Community, 40◦ 460 N, 00◦ 100 W, N = 31. (B) P. montsicciana: PM-CAM, Camarasa gorges, Catalonia, 41°53′ N, 00°53′◦E, N = 31; PM-MCO, La Móra Comdal rocks, Catalonia, PM-CAM, Camarasa gorges, Catalonia, 41 530 N, 00◦ 530 E, N = 31; PM-MCO, La Móra Comdal rocks, Catalonia, 42◦ 070 N, 01◦ 210 E, N = 42; PM-MRB, Mont-rebei pass, Catalonia, 42◦ 040 N, 00◦ 400 E, N = 46; PM-TER, Terradets gorges, Catalonia, 42◦ 020 N, 00◦ 530 E, N = 21. 2.2. Niche Analysis Ecological niche differences between P. montsicciana and P. pardoi have been assessed both in geographical (G) and environmental (E) space. In G-space we performed an Ecological Niche Modeling (ENM) using the maximum entropy algorithm implemented Diversity 2021, 13, 205 4 of 14 in MaxEnt v.3.3 [31], which requires presence data only. In E-space we used the approach developed by [32], in which background climatic areas around occurrence records are used to perform a principal component analysis (hereafter PCA-env). The georeferenced occurrences were compiled from: (1) the maps of P. montsicciana and P. pardoi included in [9] and which were drawn using personal, literature and herbarium records; (2) the list of occurrences of P. montsicciana obtained during the elaboration of the Spanish Red Book of threatened vascular flora [33]; (3) the map of P. montsicciana elaborated during the project LIFE RESECOM [34]; (4) the list of occurrences of P. pardoi included in the Banc de Dades de Biodiversitat de la Comunitat Valenciana (BDBV, http://www.bdb.gva.es/; accessed on 29 December 2020); and (5) the records of P. pardoi included in [35]. All occurrences were originally UTM (Universal Transverse Mercator) squares of 1 × 1 km and have been transformed to geographic latitude and longitude coordinates. As we do not know the datum of the source data, we transformed them using the midpoint of the UTM square to minimize the datum effect. In total, after removing duplicate records within each pixel [30 arc-sec (ca. 1 km)], 60 occurrences of P. montsicciana and 31 of P. pardoi were used in the niche analyses (see Table S1 for the list of all occurrences and Figure 1 for the map). A set of 19 bioclimatic variables was downloaded from the WorldClim website (www.worldclim.org; accessed on 15 January 2021) under current climatic conditions (1970– 2000). Despite that both taxa are strictly rocky plants, this variable (rocky vs. non-rocky substrates) has not been considered due to technical limitations; at the geographic scale we are working (ca. 1 km), and given the nature of rocky outcrops (sometimes they are just occupying a few dozen meters), it is not reliable to assign cells as of ‘rocky substrate’ or ‘non-rocky substrate’. Therefore, our estimation of the ecological niche of the two taxa should be equated to their climate niches. A pairwise Pearson correlation analysis using the “SDM Toolbox” extension for Ar- cGIS [36] was conducted to retain only seven relatively uncorrelated climatic variables (r < |0.85|): isothermality (bio3), minimum temperature of coldest month (bio6), tempera- ture annual range (bio7), mean temperature of wettest quarter (bio8), mean temperature of driest quarter (bio9), precipitation seasonality (bio15), and precipitation of coldest quarter (bio19). For the variables selection we made a consensus to include the most influential variables for both taxa based on their relative contribution to the models (percentage contribution, permutation importance and jackknife of regularized gaining train), and an exploration of curve response shape. For the final model (i.e., using only seven variables), MaxEnt was run 100 times using the “subsample” method (with 25% of the localities randomly selected to test the model). The AUC (area under the curve) values were used as a metric to assess models’ performance. The resulting maps were modified and ex- ported with ArcGIS v.10.2 [37] considering as cut-off value the maximum sensitivity plus specificity threshold (MSS) as recommended by [38]. To perform the PCA-env we followed [39], where the background for each species was selected from a minimum convex polygon (convex hull) with a buffer size of 0.1 degrees. Original occurrences were smoothed using a kernel density function and then projected in the 500 × 500 gridded environmental space. Niche differences between P. montsicciana and P. pardoi were measured using the niche overlap metric of Schoener’s D [40,41], which ranges from 0 (representing no overlap) to 1 (complete overlap). The E-space analyses were performed using the original R code reported in [32] and later modified by [39], and incorporating last modifications reported in [42]. Additionally, we extracted individual values of the selected variables for the occurrences with ArcGIS, in order to compare the species niche preferences through a boxplot representation and statistical test. First, the Shapiro–Wilk test was performed to check the normality of variables. As all were not normally distributed (p-value is less than alpha level, 0.05), the non-parametric analysis Wilcoxon-Mann-Whitney was finally conducted. Both analyses were performed in R v.4.0.3 [43] using the package stats v.4.0.3. the species niche preferences through a boxplot representation and statistical test. First, the Shapiro–Wilk test was performed to check the normality of variables. As all were not normally distributed (p-value is less than alpha level, 0.05), the non-parametric analysis Wilcoxon-Mann-Whitney was finally conducted. Both analyses were performed in R Diversity 2021, 13, 205 5 of 14 v.4.0.3 [43] using the package stats v.4.0.3. 3. Results 3.3.1. Results Genetic Structure and Divergence between Taxa 3.1. Genetic Based Structure on the ΔKand Divergence approach, K =between Taxa 3 was the most likely number of genetic clusters after Based running the the ∆K approach, onsoftware Structure,Kalthough = 3 was K the = most 2 alsolikely number provided good ofresults genetic(Figure clusters2A). afterIn running the software Structure, although K = 2 also provided good the K = 2 grouping scheme, all populations of P. montsicciana (PM-CAM, PM-TER, PM-results (Figure 2A). In the MRB K= 2 grouping and PM-MCO;scheme, all populations for population of P. codes, see montsicciana Figure 1) were (PM-CAM, merged in aPM-TER, PM- single cluster, MRB andthose whereas PM-MCO; for population of P. pardoi codes, see (PP-CBD, PP-BAR Figure and 1) were PP-AGV) merged formed in a single the second cluster, cluster (Fig- whereas ure 2B). those With K of=P.3,pardoi (PP-CBD, the cluster of P. PP-BAR and was montsicciana PP-AGV) formed split into two: the onesecond cluster containing the (Figure 2B). With K = 3, the cluster of P. montsicciana was split into southernmost populations of the species (PM-CAM, PM-TER) and the other with the two: one containing the southernmost northernmost populations ones (PM-MRB and of the species (PM-CAM, PM-MCO), whereas allPM-TER) and the other the populations of P. with pardoithere- northernmost ones (PM-MRB and PM-MCO), whereas all the populations mained in the same genetic cluster (Figure 2B). Barrier results also indicated the lack of P. pardoi of remained in the same genetic differences genetic the between cluster (Figure 2B). populations Barrier of P. pardoi,results whilealsothe indicated the lack largest barriers of were genetic differences put between the twobetween taxa (the the populations barrier P. pardoi, and betweenofPM-CAM while the largest PP-AGV barriers received up towere 92% put bootstrap support; i.e. this barrier appeared in 920 out of 1000 matrices, Figure to between the two taxa (the barrier between PM-CAM and PP-AGV received up 3).92% The bootstrap support; i.e. this barrier appeared in 920 out of 1000 matrices, Figure 3). The PCoA analysis had relatively high values of percent of variance explained by first two PCoA analysis had relatively high values of percent of variance explained by first two axes (68.63% and 56.12% at the population and individual level, respectively). The PCoA axes (68.63% and 56.12% at the population and individual level, respectively). The PCoA plot at the population level showed a certain level of genetic divergence between the two plot at the population level showed a certain level of genetic divergence between the two species, as populations of P. montiscciana and those of P. pardoi appeared in the right and species, as populations of P. montiscciana and those of P. pardoi appeared in the right and left quadrants, respectively (Figure 4A). Such pattern, however, was less clear when we left quadrants, respectively (Figure 4A). Such pattern, however, was less clear when we plotted all the studied individuals (because a few individuals of PM-MCO appeared at plotted all the studied individuals (because a few individuals of PM-MCO appeared at the left side of the first PCoA axis; Figure 4B). Both assays of AMOVA (Table 1) revealed the left side of the first PCoA axis; Figure 4B). Both assays of AMOVA (Table 1) revealed that the maximum percentage of variation was found within populations irrespective if that the maximum percentage of variation was found within populations irrespective if considered together (55.49%) or as belonging to two species (50.87%). When sorting the considered together (55.49%) or as belonging to two species (50.87%). When sorting the populations by species, the percentage of variation between P. montsicciana and P. pardoi populations by species, the percentage of variation between P. montsicciana and P. pardoi was18.60% was 18.60%(Table (Table1). 1). Figure2.2.Results Figure Resultsof ofStructure Structurev.2.3.4 v.2.3.4analysis analysisfor forthe the pair pair Petrocoptis Petrocoptis montsicciana/P. montsicciana/P. pardoi. pardoi. (A)(A) graph- graphical ical representation of ∆Kofvalues representation ΔK values obtained obtained fromfrom Structure Structure Harvester Harvester v.0.6.94; v.0.6.94; (B)(B) assignationof assignation individuals to genetic clusters at K = 2 and K = 3 after running Clumpp v.1.1.2 and Distruct v.1.1. Diversity 2021, 13, x 6 of 14 Diversity 2021, 13, x 6 of 14 Diversity 2021, 13, 205 of individuals to genetic clusters at K = 2 and K = 3 after running Clumpp v.1.1.2 and Distruct 6 of 14 of individuals to genetic clusters at K = 2 and K = 3 after running Clumpp v.1.1.2 and Distruct v.1.1. v.1.1. Figure 3.3.Analysis with Barrier Analysis v.2.2, with the the firstfirst barrier checked. Significance was tested by 1000 Figure Figure 3. Analysis withwith Barrier Barrier v.2.2,v.2.2, withwith the first barrierbarrier checked. checked. Significance Significance was by was tested tested 1000by bootstrap matrices 1000 bootstrap of [27] genetic distance (D A). Numbers indicate bootstrap support. Codes PM- bootstrap matricesmatrices of [27] distance of [27] genetic genetic distance (DA ). Numbers (DA). Numbers indicate bootstrap indicate bootstrap support. support. Codes PM- Codes CAM, PM-MCO, PM-CAM, PM-MRB PM-MCO, PM-MRB and PM-TER and PM-TER are of populations of P.ofmontsicciana, and codes P. montsicciana, PP-AGV, CAM, PM-MCO, PM-MRB and PM-TER are ofare of populations populations of P. montsicciana, and and codes codes PP-AGV, PP-AGV, PP-BAR and PP-CBD of P. pardoi (see Figure 1). PP-BAR PP-BAR andand PP-CBD PP-CBD P. P. of of pardoi pardoi (see (see Figure Figure 1).1). Principal Figure4.4.Principal Figure Coordinate Coordinate Analysis Analysis (PCoA) (PCoA) performed performed fromfrom pairwise pairwise genotypic genotypic distances distances Figure 4. Principal betweenthe Coordinate thestudied studiedpopulations Analysis populationsofofthe (PCoA) thepair performed from pairwise Petrocoptismontsicciana/P. pairPetrocoptis genotypic montsicciana/P.pardoi, distances pardoi,atatpopulation population(A)(A) and between between the studied populations of the pair Petrocoptis montsicciana/P. pardoi, at population (A) at individual and at individual(B)(B) level. Codes level. Codes PM-CAM PM-CAM(red),(red),PM-MCO PM-MCO(yellow), (yellow),PM-MRB PM-MRB (blue)(blue) and PM- PM-TER and at individual (B) level. Codes PM-CAM (red), PM-MCO (yellow), PM-MRB (blue) and PM- TER (green) (green) arepopulations are of of populations of P.of P. montsicciana, montsicciana, and codes and codes PP-AGV PP-AGV (dark(dark red), red), PP-BARPP-BAR (light(light blue) and TER (green) are of populations of P. montsicciana, and codes PP-AGV (dark red), PP-BAR (light blue) and (pink) PP-CBD PP-CBD of (pink) P. of (see pardoi P. pardoi (see Figure 1).Figure 1). blue) and PP-CBD (pink) of P. pardoi (see Figure 1). Table 1.1.Analysis of of molecular variance (AMOVA) in the pair Petrocoptis montsicciana/P. Petrocoptis pardoi montsicciana/P. Table Table Analysis 1. Analysis molecular of molecular variance variance (AMOVA) (AMOVA) in theinpair thePetrocoptis pair montsicciana/P. pardoi pardoi considering considering (1) (1)all allpopulation populationtogether together(one (one group) group)and and(2) that (2) the that populations the populations belong belongtotodifferent different considering (1) all population together (one group) and (2) that the populations belong to different species (two groups). species species (two(two groups). groups). Model Partitioning Variance (%) F-Statistic p Model Model Partitioning Partitioning Variance (%) (%) Variance F-Statistic F-Statistic p p Among populations 44.51 FST = 0.445 <0.001 1 group Among populations Among populations 44.51 44.51 FST = 0.445 <0.001 <0.001 FST = 0.445 1 group 1 group Within populations 55.49 WithinWithin populations populations 55.49 55.49 2 groups: Among species 18.60 FCT = 0.186 <0.01 P. montsicciana— Among populations within species 30.54 FSC = 0.375 <0.001 P. pardoi Within populations 50.87 FST = 0.491 <0.001 Diversity 2021, 13, x 7 of 14 Diversity 2021, 13, 205 7 of 14 Among species 18.60 FCT = 0.186 <0.01 2 groups: Among populations within spe- P. montsicciana— 30.54 FSC = 0.375 <0.001 3.2. Climatic Niches of P. montisicciana cies and P. pardoi P. pardoi Within populations Niche models produced with MaxEnt for both50.87 species hadFSTvery = 0.491 <0.001 high AUC (>0.99), indicating good performance (Table 2). According to jackknife test and percentage contri- 3.2. Climatic bution Nichesbio3, to models, of P.bio7 montisicciana and P.important and bio9 were pardoi for the two Petrocoptis taxa, although for P.Niche pardoimodels bio19 ranked producedfirstwith according MaxEnt to for jackknife test (Table both species had 2). veryThe highpredicted areas AUC (>0.99), of suitablegood indicating habitat for both taxa performance (Figure (Table 5B) included 2). According their respective to jackknife test and occurrence percentage points, contri- although bution to new models,suitable bio3,areas appeared: bio7 and withimportant bio9 were few exceptions, for thethese were located two Petrocoptis in although taxa, the same latitude for thatbio19 P. pardoi the occurrence ranked firstrecords according P. jackknife for to montsicciana, while areas test (Table 2). Thefurther south predicted were areas of evidenced suitable for P. for habitat pardoi (Figure both 5). Therefore, taxa (Figure neithertheir 5B) included of the distribution respective models was occurrence ablealt- points, to predictnew hough the occurrences suitable areas of the other species. appeared: with few exceptions, these were located in the same latitude that the occurrence records for P. montsicciana, while areas further south were Table 2. Occurrences evidenced of Petrocoptis for P. pardoi (Figuremontsicciana 5). Therefore, P. pardoi, of and neither model the performance distributionassessed modelsby the able was area under the curve (AUC ± standard deviation), to predict the occurrences of the other species. threshold values as the maximum training sensitivity plus specificity threshold (MSS), and importance of all variables based on jackknife analyses (in parenthesis Table the percentage 2. Occurrences contribution of Petrocoptis to models). montsicciana and P. pardoi, model performance assessed by the area under the curve (AUC ± standard deviation), threshold values as the maximum training sen- sitivityModel plus specificityOccurrences threshold (MSS), and ± SD AUCimportance ofThreshold all variables based onVariables jackknife anal- yses (in parenthesis the percentage contribution to models). bio7 (24.9) > bio9 (35.2) > bio3 (24.9) > bio15 P. Model montsiccianaOccurrences SD ±Threshold 60 AUC ± 0.994 0.002 0.1355 Variables (0.3) > bio19 (5.2) > bio6 bio7 (24.9) > bio9 (2.7) (35.2)>>bio8 bio3(3.5) (24.9) > P. montsicciana 60 0.994 ± 0.002 0.1355 bio15 (0.3) > bio19 (5.2) > bio6 (2.7) > bio19 (13.0) > bio9 (30.6) bio8 (3.5) > bio3 (31.7) > bio7 (17) P. pardoi 31 0.998 ± 0.001 0.4936(13.0) > bio9 (30.6) > bio3 (31.7) > bio19 > bio6 (3.0) > bio15 (4.4) P. pardoi 31 0.998 ± 0.001 0.4936 bio7 (17) > bio6 (3.0) >> bio15 (4.4) > bio8 bio8 (0.4) (0.4) 5. (A) Distribution Figure 5. Distribution of of occurrence occurrence records records for for P. montsicciana and P.P. pardoi pardoi used used in in the the ecological ecological niche niche modeling modeling analyses and analyses and background background climates climates represented represented by by polygon polygon areas areas used used for for the the PCA-env PCA-env analysis. analysis. (B) (B) MaxEnt MaxEnt distribution distribution models delimitated models delimitated by by the the maximum maximum training training sensitivity sensitivity plus plus specificity specificity (MSS) (MSS) threshold threshold obtained obtained for for the the both both Petrocoptis Petrocoptis taxa for the present time. taxa for the present time. The The first first two twocomponents componentsofofthe PCA-env the PCA-envexplained explained76.54% 76.54%of the totaltotal of the climatic var- climatic iation (Figure 6). For both axes (PC1 = 55.47% and PC2 = 21.07%) the precipitation variation (Figure 6). For both axes (PC1 = 55.47% and PC2 = 21.07%) the precipitation varia- bles bio15 bio15 variables and bio19 and were bio19the most were informative the factors (Figure most informative factorsS2). The 100% (Figure S2). of occurrence The 100% of occurrence density showed non-overlap ranges between P. montisicciana and P. pardoi (Figure 6), indicating divergence in climatic niches. Evidence also came from values of Diversity 2021, 13, x 8 of 14 Diversity 2021, 13, x 8 of 14 Diversity 2021, 13, 205 8 of 14 density showed non-overlap ranges between P. montisicciana and P. pardoi (Figure 6), in- density showed dicatingnon-overlap divergenceranges between in climatic P. montisicciana niches. Evidence also andcame P. pardoi from(Figure values 6), in- of Schoener’s D dicating divergence index (D =in0, Schoener’s Dclimatic no index (Dniches. overlap) and = 0, noEvidence niche also andcame metrics overlap) nichefrom (niches werevalues metrics of Schoener’s non-stabilized, (niches were totallyD unfilled, non-stabilized, index (D =totally 0, no and inoverlap) expansion) unfilled, and niche and(Table metrics 3). When in expansion) (niches were analyzing (Table non-stabilized, the 3). When bioclimatic totally analyzing variables unfilled, individually, the bioclimatic all variables and in expansion) (Table showed significant individually, 3). When analyzing differences all showed the< bioclimatic (p-value significant 0.05)(p-value differences variables < 0.05) inindividually, in Wilcoxon-Mann-Whitney all although test, Wilcoxon-Mann-Whitney showed significant test values test, differences althoughfortest (p-value biovalues 8 and bio9 < 0.05) are for bio andinto 8near Wilcoxon-Mann-Whitney p-value bio9 (Table are near S2 and(Table to p-value test,S2 Figure although 7). and Figure 7). test values for bio 8 and bio9 are near to p-value (Table S2 and Figure 7). Figure Figure6. 6. Global Global climatic climatic space space (PCA-env), (PCA-env), projecting projecting the the realized realized niches of Petrocoptis niches of Petrocoptis montsicciana montsicciana andP.P.climatic and Figure 6. Global pardoi pardoiinin aashading shading space gradation gradation (PCA-env), scale scaleaccording projecting according tototheir the realized theirdensity niches density ofofoccurrence occurrence of Petrocoptis per percell. montsicciana cell.The Thesolid and P. pardoi solid lines lines illustrate in aillustrate shading the the of gradation 100% 100% scaleof occurrence according occurrence density to their density and dashed anddensity dashed oflines lines occurrence the the 100% 100% per of cell. of Thebackground background areas. ar- eas. solid lines illustrate the 100% of occurrence density and dashed lines the 100% of background ar- eas. Table 3. Metrics of niche dynamics comparing species pairs (1 → 2), according to definitions from [32]. Table 3. Metrics of niche dynamics comparing species pairs (1  2), according to definitions from Metrics are interpreted as: overlap (Schoener’s D; 0 = no overlap, 1 = complete overlap), unfilling [32]. Metrics Table 3. Metrics of nicheare interpreted dynamics as: overlap comparing (Schoener’s species pairs (1 D;2), 0 =according no overlap,to 1definitions = completefrom overlap), unfil- (0 = totally filled, 1 = totally unfilled), stability (0 = non-stabilized, 1 = totally stabilized), expansion [32]. Metricsling are(0interpreted = totally filled, 1 = totally as: overlap unfilled),D;stability (Schoener’s (0 = non-stabilized, 0 = no overlap, 1 = complete1 overlap), = totally stabilized), unfil- ex- (0 = in equilibrium, 1 = fully in expansion). Abbreviations used: PM = Petrocoptis montsicciana, pansion ling (0 = totally (0 1= =intotally filled, equilibrium, 1 =stability unfilled), fully in (0 expansion). Abbreviations = non-stabilized, used: 1 = totally PM = Petrocoptis stabilized), ex- montsicci- pansion (0 PP=ana, P. pardoi. in=equilibrium, PP = P. pardoi. 1 = fully in expansion). Abbreviations used: PM = Petrocoptis montsicci- ana, PP = P. pardoi. Species Comparisons Species Niche Overlap Niche Overlap (D) (D) Niche NicheDynamics Dynamics Species Comparisons Sp. Sp. 1 1 Sp.Niche 2 Sp. 2 Overlap (D) Niche Unfilling Unfilling DynamicsStability Stability Expansion Expansion Sp. 1 PMSp. 2 PP 0 Unfilling Stability 1 Expansion 0 1 PM PP 0 1 0 1 PM PPPP PM 0 0 1 10 0 1 1 PP PPPM PM 0 0 1 10 0 1 1 Figure 7. Cont. iversity 2021, 13, x 2021, 13, 205 Diversity 9 of 14 9 of 14 Figure 7. Boxplots of Boxplots Figure 7. the sevenofselected variables the seven related selected with related variables temperature (bio3, bio6, bio7, with temperature bio8, (bio3, bio6, bio7, bio8, bio9) and precipitation (bio15, bio19) where the Petrocoptis montsicciana and P. pardoi currently bio9) and precipitation (bio15, bio19) where the Petrocoptis montsicciana and P. pardoi currently occur. occur. 4. Discussion and Conclusions 4. Discussion 4.1. and Both Conclusions Genetic and Ecological Data Confirm the Distinction between P. montsicciana and P. pardoi 4.1. Both Genetic and Ecological Data Confirm the Distinction between P. montsicciana and P. pardoi The re-analysis of the raw allozyme data of [10] with the four genetic structure approx- imations (Structure, Barrier, PCoA, and AMOVA) confirms the conclusions presented in the The re-analysis of the raw allozyme data of [10] with the four genetic structure ap- previous paper: P. montsicciana and P. pardoi are assigned to two genetically differentiated proximations (Structure, Barrier, PCoA, and AMOVA) confirms the conclusions pre- clusters, although this differentiation is not large. In [10], the authors noted that their sented in the previous paper: P. montsicciana and P. pardoi are assigned to two genetically genetic results were compatible with a scenario of a progenitor–derivative (P–D) species differentiated clusters, although this differentiation is not large. In [10], the authors noted pair, based on both the pattern of genetic structure and allele profiles. In contrast to the that their genetic results were compatible with a scenario of a progenitor–derivative (P– classic allopatric mode of speciation, in a P–D species pair the derived (D) species buds off D) species pair, based on both the pattern of genetic structure and allele profiles. In con- and acquires new traits—or becomes fixed for features within the polymorphism of the trast to the classic allopatric mode of speciation, in a P–D species pair the derived (D) parental (P) species—while the P species remains largely unchanged [44–48]. Although the species buds off and acquires new traits—or becomes fixed for features within the poly- D species may show a series of changes with respect to P that may include new morpho- morphism of the parental (P) species—while the P species remains largely unchanged [44– logical characters, a new breeding system, selection for a new habitat, and chromosomal 48]. Although the D species restructuring [45], may theseshow couldanot series occur of changes withcases, [48]. In such respect to P that genetic may can markers in- be useful clude newtools morphological characters, a new for detecting P–D species pairs [48–50]. breeding system, selection for a new habitat, and chromosomal In therestructuring study of [10],[45], these could ecological not occur [48]. and morphological In such cases, differences between genetic P. montsicciana markers can be useful tools for detecting P–D species pairs [48–50]. and P. pardoi were not explored, as they were assumed to be minimal, if any. Because In thethey study of [10], have alsoecological the same and morphological chromosome number, differences therefore, between P. montsicciana the allozyme patterns were the and P. pardoi only evidence for these authors to support the assignment as a P–DBecause were not explored, as they were assumed to be minimal, if any. species they pair. Certainly, have also the the same allozymechromosome data closely number, matched therefore, the allozyme the expectations: (1) patterns the spectrum were of thealleles only observed evidence forin these authors P. pardoi to support (D species) wasthe assignment a subset of those as afound P–D species pair. Certainly, in P. montsicciana the (P species), as the allozyme data latterclosely had amatched higher numberthe expectations: of alleles (1)than thethe spectrum former of (49alleles vs. 39) observed and very in P. few unique pardoi (D species) alleles werewas adetected subset of in those the Dfound species in(only P. montsicciana four, while (Pthespecies), P speciesas the hadlatter 14); (2) levels had a higher number of alleles than the former (49 vs. 39) and of within-population genetic variation observed in the D species were lower very few unique allelesthan in the were detected in the[P P species D (percentage species (onlyoffour, while the P polymorphic species loci) = 56.3had vs.14); (2)Alevels 70.3, (mean of number within- of alleles populationper genetic locus)variation = 1.9 vs.observed 2.2; andinHthe D species expected e (unbiased were lower than in the P species heterozygosity) = 0.192[Pvs. 0.239]; (percentage andof polymorphic (3); the two speciesloci) = 56.3 vs. 70.3, showed highA genetic (mean number of alleles similarities, as theperinterspecific locus) = 1.9 values of vs. 2.2; andgenetic He (unbiased identityexpected heterozygosity) (I) were only slightly lower = 0.192 vs. 0.239]; and than intraspecific ones(3);(mean the two I = spe- 0.724 between cies showed high genetic populations ofsimilarities, P. montsicciana as the andinterspecific P. pardoi vs. values meanofI genetic = 0. 809 identity between (I) were populations of only slightly lower thanand P. montsicciana intraspecific ones between mean I = 0.870 (mean I populations = 0.724 between populations of P. of P. pardoi). montsicciana and TheP. results pardoi vs. mean of the I = 0. 809niche ecological between populations analyses performed of P.here montsicciana are in close andagreement mean I = 0.870 withbetween populations the diversity and geneticof P.structure pardoi). patterns; the two Petrocoptis species are showing considerable The results differencesniche of the ecological regarding analysestheirperformed climate niche. here These are in differences close agreement are found both in the G-space with the diversity and genetic(as the niche models structure patterns;builttheusing two the occurrences Petrocoptis species P. montsicciana of are showing are not considerablecapable to predict differences the occurrences regarding their climate P. pardoi, of niche. anddifferences These viceversa; are Figurefound 5B)both and in the E- space (there is no overlap in climatic space between the in the G-space (as the niche models built using the occurrences of P. montsicciana are not 100% of occurrence densities of Diversity 2021, 13, 205 10 of 14 P. montsicciana and P. pardoi; Figure 6). Mostly consequence of its much smaller range compared to P. montsicciana, the climatic conditions in which P. pardoi grows are much narrower with respect to the former regarding the studied variables (Figure 7); nevertheless, signs of climatic divergence are clear, as P. pardoi populations occur either (1) in one of the extremes of the range of variability of P. montsicciana (bio6, bio8, and bio9) or (2) outside this range (bio3, bio7, bio15, and bio19). A straightforward explanation for these signs of climatic divergence is that the two taxa are distributed at different latitudes and therefore there is no overlap in the climatic background space (Figure 6). In other words, the pattern of climate differentiation between the two taxa of Petrocoptis would have been driven, among other factors, by type and abun- dance of the substrate where the species are dwelling. Geologically, both taxa occur roughly on the same kind of rocks (P. pardoi on calcareous conglomerates, P. montsicciana on either calcareous conglomerates or limestone rocks). These kinds of habitats, though frequent in NE Spain [51] are, however, isolated and scatteredly distributed—generally occurring as outcrops or small ranges (Figure S1). Given that Petrocoptis spp. are edaphic specialists instead of generalists, the availability of suitable substrates (rocks with cracks or fissures with enough water supply) where the seeds of Petrocoptis can germinate and develop into mature plants would have, thereby, contributed to delineate their current climate niches. It should be noted that azonal vegetation as the one of rocky outcrops formed by plants like Petrocoptis spp. depend on the microclimate instead of the regional climate. Rocky habitats are presumed to be long-term persistent habitats that have offered opportunities for conserving specific ecological niches [52]. In this sense, it would be highly advisable to characterize the climatic niche of the two Petrocoptis taxa at the microclimatic scale (with several data loggers monitoring the local conditions of the rocky outcrops for 1–2 years, e.g., [53]). Such a design would be able to discern whether the niche differentiation at the macroclimatic scale can be transferred to the microclimatic one. Certainly, the isolated nature of rocky outcrops is playing a major role in population differentiation in the Petrocoptis taxa studied here; populations would rarely exchange genes because are geographically isolated and the species have very low dispersal abilities (seeds are dispersed by gravity or by ants [54], with spider webs behaving as seed receptacles [55]), as stated by many authors in rocky habitats in other European countries [56–58]. Many rocky species show very high genetic divergence among populations (with values of FST for nuclear markers typically around or over 0.5, as found in the present study), including Aquilegia spp. from Sardinia [59], Begonia luzhaiensis from Guangxi Province of SW China [60], Hypochaeris leontodontoides from Morocco [61], or Pitcairnia geyskesii from French Guiana [62]. Petrocoptis montsicciana and P. pardoi are separated by about 150 km and, in addition to this large distance (well enough to interrupt gene exchange), there is also a widely-recognized biogeographical barrier between them, the Ebro River basin (e.g., [63,64]). Speciation at the genus level—which is likely also by geographic isolation [2] as all taxa of Petrocoptis are chasmophytic—seems, however, to be incomplete given the small morphological differences between the alleged species, although other factors might be involved, such as morphological stasis or morphological plasticity. Despite many attempts, the lack of solid diagnostic morphological characters for distinguishing taxa of Petrocoptis (with polymorphism within taxa often exceeding interspecific boundaries [2]) would reflect some degree of morphological stability in a genus that originated long time ago (Late Miocene [4]). It is well known that long-term stability of habitat and climatic conditions (i.e. homogenizing selective pressures) can lead to morphological stasis [65]. As some authors recall (e.g., [66,67]), several Mediterranean disjunct rocky endemics show low levels of morphological variation thanks to their persistence in stable habitats (cliffs and narrow gorges served as glacial refugia in the Mediterranean Basin as they were buffered against the harsh Pleistocene climatic conditions [68]) despite a long history of isolation. Further studies covering all taxa of Petrocoptis and combining high-resolution molecular markers (ideally including those revealing adaptive genetic diversity) will help to shed light on Diversity 2021, 13, 205 11 of 14 the actual reasons for the apparent delay in morphological differentiation with respect to genetic and ecological niche data. 4.2. Conservation Implications Whether P. monticciana/P. pardoi constitute different species cannot be stated with the present results. The reusing of formerly published genetic data (allozyme banding patterns) combined with the generation of new ones (climatic data) allows, however, to suggest that the two entities should be managed independently regarding their conservation, an important issue as both are of conservation concern: P. montisicciana is listed as ‘Near Threatened’ (NT) and P. pardoi as ‘Vulnerable’ (VU) in the most recent version of the Spanish Red List [69]. More importantly, they are protected by law in Aragon (both P. montsicciana and P. pardoi [70]), in Catalonia (P. montisicciana [71]), and in the Valencian Community (P. pardoi [72]). While whether the two taxa should be regarded as ‘Evolutionarily Significant Units’ (ESUs) can be only confirmed by further genetic studies, the significant genetic differences already detected in [10] and confirmed here clearly indicate that P. montsicciana and P. pardoi must be regarded at least as separate ‘Management Units’ (MUs). For example, translocation between MUs could be done only under special circumstances [73], as this can lead to the breakdown of the co-adapted gene complexes and, subsequently, to outbreeding depression [74,75]. A final conservation issue to be addressed is the usefulness of niche modeling to iden- tify new populations of rare or endangered plant species, with many successful examples around the world (e.g., [76–78]), without neglecting the crucial role of in situ and ex situ conservation [79,80] especially for endemic and threatened taxa already known [81–83]). Although the Iberian Peninsula is floristically relatively well known, new field work aimed to discover new populations of Petrocoptis montsicciana/P. pardoi should be directed towards the highest suitability areas (those painted with the darkest colors in Figure 5B) where rocky outcrops are present. Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/d13050205/s1, Figure S1: General aspect of Petrocoptis montsicciana, P. pardoi and their habitats. Figure S2 Contribution of each environmental variable to the spatial distribution of the PCA- env and direction of the seven climatic variables to the first to PCA-env axes. Table S1. Occurrences of Petrocoptis montsicciana and P. pardoi used in the niche analyses. Table S2. Comparison of average (SD) of individual climatic values between Petrocoptis montsicciana (60 occurrences) and P. pardoi (31 occurrences). Author Contributions: Conceptualization, J.L.-P. and N.N.; analysis, S.H.-M., J.L.-P. and N.N.; data supplier, M.G., E.L., J.L.-P., L.S.; maps D.P.-P., writing—original draft preparation, J.L.-P., S.H.-M. and N.N.; writing—review, E.C., M.G., E.L., D.P.-P. and L.S. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by Generalitat de Catalunya, grant number 2017-SGR1116. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Data are available from the authors by direct request. Acknowledgments: We are grateful to Sara López-Vinyallonga for her technical assistance and to José Antonio Fernández Prieto for sharing his profound knowledge on the genus Petrocoptis in the Iberian Peninsula. Conflicts of Interest: The authors declare no conflict of interest. References 1. Buira, A.; Aedo, C.; Medina, L. Spatial patterns of the Iberian and Balearic endemic vascular flora. Biodivers. Conserv. 2017, 26, 479–508. [CrossRef] 2. Mayol, M.; Rosselló, J.A. A synopsis of Silene subgenus Petrocoptis (Caryophyllaceae). Taxon 1999, 48, 471–482. [CrossRef] Diversity 2021, 13, 205 12 of 14 3. Jafari, F.; Zarre, S.; Gholipour, A.; Eggens, F.; Rabeler, R.K.; Oxelman, B. A new taxonomic backbone for the infrageneric classification of the species-rich genus Silene (Caryophyllaceae). Taxon 2020, 69, 337–368. [CrossRef] 4. Cires, E.; Fernández-Prieto, J.A. Phylogenetic relationships of Petrocoptis A. Braun ex Endl. (Caryophyllaceae), a discussed genus from the Iberian Peninsula. J. Plant Res. 2015, 128, 223–238. [CrossRef] [PubMed] 5. De Castro-Arrazola, I.; March-Salas, M.; Lorite, J. Assessment of the potential risk of rock-climbing for cliff plant species and natural protected areas of Spain. Front. Ecol. Evol. 2021, 9, 229. 6. De Bolòs, O.; Rivas-Martínez, S. Comentarios sobre el género Petrocoptis. Petrocoptis montsicciana sp. nova. An. Inst. Bot. A J. Cavanilles 1970, 26, 53–60. 7. Montserrat, P.; Fernández-Casas, J. Petrocoptis. In Flora Iberica; Castroviejo, S., Laínz, M., López-González, G., Montserrat, P., Muñoz-Garmendia, F., Paiva, J., Villar, L., Eds.; Real Jardín Botánico, Consejo Superior de Investigaciones Científicas (CSIC): Madrid, Spain, 1990; Volume 2, pp. 304–312. 8. De Bolòs, O.; Vigo, J. Flora dels Països Catalans. Volum II (Crucíferes-Amarantàcies); Barcino: Barcelona, Spain, 1990. 9. Walters, S.M. Petrocoptis A. Braun. In Flora Europaea, 1, 2nd ed.; Tutin, T.G., Burges, N.A., Chater, A.O., Edmonson, J.R., Heywood, V.H., Moore, D.M., Valentine, D.H., Walters, S.M., Webb, D.A., Eds.; Cambridge University Press: Cambridge, UK, 1993; pp. 190–191. 10. López-Pujol, J.; Bosch, M.; Simon, J.; Blanché, C. Allozyme diversity of two endemic Petrocoptis species: P. montsicciana and its close relative P. pardoi (Caryophyllaceae). Can. J. Bot. 2001, 79, 1379–1389. 11. Dayrat, B. Towards integrative taxonomy. Biol. J. Linn. Soc. 2005, 85, 407–415. [CrossRef] 12. Padial, J.M.; Miralles, A.; de la Riva, I.; Vences, M. The integrative future of taxonomy. Front. Zool. 2010, 7, 16. [CrossRef] 13. Su, X.; Wu, G.; Li, L.; Liu, J. Species delimitation in plants using the Qinghai-Tibet Plateau endemic Orinus (Poaceae: Tridentinae) as an example. Ann. Bot. 2015, 116, 35–48. [CrossRef] 14. Massó, S.; López-Pujol, J.; Vilatersana, R. Reinterpretation of an endangered taxon based on integrative taxonomy: The case of Cynara baetica (Compositae). PLoS ONE 2018, 13, e0207094. [CrossRef] [PubMed] 15. Safner, T.; Miller, M.P.; McRae, B.H.; Fortin, M.-J.; Manel, S. Comparison of Bayesian clustering and edge detection methods for inferring boundaries in landscape genetics. Int. J. Mol. Sci. 2011, 12, 865–889. [CrossRef] 16. Wilkinson, M.D.; Dumontier, M.; Aalbersberg, I.J.; Appleton, G.; Axton, M.; Baak, A.; Blomberg, N.; Boiten, J.-W.; da Silva Santos, L.B.; Bourne, P.E.; et al. The FAIR Guiding Principles for scientific data management and stewardship. Sci. Data 2016, 3, 160018. [CrossRef] [PubMed] 17. Pasquetto, I.V.; Randles, B.M.; Borgman, C.L. On the reuse of scientific data. Data Sci. J. 2017, 16, 8. [CrossRef] 18. Scerri, E.M.L.; Kühnert, D.; Blinkhorn, J.; Groucutt, H.S.; Roberts, P.; Nicoll, K.; Zerboni, A.; Orijemie, E.A.; Barton, H.; Candy, I.; et al. Field-Based sciences must transform in response to COVID-19. Nat. Ecol. Evol. 2020, 4, 1571–1574. [CrossRef] 19. Wagensommer, R.P.; Medagli, P.; Turco, A.; Perrino, E.V. IUCN Red List evaluation of the Orchidaceae endemic to Apulia (Italy) and considerations on the application of the IUCN protocol to rare species. Nat. Conserv. Res. 2020, 5 (Suppl. S1), 90–101. [CrossRef] 20. Kuzmanović, N.; Lakušić, D.; Frajman, B.; Stevanoski, I.; Conti, F.; Schönswetter, P. Long neglected diversity in the Accursed Mountains (western Balkan Peninsula): Ranunculus bertisceus is a genetically and morphologically divergent new species. Bot. J. Linn. Soc. 2021. [CrossRef] 21. Pritchard, J.K.; Stephens, M.; Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 2000, 155, 945–959. [CrossRef] 22. Evanno, G.; Regnaut, S.; Goudet, J. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol. Ecol. 2005, 14, 2611–2620. [CrossRef] 23. Earl, D.A.; von Holdt, B.M. STRUCTURE HARVESTER: A website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 2012, 4, 359–361. [CrossRef] 24. Jakobsson, M.; Rosenberg, N.A. CLUMPP: A cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 2007, 23, 1801–1806. [CrossRef] 25. Rosenberg, N.A. DISTRUCT: A program for the graphical display of population structure. Mol. Ecol. Notes 2004, 4, 137–138. [CrossRef] 26. Manni, F.; Guérard, E.; Heyer, E. Geographic patterns of (genetic, morphologic, linguistic) variation: How barriers can be detected by using Monmonier’s algorithm. Hum. Biol. 2004, 76, 173–190. [CrossRef] [PubMed] 27. Nei, M.; Tajima, F.; Tateno, Y. Accuracy of estimated phylogenetic trees from molecular data. J. Mol. Evol. 1983, 19, 153–170. [CrossRef] [PubMed] 28. Dieringer, D.; Schlötterer, C. Microsatellite analyser (MSA): A platform independent analysis tool for large microsatellie data sets. Mol. Ecol. Notes 2003, 3, 167–169. [CrossRef] 29. Peakall, R.; Smouse, P.E. GENALEX 6: Genetic analysis in Excel. Population genetic software for teaching and research. Mol. Ecol. Notes 2006, 6, 288–295. [CrossRef] 30. Excoffier, L.; Laval, G.; Schneider, S. Arlequin (version 3.0): An integrated software package for population genetics data analysis. Evol. Bioinform. Online 2005, 1, 47–50. [CrossRef] 31. Phillips, S.J.; Anderson, R.P.; Schapire, R.E. Maximum entropy modeling of species geographic distributions. Ecol. Model. 2006, 190, 231–259. [CrossRef] Diversity 2021, 13, 205 13 of 14 32. Broennimann, O.; Fitzpatrick, M.C.; Pearman, P.B.; Petitpierre, B.; Pellissier, L.; Yoccoz, N.G.; Thuiller, W.; Fortin, M.-J.; Randin, C.; Zimmermann, N.E.; et al. Measuring ecological niche overlap from occurrence and spatial environmental data. Glob. Ecol. Biogeogr. 2012, 21, 481–497. [CrossRef] 33. Guardiola, M.; Sáez, L. Petrocoptis montsicciana O. Bolòs & Rivas Mart. In Atlas y Libro Rojo de la Flora Vascular Amenazada de España—Adenda 2017; Moreno Saiz, J.C., Iriondo Alegría, J.M., Martínez García, F., Martínez Rodríguez, J., Salazar Mendías, C., Eds.; Ministerio para la Transición Ecológica-Sociedad Española de Biología de la Conservación de Plantas: Madrid, Spain, 2019; pp. 134–135. 34. Red de Seguimiento para Especies de Flora y Hábitats de Interés Comunitario en Aragón, LIFE12 NAT/ES/000180 RESECOM. Petrocoptis montsicciana O. Bolòs & Rivas Mart. 2014. Available online: http://www.liferesecom.ipe.csic.es/eic/Petrocoptis%20 montsicciana%20_LIFE.doc (accessed on 27 December 2020). 35. López, S.; Fabregat, C. Nuevos datos para la flora de Aragón. Flora Montiber. 2011, 49, 85–95. 36. Brown, J.L. SDM toolbox: A python-based GIS toolkit for landscape genetic, biogeographic and species distribution model analyses. Methods Ecol. Evol. 2014, 5, 694–700. [CrossRef] 37. ESRI. ArcGIS: Release 10.2; Esri Inc.: Redmond, WA, USA, 2013. Available online: http://www.esri.com/ (accessed on 29 December 2020). 38. Liu, C.; Newell, G.; White, M. On the selection of thresholds for predicting species occurrence with presence-only data. Ecol. Evol. 2016, 6, 337–348. [CrossRef] 39. Silva, D.P.; Vilela, B.; Buzatto, B.A.; Moczek, A.P.; Hortal, J. Contextualized niche shifts upon independent invasions by the dung beetle Onthophagus taurus. Biol. Invasions 2016, 18, 3137–3148. [CrossRef] 40. Schoener, T.W. Nonsynchronous spatial overlap of lizards in patchy habitats. Ecology 1970, 51, 408–418. [CrossRef] 41. Warren, D.L.; Glor, R.E.; Turelli, M. Environmental niche equivalency versus conservatism: Quantitative approaches to niche evolution. Evolution 2008, 62, 2868–2883. [CrossRef] [PubMed] 42. Herrando-Moraira, S.; Nualart, N.; Herrando-Moraira, A.; Chung, M.G.; Chung, M.Y.; López-Pujol, J. Climatic niche characteristics of native and invasive Lilium lancifolium. Sci. Rep. 2019, 9, 14334. [CrossRef] 43. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2021. Available online: https://www.R-project.org/ (accessed on 29 December 2020). 44. Gottlieb, L.D. Isozyme evidence and problem solving in plant systematics. In Plant Biosystematics; Grant, W.F., Ed.; Academic Press: Orlando, FL, USA, 1984; pp. 343–357. 45. Gottlieb, L.D. Rethinking classic examples of recent speciation in plants. New Phytol. 2003, 161, 71–82. [CrossRef] 46. Crawford, D.J. Electrophoretic data and plant speciation. Syst. Bot. 1985, 10, 405–416. [CrossRef] 47. Crawford, D.J. Genetics of plant isozymes. In Isozymes in Plant Biology; Soltis, D.E., Soltis, P.S., Eds.; Dioscorides Press: Portland, OR, USA, 1990; pp. 146–164. 48. Crawford, D.J. Progenitor-derivative species pairs and plant speciation. Taxon 2010, 59, 1413–1423. [CrossRef] 49. Gottlieb, L.D. Genetic differentiation, sympatric speciation and the origin of a diploid species of Stephanomeria. Am. J. Bot. 1973, 60, 545–553. [CrossRef] 50. Crawford, D.J. Phylogenetic and systematic inferences from electrophoretic studies. In Isozymes in Plant Genetics and Breeding, Part A; Tanksley, S.D., Orton, T.J., Eds.; Elsevier Science Publishers: Amsterdam, The Netherlands, 1983; pp. 257–287. 51. Gutiérrez, F.; Gutiérrez, M. (Eds.) Landscapes and Landforms of Spain; Springer: Dordrecht, The Netherlands, 2014. 52. Lososová, Z.; Divíšek, J.; Chytrý, M.; Götzenberger, L.; Těšitel, J.; Mucina, L. Macroevolutionary patterns in European vegetation. J. Veg. Sci. 2021, 32, e12942. [CrossRef] 53. Fernández-Mazuecos, M.; Jiménez-Mejías, P.; Rotllan-Puig, X.; Vargas, P. Narrow endemics to Mediterranean islands: Moderate genetic diversity but narrow climatic niche of the ancient, critically endangered Naufraga (Apiaceae). Perspect. Plant. Ecol. Evol. Syst. 2014, 16, 190–202. [CrossRef] 54. Prieto-Mossi, J.; García-Mut, L.; Estrelles, E.; Ibars, A.M. Protocolo de germinación y cultivo de Petrocoptis pardoi Pau (Caryophyl- laceae). Bot. Asppects 2017, 3.2, 7–12. 55. Sainz, H.; Franco, F.; Arias, J. Estrategias para la Conservación de la Flora Amenazada de Aragón; Diputación General de Aragón: Zaragoza, Spain, 1996. 56. Wagensommer, R.P.; Fröhlich, T.; Fröhlich, M. First record of the southeast European species Cerinthe retorta Sibth. & Sm. (Boraginaceae) in Italy and considerations on its distribution and conservation status. Acta Bot. Gall. Bot. Lett. 2014, 161, 111–115. 57. Panitsa, M.; Kontopanou, A. Diversity of chasmophytes in the vascular flora of Greece: Floristic analysis and phytogeographical patterns. Bot. Serb. 2017, 41, 199–211. 58. Liber, Z.; Surina, B.; Nikolić, T.; Škrtić, D.; Šatović, Z. Spatial distribution, niche ecology and conservation genetics of Degenia velebitica (Brassicaceae), a narrow endemic species of the north-western Dinaric Alps. Plant. Syst. Evol. 2020, 306, 64. [CrossRef] 59. Garrido, J.L.; Fenu, G.; Mattana, E.; Bacchetta, G. Spatial genetic structure of Aquilegia taxa endemic to the island of Sardinia. Ann. Bot. 2012, 109, 953–964. [CrossRef] 60. Tseng, Y.-H.; Huang, H.-Y.; Xu, W.-B.; Yang, H.-A.; Peng, C.-I.; Liu, Y.; Chung, K.-F. Phylogeography of Begonia luzhaiensis suggests both natural and anthropogenic causes for the marked population genetic structure. Bot. Stud. 2019, 60, 20. [CrossRef] Diversity 2021, 13, 205 14 of 14 61. Jiménez-López, F.J.; Ortiz, M.A.; Berjano, R.; Talavera, S.; Terrab, A. High population genetic substructure in Hypochaeris leontodontoides (Asteraceae), an endemic rupicolous species of the Atlas Mountains in NW Africa. Alp. Bot. 2016, 126, 73–85. [CrossRef] 62. Boisselier-Dubayle, M.-C.; Leblois, R.; Samadi, S.; Lambourdière, J.; Sarthou, C. Genetic structure of the xerophilous bromeliad Pitcairnia geyskesii on inselbergs in French Guiana—A test of the forest refuge hypothesis. Ecography 2010, 33, 175–184. [CrossRef] 63. Orellana, M.R.; López-Pujol, J.; Blanché, C.; Bosch, M. Genetic diversity in the endangered dysploid larkspur Delphinium bolosii and its close diploid relatives in the series Fissa of the Western Mediterranean area. Biol. J. Linn. Soc. 2007, 92, 773–784. [CrossRef] 64. Vergara, M.; Basto, M.P.; Madeira, M.J.; Gómez-Moliner, B.J.; Santos-Reis, M.; Fernandes, C.; Ruiz-González, A. Inferring population genetic structure in widely and continuously distributed carnivores: The stone marten (Martes foina) as a case study. PLoS ONE 2015, 10, e0134257. [CrossRef] 65. Milne, R.I.; Abbott, R.J. The origin and evolution of Tertiary relict floras. Adv. Bot. Res. 2002, 38, 281–314. 66. Molins, A.; Bacchetta, G.; Rosato, M.; Rosselló, J.A.; Mayol, M. Molecular phylogeography of Thymus herba-barona (Lamiaceae): Insight into the evolutionary history of the flora of the western Mediterranean islands. Taxon 2011, 60, 1295–1305. [CrossRef] 67. Bobo-Pinilla, J.; Barrios de León, S.B.; Seguí Colomar, J.; Fenu, G.; Bacchetta, G.; Peñas de Giles, J.; Martínez-Ortega, M.M. Phylogeography of Arenaria balearica L. (Caryophyllaceae): Evolutionary history of a disjunct endemic from the Western Mediterranean continental islands. PeerJ 2016, 4, e2618. [CrossRef] 68. Médail, F.; Diadema, K. Glacial refugia influence plant diversity patterns in the Mediterranean Basin. J. Biogeogr. 2009, 36, 1333–1345. [CrossRef] 69. Moreno, J.C. (Ed.) Lista Roja de la Flora Vascular Española. Actualización con los Datos de la Adenda 2010 al Atlas y Libro Rojo de la Flora Vascular Amenazada; Sociedad Española de Biología de la Conservación de las Plantas: Madrid, Spain, 2011. 70. Diputación General de Aragón. Decreto 181/2005, de 6 de septiembre, del Gobierno de Aragón, por el que se modifica parcialmente el Decreto 49/1995, de 28 de marzo, de la Diputación General de Aragón, por el que se regula el Catálogo de Especies Amenazadas de Aragón. Bol. Of. Aragón 2005, 114, 11527–11532. 71. Generalitat de Catalunya. RESOLUCIÓ AAM/732/2015, de 9 d’abril, per la qual s’aprova la catalogació, descatalogació i canvi de categoria d’espècies i subespècies del Catàleg de flora amenaçada de Catalunya. D. Of. General. Catalunya 2015, 6854. CVE-DOGC-A-15106031-2015. 72. Generalitat Valenciana. Ordre 6/2013, de 25 de març, de la Conselleria d’Infraestructures, Territori i Medi Ambient, per la qual es modifiquen les llistes valencianes d’espècies protegides de flora i fauna. D Of. Comunitat Valencia. 2013, 6996, 8682–8690. 73. Moritz, C. Conservation units and translocations: Strategies for conserving evolutionary processes. Hereditas 1999, 130, 217–228. [CrossRef] 74. Fenster, C.B.; Dudash, M.R. Genetic considerations for plant population restoration and conservation. In Restoration of Endangered Species: Conceptual Issues, Planning and Implementation; Bowles, M.L., Whelan, C.J., Eds.; Cambridge University Press: Cambridge, UK, 1994; pp. 34–62. 75. Frankham, R.; Ballou, J.D.; Eldridge, M.D.B.; Lacy, R.C.; Ralls, K.; Dudash, M.R.; Fenster, C.B. Predicting the probability of outbreeding depression. Conserv. Biol. 2011, 25, 465–475. [CrossRef] [PubMed] 76. Guisan, A.; Broennimann, O.; Engler, R.; Vust, M.; Yoccoz, N.G.; Lehmann, A.; Zimmermann, N.E. Using niche-based models to improve the sampling of rare species. Conserv. Biol. 2006, 20, 501–511. [CrossRef] [PubMed] 77. Menon, S.; Choudhury, B.I.; Khan, M.L.; Peterson, A.T. Ecological niche modeling and local knowledge predict new populations of Gymnocladus assamicus a critically endangered tree species. Endang. Species Res. 2010, 11, 175–181. [CrossRef] 78. Särkinen, T.; Gonzáles, P.; Knapp, S. Distribution models and species discovery: The story of a new Solanum species from the Peruvian Andes. PhytoKeys 2013, 31, 1–20. [CrossRef] 79. Wagensommer, R.P.; Perrino, E.V.; Silletti, G.N. Carex phyllostachys A. Mey. (Cyperaceae) new for Italy and phytogeographical considerations. Phyton 2014, 54, 215–222. 80. Schwartz, K.R.; Parsons, E.C.M.; Rockwood, L.; Wood, T.C. Integrating in-situ and ex-situ data management processes for biodiversity conservation. Ecol. Evol. 2017, 120, 1–7. [CrossRef] 81. Işik, K. Rare and endemic species: Why are they prone to extinction? J. Bot. 2011, 35, 411–417. 82. Perrino, E.V.; Silletti, G.N.; Erben, M.; Wagensommer, R.P. Viola cassinensis subsp. Lucana (Violaceae), a new subspecies from Lucanian Apennine, southern Italy. Phyton 2018, 58, 109–115. 83. Coelho, N.; Gonçalves, S.; Romano, A. Endemic plant species conservation: Biotechnological approaches. Plants 2020, 9, 345. [CrossRef]