The Natural Oriented Reserve “Saline di Trapani e Paceco” is located along the western coast of Sicily and covers approximately 1,000 ha. The reserve is characterized by regulated salinas that comprise a number of interconnected ponds, built for the production of salt. Nonetheless, the area is recognized as a Special Protection Area (SPA) as part of the NATURA 2000 network (SIC ITA010007). Different pond types can be distinguished by their positioning and the salinity of the water. At the beginning of the salt process, the first pond receives water directly from the sea or via a channel connected to the sea. This pond, referred to as ‘salina fridda’, has a salinity that is similar or slightly higher than that of seawater. In later phases, water progressively flows to the next ponds, from where it will start to evaporate and salinity levels will increase. Ruppia spiralis could be collected from the first two kinds of ponds in the salt process (‘salina fridda’ and ‘salina caure d’acqua frisca’) and the accompanying channels.

Along the coast, sea currents follow a seasonal pattern. Throughout the summer, the reserve is subject to the Tyrrhenian Sicilian Current that directs relatively high saline water from north to south (Thompson et al., 2016; Menna et al., 2019). Because of the irregular shape of the coastline, different areas experience this current at different velocities. As such, the zone between the mainland and Isola Grande becomes almost completely sheltered, forming coastal lagoons at shallow depths ( Figure 1 ). Here, additional populations of R. spiralis were sampled.

Ruppia spiralis plants were collected in 10 water bodies (identified by a code from A to J) during the first week of August 2020 ( Table 1 ). At each locality, 16 to 30 individual shoots were collected following a 30-m linear transect. For sites with a low density of stands, the linear transect was extended until 30 individual shoots were collected. For population C, we sampled without transect as stands were so small and scattered. Each plant sample was individually labeled and dried, and for each population, a reference herbarium was deposited at BRVU (Herbarium of Vrije Universiteit Brussel, formally transferred in 2018 to BR, Botanic Garden, Meise, Belgium).

Genomic DNA was isolated from dry plant material using the E.Z.N.A. SP Plant DNA Mini Kit (Omega Bio‐tek, Norcross, GA, USA). Polymerase chain reaction (PCR) was used to amplify 15 microsatellite loci ( Supplementary Table S1 ). For the total reaction mixture, 2.5 μl genomic DNA was combined with 1.25 μl primer mix (forward and reversed fluorescent labeled primers), 6.25 μl QIAGEN Multiplex PCR Master Mix, and 2.5 μl H₂O. The PCR reactions were performed in a thermal cycler (Bio-Rad T100) as follows: 95°C for 15 min, 35 cycles of 95°C for 30 s (denaturation), 57°C for 1 min 30 s (annealing), 72°C for 40 s (elongation), and a final extension step of 60°C for 30 min. PCR products were run on an ABI3730XL sequencer (Macrogen, Seoul, Korea), and fragments were analyzed with GeneMarker V2.6.0 (SoftGenetics LLC^®). We encountered 10 ramets with five alleles for one locus during scoring. To allow data analyses as a tetraploid taxon, the least common allele was omitted from those few cases.

To detect haplotype variation and rule out any different species identity, one sample of each population was used for a comparative analysis of chloroplast sequences. After passing quality inspection (DNA concentration between 5 and 15 ng/µl), the constructed library (TruSeq Nano DNA Kit) was sequenced by 300-bp × 2 paired-end sequencing in an Illumina MiSeq platform (Macrogen, Seoul, South Korea). Raw data were used for de novo assembling of several genes/introns from a nearly complete chloroplast genome. All assemblies used GenBank accession MN233650 of Ruppia sinensis Yu and den Hartog (Yu et al., 2019) as a reference genome (similar method as in Triest et al., 2021b), and the mapping of reads of Sicilian samples was performed in Geneious Prime v 2019.2.1 (^©2005–2019 Biomatters Ltd.) (https://www.geneious.com). The assemblages averaged 12,913–57,881 reads with a mean depth of reads ranging from 24 to 109 coverage. The de novo targeted chloroplast assemblies were processed against R. spiralis sequences (under R. cirrhosa in GenBank) of the intergenic spacer trnH-psbA (JN113267; 452 bp) and three mononucleotide microsatellite regions, ccmp2 (JN113250; 202 bp), ccmp3 (JN113258; 157 bp), and ccmp10 (JN113260; 166 bp), representing haplotype B from the Mediterranean as defined by Triest and Sierens (2014). The consensus sequences were aligned with ClustalW (Thompson et al., 1994). Lastly, chloroplast microsatellite regions that showed variation among populations were amplified for individuals, to check if haplotypes occur homogeneously within a population.

The nuclear ribosomal cistron (18S, ITS1, 5.8S, ITS2, and 26S) was assembled from a Ruppia spiralis sequence of a partial nuclear ribosomal cistron (GenBank JQ937105; submitted under Ruppia maritima) containing an internal transcribed spacer 1 (partial sequence), 5.8S ribosomal RNA gene (complete sequence), and internal transcribed spacer 2 (partial sequence) and was subsequently used as a seed in Geneious Prime v 2019.2.1 to progressively enlarge the flanking regions. A base pair of 6,545 was obtained for one selected sample and further used as a reference to map every other sample, averaging 11,570–20,436 reads, with a mean depth of reads ranging from 380 to 652 coverage. The generated consensus sequences of 7,338-bp length of 10 samples (same as used for chloroplast sequences) were aligned for a comparative description of mutated positions using MAFFT v7.388 (Katoh et al., 2002; Katoh and Standley, 2013).

Prior to all analyses, the package Poppr in R v1.2.5001 (Kamvar et al., 2014; Kamvar et al., 2015; R Core Team, 2019) was used to detect repeated multilocus genotypes (MLGs) within populations. Because repeated MLGs are likely the result of clonal growth, only one representative sample per site was retained. For the calculation of population genetic variables, the software of SPAGeDi v1.5 was used, which is suitable for tetraploids (Hardy and Vekemans, 2002; Vekemans and Hardy, 2004). At each site, the number of alleles (A), effective number of alleles (A_e), allelic richness (A_r), observed heterozygosity (H_o), expected heterozygosity (H_e), and individual inbreeding coefficient (F_i) were calculated. For comparison purposes, allelic richness was estimated for subsamples of the smallest sample size (A_r at k = 22). Clonal richness (R) was calculated with the formula R = (G−1)/(N−1) (Dorken and Eckert, 2001).

To detect isolation by distance (IBD) patterns, a one-sided Mantel test was performed on a matrix of the pairwise genotypic differentiation (pairwise F_st) between populations and a matrix of pairwise geographical distances, randomizing with 1,000 permutations using SPAGeDi v1.5. This test is based on an analysis of variance (ANOVA) and was performed separately for two different geographic distance matrices. The first included Euclidean distances between populations, the second was based on calculated geographic distances between populations following hydrological connections in salinas or coastal water flow. This was complemented with a more detailed ANOVA at the population level using a similar pairwise comparison and 1,000 permutations but dividing distance among four classes. At first, distance classes were defined by obtaining an equal number of pairwise comparisons of populations (3.42–7.79–12.88–22.58 km). Next, distance classes (4–9–14–25 km) were proposed based on the population’s respective distances from each other. The variance of the pairwise kinship coefficient (F_ij) (Loiselle et al., 1995) for all populations was in the same way calculated over complete distance and divided among distance classes. The genetic structure among populations (F_st), inbreeding within populations (F_is), and overall inbreeding (F_it) were estimated with ANOVA, and significance was tested with 1,000 random permutations. To identify fine-scale spatial genetic structures (FSGS), i.e., the non-random spatial distribution of genotypes (Vekemans and Hardy, 2004), pairwise kinship coefficient F_ij (Loiselle et al., 1995) was assessed between genets at the transect level. Population C was omitted from this analysis because the low coverage and fragmented distribution did not allow sufficient pairwise comparisons of individual genets within given distance classes to obtain statistical support. Pairwise kinship within sites collectively and individually was analyzed (one-sided ANOVA, 1,000 permutations) using the whole sample as a reference for allele frequencies. Five distance classes were tested, defined by an equal number of pairs (4–8–12.5–19–100 m) or based on genetic structures that have previously been found in seagrasses (1–5–10–15–100 m) (Orth et al., 1994; Phan et al., 2017; Dierick et al., 2021).

Patterns of population genetic structure, such as groups with genetically related individuals, can be visualized in clusters. Hence, a principal component analysis (PCA) was performed at the individual level using the R package ade4 (Dray and Dufour, 2007; R Core Team, 2019). Because PCA focuses on loci containing the most variation and we wanted to maximize separability between groups, complementary discriminant analysis of principal components (DAPCs) was computed using the Adegenet package in R (Jombart, 2008; R Core Team, 2019). Both are acceptable multivariate methods to analyze genetic patterns; however, the latter is designed to identify complex genetic structures in populations (Jombart et al., 2010). Using DAPC, genetic structure was inferred from 40 principal components and populations were partitioned into clusters (K = 3). Thereafter, we calculated the genetic differentiation (ρ_ST) among populations and clusters using an analysis of molecular variance (AMOVA) with the software GenoDive 3.05 (Meirmans and van Tienderen, 2004; Meirmans and Liu, 2018).

To interpret the observed genetic variation within populations, potential patterns of historical gene flow were investigated. Therefore, coalescent models were run using Migrate-n (Beerli et al., 2019) with two sets of populations based on the DAPC results. First, gene flow directionalities were explored between populations D and F, which belong to the same cluster. In the second analysis, the remaining eight populations were partitioned into a northern and southern group. Because Migrate-n is designed for diploid species (2×) and R. spiralis is polyploid (4×), our datasets first had to be restructured. As suggested by Beerli (2012), alleles of one individual were redistributed into two new diploid individuals, making sure to maximize heterozygosity and minimize effects on allele frequencies. While deduplicating is evident for diploid and tetraploid loci, it brings the challenge of determining which allele has to be duplicated from loci with three alleles. We rewrote them as two diploids using the allelic dosages. For instance, if a locus had three alleles (ABC) and we could infer which allele had the highest allelic dosage (e.g., B), this allele was duplicated (AB and BC). To analyze allelic dosages, GeneMarker V2.6.0 (SoftGenetics LLC^®) was used. When the allelic dosages were unclear, the allele with the highest frequency for this population was duplicated. To ensure consistency with the original data, chi-square tests were run on every locus for each group separately ( Supplementary Table S2 ). In addition, individuals containing missing allele values were discarded (∼28%), since using only completely genotyped individuals provided a sufficient large sample size.

Migrate-n allows to estimate the parameters θ (population size) and M (migration rate) under the assumption that they are constant through time. Uniform priors for these parameters were set (for θ, minimum = 0, maximum = 10, and delta = 0.1; for M, minimum = 0, maximum = 100, and delta = 10), and the program was run on a subsample of 24 (for populations D and F) or 30 (for north–south populations) individuals per population with the Brownian motion microsatellite model using Bayesian inference. Parameter values were searched using Markov chain Monte Carlo (MCMC) for 1 million generations with a burn-in of 10,000. Three models were explored to test different hypotheses of gene flow. The first model represented a full migration model with bidirectional migration rates. The additional two models tested unidirectional migration, each in the opposite direction. The Bezier approximation score for all loci was compared between models to find the most likely direction of gene flow. Using the estimated θ and M parameter values, the number of migrants per generation was calculated from that model by the following equation: Nm = (θ*M)/4.

The software DIYABC Random Forest (RF) v1.0.14-beta (Collin et al., 2021) was used to gain insight into the demographic history of R. spiralis populations. Seven scenarios of historical divergence were investigated in a model that was designed using the three groups adopted from the DAPC results (ABCE, GHIJ, and DF). In doing so, parameter priors were assumed uniform (for N_A, min = 10, max = 10,000; for N₁, min = 10 max = 10,000; for N₂, min = 10, max = 10,000; for N₃, min = 10 max = 10,000; for t_d, min = 0 max = 5,000; for t₁, min = 0 max = 5000) and all loci were considered dinucleotidic. A training dataset was generated by running 70,000 simulations. Hereby, all summary statistics available in DIYABC RF v1.0.14-beta were used. Thereafter, RF algorithms implemented in the program were used to predict the best scenario. The full training dataset was considered, and RF analyses were ran using 500 trees. The scenario with the highest classification vote was selected as the best-suited scenario for the target dataset. For this scenario, parameters (N_A, N₁, N₂, N₃, t₁, t_d) were estimated.

Since all populations contained several individuals with more than two alleles for several loci, there is an indication of tetraploidy. In two populations (A and C), we even detected a few individuals with up to five alleles for the same locus ( Supplementary Table S1 ). Although this might suggest the existence of hexaploid individuals, for the purpose of this study these samples were considered tetraploids, by considering the most common alleles. We observed plenty of multiallelic (three or four alleles per locus) situations in populations D and F when compared to all other populations. In total, we identified 147 alleles from 15 loci in 193 individuals (mean number A = 9.7, effective number A_e = 4.2, allelic richness A_r = 5.5), with minor differences among populations ( Table 1 ). Population J had the lowest number of alleles (A = 2.9 and A_e = 2.6). Averaged across all populations and loci ( Table 1 ), observed heterozygosity was lower than expected heterozygosity (H_o = 0.599, H_e = 0.671). The individual inbreeding Fi was significantly positive for seven out of 10 populations and significantly negative for the remaining ( Table 1 ). Population G showed the highest level of heterozygote shortage (F_i = 0.361), whereas population A showed the highest level of heterozygote access (F_i = -0.221). We identified 57 repeated MLGs from the original dataset. The largest clonal richness (R = 0.86-1 was observed in the ‘salina fridda’ water body) ( Table 1 ). Populations A and J featured genets that were extended along the majority of the transect, although they were intermingled with other genets.

Chloroplast sequences of samples from each population resulted had exactly similar lengths for the trnH-psbA intergenic spacer (452 bp) and of two mononucleotide microsatellite regions ccmp3 (157 bp) and ccmp10 (166 bp). For ccmp2, variants of either 202 and 203 bp were found in population C and populations E, F, and G, respectively, or were encountered mixed among individuals within the same population, namely, populations A, B, D, H, I, and J (GenBank accessions were provided in the Methods section). These sequences respectively correspond to haplotypes B1 and B3 of R. spiralis as defined in Triest and Sierens (2014).

The nuclear ribosomal cistron sequences corresponded to ITS-A of R. spiralis as defined in Triest and Sierens (2014). However, the complete cistron showed three additional variants due to mutations in position 690 (C or G)– and position 7,064 (C or T) ( Table 1 ). An rRna-cistron variant (C–C, respectively at positions 690 and 7,064) was observed in populations A and H, a C–T variant in populations B, D, and F, and a G–C variant in populations C, E, I, and J (GenBank accession ON527176, ON527177, and ON527178). The overlapping reads were consistent in most samples, although populations D and F showed an incomplete homogenization at nucleotide position 7,064 (61%–67% of T).

Global ANOVA-based F-statistics were calculated to obtain a general imprint of the genetic structure of R. spiralis along the analyzed coast of Sicily. Overall, we found moderate population differentiation (F_st = 0.185), moderate within-population inbreeding (F_is = 0.193), and higher overall inbreeding (F_it = 0.343). Mantel tests confirmed significant patterns of IBD (p < 0.01) among populations ( Figure 2 ). Firstly, populations are significantly more differentiated (F_st) over larger geographic distance (y = 0.004x + 0.14 with R² = 0.24). Secondly, Euclidean distances reflected genetic structure slightly better (R² = 0.24) than distances following hydrological connections (R² = 0.20). More precisely, pairwise F_st permutation testing of all populations within four distance classes revealed a significantly lower differentiation (p < 0.01) below 4 km and a significantly higher differentiation (p < 0.05) beyond 14 km ( Figure 2A ). As expected, an opposite pattern was obtained for the relatedness among individuals between populations F_ij (y = -0.0023x + 0.0019 with R² = 0.23 and p < 0.01). We found higher kinship below 4 km (p < 0.01).

A fine-scale genetic structure was detected among genets within populations. Across all populations, individual kinship over transect distance resulted in significant slopes of the regression with ln(distance) (y = -0.0091 * ln(x) + 0.162, with R² = 0.004 and p < 0.01). Divided among distance classes, genets that were less than 1 m apart had a significantly higher kinship than average (p < 0.001). However, a significant FSGS on genets was only observed in five out of nine populations (channel and lagoon populations, namely, populations A, B, G, H, J) with elevated values of kinship mostly within 1 m ( Figure 2B ). In populations A and J, a significantly higher value of kinship was found at higher distances (15–100 m). This is likely the result of gaps in the dataset that follow from the exclusion of pairs between ramets of the same genet. An FSGS was absent in the ‘salina fridda’ habitat.

The DAPC showed three strongly separated groups (K = 3) ( Figure 3A ) and is reflected by the PCA ( Supplementary Figure S3 ). Clusters 1 and 2, which are only separated along the y-axis of the DA, discriminate between the northern populations and southern populations. Cluster 3 exclusively contains samples of populations D and F and are clearly separated along the x-axis. Almost all populations contain low levels of admixture between clusters ( Figure 3B ). Individuals that were assigned to a different cluster than the majority of its population are characterized by alleles that can be considered common within the complete region of this study. The three clusters are differentiated by an AMOVA-based ρ_st value of 0.251 (ρst among populations = 0.279).

The Migrate-n analyses revealed the north to south directional models to be the most likely models to predict gene flow along the western coast of Sicily, when compared to the south to north unidirectional and bidirectional models ( Table 2 ). These models estimated a gene flow of Nem = 9.9 from the source population D to sink population F and a gene flow of Nem = 2.7 from the source populations A, B, C, and E to the sink populations G, H, I, and J.

The ABCDIY-RF analysis suggested that northern and southern populations have diverged more recently than, but not long after, the divergence between northern populations and populations D and F (Scenario 7 in Figure 4 , Table 3 ). Considering a generation time of 1 year, with potential overlap between generations, the divergence along this coastal stretch of the Western Sicily coastline appeared fairly recent and was estimated between 430 (lower CI at t₁) and 4615 (upper CI at t₁) generations.

Overall, the observed genetic structure along the western coast of Sicily contains patterns of distance-related differentiation between populations. A Mantel test resulted in a significant pattern of IBD. Moreover, pairwise F_st permutation testing within distance classes revealed a significantly lower differentiation between sites that were less than 4 km apart and a significantly higher differentiation between sites that were over 14 km apart. Nevertheless, the three clusters as suggested by the DAPC indicate that genetic structure is not solely reflecting an overall IBD pattern. One of the gene pools (cluster 3) that comprise individuals of populations D and F is characterized by a greater number of multiallelic loci and many unique alleles. Although geographically closely located to cluster 1, populations D and F are genetically relatively strongly diverged from others. This high differentiation can be explained by the historical divergence and by the relatively lower effective population size (N_e) of populations D and F. According to the DIYABC-RF analysis, D and F most likely diverged from (or alternatively, currently are more related to) the northern populations before the divergence between north and south over this regional distance. Since these time events were estimated relatively close to each other, the severity of differentiation is possibly driven further by the lower N_e of populations D and F. This is possible when restricted gene flow contributed to the maintenance of only limited alleles and/or alleles with deviant lengths. The main contributing alleles of the DAPC do tie in with this hypothesis. It could be discussed that model choice was not fully unambiguous. However, this is likely the result of gene flow allowing migrants from other groups to enter the gene pool, thus not having fully separated groups in the analysis. Indeed, individuals from cluster 1 were predominantly present in the north and individuals from cluster 2 in the south, although not exclusively. Based on the Migrate-n analysis between northern and southern populations, we can exclude that their differentiation is due to effects of drift only and can provide evidence of historical gene flow (i.e., Nem >1) (Wright, 1943). However, the estimated Nem values are low, indicating limited overall connectivity between north and south. In fact, the degree of migration may still be overestimated in the calculation of Nem. This is because a tetraploid migrant of R. spiralis carries more allele copies than a diploid migrant, hence having a greater effect on the evening out of allele frequencies (Meirmans et al., 2018). Thus, the ploidy level of R. spiralis ultimately increases the effect of migration to buffer population differentiation. Migrate-n additionally revealed that gene flow most likely followed a north to south direction. Interestingly, this directionality corresponds with the direction of sea currents along the western coast of Sicily and thus supports our hypotheses that hydrological connectivity may play an important role in the colonialization of R. spiralis in new habitats. The main sea current along the western coast of Sicily moves from north to south but differs in velocity at the local scale (Pisacane et al., 2018; Menna et al., 2019; Palma et al., 2020). Local conditions such as reduced velocity between the eastern coast of Isola Grande and the mainland might partially retain dispersal and restrict gene flow even further. In fact, isolating effects of currents in genetic structure have been observed in seagrasses (Ruiz-Montoya et al., 2015; Jahnke et al., 2019).

In contrast to our results, a study of the Southern Iberian Peninsula found a stronger influence of dispersal by birds than by sea currents (Martínez-Garrido et al., 2017). Previously, it has been discussed that these findings might be an effect of repeated sampling sites within the same lagoon (Triest et al., 2017). It is important to note that birds can only be considered a main disperser when they are abundant during the period of seed production, which for R. spiralis would be around August (Verhoeven, 1979). Although the natural reserve in Trapani is the main corridor in the central-western Mediterranean for migratory birds (Sorci et al., 1991), we expect that this will not have a considerable contribution to the dispersal of R. spiralis seeds in the Mediterranean. In studies on Ruppia where bird-mediated dispersal is evidenced in isolated lagoons, a scattered genetic structure is observed (Triest and Sierens, 2015). If waterfowl had a strong impact on gene dispersal of R. spiralis, we would expect more pronounced homogenization between populations. Although microsatellite alleles were substantially different and indicated three major clusters, the salina populations of western Sicily feature mixed populations of individuals assigned to different genetic clusters. The chloroplast haplotypes B1 and B3 indicated additional traces of mixed gene pools in most populations situated most closely and connected to the coastal zone, namely, A, B, D, H, I, and J. The more isolated inland locations appeared to contain a single haplotype. Three nuclear ribosomal cistron variants further subdivided the few studied individuals of the three major clusters whereas the maternal chloroplast markers suggest mixing of variants after seed dispersal between populations. Populations D and F remain most divergent for all markers and had incomplete homogenization of the ribosomal cistron suggesting historical introgression events. This idea is also supported by their higher level of allelic diversity and larger occurrence of individuals with up to four alleles in a locus. Together with the potential occurrence of hexaploid individuals in population A, the genetic clusters and mixing of variants suggest a dynamic situation as a result of occasional dispersal between populations. Even if tests of directionality point to a historical north to south migration for the region, one should consider as well random dispersal patterns and occasions of locally mixed gene pools comparable to observations of multiple hybridization events in the Camargue region (Beirinckx et al., 2020). Therefore, it is important to note that birds will certainly have their role in seed dispersal of R. spiralis. However, our results imply that the influence of currents could be more important in leaving traces of accumulated historical gene flow. In this aspect, R. spiralis resembles largely the survival strategy of seagrasses that depend on coastal currents than of, e.g., brackish water-submerged macrophytes such as Stuckenia pectinata (L.) Böerner (syn. Potamogeton pectinatus L.) that largely depend on bird-mediated seed dispersal. No doubt that R. spiralis undergoes both ways of dispersal, although hydrodynamics appear to be explanatory.

Spatial autocorrelation analyses on genets showed different FSGS patterns between populations. For example, significantly higher values of kinship at a close distance (<1 m) were detected in populations B, G, H, and J, whereas no significant structures were observed in four other populations (D, E, F, I). We can assume that gene dispersal remained stable over several generations in populations where FSGS was detected. Under these conditions, local seed rains might form, whereby seedlings settle in the proximity of the parental species (Triest et al., 2014). However, patterns on FSGS among genets do not allow to distinguish between the dispersal of seeds or pollen. In contrast, the absence of FSGS might be related to effects of physical disturbance that may upset local seed rains and/or may cause a high mortality rate and decrease the reproductive output. Consequently, the frequency of disturbances will influence the pattern of fine-scale gene dispersal. The fact that we could not find significant FSGS in the ‘salina fridda’ indicated that this habitat presents more disturbances for R. spiralis. It has been demonstrated that the waterbird communities in a Mediterranean salina differ between ponds and from natural salt marshes (Birtsas et al., 2011). Moreover, it was shown that flamingos (Phoenicopterus roseus P.) specifically prefer the ‘salina fridda’ (or prebasin in Birtsas et al., 2011) and were the most numerous species in that pond. The feeding behavior of flamingos has previously been brought up as a disturbance factor in Ruppia ponds (Rodriguez-Pérez and Green, 2006). Therefore, it can be hypothesized that flamingos negatively impact the formation of FSGS through random mixing of propagules in R. spiralis populations. Likewise, hydrological processes such as the opening and closing of water supply may be considered a disturbance factor.

Ruppia spiralis can reproduce sexually through outcrossing or clonally through rhizomes (Verhoeven, 1979; Gesti et al., 2005). It is not known which factors either affect clonal growth or induce flower production. It has been suggested that the stability of the habitat could explain the observed variation in clonal structure (Martínez-Garrido et al., 2017), as was observed in seagrasses (Eriksson, 1993; McMahon et al., 2017). Nearly all populations sustained a high genotypic diversity and a low frequency of clones, although much variation in clonal richness was found between populations (R = 0.48–1). In addition, we found high allelic diversity and high values of heterozygosity, suggesting that outcrossing is the preferred mode of reproduction. The flowers of R. spiralis are well adapted for outcrossing (Ackerman, 1995; Triest et al., 2018). Pollen is released at the water surface and floats upon the surface toward female flowers (Verhoeven, 1979). Because of this different reproductive ecology of R. spiralis, high inbreeding levels have not been detected within the species, whereas it has been in sister taxa R. maritima and R. brevipedunculata (Triest and Sierens, 2015; Triest et al., 2021a). Next to reproduction strategy, habitat preferences and ploidy level may be considered genetic drivers in Ruppia species. Overall, R. spiralis prefers more permanent waterbodies than R. maritima (Verhoeven, 1979; Triest et al., 2015; Beirinckx et al., 2020). Therefore, R. spiralis can form perennial populations that provide more time for mutations to accumulate within the clones. Similarly, more mutations per generation are possible in the tetraploid R. spiralis compared to the diploid R. maritima, because the amount of chromosome copies is doubled (Meirmans et al., 2018). The combination of these characteristics makes that R. spiralis can maintain a large amount of genetic diversity. Notwithstanding the overall high diversity, a higher proportion of clones was sampled in some populations (A, J, B, H, and G in decreasing order). These populations correspond to those with significant FSGS and consequently more stable habitats. It has been hypothesized that stable habitats increase the probability of rhizome survival, allowing clones to spread and colonize over a long period (Martínez-Garrido et al., 2017). Our results might support this hypothesis, although the observed high diversity is in favor of supporting a survival strategy of limited local (patchy) regrowth.

Anthropogenic disturbances are threatening wetland ecosystems all over the world. Because R. spiralis has a wide tolerance from brackish water to saltwater and can persist under many disturbances, it can provide habitat for marine fish and invertebrates and resting areas for many species of migratory birds in a variety of aquatic systems (Verhoeven, 1980). Salinas have been identified as valuable surrogate wetlands for shore and migratory birds, which may buffer against the impact of natural habitat loss (Weber et al., 1999; Masero, 2003). In the Mediterranean Basin, where salinas extend over 1,000 km² in area (Sadoul et al., 1998), their importance for conservation is without a doubt. However, there has been a trend to cease salt production in many parts of the Mediterranean. As a result, many salinas are inactive and their hydrological infrastructures are collapsing (Sadoul et al., 1998; Petanidou and Dalaka, 2009). Present data on R. spiralis highlight the importance of water flow in salinas. We therefore recommend active management of salinas, i.e., ponds have to be (temporarily) connected to and supplied with seawater, for maintaining long-term viable populations. A similar management plan has been advised for the conservation of salina water-bird species with regard to the protection of their nests (Birtsas et al., 2011). Since R. spiralis is the most dominant submerged macrophyte species that can persist in these salinas and deliver conditions for high-quality food in sufficient quantities to support waterfowl (Triest and Sierens, 2009), this will be beneficial for a wide range of wildlife.

The taxonomy of the Ruppia genus is complicated and has undergone many revisions over the past decades (den Hartog and Kuo, 2006; Triest and Sierens, 2014; Ito et al., 2017; den Hartog and Triest, 2020). In general, three species are recognized in the Mediterranean: R. spiralis, R. maritima, and the more restricted R. drepanensis Tineo. Several genome duplication events have characterized these species, causing different numbers of chromosomes to be observed (Triest and Symoens, 1991; Talavera et al., 1993; Mannino and Geraci, 2016; Beirinckx et al., 2020). Traditionally, R. spiralis is considered to be a tetraploid taxon. However, an earlier study on Sicilian Ruppia populations revealed the presence of a diploid R. spiralis (Mannino and Geraci, 2016). In the present study, we found, next to tetraploid R. spiralis, indications of hexaploid individuals of R. spiralis. A hexaploid R. spiralis has only once before been detected in northern Germany (Reese, 1962). In western Sicily, recent evidence also revealed the presence of R. drepanensis (De Castro et al., 2021). The diversity of Ruppia that coexists in this area may follow from the diversity gradient, observed over the whole Mediterranean Basin (Triest and Sierens, 2014). Hereby, the strait of Sicily could function as a funnel, where the high diversity of the western basin passes and accumulates. Further investigation is certainly needed to confirm this hypothesis and to confirm if hexaploid R. spiralis truly exists in Sicily.

In conclusion, high levels of allelic diversity were detected in R. spiralis at the regional scale. Distinct genetic clusters indicate that its distribution along western Sicily is shaped by different historical colonization events and geographic distance. The observed genetic structure follows an overall north-to-south migration pattern, confirming expectations of gene dispersal by currents. To allow gene dispersal between populations, managers should be aware of the importance of hydrological connectivity. Despite significant structures at a regional scale, significant FSGS were only detected in some populations. Our results suggest that salinas are valuable habitats for harboring diverse R. spiralis populations at a regional scale, although these may be subject to a local physical disturbance.