Eokochia saxicola is a multi-branched perennial evergreen shrub found exclusively on maritime rocks (both calcareous and volcanic) close to the sea level in the area occupied by halophilic vegetation (Strumia et al., 2020b). Characteristic features of this species include protogyny, papillose stigma, and recalcitrant pollen, typical of wind pollinated plants (Barone Lumaga et al., 2016). Like other halophytes, the main seed dispersal mechanism is hypothesized to be hydrochory, owing mainly to the presence of diaspores able to float for several days. Moreover, a decrease in the germination rate of aging seeds and the inability for seedlings to anchor at suitable establishment sites have been assessed (Strumia et al., 2020a). According to IUCN, E. saxicola is threatened with extinction as it falls in the category Endangered (EN) (Santangelo et al., 2012; Orsenigo et al., 2018). The main threats are landslides in the few sites of occurrence and supposed reduction of fertility and reproduction success due to the small number of living individuals. We sampled E. saxicola in three locations representing the only known populations of its current distribution: Capri Island near Naples (Campania), Strombolicchio islet (Aeolian Archipelago, Sicily), and Palinuro, along the Cilento and Vallo di Diano and Alburni Nature Park coastline (Campania). Field collections were conducted from 2015 to 2017. Due to the inaccessibility of the cliffs (Strumia et al., 2020a), only the most accessible individuals were sampled. In each sampling site, the maximum estimate number of individuals was recorded through visual census (i.e., by binoculars from a boat) (Table 1).

Genomic DNAs of 51 Eokochia saxicola individuals were extracted from preserved leaf material using the DNeasy Plant Mini Kit (Qiagen, Hilden, Germany), following the protocol of the manufacturer. The quantity and the purity of DNA extractions were checked using the NanoDrop ND 1000 Spectrophotometer (Thermo Fisher Scientific, Delaware), whereas the concentration was assessed with the Qubit fluorometer system (Invitrogen, United States) and the Quant-IT ds-DNA BR Assay kit (Invitrogen). Samples were sequenced using a slight modification to the standard ddRAD protocol (Peterson et al., 2012), with the restriction enzymes EcoRI and TaqI. A single library with barcoded individuals was sequenced in a single Illumina HiSeq 2500 lane for 150-base single-end reads. From a subset of samples (16 E. saxicola + S. hirsuta), a “genome skim” was performed (Straub et al., 2012) to gather sufficient data for complete plastid assembly without prior enrichment or isolation of plastid DNA (Coissac et al., 2016; Supplementary Table SM1). Libraries were constructed according to the Nextera DNA Library Prep Kit (Illumina) protocol and Nextera XT indexes (Kit v2 Set A), which were used to multiplex the individual samples. Libraries were sequenced in a single Illumina HiSeq 2500 lane for 150-base single-end reads.

The Eokochia saxicola and Spirobassia hirsuta raw sequence reads were quality checked with the FastQC v0.11.9¹, and low-quality reads and adapters were subsequently trimmed with Trimmomatic v0.39 (Bolger et al., 2014). We implemented a de novo assembly of S. hirsuta plastome by setting kmer = 121 in Velvet v1.2.10 software (Zerbino and Birney, 2008), and trimmed reads were then further assembled using Pilon v1.23 (Walker et al., 2014). We implemented FASTQ Screen v0.14 (Wingett and Andrews, 2018) for filtering the plastid reads from each E. saxicola individual by using a de novo-assembled S. hirsuta plastome (unpublished), and orphan reads were removed through the python script FASTQ Combine Paired End.py². The E. saxicola filtered reads were mapped to S. hirsuta plastome by using BOWTIE2 v2.3.5.1 (Li and Durbin, 2009). The alignments (BAM files) produced a single-variant call format (vcf) file for all samples, where the SNPs were identified using bcftools v.1.9 with the setting “mpileup -Ou” and called via bcftools using the −mv and –ploidy 1 functions. The vcf file was converted to a Nexus format, and a haplotype network (also including available Amarantaceae plastome sequences) was built by using the TCS method implemented in PopArt v1.7 (Leigh and Bryant, 2015).

Raw reads were demultiplexed in individual paired-end libraries and filtered for low-quality reads with the process_radtags pipeline in STACKS v2.4 (Rochette et al., 2019). Quality control of demultiplexed reads was performed with FastQC v0.11.9. The catalog RAD-loci was built with de novo assembly using the de novo_map.pl wrapper in STACKS v2.4 by using parameters m = 1, M = 8, and n = 1, following the optimization procedure described in Paris et al. (2017). The population program in a Denovo_map.pl wrapper was used for collecting the polymorphic SNPs in a vcf file. The vcf files (with 48 individuals) were analyzed with the –missing function of PLINK v1.07 (Purcell et al., 2007) for calculating the percentage of missing data (approximately ∼81% of missing data). Lastly, to avoid the risk of calling false SNPs, we reduced the dataset to include all sites with a minor allele frequency (MAF) greater than or equal to 0.05 and with a (–max-obs-het) setting of 0.70 (as also suggested in Gargiulo et al., 2021). Finally, to account for the effects of missing data per individual, VCFTOOLS v0.1.14 (Danecek et al., 2011) was implemented to filter the minimum number of ddRAD loci shared by at least 70% of individuals (i.e., the full dataset). To minimize the potential effect of missing data in altering the results, some population demographic analyses were also or only performed on a subsample of individuals by maximizing the number of reads per sample and minimizing the amount of missing data. Therefore, a reduced dataset was constructed with 25 individuals (i.e., only those ranging between 2,000,000 and 500,000 reads, namely 10 individuals from Palinuro, eight from Capri, seven from Strombolicchio, see Supplementary Text 1.1), which were filtered so that missing data were present in two or less individuals for each locus (i.e., ddRAD loci present in 95% of the individuals).

Analysis of full plastid genome variation in E. saxicola revealed a surprisingly low level of intraspecific variation (only three haplotypes differing for a maximum of two mutation steps) among the three isolated populations with a considerable number of differences (in total plastome length too) with the related S. hirsuta. Since large interspecific differences can be expected between paleoendemic species separated approximately 10 Mya (cfr. Kadereit and Freitag, 2011), both haplotype admixture in populations and very low haplotype variation at the intraspecific level support a recent separation of current E. saxicola populations. The Palinuro and Capri populations share two haplotypes (H1 and H2), differing by a single-base mutation (Figure 1). The most common haplotype (H2) was found in the Strombolicchio population, but, here, the “distant” haplotype H3 (i.e., differing by two mutational steps) was also detected (Figure 1). Even if the scenario of dispersal from Palinuro is the best supported, actual co-occurrence of two related haplotypes in both insular Capri and mainland Palinuro does not rule out an alternative scenario where these populations were likely part of a larger, ancestral population that later split into today isolated populations without a bottleneck phase. Interestingly, this scenario is the second-best supported scenario (i.e., Scenario 3) in the DYABC analysis based on nuclear markers (discussed below).

According to the morpho-bathymetric reconstruction of the Southern Tyrrhenian Sea (Aiello et al., 2015) and considering that, during the last glacial phases, the sea level decreased up to 100 meters, Capri was once connected to the mainland cliffs of Southern Italy. Moreover, the coastline shape between Capri and Palinuro drastically changed during the Quaternary due to the combination of tectonics, sedimentary inputs, volcanism, and sea-level oscillations (Santangelo et al., 2017). In this scenario, which occurred in the last 1.8 Mya, both frequency and surface extension of unsuitable (defined as sandy coasts related to alluvial plains) and suitable (defined as rocky coasts related to pre-quaternary geological substrates) habitats for E. saxicola were transformed from the sea-level oscillations. Moreover, these oscillations most likely reduced the “jumping” distance from the Palinuro cliffs to Capri, favoring a linear “hopping” migration along the coastline, as also reported for several other Mediterranean halophytes (Clausing et al., 2000; Kadereit et al., 2005; Weising and Freitag, 2007). The current distribution of E. saxicola along the Palinuro cliffs (ca. < 15 km) further supports the evidence of a “hopping” dispersal ability along with these coastal rocky habitats.

In contrast, the relatively young volcanic islet of Strombolicchio (0.2 Mya) has been completed isolated from any mainland, even during the Quaternary seal-level oscillations. Accordingly, long distance dispersal or “jumping” is the only plausible scenario for the colonization of E. saxicola on Strombolicchio. Surprisingly, despite the rarity of E. saxicola populations and the small size of the islet, at least two occasional events (two different plastid haplotypes) of colonization likely occurred as founder events for the Strombolicchio population (Figure 1). Hydrochory, as for many other halophytes of the Mediterranean Sea (Kadereit et al., 2005; Arafeh and Kadereit, 2006; Weising and Freitag, 2007; De Castro et al., 2020), may have contributed to the long-distance dispersal of E. saxicola seeds. Indeed, sea-drifted seeds depend entirely on currents for their dispersal, and there is some evidence from the surface circulation observed in the southern Tyrrhenian Sea (Vetrano et al., 2004; Rinaldi et al., 2010) that is consistent with the known distribution of E. saxicola. Changes in sea currents due to the presence of a vortex north of Capri renders the north and south Tyrrhenian Sea “isolated” from each other (Procaccini et al., 2001) and impedes seeds dispersal toward the north.

According to the average reported speed of south Tyrrhenian surface circulation (Rinaldi et al., 2010), the floating time ability of the E. saxicola (Barone Lumaga et al., 2016) would have been sufficient to reach Strombolicchio islet from the Palinuro/Capri coastlines. Haplotype admixture between Capri and Palinuro does not allow detecting the source population, but nuclear data point to the Capri population as more proximate to the Strombolicchio one. Nevertheless, the presence of a second haplotype in the young islet, assuming that this haplotype may not have locally evolved or that it was present in Palinuro/Capri and gone undetected/extinct, points to a second colonization event from an unknown source population. This latter hypothesis is supported by the distinctiveness of exclusive Strombolicchio haplotype compared with those found in Palinuro/Capri (i.e., two mutation steps). The presumed source population, likely distributed along the south Tyrrhenian coastline, could be extinct today or even still existing. Despite the long history of floristic investigations in the Mediterranean Basin, the existence of unknown populations of E. saxicola could not be excluded due to the difficulty in recognizing the species (Strumia et al., 2015). Implementation of botanical surveys along the South Tyrrhenian coastlines, considering the Sea surface circulation and suitable rocky habitats, could fill this gap. Furthermore, a previous study on the genus Limonium showed that rising sea levels during the warm phases of the glaciation had led to the steady decline of coastal populations (Koutroumpa et al., 2021). We, therefore, cannot exclude that this flooding may have played a similar role in the extinction of some E. saxicola populations along the coastline.

The complete analysis of plastid genome variation revealed partial admixture and a very low level of between population differentiation (only 1–2 mutation steps). The different analyses of nuclear markers, performed with different settings, concordantly shows three main genetic clusters matching the three sampling localities of Palinuro, Capri, and Strombolicchio (Figure 2C). In contrast to plastidial admixture, the nuclear markers revealed a present story of strong genetic isolation for the three living populations (Figures 2–4; Table 3). Indeed, all sampled individuals have been correctly assigned to their source population in all analyses, thus almost indicating a complete absence of genetic admixture and ongoing gene flow (with corresponding high Fst values).

MIGRATE and BayesAss estimates revealed that historical (mh) and contemporary (mc) levels of migration between the three genetically distinct populations are very small and similar in magnitude (Supplementary Table 2), as also supported by the large generation times estimated by the simulation performed in DIYABC. The low levels of contemporary migration are not surprising because geographic distance and cliff habitat make dispersal between populations highly unlikely. However, as historical migration rates are also very low, concordant patterns strongly imply that the high levels of genetic structure currently observed among extant populations stem from a low colonization ability, which thus seems to be a long-standing life-history trait of E. saxicola. One evolutionary implication of similar low contemporary and historical migration rates over time is that E. saxicola has a long history of living in relatively small, isolated populations (still preserving enough genetic variation). In contrast to other island narrow endemics (Molins et al., 2009; Mayol et al., 2012; López-Pujol et al., 2013), for which the sea level oscillations represented a switch-off-switch-on the barrier to gene flow, sea levels theoretically do not represent a major barrier for seed dispersal of E. saxicola, as well as for other halophytes of both sandy and rocky coasts. Nonetheless, phases of sea lowering could have led to both the increase in available habitat through land emersion and the promotion of hopping colonization by the decrease in distance between isolated populations. However, phases of the current sea rising do not totally halt the potential dispersion of E. saxicola but decrease the probability of its occurrence due to the increased distance between a few suitable sites.

The DIYABC analysis identified two introductions temporally separated from the Palinuro to Capri, and from there to Strombolicchio (DIYABC model selection, Scenario 1) as the most likely hypothesis for E. saxicola presents distribution. By assuming a generation time of 10 years, we can consider that the introduction to Capri, approximately 0.25 Mya, and the colonization of Strombolicchio (0.09 Mya) took place around the late Pleistocene. The presence of inaccessible known plant spots (on Capri) and the possible occurrence of other unknown plant spots along the Palinuro cliffs suggest we probably have underestimated the census population sizes (N) for the three populations. Still, Ne estimations are near and even larger than N (Table 3), i.e., most of the individuals are breeding individuals in the population. Overall, the effective population size of E. saxicola is historically low because climatic oscillations during the quaternary glacial-interglacial cycles periodically forced populations to track their optimum into smaller cliff areas of suitable habitat. The almost exclusive preference of E. saxicola for north-facing cliffs suggests an adaptation to cooler conditions, which further highlights the possibility of different phases of population expansion/contraction during these oscillation cycles.

While geographic isolation seems the obvious explanation for the current population genetic structure, it partly contrasts with the levels of genetic diversity still detected within each population. Such genetic diversity within populations is unlikely to have resulted from recent or contemporary gene flow but rather points at relative demographic stability (as suggested by Ne estimations) both in recent and ancient timescales. All sampled individuals, despite a clumped distribution due to the accessibility of the cliff faces, were genetically different, indicating that clonality in this species is rare. Moreover, at least at the intrapopulation level, the long lifespan of adult plants, together with their high resprouting ability and wind pollination, may contribute to maintaining a high level of outcrossing and multiple random mating among the few individuals (i.e., observed heterozygosity was comparable to the expected one). Interestingly, wind pollination is an almost unique trait among Mediterranean narrow endemics. Indeed, anemophily is a plesiomorphic trait in E. saxicola and Amarantaceae, whereas, in the few other known wind-pollinated endemics, this mechanism is rather a secondary adaptation to the low availability of biotic pollination (Traveset and Navarro, 2017).

Regardless of the extremely small number of individuals in each population, the moderate degree of intrapopulation co-ancestry (Figure 5) suggests that a consistent level of crosses among unrelated individuals occurs. Even if the reproductive success of individual plants is unknown, the pattern of genetic relatedness among individuals of these small populations suggests that most of the adult plants contribute to the next generation, i.e., that random mating occurs. This reproductive behavior is consistent with the observation of Ho not different from He and a moderate excess of homozygotes in all populations compared with Hardy–Weinberg equilibrium expectations (Fis average, 0.10). The relative excess of homozygotes (higher Fis) observed in Palinuro denotes some heterogeneity in this mainland population compared to the islets, consisting of populations composed of several subgroups (the so-called Wahlund effect). Indeed, in contrast to the other two insular populations sampled on single cliffs, the Palinuro population was sampled along different cliffs. This subpopulation structure is also partially evident in the heatmap graph of Fst (Figure 4).

As reported for other extremely narrow endemics (see Médail and Baumel, 2018), E. saxicola still maintains some degree of genetic diversity even at the intrapopulation level (despite its small population sizes); this genetic trait may have represented life insurance for the long-term survivorship of this rare species. This survivorship was further aided by the low level of anthropogenic disturbance and high stability typical of cliff plant communities (Davis, 1951). Indeed, maintaining genetic diversity can reduce the probability of extinction of small populations by providing the standing genetic variation for local adaptation (Leimu and Fischer, 2008). Therefore, the present Palinuro population still holds enough genetic diversity to represent a valuable “donor population” for future reintroduction programs of this species. The high heterozygosity can be relevant in terms of conservation guidelines as genetic diversity of E. saxicola can be preserved in the long term even in its current isolated populations, confirming that narrow endemics are not necessarily “evolutionary dead-ends,” but rather may represent species that have a strong evolutionary legacy (Médail and Baumel, 2018). From a conservation perspective, we suggest that, at least in the short term, genetic factors may have little impact on the persistence of small populations of this rare paleoendemic species, but, rather, ecological factors (as persistence of suitable habitat) may play a larger role in determining whether populations survive in the long term. In this respect, both the IUCN category and the potential threats of extinction of this species (Santangelo et al., 2012; Orsenigo et al., 2018) should be reconsidered.