Genetic structure of the South European toothcarp Aphanius fasciatus (Actinopterygii: Cyprinodontidae) populations in the Mediterranean basin with a focus on the Venice lagoon

ABSTRACT The genetic structure of Aphanius fasciatus populations has been analysed using two mitochondrial DNA (mtDNA) markers (16S rRNA and D-loop) obtained from specimens collected in nine sites from the Venice lagoon, Comacchio saltworks and Corsica. Available GenBank sequences were also included, in order to extend the results on a Mediterranean scale. Genetic polymorphism within the Venice lagoon was very low, with most of the specimens analysed (66% for 16S rRNA and 83% for D-loop) sharing the same haplotype for either of the two markers. The genetic homogeneity found within the Venice lagoon may be the consequence of the northward migration of southern Adriatic populations after the Last Glacial Maximum: mismatch analysis showed indeed clear signs of a rapid demographic and spatial expansion. To explain this genetic homogeneity other hypotheses were also suggested, such as adaptation to the high variability of brackish water habitats, and artificial introductions. On a Mediterranean scale, phylogenetic analyses showed the presence of five distinct geographical lineages: Aegean Sea, Greek coast of the Ionian Sea, Adriatic Sea, Tyrrhenian Sea and Southern Sicily. Analysis of molecular variance revealed a genetic partitioning mainly due to differences between groups presumably due to late Miocene geological events, while less polymorphism was present within groups and populations.


Introduction
Species belonging to the large group of killifish (Cyprinodontiformes) are characterised by a small size and a short life cycle. These species have evolved special adaptations enabling them to survive in shallow fresh-to-hyperhaline waters, coping with marked spatio-temporal fluctuations of environmental parameters (Kroll 1984). They have been widely used in environmental monitoring programmes and in studies concerning the adaptations to environmental pressures that characterise these aquatic systems (Burnett et al. 2007). In the Mediterranean Sea, the South European toothcarp Aphanius fasciatus (Valenciennes, 1821) can be considered an estuarine resident (Franzoi et al. 2010), inhabiting mainly estuaries, lagoons, coastal ponds and saltworks. This species shows low dispersal capabilities, particularly due to the absence of a planktonic larval stage (Maltagliati 1999). The natural fragmentation of transitional water habitats and the resulting geographical isolation of A. fasciatus populations have made this species a good model for investigating interpopulation gene flow as well as species genetic structuring (Maltagliati 1998(Maltagliati , 1999Maltagliati et al. 2003;Triantafyllidis et al. 2007;Pappalardo et al. 2008;Annabi et al. 2012;Ferrito et al. 2013;Buj et al. 2015). Life-history traits and habitat constraints are known to interact in shaping the genetic structure of a species (Maltagliati 1998;Cognetti & Maltagliati 2000), particularly among brackish-water species.
The analysis of the genetic structure of A. fasciatus samples has revealed marked geographical differentiation among populations, while some aspects of its genetic structure have been used to detect environmental stress or specific adaptations with respect to certain parameters, such as temperature and salinity (Maltagliati 2002;Cimmaruta et al. 2003;Angeletti et al. 2010). Apart from the effect of the naturally fragmented distribution of the typical habitats of the species on its present-day structure, past geological events may also have contributed. The Messinian Salinity Crisis is well known to have shaped species composition and distribution of aquatic communities within the Mediterranean basin (Huyse et al. 2004;Carnevale et al. 2006). In particular for A. fasciatus, the findings of Hrbek and Meyer (2003) and Triantafyllidis et al. (2007) highlighted the critical role of the Messinian Crisis in the evolutionary history of the species, although subsequent gene flow seems to have erased traces of this event to a great extent. In a more recent past, sea level oscillations during glacial-interglacial periods have had an effect, perhaps minor, on the geographic distribution and the genetic structuring of aquatic species (Stefanni & Thorley 2003).
Even though the phylogeographic structuring of the species has been studied by many authors, considering different habitats along Mediterranean coasts (i.e. on a macro-geographic scale), information on small-scale spatial distribution is limited (Maltagliati et al. 2003).
For the Adriatic Sea, information was recently provided for the eastern part by Buj et al. (2015), while for the Italian coast only samples from the Lesina (Maltagliati 1999(Maltagliati , 2002Cimmaruta et al. 2003;Ferrito et al. 2013) and Grado (Ferrito et al. 2013) lagoons were analysed. Situated in the northern part of the Adriatic basin, the Venice lagoon is one of the biggest and most impacted transitional water systems of the Mediterranean, structured in a complex network of different habitats. Within the lagoon, previous surveys on the fish community have highlighted a close link of the species with the salt marsh habitats (Franco et al. 2006a, b). Small intertidal creeks and artificial ditches, within lagoon islands, host abundant populations of A. fasciatus showing significant differentiation in demographic and life-history traits, even over small geographical distances (Cavraro et al. 2013(Cavraro et al. , 2014. This is mainly due to geographical isolation and/or different environmental pressures among habitats (i.e. food availability and mortality rates). Nevertheless, no data are available on a possible genetic differentiation, in particular for nearby populations, as found for example in Sardinia by Maltagliati et al. (2003).
In this study, the genetic structure of A. fasciatus was investigated through the use of two mitochondrial DNA (mtDNA) markers. Until now, most of the studies on A. fasciatus genetics using mitochondrial markers have mainly focused separately either on the D-loop region (Tigano et al. 2004(Tigano et al. , 2006Rocco et al. 2007;Pappalardo et al. 2008;Ferrito et al. 2013;Buj et al. 2015) or on a portion of the gene encoding the 16S ribosomal subunit (Triantafyllidis et al. 2007). In the first case, populations from the Central Mediterranean were mainly analysed, while in the second Greek/Eastern Mediterranean populations were investigated. Only Annabi et al. (2013) considered both markers, but limited their analysis to 16 sequences from southeast Tunisia. In the present study, we use the two mtDNA fragments in order to (1) analyse the genetic structuring of A. fasciatus Northern Adriatic populations, focusing for the first time on the micro-geographic scale of the Venice lagoon; and (2) put this on a wider scale by analysing the obtained data in the framework of the full Mediterranean distribution of the species.

Materials and methods
Nine A. fasciatus samples were collected ( Figure 1): seven were from the lagoon of Venice, while the other two were from the saltworks of Comacchio (about 80 km south of Venice) and the lagoon of Diana (Corsica). Sampled fish were preserved in absolute ethanol and then transferred to the laboratory (University of Ioannina and University of Thessaloniki, Greece) for DNA extraction/ amplification and statistical analyses. For each station, seven to 10 specimens were processed. DNA was extracted with the CTAB (cetyltrimethylammonium bromide) method described by Hillis et al. (1996). Two segments of mtDNA were then amplified for each individual: the first (780 bp) comprised a fragment of the 16S rRNA, while the second (339 bp) was a section of the mitochondrial control region (D-loop). The conditions for amplification were as follows: for 16S, 30 cycles of amplification were performed with denaturation at 94°C for 35 s, annealing at 52°C for 35 s and extension at 72°C for 120 s; for D-loop, 35 cycles of amplification were used with denaturation at 94°C for 30 s, annealing at 51°C for 30 s and extension at 72°C for 40 s. For the 16S segment primers L3002 and H4280 were used (Hrbek & Meyer 2003), whereas for the D-loop segment the forward primer designed by Tigano et al. (2004) and the reverse primer designed by Meyer et al. (1990) were used. The amplified fragments were then sent out for sequencing at the VBC Biotech Service GmbH (Austria). The obtained sequences were aligned with Geneious 5.6 (http:// www.geneious.com/) and subsequently edited manually. The best substitution model for each region was chosen using the software jModelTest v. 2.1.4 (Guindon & Gascuel 2003;Darriba et al. 2012). Phylogenetic analyses were carried out in MEGA 5 (Tamura et al. 2011), using maximum likelihood (ML, using the nearest-neighbour-interchange heuristic method with complete deletion of gaps and missing data), neighbour joining (NJ, using the Tajima-Nei model with complete deletion of gaps and missing data) and maximum parsimony (MP, using the close-neighbour-interchange search method). The robustness of resulting trees was assessed by means of bootstrap analysis (100 pseudoreplicates for ML and MP, 1000 for NJ). The software Arlequin 3.5 (Excoffier & Lischer 2010) was used to calculate haplotype (h) and nucleotide (π) diversity of the different populations and to assess the degree of genetic structuring of the species through analysis of molecular variance (AMOVA, Excoffier et al. 1992), using 1000 permutations and partitioning the sequences into five geographical groups (see Results section for details), according to the results of phylogeny reconstruction. Genetic distances between groups were estimated in MEGA 5 using the Kimura 2-parameter model with gamma distributed rates. Furthermore, for the samples collected within the Venice lagoon, demographic population history, in terms of recent expansion or decline, was examined using a mismatch analysis coupled with the Fs statistics (Fu 1997). Network analysis (Network 4.6.1.1, Fluxus Technology Ltd) was used to investigate haplotype relationships on a Mediterranean scale using available sequences of A. fasciatus deposited in GenBank (see Appendix).

Venice lagoon
Specimens of A. fasciatus from the seven sites of Venice lagoon showed low haplotype and nucleotide diversity (Table I). Analyses based on the two different segments of mtDNA allowed the identification of 11 haplotypes for the 16S rRNA segment and seven for the D-loop. In total, 49 out of 63 amplified specimens for 16S rRNA and 54 out of 65 for D-loop shared the same haplotype (16-AD-12 and DL-AD-06, respectively). Except for a few cases, all other haplotypes were represented by single specimens for both markers (Table I). Phylogenetic reconstruction confirmed the genetic homogeneity of A. fasciatus within the Venice lagoon. Pairwise population differentiation tests (Raymond & Rousset 1995) showed no significant differences among sites (P > 0.05, after Bonferroni correction), with all samples grouping in the same cluster ( Figure 2). This is also evident in the network analysis for both markers (

Mediterranean Sea
Sequences from the Venice lagoon were analysed together with those of A. fasciatus sampled from Comacchio saltworks and Diana lagoon (Corsica) as well as with available sequences from GenBank. This led to a consistent reconstruction of the phylogeography of the species within the Mediterranean Sea. Overall, 49 haplotypes were found for 16S rRNA segment and 82 haplotypes for D-loop. The best-fit substitution models found were GTR+Γ for 16S (gamma-shape parameter α = 0.499; base frequencies: A = 0.249, C = 0.275, G = 0.178, T = 0.298) and HKY+ Γ for D-loop (gamma-shape parameter α = 0.252; base frequencies: A = 0.367, C = 0.184, G = 0.149, T = 0.300). Phylogenetic trees constructed from these haplotypes with the three different methods (ML, NJ, MP), using the closest available models with the software used, Table I Sicily. Although the distinction into groups appeared quite clear, some haplotypes from Southern Sicily (DL-SS-18 to 21) had an ambiguous placement. Furthermore, it should be pointed out that, for both markers, a single specimen from the saltworks  The phylogenetic reconstruction of A. fasciatus presented above was confirmed by network analysis (Figure 3), which underlined a greater distance for the Aegean group. This was confirmed by the estimates of genetic distances. For 16S rRNA, genetic distance of Aegean populations from the other groups ranged between 0.028 and 0.032, while distances among Adriatic, Tyrrhenian and Ionian populations were in the range of 0.012-0.016 (there are no sequences available for the Aegean Sea for D-loop). As pointed out for the phylogenetic trees, also in the network analysis Southern Sicily haplotypes showed an unclear position, with few haplotypes (SS-18 to SS-21) clustering with the Tyrrhenian Sea group. Thus, considering the genetic distances, Southern Sicily haplotypes were split into two groups. As in the case of the Aegean populations, the most divergent haplotypes from Southern Sicily (SS-01 to SS-17) showed a high genetic distance (0.039-0.050), nearly two-fold higher than that among the other groups (0.021-0.027).
Analysis of molecular variance, together with statistically significant fixation indices, showed a high degree of divergence of A. fasciatus populations (Table II), which were characterised by a low variance of genetic polymorphism (3.64% for 16S and 12.61% for the D-loop). Most of the genetic polymorphism derived from differences among geographical groups. In particular, for 16S rRNA this difference explained 91.47% of the total variance, while the percentage of variance dropped to 54.54% for the D-loop. Only in the latter fragment did the differences among populations within groups contribute significantly to the genetic polymorphism of the species (32.84% of variance for D-loop and 4.90% of variance for 16S).

Discussion
Data collected in the present study allowed us to obtain information from the northern part of the A. fasciatus distribution area along the Mediterranean coasts. First, the situation of Venice lagoon is discussed, followed by the contextualisation of the new sequences obtained within the phylogeography of the species on a Mediterranean scale, by a comparison with sequences available from the literature.

Venice lagoon
In agreement with the findings of previous studies (Triantafyllidis et al. 2007;Annabi et al. 2012;Ferrito et al. 2013), the analysis carried out on a small spatial scale, within the Venice lagoon, showed a substantial genetic homogeneity. This homogeneity could be interpreted as the result of a rapid spatial and demographic expansion that occurred during the Pleistocene. About 18,000 years ago, during the Last Glacial Maximum (Emiliani 1955), mean sea level was 160-180 m lower than at the present day (Van Straaten 1965). Regarding the Adriatic Sea, during that period the north-central portion emerged and the northern coast was about 300 km southward. When the sea level began to rise at the end of the glacial period, fish species living in the southern Adriatic would have moved northward as a consequence. Also, A. fasciatus would have begun to expand northwards, experiencing a rapid spatial and demographic expansion, as confirmed by the mismatch analysis of the Venice lagoon samples. This hypothesis would be in accordance with the findings of Buj et al. (2015), who analysed A. fasciatus populations along the eastern Adriatic coast, finding traces of recent demographic expansion in the new habitats available after the Last Glacial Maximum. A similar pattern has been also found for other Mediterranean fish species such as Atherina boyeri (Francisco et al. 2006;Milana et al. 2012), Pomatoschistus minutus (Stefanni & Thorley 2003) and Dicentrarchus labrax (Bahri-Sfar et al. 2000). This picture would be consistent with the low haplotype and nucleotide diversity (Grant & Bowen 1998) and the grouping into a single cluster of the current Northern Adriatic populations with those of Apulia (Ferrito et al. 2013).
Another possible explanation of the genetic homogeneity within lagoon waters may lie in the migration of individuals among geographically close habitats, contrary to the significantly restricted home range found for North American killifish within salt marsh systems (Lotrich 1975;Able et al. 2006). Such small-scale migrations could occur particularly in winter, when fish may move and concentrate in deeper and warmer waters within the lagoon, preventing further differentiation of haplotypes. Only Maltagliati et al. (2003) reported substantial genetic differences between two populations of A. fasciatus living in two adjacent coastal ponds, based on RAPD analysis (Random Amplification of Polymorphic DNA). This difference may arise from the different markers used or, alternatively, it may be the result of different processes, of both natural and anthropogenic origin. The future use of additional genetic markers, such as microsatellites, could provide further information about this point.
Within the Adriatic context, the most peculiar situation was that of a specimen sampled in the saltworks of Comacchio, which showed the same haplotype expressed by individuals collected within the Greek lagoon of Messolonghi, in the Ionian Sea. Even if Maltagliati (1999) excluded the translocation of A. fasciatus specimens among different sites, A. fasciatus may have been moved together with fish species of commercial importance for aquaculture purposes, as discussed previously. The same hypothesis may be advocated in order to explain the clustering of the north-eastern Sicily sample of Ganzirri within the Adriatic group: in this site, A. fasciatus was introduced at the end of the 20 th century . Results of the present study could also be discussed taking into account the hypothesis of human manipulation and translocation of A. fasciatus specimens in the Venice lagoon. Giandomenico Nardo, who described the genus Aphanius, in his work of 1847 stated that A. fasciatus was firstly recorded within the Venice lagoon at the beginning of the 19 th century, after the accidental introduction with juvenile mullets used in traditional fish farming, and quickly became one of the most abundant species in the most confined area of the lagoon. In support of the fact that A. fasciatus could not have been present within the Venice lagoon before the beginning of the 19 th century, there is the absence of this species from the exhaustive work on the North Adriatic fish fauna written by Stefano Chiereghin (1816) some decades before the work of Nardo. As discussed below, the haplotype found near the Po river delta (16S IO 06), originating from the Messolongi lagoon (Greece), could support the hypothesis of past translocations of A. fasciatus specimens along the Mediterranean coasts, even into the Venice lagoon. Further studies, also taking into account samples in historic collections and museum archives, are needed to clarify this point.

Mediterranean Sea
The phylogenetic analysis based on the two segments of mitochondrial DNA scored in this study made it possible to describe the genetic structuring of A. fasciatus on a Mediterranean scale. Although the two segments used exhibit different mutation rates, with 16S rRNA being more conserved than D-loop, the results obtained were consistent between the two markers and with the existing literature (Tigano et al. 2004(Tigano et al. , 2006Rocco et al. 2007;Triantafyllidis et al. 2007;Pappalardo et al. 2008;Annabi et al. 2012Annabi et al. , 2013Ferrito et al. 2013;Buj et al. 2015). Only Hrbek and Meyer (2003) maintained that signs of genetic structuring of A. fasciatus across the Mediterranean Sea would have been diluted by subsequent gene flow. On the contrary, results of the present study, through the analysis of a larger number of samples and the comparison of two different markers, allowed the identification of a clear geographical structuring of the species in the Mediterranean. The analysis of molecular variance identified five distinct geographical groups, formed by quite homogeneous populations. The F st indices calculated for both markers were in accordance with those found for other species of killifish from both  (Ashbaugh et al. 1994;Doadrio et al. 1996;Dunham & Minckley 1998). These results are consistent not only with the biogeographical history of A. fasciatus, but also with its ecological characteristics. The close link with transitional environments and the limited dispersal potential, even during the early stages of the life cycle, determine a high level of isolation of populations. However, the phylogeographic structure resulting from the analysis of mtDNA should be mainly influenced by past events on a geological scale (Milana et al. 2012). As proposed by Hrbek and Meyer (2003) and Reichenbacher and Kowalke (2009), this species originated about 5 million years ago, during the Messinian Salinity Crisis of the late Miocene. In that period, the partial draining of the Mediterranean basin would have promoted an east-west genetic differentiation of Aphanius populations. Within the five-group configuration observed in this study, the high genetic distances of Aegean and South Sicily populations with respect to the other groups could reflect a taxonomic difference. As mentioned above, the differentiation of Aegean populations could be dated to the late Miocene, after the Messinian Salinity Crisis (Hsü et al. 1977;Triantafyllidis et al. 2007). Kottelat et al. (2007) described in the Aegean Sea, on the basis of phenotypic traits, a new species for the genus Aphanius very similar to A. fasciatus, Aphanius almiriensis. At present, its description is not confirmed by targeted genetic analyses, but sequences published in Triantafyllidis et al. (2007) may belong to this recently described species. In accordance with the hypothesis formulated in Ferrito et al. (2013), the presence of a new (sub-)species could explain also the high genetic differentiation of South Siciliy samples, also reported by Tigano et al. (2006) and Rocco et al. (2007).
In conclusion, Aphanius fasciatus within the Venice lagoon show genetic homogeneity, despite the geographic complexity and the wide extension of this site. The signs of a recent spatial and demographic expansion, ascribed to a post-glacial northward spreading of South Adriatic populations or to a more recent founder effect, has been discussed, taking into account the ecology and life history of the species. On the contrary, at a Mediterranean scale, a well-defined genetic structuring of the species has been confirmed. However, further analyses with additional molecular markers are needed in order to resolve some of the undefined nodes, particularly those related to the phylogenetic position of South