A total of 112 samples of A. fasciatus collected from eight localities along the eastern Adriatic coast were included in the study (Fig. 1, Table 1). Only one locality, Pantan, is an inshore brackish lagoon, whereas the remaining locations are hypersaline water bodies, entirely or partially altered into salt pans. Each population/locality was represented by 11-15 samples. Sequence sets for phylogenetic reconstructions were supplemented with available GenBank sequences (Supplementary material Table S1).

DNA extraction was performed using a standard extraction kit (DNeasy tissue kit, Qiagen). Two mitochondrial genetic markers with different mutation rates were analysed: the gene for cytochrome b (cyt b) and mitochondrial control region (mtDNA CR). The reaction volume for polymerase chain reaction (PCR) amplifications contained 25 µl AmpliTaq Gold Master Mix (Applied Biosystems Inc.), 2 µl of each primer (10 pmol), 4 µl template DNA and 17 µl water. The primers used for the amplification of the cyt b gene were H15915 and L14724 (Schmidt and Gold 1993), whereas the primers A and E from Lee et al. (1995) were used for the mtDNA CR. The thermal profile for PCR amplification of cyt b was 10 min at 95°C, 35 cycles of 60 s at 92°C, 90 s at 48°C, 180 s at 72°C and 4 min at 72°C. PCR conditions for the amplification of mtDNA CR were 10 min at 95°C, 35 cycles of 45 s at 94°C, 45 s at 48°C, 60 s at 72°C and 7 min at 72°C. Sequencing with both forward and reverse primers was carried out by the Macrogen Service Center using the same primers as for the PCR amplifications.

Homologous regions of cyt b and mtDNA CR were aligned manually against previously published sequences using BioEdit Sequence Alignment Editor (Hall 1999). Chromatograms and alignments were checked visually.

The measures of DNA polymorphism obtained to describe the intrapopulational genetic diversity of the eastern Adriatic populations were haplotype diversity (Hd) and nucleotide diversity (Π). The level of genetic differentiation between populations was assumed by calculating statistics based on haplotypes (HST) and the one based on nucleotide sequences (KST). The decision to accept or reject the null hypothesis (there is no genetic differentiation between populations) was made based on the permutation test. Calculations were made using DNAsp software (Librado and Rozas 2009) and were conducted on both data sets (cyt b and mtDNA CR)

In order to obtain insight into the intraspecific and population genetic structure, cyt b and mtDNA CR sequences from the eastern Adriatic populations were analysed separately using the statistical parsimony method (Templeton et al. 1992) under a 95% parsimony connection limit, using TCS software (version 1.21; Clement et al. 2000). With the aim of identifying most likely evolutionary processes responsible for recent phylogeographic structure, nested clade analyses (NCAs; Templeton 1998) were conducted on statistical parsimony networks based on both genetic markers. NCAs were performed using ANeCa software (Panchal 2007), which includes both TCS v1.18 (Clement et al. 2000) and GeoDis v2.2 (Posada et al. 2000), thus automating the nesting algorithm and the inference key that interprets the population structure and history at each nesting level.

To ascertain whether the observed pattern of genetic diversity is in accordance with the neutral evolution model, tests of neutrality were performed: Fu and Li’s tests, and Tajima’s test as implemented in DNAsp ver. 5.10.01. Since significant departures from the null hypothesis suggest either selective pressure on the locus under study or changes in population size (Ramos-Onsins and Rozas 2002), the results of the neutrality tests were used to detect traces of population growth. In particular, Fu’s FS statistic is expected to have negative values in expanded populations and it is considered one of the strongest tests in detecting population enlargements (Ramos-Onsins and Rozas 2002). For populations where genetic structures expressed indications of population expansion, their onsets were calculated on the basis of the peaks of theoretical mismatch curves, which were estimated using the same software. The mismatch distribution (the distribution of the number of pairwise mutation differences between sequences) is expected to be unimodal in recently expanded populations and irregular in populations of constant size (Rogers and Harpending 1992). The time of the population growth measured in the units of mutational time is expressed as the parameter τ (τ = 2ut; u is the mutation rate per sequence, t is the absolute time) (Rogers and Harpending 1992). Thereafter, the time of population expansion was calculated from the parameter τ.

Phylogeographic analyses were conducted on two data sets. The cyt b haplotypes of A. fasciatus obtained in this study were pooled with cyt b sequences of several other Aphanius species (A. sophiae [Heckel, 1847], A. mesopotamicus Coad, 2009, A. farsicus Teimori, Esmaeili and Reichenbacher, 2011, A. isfahanensis Hrbek, Keivany and Coad, 2006, A. vladykovi Coad, 1988, A. iberus [Valenciennes, 1846], A. baeticus Doadrio, Carmona and Fernandez-Delgado, 2002, A. saourensis Blanco, Hrbek and Doadrio, 2006 and A. dispar [Rüppell, 1829]) retrieved from GenBank (S1) to assemble the first data set. The purpose of phylogenetic reconstruction of this data set was to ascertain the position of eastern Adriatic Aphanius populations and their relationships with other Aphanius species, since they were not included in any of the previous analyses. The phylogenetic reconstruction was conducted using two methods of phylogenetic inference: maximum parsimony (MP), as implemented in PAUP (version 4.0b10; Swofford 2002), and Bayesian inference (BAY), implemented in MrBayes (version 3.1.2, Ronquist and Huelsenbeck 2003). For MP analyses, the heuristic search mode with 100 replicates was used, with randomized input orders of taxa, and tree bisection-reconnection branch swapping with all codon sites and nucleotide substitutions types weighted equally. Nonparametric bootstrapping (100 pseudoreplicates, 10 addition-sequence replicates) was used to assess the branch support. Each BAY analysis consisted of two simultaneous runs. For each, four Markov Chain Monte Carlo chains were run for three million generations with trees sampled every 100 generations. The first 20% of the sampled trees were discarded and Bayesian posterior probabilities were estimated from the 50% majority-rule consensus tree of the retained trees. The sequence of Poeciliopsis monacha Miller, 1960 (GenBank accession number AF047346), belonging to the same order but a different family than Aphanius, was used as an outgroup.

The second data set was composed of the mitochondrial control region sequences obtained in this study, and mtDNA CR haplotypes available from the GenBank (S1, localities marked in Fig. 1). This data set contained haplotypes of A. fasciatus only; however, more localities were included to enable more detailed analysis of the intraspecific phylogeographic structure. Since the mitochondrial control region accumulates mutations up to ten times faster than the rest of the mitochondrial genome, it is considered suitable for investigations of rapidly evolving populations (Meyer 1993) and is therefore the genetic marker of choice for investigating intraspecific relationships. In order to visualize genealogical relationships among closely related haplotypes (belonging to the same species) and to enable the possibility of presence of ancestral haplotypes in the recent genetic structure, a phylogenetic network was constructed. Network construction was based on the median-joining algorithm using Network 4.5.1.6. software (Fluxus Technology Ltd.). In addition to the sequences of A. fasciatus, the sequence of A. dispar (GenBank accession number UO6052) was included to enable identification of ancestral haplotypes.

The molecular clock calibration, namely the estimation of divergence per pairwise comparisons per million years for both genetic markers, was conducted using the Bayesian MCMC coalescent method implemented in BEAST v1.7.0. (Drummond et al. 2012). Since only one evolutionary event has been proposed with a time frame (separation of A. baeticus from A. iberus was associated with the opening of the Strait of Gibraltar 5.5 million years ago; Perdices et al. 2001), this was used to calibrate the phylogenetic tree based on cyt b haplotypes. A relaxed clock model was employed (since general clock-like behaviour was rejected) with branch rates drawn from an uncorrelated lognormal distribution and a Yule speciation prior. The time of divergence of A. fasciatus estimated from this analysis was used to infer the divergence rate of A. fasciatus mtDNA CR, in addition to the number of mutations per million years. This was subsequently used to calculate divergence times of different populations using the Time estimates sub-program of the Network 4.5.1.6. software.

The two genetic markers showed similar results regarding intrapopulational genetic polymorphism and revealed differences among the populations (Table 1). Based on the cyt b marker, the range of the number of haplotypes detected per population was 2-6, the haplotype diversity 0.133-0.571, and nucleotide diversity 0.00012-0.00129, whereas for mtDNA CR the range of haplotype number per population was 1-7, for haplotype diversity 0-0.773 and for nucleotide diversity 0-0.00294. Both analysed genetic markers indicated the highest measures of genetic polymorphism in the Pag, Nin and Narta populations, while populations from Pantan and Ulcinj displayed the smallest amount of genetic diversity.

Both statistics calculated to infer genetic divergence among populations (HST, KST; Table 2) revealed the same result: the Adriatic populations form two genetically distinct groups. The Northern group consists of the populations from Sečovlje (Slovenia), Pag, Dinjiška, Nin and Pantan (Croatia). The Southern group includes the populations from Ston (Croatia), Ulcinj (Montenegro) and Narta (Albania). Based on control region sequences, no significant genetic divergence was detected among the populations belonging to the same group, whereas populations from different groups were genetically different. However, the cyt b sequences showed genetic similarity for the majority of populations from the same group, with a single exception in each group (Dinjiška and Narta). When the two groups were compared, genetic polymorphism was higher in the Southern group.

The demographic history of the two population groups from the eastern Adriatic coast is based on different events (Table 3), as inferred by the NCA based on cyt b sequences (Fig. 2b). The NCA performed on mtDNA CR haplotypes (Fig. 2a) implied that allopatric fragmentation was responsible for the observed structure of the whole cladogram, while no conclusive outcomes could be reached for any of the internal clades.

Fu’s F_S and Tajima’s D (Table 2) had significantly negative values for both genetic markers only for the Narta population, implying its demographic expansion and also suggesting very little overall deviation from the mutation-drift equilibrium. Nevertheless, traces of population growth were also noted for the Nin population (Fu’s F_S statistic was significantly negative for both markers and Tajima’s D was significantly negative for cyt b sequences) and probably for the Pag population (Fu’s F_S was significantly negative for mtDNA CR and Tajima’s D for cyt b sequences) Table 2). Based on the parameter τ, the time of the onset of population growth dates back to the Pleistocene (324000-344000 years ago (YA) for the Narta population; 122000-444000 YA for the Nin population based on cyt b and control region sequences, respectively; and 177000 YA for the Pag population).

Both phylogenetic reconstruction methods corroborated the monophyletic position of A. fasciatus and the inclusion of samples from the eastern Adriatic coast into the A. fasciatus phylogenetic lineage (Fig. 3). This is a sister lineage to the lineage comprising A. sophiae, A. mesopotamicus, A. farsicus, A. isfahanensis and A. vladykovi. The third phylogenetic group is composed of A. iberus, A. baeticus and A. saourensis, whereas A. dispar is more distantly related to the remaining investigated Aphanius species.

The median-joining network comprising 45g the newly obtained mtDNA CR haplotypes and haplotypes from several localities in Italy, Malta and Tunisia (Fig. 4) shows the interesting intraspecific genetic configuration of A. fasciatus. Structuring correlated with geographic position is clearly visible, with haplotypes belonging to the same population almost always clustering together. However, populations with haplotypes separated by a reduced number of mutational steps are not necessarily geographically close and vice versa. Separation of the eastern Adriatic populations into two evolutionary groups (Northern and Southern) is clearly visible in the median-joining network (Fig. 4) and in the statistical parsimony networks (Fig. 2).

Based on molecular clock calibration, the evolutionary rate of cyt b was estimated as 1.2% per million years, whereas for the mitochondrial control region, the estimate was 3.6%. The divergence times of different A. fasciatus populations are marked on the median-joining network (Fig. 4). The last common ancestor of northern and southern Adriatic populations is dated back to around 0.58 million years ago.

Based on the results obtained, it can be concluded that the populations from Pag and Nin within the Northern group and the Ston and Narta within the Southern group are the most important populations for the further survival of Aphanius fasciatus in the Adriatic Basin. Although Pantan is under a special protection regime precisely because of the toothcarp population, it seems that salt pans have greater significance for the survival of this species, though inshore salt marshes are very important for overall biodiversity. Though those localities (especially Pag, Nin, Ston and Narta) are under substantial human impact, they neverthless support populations of Aphanius fasciatus with high genetic diversity and that seem to be in expansion. For the existence of a species, maintenance of the highest genetic diversity is extremely important because it will allow adaptation to possible environmental changes. Therefore, populations from Pag, Nin, Ston and Narta, in addition to the Pantan population, should receive special conservation measures. Salt pans, as artificial anthropogenic formations and extreme habitats, are unfavourable places for most fish species. On the other hand, reduced competitive pressure makes them a convenient place for the toothcarp. Despite their recorded lower genetic diversity, natural, brackish wetland habitats are also very important for the conservation of populations of this species.