Reusing Old and Producing New Data Is Useful for Species Delimitation in the Taxonomically Controversial Iberian Endemic Pair Petrocoptis montsicciana / P. pardoi (Caryophyllaceae)

: Petrocoptis montsicciana and P. pardoi are two Iberian endemic taxa of Caryophyllaceae family with an unclear taxonomic delimitation, being variously treated as independent species, subspecies or even synonyms. In the present study, allozyme raw data obtained in the early 2000s have been reused with improved tools to survey genetic structure, and complemented with modeling and niche comparative analyses to shed light on species delimitation. Genetic structure was investigated using four approaches: Bayesian clustering, Monmonier’s algorithm, Principal Coordinate Analysis (PCoA), and Analysis of Molecular Variance (AMOVA). Ecological niche differences have been assessed through Ecological Niche Modeling (ENM) using MaxEnt, and Principal Component Analysis using 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 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 abundance of rocks where species live although differences at the microclimate instead of the regional climate should be explored in future research. Given the genetic distinction between P. montsicciana and P. pardoi , both taxa should be regarded as separate ‘Management Units’ (MUs).


Introduction
Petrocoptis A. Braun ex Endl. (Caryophyllaceae) is one of the flagship genera of higher plants of the Iberian Peninsula, given to a series of reasons that include: (1) its alleged endemic status; (2) the conservation concerns of several of its members (some are narrow endemics and highly threatened); and (3) the beauty of the plants and the habitats where they grow ( Figure S1). Petrocoptis is one of the 27 endemic genera in the Iberian flora, with nine species [1], although it should be included within Silene L. according to some authors [2]. However, recent molecular studies (based both on nuclear and plastid sequences) have confirmed its independence from Silene [3], having probably diverged approximately in Late Miocene, ca. 10.6 Mya [4]. Most of Petrocoptis taxa have restricted 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. pardoi 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 policies 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].

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 10 5 iterations and 10 6 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]. Programs Clumpp v.1.1.2 [24] and Distruct v.1.1 [25] were used for the post-processing of the results obtained from Structure. Second, the location of potential genetic barriers between populations was explored through the Monmonier's algorithm implemented in Barrier v.2.2 [26]. The significance of barriers was tested for nSSR data by bootstrapping 1000 Nei's DA genetic distances [27] matrices that were previously obtained with Microsatellite Analyzer (MSA) v.4.05 [28]. Third, GenAlEx v.6.1 [29] was used to perform a Principal Coordinate Analysis (PCoA) based on codominant genotypic distances both at the individual and population levels. Fourth, and last, an Analysis of Molecular Variance (AMOVA) using the program ARLEQUIN v.3.11 [30] was performed. AMOVA was carried under two hypotheses: (1) all populations belong to a single species; and (2) taking into account the existence of two different species and therefore establishing three hierarchical levels of analysis: among taxonomic groups, among populations within species and within populations.

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 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/); 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) under current climatic conditions . 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 ArcGIS [36] was conducted to retain only seven relatively uncorrelated climatic variables (r < |0.85|): isothermality (bio3), minimum temperature of coldest month (bio6), temperature 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 exported 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.

Genetic Structure and Divergence Between Taxa
Based on the ΔK approach, K = 3 was the most likely number of genetic clusters after running the software Structure, although K = 2 also provided good results (Figure 2A). In the K = 2 grouping scheme, all populations of P. montsicciana (PM-CAM, PM-TER, PM-MRB and PM-MCO; for population codes, see Figure 1) were merged in a single cluster, whereas those of P. pardoi (PP-CBD, PP-BAR and PP-AGV) formed the second cluster (Figure 2B). With K = 3, the cluster of P. montsicciana was split into two: one containing the southernmost populations of the species (PM-CAM, PM-TER) and the other with the northernmost ones (PM-MRB and PM-MCO), whereas all the populations of P. pardoi remained in the same genetic cluster ( Figure 2B). Barrier results also indicated the lack of genetic differences between the populations of P. pardoi, while the largest barriers were put between the two taxa (the barrier between PM-CAM and PP-AGV received up to 92% 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 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 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 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 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 populations by species, the percentage of variation between P. montsicciana and P. pardoi was 18.60% (Table 1).

Climatic Niches of P. montisicciana and P. pardoi
Niche models produced with MaxEnt for both species had very high AUC (> 0.99), indicating good performance (Table 2). According to jackknife test and percentage contribution to models, bio3, bio7 and bio9 were important for the two Petrocoptis taxa, although for P. pardoi bio19 ranked first according to jackknife test ( Table 2). The predicted areas of suitable habitat for both taxa ( Figure 5B) included their respective occurrence points, although new suitable areas appeared: with few exceptions, these were located in the same latitude that the occurrence records for P. montsicciana, while areas further south were evidenced for P. pardoi ( Figure 5). Therefore, neither of the distribution models was able to predict the occurrences of the other species.  The first two components of the PCA-env explained 76.54% of the total climatic variation ( Figure 6). For both axes (PC1 = 55.47% and PC2 = 21.07%) the precipitation variables bio15 and bio19 were the most informative factors ( Figure S2). 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 Schoener's D index (D = 0, no overlap) and niche metrics (niches were non-stabilized, totally unfilled, and in expansion) ( Table 3). When analyzing the bioclimatic variables individually, all showed significant differences (p-value < 0.05) in Wilcoxon-Mann-Whitney test, although test values for bio 8 and bio9 are near to p-value (Table S2 and Figure 7).

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 approximations (Structure, Barrier, PCoA, and AMOVA) confirms the conclusions presented in the previous paper: P. montsicciana and P. pardoi are assigned to two genetically differentiated clusters, although this differentiation is not large. In [10], the authors noted that their genetic results were compatible with a scenario of a progenitor-derivative (P-D) species pair, based on both the pattern of genetic structure and allele profiles. In contrast to the classic allopatric mode of speciation, in a P-D species pair the derived (D) species buds off and acquires new traits-or becomes fixed for features within the polymorphism of the parental (P) species-while the P species remains largely unchanged [44][45][46][47][48]. Although the D species may show a series of changes with respect to P that may include new morphological characters, a new breeding system, selection for a new habitat, and chromosomal restructuring [45], these could not occur [48]. In such cases, genetic markers can be useful tools for detecting P-D species pairs [48][49][50].
In the study of [10], ecological and morphological differences between P. montsicciana and P. pardoi were not explored, as they were assumed to be minimal, if any. Because they have also the same chromosome number, therefore, the allozyme patterns were the only evidence for these authors to support the assignment as a P-D species pair. Certainly, the allozyme data closely matched the expectations: (1) the spectrum of alleles observed in P. pardoi (D species) was a subset of those found in P. montsicciana (P species), as the latter had a higher number of alleles than the former (49 vs. 39) and very few unique alleles were detected in the D species (only four, while the P species had 14); (2) levels of withinpopulation genetic variation observed in the D species were lower than in the P species [P (percentage of polymorphic loci) = 56. The results of the ecological niche analyses performed here are in close agreement with the diversity and genetic structure patterns; the two Petrocoptis species are showing considerable differences regarding their climate niche. These differences are found both in the G-space (as the niche models built using the occurrences of P. montsicciana are not capable to predict the occurrences of P. pardoi, and viceversa; Figure 5B) and in the E-space (there is no overlap in climatic space between the 100% of occurrence densities of 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 abundance 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][57][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 the actual reasons for the apparent delay in morphological differentiation with respect to genetic and ecological niche data.

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 identify new populations of rare or endangered plant species, with many successful examples around the world (e.g. [76][77][78]), without neglecting the crucial role of in situ and ex situ conservation [79,80] especially for endemic and threatened taxa already known [81][82][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 www.mdpi.com/1424-2818/13/5/205/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