Phylogeography of Aphanius fasciatus (Osteichthyes: Aphaniidae) in the Mediterranean Sea, with a focus on its conservation in Cyprus

Aphanius fasciatus is a small fish occurring in Mediterranean brackish environments. In Cyprus it is known from three localities separated by long stretches of coast. The genetic diversity of these populations was evaluated using fragments of two mitochondrial genes. A comparison with the other available data showed that Cyprus populations represent a distinct lineage. The other lineages are concentrated in a relatively small area between the Strait of Sicily and the Western Ionian Sea, while all other areas include a subset of these lineages, suggesting that the aforementioned area might have acted as a glacial refugium. Landlocked North-African populations diverge from all other populations, suggesting that they might have originated in the Late Pleistocene, during transgression events of the Mediterranean Sea in North-African inland water bodies. The genetic diversity of A. fasciatus varied across different Cyprus populations, with a pattern mirroring the degree of environmental degradation, which likely affected population genetic variability through demographic reductions. The three Cyprus populations showed genetic uniqueness, suggesting the need of population-based management practices; the low genetic diversity of two populations, and the number of threats affecting them, suggest that the species should be considered endangered at national level and deserves protection measures.


Introduction
The Mediterranean killifish, Aphanius fasciatus (Valenciennes, 1821), is a small fish typically associated to coastal brackish-water environments of the whole Mediterranean region, with the exception of the westernmost part, where it is replaced by Apricaphanius iberus (Valenciennes, 1846) and Apricaphanius baeticus (Doadrio, Carmona & Fernández-Delgado, 2002), and possibly the Aegean Sea, where the closely related Aphanius almiriensis Barbieri, Kottelat & Stomboudi, 2007 may partially or totally replace it (Valdesalici et al., 2019;Freyhof & Yogurtcuoglu, 2020). In the South-Eastern Mediterranean Sea the distribution of A. fasciatus partially overlaps that of Aphaniops dispar (Rüppell, 1829), with which it is able to produce hybrid offspring (Villwock, 1987;Fouda, 1995;Lotan & Ben-Tuvia, 1996). In addition, A. fasciatus successfully colonised the Northern part of the Suez Canal (Fouda, 1995), with a few records in the Red Sea (Freyhof & Yogurtcuoglu, 2020), making it one of the few anti-Lessepsian migrants known to date.
Aphanius fasciatus is characterised by habitat fidelity in its adult phase. In addition, eggs are laid on benthic vegetation, and hatchlings show a lifestyle similar to that of adults (Peloso, 1946;Tortonese, 1986;Maltagliati, 1998a). Despite the sporadical presence in coastal marine waters (Torchio, 1967;Ç oker & Akyol, 2014), possibly related to dramatic events, such as floods or dystrophic crises (Maltagliati, 1998a), these features contribute to a very restricted gene flow that may occur only between adjacent populations, rendering this species a valuable model for the reconstruction of microevolutionary processes (Maltagliati, 1998a(Maltagliati, , b, 1999. Although this species is listed as ''Least Concern'' (LC) by the International Union for Conservation of Nature (IUCN) (Kottelat & Freyhof, 2007), as the other Mediterranean species of killifish it is affected by impacts such as habitat degradation and fragmentation, pollution and competition with introduced species (Valdesalici et al., 2015). This led to a strong reduction of populations, with negative impacts on their genetic diversity (Maltagliati, 2002;Cimmaruta et al., 2003;Angeletti et al., 2010), with some documented instances of local extinction (Changeux & Pont, 1995;Ferrito & Tigano, 1996;Valdesalici et al., 2015;Taybi et al., 2020). Moreover, the fragmentation of the species in discrete populations and the low gene flow detected between them (Maltagliati, 1998a(Maltagliati, , 1999Ferrito et al., 2013;Buj et al., 2015;Cavraro et al., 2017), even at very low spatial scale , suggest that management practices cannot disregard the geographical distribution of genetic variability. The absence of dispersal stages, and the scarce connectivity retrieved even between adjacent populations, further suggest that local extinction cannot be expected to recover naturally, and that each population might deserve specific management practices. Despite the high number of works taking into account the distribution of A. fasciatus across the Mediterranean Sea, however, several populations are partially or totally unknown, and the lack of demographic and genetic data prevents a precise assessment of their conservation status.
The island of Cyprus represents a paradigmatic case for A. fasciatus conservation. Despite the presence of several coastal water bodies that might be potentially suitable for this species, only three populations of A. fasciatus are known, with long stretches of coast separating them. At least one of these populations was discovered only recently, and two of them are facing environmental threats of anthropogenic origin (Englezou et al., 2018). These populations are unknown from the demographic and genetic point of view, as the vast majority of Eastern Mediterranean populations of A. fasciatus. The aim of this work is the genetic characterisation of Cyprus populations of A. fasciatus through mitochondrial molecular markers, in order to compare them with the other Mediterranean populations for which genetic data are available, and to reconstruct the phylogeography of this species in the Mediterranean Sea in a conservation context.

Sampling
A total of 72 individuals of Aphanius fasciatus were sampled in November 2017 by hand net in the three known localities in Cyprus (Englezou et al., 2018), namely Glapsides Beach, Famagusta Bay (35.1604°N; 33.9125°E), Morphou Bay/Güzelyurt Körfezi (35.2013°N;32.9143°E) and Akrotiri (34.6053°N; 32.9367°E). The wetland habitat at Famagusta Bay is constituted by two coastal ponds located in Glapsides Beach and Silver Beach, representing part of the Pedieos River delta. After an event of mass fish mortality that occurred in 2012 in Silver Beach pond (Zogaris, 2017), the species has not been observed in the latter site anymore (C. Englezou, pers. obs.); therefore, the occurrence of A. fasciatus in Famagusta Bay is nowadays restricted to Glapsides Pond. A. fasciatus habitats at Morphou are part of the Serakhis River flood plain, which historically experienced irregular, intermittent flooding (Ward et al., 1958). Several potentially suitable habitats were subjected to extensive draining during the 1950s-1960s, and the area is nowadays heavily impacted by agriculture, animal farming, mine pollution and illegal industrial waste dumping (Karafistan & Gemikonakli, 2019). Akrotiri Marsh is a natural wetland in Cyprus covering an area of approximately 1.5 km 2 . It is part of a Ramsar site, an Important Bird Area and a designated Special Protection Area, equivalent to the EU designation, according to the mirror law (26/2007) of the Cyprus Sovereign Base Areas (Zogaris, 2017). Since 1990s, the reduction of pond water salinity, due to human intervention aimed at ameliorating the stormwater drainage system in the Limassol/Amathus area determined the proliferation of the mosquitofish Gambusia holbrooki Girard, 1859, with an increased competition pressure on the A. fasciatus population.
Furthermore, for molecular comparison, individuals of A. fasciatus from Lesina Lake (41.8860°N; 15.3640°E) and Orbetello Lagoon (42.4356°N; 11.2070°E) were sampled in April 2016 and July 2018, respectively (Table 1). Specimens were euthanized in ice-cold water, fixed in 96% ethanol and preserved at -20°C until DNA extraction. Voucher specimens were deposited in the Museo di Storia Naturale of the University of Pisa (Italy).

Molecular analyses
DNA extraction was carried out on 5 mg of muscular tissue using a salting-out protocol modified after Aljanabi & Martinez (1997). Tissue lysis was obtained by incubation at 55°C for 2 h in 290 ll of TNE 19 buffer and 20 ll of protease K; residual RNA was degraded adding 10 ll of RNAse A and incubating for 10 min at room temperature. Lipids and proteins were precipitated by adding 100 ll of NaCl 6 M and centrifuging at 12,5009g for 18 min. Absolute ethanol at -20°C (800 ll) was added to the surnatant containing genomic DNA; the solution was incubated at -20°C overnight and then centrifuged at 12,5009g for 15 min to precipitate DNA. The resulting DNA pellet was cleaned from impurities by adding 300 ll of 70% ethanol and centrifuging at maximum speed for 2 min; this step was repeated thrice. The pellet was then dried at 55°C, resuspended in 100 ll of ultrapure water, and stored at -20°C.
We amplified a fragment of the mitochondrial control region (CR) and a fragment of the gene coding for the subunit I of the cytochrome c oxidase (COI) using the primer pairs CRAfF (5 0 -ACTATTCTTTGCCGGATTCTG-3 0 ) (Tigano et al., 2004) and H16498 (5 0 -CCTGAAGTAGGAACCAGATC-3 0 ) (Meyer et al., 1990), and FishF1 (5 0 -TCAACCAACCACAAAGA-CATTGGCAC-3 0 ) and FishR2 (5 0 -ACTTCAGGGT-GACCGAAGAATCAGAA-3 0 ) (Ward et al., 2005), respectively. Polymerase chain reaction (PCR) amplifications were carried out in 20 ll solutions using 1.5 mM of MgCl 2 , 0.2 mM of each dNTP, 0.1 lM of each primer, 1 U of DreamTaq DNA polymerase (Thermo Scientific), and * 2.5 ng of template DNA. For both markers the PCR profile was set as follows: initial denaturing step at 94°C for 3 min; 40 cycles of denaturing at 94°C for 45 s, annealing at 45°C for 1 min, and extending at 72°C for 1 min, and a final extending step at 72°C for 7 min. A negative control was included in each reaction. PCR products were precipitated with sodium acetate and absolute ethanol and sent to GATC Biotech (Eurofins Genomics) for sequencing.

Data analyses
All available sequences for the two examined mitochondrial regions were obtained from GenBank (Table 1; Supplementary Material 1; Supplementary Material 2) and haplotype frequencies were retrieved from the original publications whenever possible for comparison purposes. Sequences from each gene were aligned with ClustalX 2.1 (Larkin et al., 2007), and alignments were edited in BIOEDIT version 7.2.5 (Hall, 1999). The program jModelTest 2.1.6 (Guindon & Gascuel, 2003;Darriba et al., 2012), based on the hierarchical likelihood ratio test, was used to assess the best model of evolution for the sequences under the Akaike Information Criterion (AIC, Akaike, 1974 Excoffier et al., 1992) were calculated with Arlequin v. 3.5.2.2. (Excoffier & Lischer, 2010) considering the nucleotide substitution model identified through jModelTest 2.1.6. Two 12 0.682 ± 0.148 0.0027 ± 0.0008 ---AMOVAs (Excoffier et al., 1992) were carried out to examine the partition of genetic variance into (i) the ''within locality'' and ''among locality'' components (2-level AMOVA) and (ii) the ''within locality'', ''among localities within biogeographical sector'' and ''among geographical sectors'' components (3-level AMOVA). An additional AMOVA was performed in order to check the isolation between Cyprus populations. The significance of U-statistics parameters was assessed by permutation tests with 10,000 replicates as implemented in Arlequin v. 3.5.1.2. A median-joining haplotype network was constructed by using NET-WORK v. 10.1 (Bandelt et al., 1999). BAPS v. 6.0 (Corander et al., 2003;Corander & Marttinen, 2006) was employed to detect possible hidden population substructuring by clustering genetically similar sampled individuals into panmictic groups, hereafter called ''genetic clusters''. We applied Bayesian assignment analysis to the 317 bp CR fragment common to ours, Pappalardo et al.'s (2008), Ferrito et al.'s (2013), Buj et al.'s (2015) and Cavraro et al.'s (2017) datasets, for 896 sequences overall, and to the 491 bp COI fragment common to ours, Geiger et al.'s (2014), Landi et al.'s (2014) and Valdesalici et al.'s (2019) datasets, for 128 sequences overall. BAPS adopts a Bayesian approach with a stochastic optimization algorithm for analysing models of population structure, which greatly improves the speed of the analysis compared to traditional MCMCbased algorithms (Corander & Marttinen, 2006). The only prior information given was the sampling locality of each individual. When testing for genetic clusters, in order to optimise computing time and check the approximate number of clusters in our dataset, we ran an explorative analysis with two replicates for each value of k (the maximum number of clusters) from k =1 to n ? 3, where n = number of populations (Evanno et al., 2005). This analysis detected nine clusters for CR and five for COI; then, we carried out a final run with five replicates for each value of k up to k = 15. In addition, we used a number of reference individuals nri = 500 and repeated the admixture analysis 500 times per individual. Individuals with ambiguous assignment (associated P \ 0.05) were not visualised in the graphs.
Demographic history was inferred by the analysis of the distribution of the number of site differences between pairs of CR sequences (mismatch distribution), which was carried out on all populations using DnaSP. Demographic analyses for COI sequences were not carried out due to the uneven number of individuals across populations and low number of individuals in many populations. Expected values for a model of constant population size were calculated and plotted against the observed values. According to Rogers & Harpending (1992), populations that have experienced a rapid demographic growth in the recent past show unimodal distributions, whereas those at demographic equilibrium present multimodal distributions. Harpending's (1994) raggedness index (Hri; quantifying the smoothness of the mismatch distributions and distinguishing between population expansion and stability) and the sum of squared deviations (SSD), as implemented in ARLEQUIN, were used to test Rogers' (1995) sudden expansion model, which fits to a unimodal mismatch distribution (Rogers & Harpending, 1992). Moreover, population expansion was tested through Fu's (1997)  Geographical sectors were established considering Bianchi's (2007) biogeographical characterisation of the Mediterranean Sea. The sequences obtained in the present study, together with those downloaded from GenBank originated from five biogeographical sectors, namely the Sea of Sardinia (sector 4), the Tyrrhenian Sea (sector 3), the Ionian Sea (sector 9), the Adriatic Sea (sectors 6, 7 and 8) and the Levant Sea (sector 12). Although Bianchi (2007) argued that the Adriatic Sea should be subdivided into three different biogeographical sectors, we decided to consider it as a single biogeographical sector, due to the post-glacial recolonisation of the Northern and Central portions of this sea (Bianchi et al., 2012). Moreover, we decided to consider the Strait of Sicily as a further biogeographical sector, as it represents the transitional area between the Western and Eastern Mediterranean basins, even though its effect on genetic connectivity in marine organisms varies depending on species biological characteristics, such as developmental mode, duration and life style of larval stages, and adult dispersal capabilities (Bianchi, 2007;Villamor et al., 2014). Following the results by Ferrito et al. (2013) we decided to consider landlocked populations of Northern Africa as belonging to a separate biogeographical sector.

Control region (CR)
Molecular analyses allowed to obtain 317 bp CR sequences for the 96 individuals of A. fasciatus (GenBank accession numbers: MW147304.1-MW147321.1). Eight haplotypes were identified in samples from Cyprus. Interestingly, all of them were private to the island and only the most frequent haplotype was shared between the three populations, whereas the remaining seven haplotypes were locality-private (Supplementary Material 3). The highest values of haplotype and nucleotide diversity were observed in the Akrotiri population (h = 0.636 ± 0.054, p = 0.0024 ± 0.0003), while the lowest values were observed in the Morphou population (h = 0.222 ± 0.110, p = 0.0010 ± 0.0005), and the Glapsides population had intermediate values for both parameters (h = 0.290 ± 0.109, p = 0.0017 ± 0.0007). When the two parameters were calculated for all populations available in literature (Table 1), the average values of haplotype diversity and nucleotide diversity were h = 0.405 ± 0.240 and p = 0.0029 ± 0.0032, respectively. Both parameters showed high variation across Mediterranean populations, as testified by the high standard deviation associated to the average values.
The complete dataset included 132 haplotypes, most of which private to a single biogeographical sector ( Fig. 1). One haplotype was shared between Strait of Sicily and Ionian Sea, two between Ionian Sea and Adriatic Sea, and another one between Strait of Sicily, Sea of Sardinia and Tyrrhenian Sea. More generally, in the haplotype network Adriatic haplotypes belonged to two distant haplogroups, all haplotypes from Cyprus formed a haplogroup separated by six mutations from the closest Adriatic haplotype, and all haplotypes from internal North-African water bodies formed a haplogroup separated by 11 mutations from the closest haplotype from the Strait of Sicily. Haplotypes from the Ionian Sea were distributed into three haplogroups, two of which including Adriatic haplotypes as well, the third one including haplotypes from the Strait of Sicily; moreover, two haplotypes, each shared by two individuals, were not particularly close to any haplogroup. Lastly, no evident groupings were observed in haplotypes sampled in the Strait of Sicily, Sea of Sardinia and Tyrrhenian Sea (Fig. 1).
The 2-level AMOVA identified the ''among-populations'' component as the main source of variance, accounting for 78.8% of the variation retrieved, while the remaining 21.2% referred to the ''within population'' component ( Table 2). The 3-level AMOVA identified the higher source of variance in the ''among population within biogeographical sectors'' component (48.9%), while the ''among biogeographical sectors'' component accounted for the 31.1%, and the ''within population'' component for the 20.0% of the variation retrieved ( Table 2). The AMOVA carried out on Cyprus populations only identified the ''within population'' component as the main source of variance, accounting for 72.7% of the variation retrieved, while the remaining 27.3% referred to the ''among population '' component (Supplementary Material 4). All the associated U-statistics value were significantly greater than zero ( Table 2).
The Bayesian assignment analysis carried out on the complete dataset detected nine genetic clusters with maximum-associated probability value (P = 1)  Two levels were considered in the first analysis, whereas in the second analysis the 'biogeographical sector' level was added. Probability values were obtained after a permutation test with 10,000 replicates ( Fig. 2). Twenty-four sequences (2.7% of the complete dataset) downloaded from GenBank were of ambiguous assignment. The majority of populations belonged to a single cluster, but populations in the Sea of Sardinia, Strait of Sicily and Ionian Sea showed the occurrence of more than one cluster. In addition, genetic clusters showed a different repartition across biogeographical sectors (Fig. 2) (Table 3).

Cytochrome c oxidase -subunit I (COI)
In this study a 491 bp fragment of the COI gene was amplified and sequenced for 85 individuals of A. fasciatus (GenBank accession numbers: MW138950.1-MW138959.1). We were not able to obtain any additional sequence for Lesina population, possibly due to a mutation in the primer annealing region. Of the five haplotypes observed in Cyprus individuals, one was shared by the three populations and another one by Akrotiri and Glapsides and each (h = 0.620 ± 0.057, p = 0.0031 ± 0.0003), followed by Glapsides (h =0.338 ± 0.120, p = 0.0007 ± 0.0003); Morphou showed remarkably lower values for both parameters (h = 0.083 ± 0.075, p = 0.0002 ± 0.0002). As the majority of the localities sampled in literature included 2-6 individuals, which are insufficient to have reliable estimates, a comparison was made only with Orbetello (N =16; h = 0.517 ± 0.132, p = 0.0022 ± 0.0009) and the Sicilian population sampled by Landi et al. (2014) (N = 11; h = 0.600 ± 0.154, p = 0.0014 ± 0.0005) ( Table 1).
The complete dataset included 23 haplotypes, most of which restricted to a single biogeographical sector. Northern Tyrrhenian individuals (Orbetello Lagoon) were included in a cluster close to a haplotype shared by the Sicilian individuals and the two individuals from Cagliari; on its part, this haplotype was separated by only 1-2 mutations from the two haplotypes sampled in the south-Eastern Adriatic locality of Ston (Fig. 3). Cyprus haplotypes clustered together and were separated by this group by only one mutation. Individuals sampled in the landlocked North-African populations of Zahzah (Tunisia) and Touggourt (Algeria) represented a cluster separated by 4 mutations from the closest haplotype, and also individuals from the Western and north-Eastern Adriatic Sea (Varano and Lesina Lakes, Po River Mouth, Veneto and Pag) represented a clearly defined haplogroup, separated by 5 mutations from the closest haplotype (Fig. 3).
The 2-level AMOVA identified the ''among populations'' component as the main source of variance, accounting for 83.5% of the variation retrieved, while   Material 7). All the associated U-statistics value were significantly greater than zero. The Bayesian assignment analysis carried out on the complete dataset detected five genetic clusters with maximum-associated probability value (P = 1) (Fig. 4). Individuals were assigned to clusters with probability values of P = 0.01 to 1; none of the sequences showed ambiguous assignment. The majority of the sampled populations belonged to a single cluster, but individuals of the Akrotiri population were assigned to two clusters.
Genetic clusters showed a different repartition across biogeographical sectors if compared with the CR diversity pattern, with the most widespread cluster occurring in the populations of Orbetello (Northern Tyrrhenian Sea), Cagliari (Southern Tyrrhenian Sea), Sicily (Strait of Sicily) and Ston (South-Eastern Adriatic Sea). One cluster was restricted to Tunisia (Strait of Sicily) and the landlocked Algerian population of Touggourt, while another one included all individuals from the Western and North-Eastern Adriatic populations. Cyprus populations were assigned to two further clusters, one of which was represented in all populations, being the other one private to Akrotiri.

Discussion
Have Pleistocenic events moulded the genetic structure of A. fasciatus? Molecular data on A. fasciatus obtained in the present work together with those gathered from previous studies allowed to give a deep insight into the pattern of species' genetic diversity across its geographical range. Bayesian analyses carried out on the CR dataset highlighted the occurrence of eight out of nine genetic clusters in an area including the Strait of Sicily and the Western Ionian Sea (henceforth SS-WIS). In other terms, all lineages retrieved in the Western and Central Mediterranean Sea occur in this area, and three of them are private to the SS-WIS, while the lineage occurring in Cyprus is private to the Eastern Mediterranean Sea. COI sequences showed a slightly different pattern, with one lineage widespread from the Tyrrhenian Sea to the South-Eastern Adriatic Sea, one restricted to landlocked populations of Tunisia and Algeria, one occurring in the Western and North-Eastern Adriatic Sea, and two lineages restricted to Cyprus. The different degree of coverage of the two datasets, in term of both number of populations and individuals, could account for the differences retrieved between the outcomes from the two mitochondrial markers. In particular, the genetic cluster occurring in the Western and North-Eastern Adriatic Sea has been reported in the SS-WIS only at Ganzirri Lake, for which COI sequences were not available. The SS-WIS hosts a number of lineages that is much higher than those occurring in the remaining biogeographical sectors of the Mediterranean (one or two, none of which is private). This genetic richness supports the hypothesis that SS-WIS area may have acted as a sort of refuge during Pleistocenic glaciations (Chefaoui et al., 2017), representing nowadays a hotspot of species' genetic diversity. This pattern is consistent with Petit et al.'s (2003) theoretical considerations on the effect of glacial refugia on species genetic structure. On the other hand, different populations of A. fasciatus in the Strait of Sicily tend to host only one or two lineages, suggesting that gene flow between coastal populations is very limited (Maltagliati, 1998a(Maltagliati, , b, 1999, as indicated by AMOVA. Northern Tyrrhenian populations were assigned to a genetic cluster that occurs in coastal populations of Northern Tunisia, while populations from Sardinia were assigned to two genetic clusters occurring in both Western Sicily and Malta. The Adriatic Sea seems to have been colonised two times independently. The Western and North-Eastern parts of the basin host a genetic cluster that is shared with the Ganzirri Lake population. Possibly individuals of A. fasciatus colonised the Western Adriatic Sea from the South, along the Western coast. The South-Eastern part of the basin hosts a genetic cluster occurring in the Eastern Ionian Sea and in the Western Ionian population of Foce Marcellino, which likely reflects colonisation of the Adriatic Sea along the Eastern coast. An alternative hypothesis might take into consideration an unintentional introduction of Adriatic individuals in Ganzirri Lake, as suggested by Rocco et al. (2007), and supported by Cavraro et al. (2017), who argued that translocation of A. fasciatus individuals is a possible event. At this regard, humanmediated introduction has been considered the most likely explanation for the occurrence of the closely related Aphanius almiriensis Kottelat, Barbieri & Stoumboudi, 2007 in the Palude del Capitano (Southern Italy) (Valdesalici et al., 2019). More specifically, eggs or adults of A. almiriensis might have been accidentally introduced together with fishes for human consumption from Greece during the Roman Age. However, the Ganzirri Lake is not known to have hosted fish-rearing facilities in the Classical Age, as commercial aquaculture in the Cape Peloro lagoons dates back to the beginning of the XX Century and mostly concerns shellfish (Prioli in VV.AA., 2011). Although the introduction of some marine taxa in the Cape Peloro lagoons has been ascribed to commercial shellfish restocking (Di Natale, 1982), the biological cycle of A. fasciatus does not entail any dispersal stage (Tortonese, 1986) that could be effectively transferred with hard-bottom or soft-bottom marine molluscs. More generally, the hypothesis of an introduced status of the Venice Lagoon population (Nardo, 1847) and its possible Greek origin (Cavraro et al., 2017) are not confirmed by the available molecular data, showing a general genetic homogeneity throughout the whole Western and North-Eastern Adriatic Sea. The only exception to this general trend is represented by a single individual from Comacchio lagoon, which showed a CR haplotype belonging to the South-Eastern Adriatic-Eastern Ionian cluster (Fig. 3). As this is the only occurrence of this haplogroup in the Western Adriatic, we suggest that transfer of A. fasciatus between different localities might indeed occur, but it is not a relevant phenomenon in shaping the pattern of species' genetic diversity. Moreover, since the potentially introduced haplogroup occurs also in the Southern Adriatic Sea, the introduction from Greek lagoons is in our opinion a rather unlikely explanation. In addition, the occurrence of incomplete lineage sorting cannot be excluded. Otherwise, there is no co-existence of the two genetic clusters in the same population (unlike in other areas), suggesting the existence of a biogeographical break between the North-Eastern and the South-Eastern Adriatic Sea.
The populations of Cyprus do not belong to any of genetic clusters occurring in the SS-WIS, although they are not distant from them in terms of mutational steps. This suggests that glacial refugia for these species could have been elsewhere, possibly in an area within the Eastern Mediterranean Sea. A more precise reconstruction is hindered by the absence of data for populations of the Eastern Mediterranean Sea, where the species is locally frequent (Villwock, 1964;Lotan & Ben-Tuvia, 1996). On the other hand, the occurrence of A. fasciatus in the Aegean Sea has been debated, as it was historically confused with A. almiriensis, which on turn has been reported only a few times (Valdesalici et al., 2019). The hypothesis that the latter species might largely replace A. fasciatus in the Aegean Sea is consistent with results by Triantafyllidis et al. (2007), who identified through mitochondrial sequences a mean divergence of 3.45% between Aegean populations identified in that work as A. fasciatus and populations from the remaining parts of the Mediterranean Sea. This divergence is close to the 4.5% of divergence detected between A. fasciatus and A. almiriensis using COI fragments by Valdesalici et al. (2019), and the highest mutation rate of the COI gene if compared to the fragments analysed by Triantafyllidis et al. (2007) could account for the difference. Nevertheless, a biogeographical break in the Southern Aegean Sea has been highlighted in the phylogeographical pattern of several marine species (Nikula & Väinolä, 2003;Domingues et al., 2005;Pérez-Losada et al., 2007;Pannacciulli et al., 2017) and might account for the differentiation between A. fasciatus and A. almiriensis.
A stable difference between landlocked North-African populations and the remaining dataset is clear from both markers, although, unfortunately, the two datasets include different landlocked populations.
Nonetheless, landlocked populations are separated by 10 and 5 mutational steps from the closest nonlandlocked haplotype for CR and COI, respectively. This outcome suggests that landlocked populations have been isolated from the coastal ones since enough time to grant a significant diversification. Ferrito et al. (2013) proposed a pre-Pleistocenic isolation of these populations; however, this hypothesis is based on a very conservative estimate of the mutation rate of the control region if compared to the values retrieved in other fish species (McMillan & Palumbi, 1997;Burridge et al., 2008). Geological data showed that stable connections of brackish or hyperhaline inland water bodies of North Africa with the sea occurred several times in the Late Pleistocene, due to sea-level rise and consequent transgression of the sea in the depressed areas of the continent (Richards & Vita-Finzi, 1982;Causse et al., 1989Causse et al., , 2003. The last event of such transgression dates back to 55-18 kya (with a maximum around 30 kya) (Causse et al., 2003-see also Richards & Vita-Finzi, 1982), coarsely corresponding to the Å lesund interstadial period, that showed a partial reduction of the ice sheet, with subsequent increase of the sea level (Mangerud, 1981). At present, a more precise identification of a Late Pleistocene event that might have caused the isolation of North-African landlocked populations of A. fasciatus is impossible. However, we would discard the most recent dates, as they would imply an exceedingly high mutation rate. Instead, transgression events that occurred approximately 200 kya (Causse et al., 2003) would be compatible with the known range of mutation rates for fish control region. The number and extent of these events, and in particular of the most recent one, has been largely debated (Richard & Vita-Finzi, 1982;Causse et al., 1989;Vita-Finzi et al., 1991;Causse et al., 2003). Results from this study suggest that isolation of A. fasciatus in North-African inland water bodies dates back to the oldest of such events, supporting Causse et al.'s (1989Causse et al.'s ( , 2003 hypothesis of a more limited transgression occurred 30 kya. A unique phylogeographical pattern among Mediterranean fishes A comparison of the phylogeographical pattern retrieved with other fish species with a similar distribution shows clear differences, possibly related to the peculiar life cycle traits of A. fasciatus. The majority of studies carried out on Mediterranean marine demersal fish species retrieved a high connectivity across the basin (Viñas et al., 2004;Bargelloni et al., 2005;Domingues et al., 2005;Barbieri et al., 2014). In these species divergent lineages are often cooccurring and show similar frequencies, possibly testifying events of past vicariance, but also confirming the occurrence of secondary contact events (Angiulli et al., 2016;Barros-García et al., 2020). A similar pattern can be identified in benthic species with pelagic larvae (Suzuki et al., 2004;Francisco et al., 2014). A break between the Aegean Sea (including any population in the Black Sea) and the remaining Mediterranean Sea has been identified in some demersal species (Suzuki et al., 2004;Domingues et al., 2005), but in these cases the degree of genetic divergence is much lower than that identified between A. fasciatus and A. almiriensis, suggesting a less pronounced effect of the Southern Aegean as a phylogeographical barrier. Conversely, although A. fasciatus might thrive in internal water bodies, its ability to inhabit coastal environments and to occasionally use marine environments for dispersal (Torchio, 1967;Ç oker & Akyol, 2014) makes its phylogeographical pattern completely different from those identified in exclusive freshwater species, that are often restricted to a single river drainage (Perea & Doadrio, 2015;Lorenzoni et al., 2021), and usually show the occurrence of several inland glacial refugia (Costedoat & Gilles, 2009).
A comparison with fish species typically associated to brackish-water environments, and with an apparently similar ecology, provides more suitable comparisons. A strong structuring between Western and Eastern Mediterranean, and among different lagoon systems, has been revealed for the sand gobies Pomatoschistus marmoratus (Risso, 1810) and Pomatoschistus tortonesei Miller, 1969 (Mejri et al., 2009(Mejri et al., , 2011. The divergence between P. marmoratus populations from different lagoon systems might be enough to consider them as different species (Mejri et al., 2011). In general, deep divergence between separate populations and even cryptic speciation events are relatively frequent in euryhaline Gobiidae (Stefanni & Thorley, 2003;Neilson & Stepien, 2009). This is probably due to the absence of dispersal stages combined with the exclusively benthic adults; conversely, the potential for migration across stretches of sea allows a shallower structuring in A. fasciatus. A similar pattern of substantial isolation between populations occurring in different brackish-water bodies was retrieved in the pipefish Syngnathus abaster Risso, 1827, a species occurring mostly in brackishwater environments across the Mediterranean Sea and characterised by egg brooding in a ventral pouch and hatching of benthic juveniles (Sanna et al., 2013). However, unlike A. fasciatus, S. abaster shows as well a further structuring into three deeply divergent clades occurring in different areas of the Mediterranean Sea (Western Mediterranean, Tyrrhenian ? Adriatic Sea and Strait of Sicily), suggesting that this strictly sedentary benthic species is in fact unable to effectively disperse by rafting, and that its dispersal capabilities are closer to those of euryhaline gobies than to those of A. fasciatus. The three-spined stickleback Gasterosteus aculeatus Linnaeus, 1758 is characterised by a high tolerance towards salinity variations, although it is not as pronounced as in A. fasciatus. It is, however, a species with boreal affinity, showing a higher capability for migration across marine environments, and a rather scattered distribution in the Mediterranean basin, where it occurs mostly in freshwater environments. As a consequence, stickleback populations occurring in coastal Mediterranean environments are completely isolated from each other and originated independently from an ancient, panmictic European population (Mäkinen & Merilä, 2008). Finally, the Iberian killifish Apricaphanius iberus is very similar to A. fasciatus as regards its ecology and life history traits, despite showing a much narrower distribution. This species and A. fasciatus show similar levels of genetic divergence between populations, despite the very different spatial scale (Pappalardo et al., 2015). In particular, Pappalardo et al. (2015) retrieved the same CR haplotype shared between the populations of S'Ena Arrubia (Sardinia) and Trapani (Sicily); by comparison, all the five populations of A. iberus assayed, coming from a relatively restricted area in South-Western Spain, included only private haplotypes.
In conclusion, the phylogeographical pattern retrieved in A. fasciatus is completely different from both marine demersal fishes and freshwater Mediterranean fishes; the effect of the South Aegean break might recall a pattern identified in some marine species, but it is far more pronounced in the A. fasciatus group, where it actually leads to a differentiation at species level. However, the phylogeographical pattern retrieved in A. fasciatus does not fully match that retrieved in other small species typically associated to brackish-water environments, usually showing a more pronounced structuring and a lower connectivity even with smaller distances between populations. This difference is possibly due to the capability of A. fasciatus adults to sporadically cross marine environments (Torchio, 1967;Ç oker & Akyol, 2014) and thus move between coastal ponds. On the other hand, the comparison with the phylogeographic structure of the stickleback, a species that frequently uses marine environments for dispersal, and has been demonstrated to be able to perform transoceanic migrations, shows that migration events are nevertheless sporadic enough to ensure population differentiation.

Genetic diversity and conservation of A. fasciatus in Cyprus
Based on our molecular data, populations of A. fasciatus in Cyprus are genetically divergent from the other Mediterranean populations. The group of haplotypes from Cyprus is separated by five (CR) or two (COI) mutational steps from the closest extra-Cyprus haplotype, respectively (Figs. 2, 4). The Bayesian assignment analysis identified Cyprus haplotypes as separate genetic clusters for both mitochondrial regions, underlining the high degree of uniqueness of these populations. Moreover, the two mitochondrial regions were consistent in showing remarkable genetic divergence among the Cyprus populations. At the within-population level, the sample from Akrotiri marshes was characterised by the highest values of genetic diversity parameters; the population from Morphou Bay/Güzelyurt Korfezi showed the lowest values of genetic diversity, while the population from Glapsides exhibited intermediate values.
Comparing these values with literature data for CR, both Glapsides and Morphou have haplotype diversity values below the average, with only a few populations with values below those of Morphou. All Cyprus populations showed nucleotide diversity values below the average for CR. In the case of COI, however, the Akrotiri population exhibited the highest values for both haplotype and nucleotide diversity, although it should be noted that information for this gene is available for only five populations. The pattern of within-population genetic diversity in Cyprus mirrors the levels of environmental stress affecting biotopes at the three localities. Although the impact of drainages and the redirection of freshwater through the Akrotiri area might have fostered the success of the invasive mosquitofish Gambusia holbrooki, with negative impact on the distribution of A. fasciatus, Akrotiri marsh is a protected area. Therefore, the high within-population genetic diversity of Akrotiri population may be a consequence of the protection measures adopted. Moreover, this A. fasciatus population exhibited very high values of haplotype diversity, even if compared with other Mediterranean populations. The brackish pond in Glapsides Beach is located on the delta fan of the ephemeral Pedieos River, and together with the pond in Silver Beach represents the only suitable habitat for A. fasciatus in Famagusta Bay. The species, however, disappeared from Silver Beach after a mass mortality event in 2012 (Zogaris, 2017), possibly due to heavy habitat alterations related to urbanisation of the area (Seffer et al., 2011). Further samplings did not show a recolonisation of this area (C. Englezou, pers. obs.). The intermediate genetic diversity retrieved in this population is compatible with a moderately disturbed environment, where local mortality events can be the cause of population decline, with consequent genetic loss. The Morphou area was characterised by the presence of multiple environmental stressors, most of them of anthropogenic origin (Englezou et al., 2018), and unsurprisingly, this population showed the lowest values of genetic diversity, suggesting that it may have experienced relevant demographic decline(s) with subsequent genetic loss. Nevertheless, the highly significant value of Fu's F S might suggest a recent population expansion, but this is not confirmed by Harpending's raggedness index and Ramos-Onzins & Rozas' R 2 . A reduction of genetic diversity due to mass mortality events has already been reported in A. fasciatus (Ferrito & Tigano, 1996;Maltagliati, 2002), although several allegedly extinct populations showed clues of recovery (Lo Duca & Marrone, 2009;Valdesalici et al., 2015). However, genetic loss resulting from population decline might impair the ability of the population to withstand further environmental stress, thus increasing the likelihood of its extinction (Markert et al., 2010).
Although the three populations of A. fasciatus from Cyprus may be seen as a statistically supported single genetic cluster, each of them shows characteristics of uniqueness. In fact, the presence of locality-private haplotypes, together with results of the AMOVA, which identified the among-population component as the highest source of variation, suggests that gene flow among populations in Cyprus is restricted, or even absent, and that they have been isolated for long time. In fact, although Cyprus coastline hosts several habitats that might be potentially suitable for A. fasciatus, extensive sampling failed in finding additional populations, and the three known biotopes in which A. fasciatus is present are separated by at least 170 km of coastline (280 km in the case of Glapsides-Morphou distance).
From a strict conservationist point of view, our results suggest that each population of A. fasciatus in Cyprus should be treated as a single unit of conservation, in view of its genetic uniqueness, suggesting that recolonisation from one of the other localities should clearly be regarded as a highly unlikely event. Moreover, genetic data, together with results obtained in previous surveys by Zogaris (2017) and Englezou et al. (2018), highlighted that two out of three populations show traces of demographic and genetic deterioration, possibly related to the high degree of anthropogenic environmental stress. Although our results suggested that the Glapsides Beach population is in a genetically better condition than that of Morphou Bay, the local extinction in the close site of Silver Beach advises that its status should be carefully monitored. Recolonisation of the brackish environment of Silver Beach could occur naturally either through migration across Famagusta Bay, where A. fasciatus has already been observed in marine environments (Ç oker & Akyol, 2014), or through seasonal flooding of the Pedieos River delta. However, recent surveys showed that the brackish environment at Silver Beach is currently deteriorated, with a reduction of the area and an almost complete replacement of macrophytes, playing a relevant role in A. fasciatus biological cycle, by algal films (C. Englezou, pers. obs.). Therefore, although technically feasible, the reintroduction in the Silver Beach environment would most likely be inconsequential. The grim situation of the Morphou Bay population, affected by multiple anthropogenic stressors, might be partially relieved by the inclusion of this area in the proposed Natura 2000 site Akdeniz Special Environmental Protected Area (S.E.P.A.), as suggested by Englezou et al. (2018), or by the designation of a new, smaller protected area corresponding to biotopes where A. fasciatus occurs.
The genetic uniqueness of Cyprus populations and the low level of connectivity between them, together with the fragmentation of their habitat and the high number of anthropogenic and environmental threats affecting this species suggest that A. fasciatus deserves high-priority conservation interest in Cyprus. However, it should be noted that, according to the International Union for Conservation of Nature (IUCN) Red List categories and criteria, A. fasciatus is included in the ''Least Concern'' (LC) category (Kottelat & Freyhof, 2007). Discrepancies between the IUCN Red List and national red lists are not uncommon (Brito et al., 2010). Although the majority of such cases is due to missed taxonomic update or absence of data for a specific area, there are several instances where a species listed as LC by the IUCN can be considered relevant for conservation at the national level. This discrepancy has been observed mostly (1) for species at the edge of their natural range, where populations are by definition scarce and scattered (Brito et al., 2010;Helfman, 2013), (2) for exploited species that are object of strict management practices in order to ensure sustainability of their harvesting (Helfman, 2013), and (3) for species that in a specific geographical area are affected by a severe habitat degradation, that impairs survival of their populations. A case of discrepancy between IUCN and national red lists is that of another European vertebrate associated with wetlands, the spade-footed toad Pelobates fuscus (Laurenti, 1768), which is listed as LC by the IUCN (Dufresnes et al., 2019), but it is considered ''Endangered'' (EN) in Italy (Andreone et al., 2004;. The more conservative measures adopted in Italy were necessary due to the limited number of known populations, the relevant habitat reduction that occurred in the last two centuries and the presence of several threats to the survival of this species on the national territory. It should be stressed that the IUCN Red List and national red lists have two fundamentally different purposes, with the former list aiming at assessing the risk of a complete extinction of a species, focusing on the global pattern of threats, while national red lists are aimed at preserving the highest level of biodiversity within the national borders, focusing on specific geographical areas with their threats and processes. In our opinion, both approaches somehow neglect the biogeographical component and its effects on the genetic diversity of organisms. Commonly, the IUCN does not consider the intraspecific patterns of diversity, including the possibility of genetic loss across its geographical range, unless it is sanctioned by the taxonomy (e.g. endemic subspecies). National red lists tend to assign to the same category species that could sustain abundant populations but are affected by extrinsic processes, as well as species at the edge of their distributional range, where populations are by necessity scanty, again without taking into account genetic diversity patterns (Helfman, 2013). The case of A. fasciatus in Cyprus recalls that of P. fuscus in Italy from several points of view: (1) both species occur in populations that are fragmented and separated by large streaks of habitat made unsuitable by human activities; (2) both species are affected by multiple anthropogenic and environmental stressors across their national distribution; (3) national populations of both species show genetic uniqueness and distinctness from other populations, although the divergence is too low to justify a taxonomic recognition Litvinchuk et al., 2013;present data); (4) the conservation of both species strictly depends on the conservation of their habitats. The listing of these species in national red lists and the subsequent enforcement of protection measures would therefore represent a further step towards the protection of wetlands and the maintaining of their ecosystem services, where these taxa can be considered as umbrella species (Andreone et al., 2004;Valdesalici et al., 2015).
More generally, the naturally fragmented distribution of A. fasciatus, the low degree of gene flow between populations, and the subsequent possibility of small-scale genetic differentiation retrieved in some of them  suggest that the management of this species should take into account the conservation status of single populations, as local extinction might lead to significant diversity losses for the species. However, even though this approach is recommendable in Cyprus, it is not possible to strictly protect all A. fasciatus populations in the Mediterranean area. Therefore, further studies on the distribution and diversity of the less-known genetic lineages should be carried out, in order to design a conservation approach aimed at maximising protection of A. fasciatus genetic diversity.
Funding Open access funding provided by Università di Pisa within the CRUI-CARE Agreement. Sampling activities were funded by Freshwater Life Project (U.K. Charity No. 1172393), while laboratory work was funded by Fondi di Ateneo (FA) of the University of Pisa.
Data availability Voucher specimens have been deposited in the Museum of Natural History of the University of Pisa. Gene sequences have been made available in GenBank.
Code availability Not applicable.

Declarations
Conflict of interest The authors do not have any conflict of interest to declare.
Ethical approval All applicable international, national, and/ or institutional guidelines for the care and use of animals were followed by the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.