Population genetic structure and demographic history of Aphanius fasciatus (Cyprinodontidae: Cyprinodontiformes) from hypersaline habitats in the eastern Adriatic

Ivana Buj 1, Jure Miočić-Stošić 1, Zoran Marčić 1, Perica Mustafić 1, Davor Zanella 1, Milorad Mrakovčić 1, Tanja Mihinjač 1, Marko Ćaleta 2

1 Department of Zoology, Faculty of Science, University of Zagreb, Rooseveltov trg 6, 10000 Zagreb, Croatia.
2 Faculty of Teacher Education, University of Zagreb, Savska cesta 77, Zagreb, Croatia. E-mail: marko.caleta@ufzg.hr

Summary: In order to investigate the phylogeography and population genetic structure of the South European toothcarp (Aphanius fasciatus), we analysed gene sequences of two mitochondrial markers (cytochrome b and mtDNA control region) in samples from eight localities along the eastern Adriatic coast and combined them with sequences from other Mediterranean localities. Since the South European toothcarp primarily inhabits hypersaline water bodies, it is a good model species for understanding patterns of colonization and dispersal of species adapted to variable conditions. The eastern Adriatic populations are separated into two groups of genetically related populations. The Northern group contains the populations from Sečovlje, Pag, Dinjiška, Nin and Pantan, whereas the Southern group contains the populations from Ston, Ulcinj and Narta. The majority of divergence events date back to the Pleistocene epoch and it is likely that sea level changes during glacial cycles played a significant role in shaping the recent genetic structure of this species. Our results imply pronounced intraspecific structuring of this species, whereas great environmental variations resulted in a smaller intrapopulational genetic diversity of A. fasciatus than seen in other Mediterranean fishes.

Keywords: colonization; conservation priorities; divergence; evolution; genetic diversity; South European toothcarp.

Estructura genética poblacional e historia demográfica de Aphanius fasciatus (Cyprinodontidae: Cyprinodontiformes) en hábitats hipersalinos del Adriático oriental

Resumen: Con el objetivo de analizar la estructura genética y filogeografía de las poblaciones del fartet oriental (Aphanius fasciatus), se analizaron secuencias de genes de dos regiones mitocondriales (citocromo b y región de control del ADNmt) en individuos de ocho localidades a lo largo de la costa oriental del Adriático y de otras localidades mediterráneas. El fartet oriental habita principalmente masas de agua hipersalinas, y es una buena especie modelo para el estudio de los patrones de colonización y dispersión de especies adaptadas a condiciones variables. Los resultados de nuestro estudio muestran que las poblaciones del Adriático oriental se separan en dos grupos genéticamente relacionados. El grupo del Norte contiene las poblaciones de Sečovlje, Pag, Dinjiška, Nin y Pantan, mientras que el grupo del Sur contiene las poblaciones de Ston, Ulcinj y Narta. La mayoría de los eventos de divergencia pudiera ser que se remontaran a la época del Pleistoceno, y es probable que los cambios del nivel del mar durante los ciclos glaciales jugasen un papel importante en la estructura genética reciente de esta especie. Nuestros resultados revelan diferencias intraespecíficas pronunciadas en el fartet oriental, y por otro lado menor diversidad genética intrapoblacional como resultado de grandes variaciones ambientales en comparación con otros estudios de peces del Mediterráneo.

Palabras clave: colonización; prioridades de conservación; divergencia; evolución; diversidad genética; fartet oriental.

Citation/Como citar este artículo: Buj I., Miočić-Stošić J., Marčić Z., Mustafić P., Zanella D., Mrakovčić M., Mihinjač T., Ćaleta M. 2015. Population genetic structure and demographic history of Aphanius fasciatus (Cyprinodontidae: Cyprinodontiformes) from hypersaline habitats in the eastern Adriatic. Sci. Mar. 79(4): 399-408. doi: http://dx.doi.org/10.3989/scimar.04198.06A

Editor: M. Pascual.

Received: December 19, 2014. Accepted: July 13, 2015. Published: October 7, 2015.

Copyright: © 2015 CSIC. This is an open-access article distributed under the Creative Commons Attribution-Non Commercial Lisence (by-nc) Spain 3.0.

Contents

Summary
Resumen
Introduction
Materials and methods
Data analyses
Results
Discussion
Acknowledgements
References

INTRODUCTIONTop

The fish community of the Adriatic Sea consists of some 440 species (Lipej and Dulčić 2010Lipej L., Dulčić J. 2010. Checklist of the Adriatic Sea Fishes. Zootaxa 2589: 1-92.). Its diversity and structure were influenced by the specificities of the Adriatic Sea (Maggio et al. 2009Maggio T., Lo Brutto S., Garoia F., et al. 2009. Microsatellite analysis of red mullet Mullus barbatus (Perciformes, Mullidae) reveals the isolation of the Adriatic Basin in the Mediterranean Sea. ICES J. Mar. Sci. 66: 1883-1891.), as a semi-closed area with its own characteristics in terms of salinity, temperature and depth (Astraldi et al. 1999Astraldi M., Balopoulos S., Candela J., et al. 1999. The role of straits and channels in understanding the characteristics of Mediterranean circulation. Prog. Oceanogr. 44: 65-108.), but also by evolutionary events that have shaped the configuration of the Mediterranean ichthyofauna. Though the Adriatic fish fauna is usually considered to be well studied (Dulčić et al. 2005Dulčić J., Soldo A., Jardas I. 2005. Adriatic fish biodiversity and review of bibliography related to Croatian small-scale coastal fisheries. AdriaMed Technical Documents nº15: 103-125.), information on intraspecific structure and diversity is scarce and insufficient. This is particularly true for species inhabiting specific habitats, such as hypersaline water bodies. Investigations of demographic history and population genetic structure of species adapted to extreme conditions can help to reveal all components of Mediterranean biodiversity, and also to understand the peculiarities of evolution in specific environments.

A suitable model species for investigating evolution in such habitats is the South European toothcarp, Aphanius fasciatus (Valenciennes, 1821), a cyprinodont distributed throughout the central and eastern coastal zone of the Mediterranean. It is especially interesting since it is endemic to the Mediterranean, but inhabits habitats different from those of typical marine species.

Aphanius fasciatus usually inhabits brackish environments, such as coastal ponds and lagoons, though there are records of inland water populations (Parenti and Tigano 1993Parenti L.R., Tigano C. 1993. Polymorphic skeletal characters in Aphanius fasciatus Teleostei: Cyprinodontiformes. Copeia 4: 1132-1137.). Along the eastern Adriatic coast, it is found almost exclusively in hypersaline water bodies. It is able to tolerate a wide range of physicochemical parameters, in particular temperature (5-39°C) and salinity (0-180 psu) (Triantafyllidis et al. 2007Triantafyllidis A., Leonardos I., Bista I., et al. 2007. Phylogeography and genetic structure of the Mediterranean killifish Aphanius fasciatus (Cyprinodontidae). Mar. Biol. 152: 1159-1167.). Due to its life history traits, such as large adhesive eggs, short generation time, high reproductive rate and, especially, the absence of larval dispersal stages, this species is considered to be an ideal organism for investigating the effects of genetic drift, selection and gene flow on population divergence (Triantafyllidis et al. 2007Triantafyllidis A., Leonardos I., Bista I., et al. 2007. Phylogeography and genetic structure of the Mediterranean killifish Aphanius fasciatus (Cyprinodontidae). Mar. Biol. 152: 1159-1167.). Furthermore, A. fasciatus has no commercial value, so the distribution of genotypes among populations has likely not been affected by human translocations.

The genetic diversity of several populations of A. fasciatus has been investigated, though with opposing results. Hrbek and Meyer (2003)Hrbek T., Meyer A. 2003. Closing of the Tethys Sea and the phylogeny of Eurasian killifishes (Cyprinodontiformes: Cyprinodontidae). J. Evol. Biol. 16: 17-36. found low genetic variation among A. fasciatus populations, despite its wide distribution and potential subjection to Mediterranean vicariant events, whereas other authors documented marked genetic structuring among populations of this species (Rocco et al. 2007Rocco L., Ferrito V., Costagliola D., et al. 2007. Genetic divergence among and within four Italian populations of Aphanius fasciatus (Teleostei, Cyprinodontiformes). Ital. J. Zool. 74: 371-379., Triantafyllidis et al. 2007Triantafyllidis A., Leonardos I., Bista I., et al. 2007. Phylogeography and genetic structure of the Mediterranean killifish Aphanius fasciatus (Cyprinodontidae). Mar. Biol. 152: 1159-1167.). Hrbek and Meyer (2003)Hrbek T., Meyer A. 2003. Closing of the Tethys Sea and the phylogeny of Eurasian killifishes (Cyprinodontiformes: Cyprinodontidae). J. Evol. Biol. 16: 17-36. explained this low genetic divergence of A. fasciatus by occasional migration events that prevented differentiation. Triantafyllidis et al. (2007)Triantafyllidis A., Leonardos I., Bista I., et al. 2007. Phylogeography and genetic structure of the Mediterranean killifish Aphanius fasciatus (Cyprinodontidae). Mar. Biol. 152: 1159-1167. analysed samples from 13 localities in Greece and Turkey and concluded that the Messinian Salinity Crisis was the most important event in shaping their recent population genetic structure. Bardakci et al. (2004)Bardakci F., Tatar N., Hrbek T. 2004. Genetic relationships between Anatolian species and subspecies of Aphanius based on RAPD markers. Biologia (Bratislava) 59: 559-566. investigated Anatolian Aphanius populations and revealed limited gene flow among populations. They assumed extreme habitat conditions to have caused occasional drastic size reductions of Aphanius populations that reduced intrapopulational diversity and increased interpopulational differentiation due to genetic drift. Furthermore, they stated that the Quaternary paleoclimate changes had modified the distribution range and intraspecific genetic structure of Aphanius species. The most recent study (Ferrito et al. 2013Ferrito V., Pappalardo A.M., Canapa A., et al. 2013. Mitochondrial phylogeography of the killifish Aphanius fasciatus (Teleostei, cyprinodontidae) reveals highly divergent Mediterranean populations. Mar. Biol. 160: 3193-3208.) found strong population genetic structuring within A. fasciatus that reflects the historical geography of the Mediterranean basin from the Late Miocene to the Pliocene. They recognized the presence of multiple glacial refugia and interglacial phylogeographic breaks for A. fasciatus.

The aims of this study were to investigate the genetic structure and demographic history of several A. fasciatus populations mostly inhabiting hypersaline water bodies along the eastern Adriatic coast, and to gain insight into the evolutionary history of this species. These populations have not previously been analysed, and therefore we aimed to ascertain their phylogenetic relationships, degree of intrapopulational diversity, as well as interpopulational differentiation and possible gene flow. Furthermore, we intended to verify the evolutionary hypotheses proposed in previous papers and to estimate time frame of diversification of the Adriatic Aphanius populations. Insight into the past history, together with a better understanding of the present status of a species, enables more precise prediction of its future alterations and is thus of great importance for conservation purposes.

MATERIALS AND METHODSTop

Sample collection

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).

sm4198fig1.jpg

Full size image

Fig. 1. – Map of the investigated area with marked sampling localities. The locations on the eastern Adriatic coast that were the main subject of the present study are marked in capital letters. From the remaining localities included in the map, mtDNA CR sequences were incorporated in the second data set and analysed to reveal the intraspecific phylogeographic structure of A. fasciatus.

Table 1. – Intrapopulational genetic diversity and tests of selective neutrality for east Adriatic Aphanius populations, based on two molecular markers. N, number of sequences; h, number of haplotypes; S, number of polymorphic sites; Hd, haplotype diversity; Π, nucleotide diversity; asterisks mark statistically significant Fu’s FS and Tajima’s D values; NG and SG explain populations belonging to the Northern and Southern group, respectively.

Population Genetic marker N h S Hd Π Fu’s Fs Tajima’s D
Sečovlje (NG) cyt b 11 2 1 0.182 0.00016 -0.410 -1.129
mtDNA CR 10 2 1 0.200 0.00051 -0.339 -1.112
Pag (NG) cyt b 15 6 11 0.571 0.00129 -1.672 -2.180*
mtDNA CR 12 4 3 0.455 0.00161 -1.590* -1.179
Dinjiška (NG) cyt b 14 4 3 0.571 0.00083 -0.544 0.006
mtDNA CR 12 2 1 0.167 0.00042 -0.476 -1.141
Nin (NG) cyt b 15 5 4 0.476 0.00047 -3.228* -1.816*
mtDNA CR 12 4 4 0.455 0.00168 -1.489* -1.747
Pantan (NG) cyt b 15 2 1 0.133 0.00012 -0.649 -1.155
mtDNA CR 11 2 1 0.182 0.00046 -0.410 -1.129
Ston (SG) cyt b 12 4 3 0.455 0.00044 -2.124* -1.629
mtDNA CR 11 3 4 0.345 0.00183 -0.019 -1.712
Ulcinj (SG) cyt b 15 2 1 0.133 0.00012 -0.649 -1.155
mtDNA CR 11 1 0 0 0
Narta (SG) cyt b 15 5 6 0.476 0.0007 -2.103* -1.983*
mtDNA CR 12 7 7 0.773 0.00294 -4.446* -1.944*
Northern group cyt b 70 13 17 0.551 0.00073 -10.425* -1.270*
mtDNA CR 57 10 8 0.293 0.00097 -11.605* -2.078*
Southern group cyt b 42 10 11 0.642 0.00083 -6.201* -1.890*
mtDNA CR 34 9 11 0.421 0.00163 -7.464* -2.38*

DNA extraction and PCR amplification

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 1993Schmidt T.R., Gold, J.R. 1993. Complete sequence of the mitochondrial cytochrome b gene in the cherryfin shiner, Lytrurus roseipinnis (Teleostei: Cyprinidae). Copeia 3: 880-883.), whereas the primers A and E from Lee et al. (1995)Lee W.-J., Conroy J., Howell W.H., et al. 1995. Structure and evolution of teleost mitochondrial control regions. J. Mol. Evol. 41: 54-66. 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.

DATA ANALYSESTop

Sequence alignment

Homologous regions of cyt b and mtDNA CR were aligned manually against previously published sequences using BioEdit Sequence Alignment Editor (Hall 1999Hall T.A. 1999. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucl. Acids. Symp. Ser. 41: 95-98.). Chromatograms and alignments were checked visually.

Analyses of population genetic polymorphism and diversity

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 2009Librado P., Rozas J. 2009. DNAsp v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451-1452.) and were conducted on both data sets (cyt b and mtDNA CR)

Genetic structure and historical demography analyses

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. 1992Templeton A.R., Crandall K.A., Sing C.F. 1992. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics 132: 619-633.) under a 95% parsimony connection limit, using TCS software (version 1.21; Clement et al. 2000Clement M., Posada D., Crandall K.A. 2000. TCS: a computer program to estimate gene genealogies. Mol. Ecol. 9: 1657-1660.). With the aim of identifying most likely evolutionary processes responsible for recent phylogeographic structure, nested clade analyses (NCAs; Templeton 1998Templeton A.R. 1998. Nested clade analyses of phylogeographic data: testing hypotheses about gene flow and population history. Mol. Ecol. 7: 381-397.) were conducted on statistical parsimony networks based on both genetic markers. NCAs were performed using ANeCa software (Panchal 2007Panchal M. 2007. The automation of nested clade phylogeographic analysis. Bioinformatics 23: 509-510.), which includes both TCS v1.18 (Clement et al. 2000Clement M., Posada D., Crandall K.A. 2000. TCS: a computer program to estimate gene genealogies. Mol. Ecol. 9: 1657-1660.) and GeoDis v2.2 (Posada et al. 2000Posada D., Crandall K.A., Templeton A.R. 2000. GeoDis: a program for the cladistic nested analysis of the geographical distribution of genetic haplotypes. Mol. Ecol. 9: 487-488.), thus automating the nesting algorithm and the inference key that interprets the population structure and history at each nesting level.

Neutrality tests and mismatch distribution

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 2002Ramos-Onsins S.E., Rozas J. 2002. Statistical properties of new neutrality tests against population growth. Mol. Biol. Evol. 19: 2092-2100.), 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 2002Ramos-Onsins S.E., Rozas J. 2002. Statistical properties of new neutrality tests against population growth. Mol. Biol. Evol. 19: 2092-2100.). 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 1992Rogers A.R., Harpending H. 1992. Population growth makes waves in the distribution of pairwise genetic differences. Mol. Biol. Evol. 9: 552-569.). 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 1992Rogers A.R., Harpending H. 1992. Population growth makes waves in the distribution of pairwise genetic differences. Mol. Biol. Evol. 9: 552-569.). Thereafter, the time of population expansion was calculated from the parameter τ.

Phylogeographic analyses

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 2002Swofford D.L. 2002. PAUP*: Phylogenetic Analysis Using Parsimony (*and Other Methods), Version 4. Sinauer Associates.), and Bayesian inference (BAY), implemented in MrBayes (version 3.1.2, Ronquist and Huelsenbeck 2003Ronquist F., Huelsenbeck J.P. 2003. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19: 1572-1574.). 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 1993Meyer A. 1993. Evolution of mitochondrial DNA in fishes. In: Hochachka P.W., Mommsen T.P. (eds), Biochemistry and Molecular Biology of Fishes, vol 2. Elsevier Science Publisher, Amsterdam, pp 1-38.) 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.

Divergence time estimation

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. 2012Drummond A.J., Suchard M.A., Xie D., et al. 2012. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29: 751-761.). 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. 2001Perdices A., Carmona J.A., Fernández-Delgado C., et al. 2001. Nuclear and mitochondrial data reveal high genetic divergence among Atlantic and Mediterranean populations of the Iberian killifish Aphanius iberus (Teleostei: Cyprinodontidae). Heredity 87: 314-324.), 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.

RESULTSTop

The phylogenetic reconstruction and analyses of population genetics included 112 newly obtained cyt b sequences that were 1140 base pairs long (GenBank accession numbers: KT280012-KT280033), and 91 new sequences of the mitochondrial control region, with lengths varying between 396 and 398 base pairs (GenBank accession numbers: KT280034-KT280057). For cyt b, 22 different haplotypes were found among the specimens from the eastern Adriatic coast, while for mtDNA CR 24 unique haplotypes were obtained in this study. Population structure and haplotype frequencies are presented in Supplementary Table S2.

Population genetic structure and diversity

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.

Table 2. – Results of the genetic differentiation tests (HST, KST) based on cyt b data set (below diagonal) and mtDNA CR data set (above diagonal). Numbers in brackets are p-values based on permutation test. P-values lower than 0.02 are considered statistically significant after false discovery rate correction. Values in bold denote populations belonging to the Northern Adriatic group, plain text correspond to the Southern group.

HST Sečovlje Pag Dinjiška Nin Pantan Ston Ulcinj Narta
Sečovlje -0.01 (0.59) -0.02 (1.00) -0.01 (0.61) -0.03 (1.00) 0.58 (0.00) 0.83 (0.00) 0.34 (0.00)
Pag 0.01 (0.35) 0.00 (0.60) -0.02 (1.00) 0.00 (0.57) 0.44 (0.00) 0.63 (0.00) 0.25 (0.00)
Dinjiška 0.41 (0.00) 0.26 (0.00) 0.00 (0.58) -0.02 (1.00) 0.61 (0.00) 0.85 (0.00) 0.37 (0.00)
Nin -0.01 (0.61) -0.01 (0.97) 0.27 (0.00) 0.00 (0.61) 0.44 (0.00) 0.63 (0.00) 0.25 (0.00)
Pantan -0.02 (1.00) 0.03 (0.19) 0.47 (0.00) 0.02 (0.33) 0.59 (0.00) 0.84 (0.00) 0.35 (0.00)
Ston 0.52 (0.00) 0.30 (0.00) 0.33 (0.00) 0.37 (0.00) 0.57 (0.00) 0.03 (0.49) 0.03 (0.20)
Ulcinj 0.74 (0.00) 0.46 (0.00) 0.50 (0.00) 0.54 (0.00) 0.77 (0.00) 0.02 (0.27) 0.12 (0.02)
Narta 0.49 (0.00) 0.32 (0.00) 0.32 (0.00) 0.36 (0.00) 0.54 (0.00) 0.37 (0.00) 0.54 (0.00)
KST
Sečovlje -0.00 (0.73) 0.00 (1.00) -0.02 (0.74) 0.00 (0.73) 0.91 (0.00) 0.97 (0.00) 0.82 (0.00)
Pag 0.00 (0.28) 0.02 (0.48) 0.00 (0.67) 0.02 (0.39) 0.84 (0.00) 0.92 (0.00) 0.79 (0.00)
Dinjiška 0.31 (0.00) 0.16 (0.00) 0.00 (0.91) 0.00 (0.78) 0.89 (0.00) 0.98 (0.00) 0.84 (0.00)
Nin 0.00 (0.94) -0.01 (1.00) 0.23 (0.00) 0.00 (0.52) 0.86 (0.00) 0.91 (0.00) 0.79 (0.00)
Pantan 0.00 (1.00) -0.01 (0.75) 0.34 (0.00) 0.00 (0.35) 0.89 (0.00) 0.98 (0.00) 0.83 (0.00)
Ston 0.91 (0.00) 0.75 (0.00) 0.81 (0.00) 0.87 (0.00) 0.92 (0.00) 0.00 (1.00) 0.00 (1.00)
Ulcinj 0.96 (0.00) 0.80 (0.00) 0.86 (0.00) 0.91 (0.00) 0.97 (0.00) 0.00 (1.00) 0.00 (1.00)
Narta 0.85 (0.00) 0.70 (0.00) 0.75 (0.00) 0.82 (0.00) 0.87 (0.00) 0.43 (0.00) 0.53 (0.00)

Demographic history

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.

Table 3. – Interpretation of the results of the nested clade analysis of cyt b and mtDNA CR sequences of eastern Adriatic populations, based on the inference key of Templeton (2005)Templeton A.R. 2005. Inference Key for the Nested Haplotype Tree Analysis of geographical Distances. Available from URL: http://darwin.uvigo.es/. Only clades with significant chi-square statistics (χ2) are included.

Clade χ2 Chain of inference Demographic event inferred
cyt b data set
1-1 67.79 1-2-3-4-NO Restricted gene flow with isolation by distance
2-1 40.71 1-2-11-12-NO Contiguous range expansion
2-2 43.00 1-2-11-12-13-14-NO Long distance colonization and/or past fragmentation (not necessarily mutually exclusive)
Total cladogram 108.05 1-2 IOI-T Inconclusive outcome
mtDNA CR data set
Total cladogram 91.00 1-19-NO Allopatric fragmentation

sm4198fig2.jpg

Full size image

Fig. 2. – Statistical parsimony networks of the eastern Adriatic Aphanius populations based on mtDNA CR sequences (a) and nested clade design for the network containing cyt b haplotypes (b). Each line in the network represents a single mutational change; circle size of haplotypes is proportional to their number; small black circles represent unobserved haplotypes. Se-Sečovlje, Pg-Pag, Ni-Nin, Di-Dinjiška, Pn-Pantan, St-Ston, Ul-Ulcinj, Na-Narta. NG and SG explain populations belonging to Northern and Southern group, respectively.

Fu’s FS 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 FS 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 FS 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).

Phylogenetic position of A. fasciatus

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.

sm4198fig3.jpg

Full size image

Fig. 3. – BAY phylogram of the cyt b haplotypes obtained in this study (symbols in capital letters), together with those retrieved from GenBank (first data set). Numbers at nodes represent Bayesian posterior probabilities and maximum parsimony branch support values, respectively.

Phylogeographic structure of A. fasciatus and divergence time estimates

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).

sm4198fig4.jpg

Full size image

Fig. 4. – Median-joining network of mitochondrial control region haplotypes. Circle size indicates haplotype frequency, whereas branch lengths are proportional to the number of mutations. Where there are more than two mutational steps, the number of mutations is displayed along that branch. Roman numerals indicate clusters with haplotypes from Contrada Manzonara (ITA) (I); Pantano Viruca, Pantano Longarini and Foce Marcellino (ITA) (II); Marina di Modica (ITA) (III); Ston (CRO), Ulcinj (MNE), Narta (ALB) (e.g. Southern Adriatic group) and Foce Marcellino (ITA) (IV); Salina Curto Marsala, Salina Chiusa Trapani (ITA), Malta and Tunisia (V); Sečovlje (SLO), Nin, Pag, Dinjiška and Pantan (CRO) (e.g. Northern Adriatic group) (VI); Ganzirri (ITA) (VII); and Lesina (ITA) (VIII). Numbers by arrows represent estimated divergence time of population or cluster.

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.

DISCUSSIONTop

The A. fasciatus populations along the eastern Adriatic coast form two groups of genetically related populations: Northern and Southern. These groups are genetically distinct, whereas populations within each group are genetically similar. Observed pattern is concordant with the pattern found by Hrbek and Meyer (2003)Hrbek T., Meyer A. 2003. Closing of the Tethys Sea and the phylogeny of Eurasian killifishes (Cyprinodontiformes: Cyprinodontidae). J. Evol. Biol. 16: 17-36. in their sample, and explained by occasional migrations that prevent differentiation. Moreover, in this study, all haplotypes were specific to the group, with only one exception (haplotype APH1 of the cyt b gene). Based on phylogeographic analyses and the fact that one observed haplotype was shared by both groups (cyt b haplotype APH1, which was found in multiple individuals from Ston and Ulcinj, and also in one fish from Pag), it could be concluded that though the two groups are separated, sporadic migrations of individuals between them are possible. Differences observed in the present genetic structure and polymorphism of the two groups probably better reflect the differences in their demographic history rather than differences in recent habitats. Previous studies also mostly explained that interspecific differences are due to evolutionary events (Rocco et al. 2007Rocco L., Ferrito V., Costagliola D., et al. 2007. Genetic divergence among and within four Italian populations of Aphanius fasciatus (Teleostei, Cyprinodontiformes). Ital. J. Zool. 74: 371-379., Triantafyllidis et al. 2007Triantafyllidis A., Leonardos I., Bista I., et al. 2007. Phylogeography and genetic structure of the Mediterranean killifish Aphanius fasciatus (Cyprinodontidae). Mar. Biol. 152: 1159-1167.). Given the geographic position of the Northern and Southern groups, it is possible that the historical Neretva River presented a barrier to A. fasciatus dispersal and promoted the diversification of these two groups. Namely, geological evidence has demonstrated that the Neretva R. was significantly longer during the Pliocene, extending along the Pelješac Peninsula (Alfirević 1965Alfirević S. 1965. Geology of the Adriatic. Matica Hrvatska, Split. (in Croatian)). It is well known that the Neretva R. represented a barrier to the dispersal of several terrestrial taxa (Podnar et al. 2004Podnar M., Mayer W., Tvrtković N. 2004. Mitochondrial phylogeography of the Dalmatian wall lizard, Podarcis melisellensis (Lacertidae). Organ. Div. Evol. 4: 307-317., Kryštufek et al. 2007Kryštufek B., Buzan E.V., Hutchinson W.F., et al. 2007. Phylogeography of the rare Balkan endemic Martino’s vole, Dinaromys bogdanovi, reveals strong differentiation within the western Balkan Peninsula. Mol. Ecol. 16: 1221-1232., Ljubisavljević et al. 2007Ljubisavljević K., Arribas O., Džukić G., et al. 2007. Genetic and morphological differentiation of mosor rock lizards, Dinarolacerta mosorensis (Kolombatovic, 1886), with the description of new species from the Prokletije Mountain Massif (Montenegro) (Squamata: Lacertidae). Zootaxa 1613: 1-22.), though this is the first indication that it also represented an obstacle to the dispersal of marine organisms.

Although the Northern group is distributed over a greater number of localities, the Southern group is genetically more diverse. Furthermore, both investigated genetic markers indicated a greater variability of the populations from Pag and Nin than the remaining populations within the Northern group, and the greatest genetic polymorphism of Narta among southern populations. Environmental factors are highly variable at all localities, particularly salinity, which can often reach extreme values (37-330 psu), while temperature fluctuations (10-32°C) are also pronounced (Bakran-Petricioli and Jakl 2007Bakran-Petricioli T., Jakl Z. 2007. Marine habitats. Handbook for inventarization and monitoring. Ministry of Culture, State Department for Nature Protection, Croatia, Zagreb. (in Croatian)). However, fish diversity is the richest at Pantan (Matić-Skoko et al. 2005Matić-Skoko S., Peharda M., Pallaoro A., et al. 2005. Species composition, seasonal fluctuations, and residency of inshore fish assemblages in the Pantan estuary of the eastern middle Adriatic. Acta Adriat. 46 (2): 201-212.) and Ulcinj (Schneider-Jacoby et al. 2006Schneider-Jacoby M., Schwarz U., Sackl P., et al. 2006. Rapid assessment of the Ecological Value of the Bojana-Buna Delta (Albania/Montenegro). Euronatur, Radolfzell.), localities where the smallest level of intrapopulational genetic polymorphism of A. fasciatus was recorded. Pantan is the only location that has not been changed anthropologically, and it is still an inshore lagoon located at the mouth of the small Pantana River. All other localities have completely or partially been transformed into salt pans where the toothcarp is the only fish species present. Therefore, it is possible that the higher genetic diversity observed in the Pag, Nin and Narta populations is a consequence of lower competitive pressure, whereas extreme conditions do not limit the expansion of Aphanius populations.

The intrapopulational genetic diversity of all investigated A. fasciatus populations, even those expressing the highest levels of genetic polymorphism, was lower than the diversity recorded for other Mediterranean fish species (e.g. haplotype diversity of mtDNA CR of Chromis chromis [L., 1758] ranged between 0.792 and 1 [Domingues et al. 2005Domingues V.S., Bucciarelli G., Almada V.C., et al. 2005. Historical colonization and demography of the Mediterranean damselfish, Chromis chromis. Mol. Ecol. 14: 4051-4063.], while haplotype diversity of cyt b of Pomatoschistus minutus [Pallas, 1770] in the Mediterranean ranged between 0.822 and 0.844 [Gysels et al. 2004Gysels E.S., Helemans B., Patarnello T., et al. 2004. Current and historic gene flow of the sand goby Pomatoschistus minutus on the European Continental Shelf and in the Mediteranean Sea. Biol. J. Linn. Soc. 83: 561-576.]). Less polymorphism was also observed for A. iberus (Perdices et al. 2001Perdices A., Carmona J.A., Fernández-Delgado C., et al. 2001. Nuclear and mitochondrial data reveal high genetic divergence among Atlantic and Mediterranean populations of the Iberian killifish Aphanius iberus (Teleostei: Cyprinodontidae). Heredity 87: 314-324.) and it has been hypothesized that rapid and pronounced changes of physicochemical and biological features may reduce population size, which, in combination with genetic drift, inbreeding and founder effects, could result in a loss of genetic diversity (Cognetti and Maltagliati 2000Cognetti G., Maltagliati F. 2000. Biodiversity and adaptive mechanisms in brackish water fauna. Mar. Pollut. Bull. 40: 7-14.).

The phylogenetic position of A. fasciatus indicates its closer relatedness to several Iranian Aphanius species (A. sophiae, A. mesopotamicus, A. farsicus, A. isfahanensis and A. vladykovi) than to other Mediterranean species (A. iberus, A. baeticus and A. saourensis). That conclusion is in accordance with the results of previous studies (Hrbek and Meyer 2003Hrbek T., Meyer A. 2003. Closing of the Tethys Sea and the phylogeny of Eurasian killifishes (Cyprinodontiformes: Cyprinodontidae). J. Evol. Biol. 16: 17-36.), though samples were not composed of the same species and populations. Populations from the eastern Adriatic coast, analysed here for the first time, clustered into the A. fasciatus lineage.

Better resolution of intraspecific phylogenetic relationships among A. fasciatus populations was obtained through the phylogenetic reconstruction of the mitochondrial control region. The long-term isolation and abruption of gene flow among geographically close populations on the one hand, and the homogenity of some distinct populations with clearly evidence of gene flow among them on the other, are the most prominent though contradictory features of the genetic structure of A. fasciatus. For a more precise explanation of the observed phenomenon, further investigations of the evolutionary history of this species are needed throughout its entire distribution range. Nevertheless, the results of this study also point out significant intraspecific genetic structuring of A. fasciatus (Rocco et al. 2007Rocco L., Ferrito V., Costagliola D., et al. 2007. Genetic divergence among and within four Italian populations of Aphanius fasciatus (Teleostei, Cyprinodontiformes). Ital. J. Zool. 74: 371-379., Triantafyllidis et al. 2007Triantafyllidis A., Leonardos I., Bista I., et al. 2007. Phylogeography and genetic structure of the Mediterranean killifish Aphanius fasciatus (Cyprinodontidae). Mar. Biol. 152: 1159-1167., Ferrito et al. 2013Ferrito V., Pappalardo A.M., Canapa A., et al. 2013. Mitochondrial phylogeography of the killifish Aphanius fasciatus (Teleostei, cyprinodontidae) reveals highly divergent Mediterranean populations. Mar. Biol. 160: 3193-3208.).

Our results imply the importance of the Pleistocene epoch in shaping the recent population genetic structure and diversity of eastern Adriatic toothcarp populations. During the glacial maxima (in particular the Last Glacial Maximum, 25000-15000 YA), the entire northern part of the Adriatic Basin, down to the Zadar-Pescara line, was a wide fluvio-lacustrine plain with a coastline located at the northern edge of the Meso-Adriatic depression (Pirazzoli 2000Pirazzoli P.A. 2000. Sea-Level Changes - the Last 20,000 Years. Wiley and Sons Ltd. Chichester, pp 1-211.). The northern border of the Adriatic Sea during the Last Glacial Maximum was located just north of the line connecting the Pantan and Lesina lagoons. The intrapopulational genetic divergence of the majority of the investigated Adriatic populations is of Pleistocene origin, implying that the colonization of the northern and southern parts of the eastern Adriatic localities did not occur from a single locality after glaciations, but rather that populations remained in isolated patches of residual sea water or in temporary lakes that formed in the northern part of the Adriatic Basin. When connections between those localities were reestablished, gene flow among populations was enabled. Rocco et al. (2007)Rocco L., Ferrito V., Costagliola D., et al. 2007. Genetic divergence among and within four Italian populations of Aphanius fasciatus (Teleostei, Cyprinodontiformes). Ital. J. Zool. 74: 371-379. indicated that paleoclimatic fluctuations in the Quatarnary influenced the genetic structure of A. fasciatus, though the detailed timing of population divergence has not been determined before. The conclusion that these two population groups have distinct evolutionary courses, despite inhabiting similar habitats located in geographic proximity, is corroborated by the NCP analysis results (see Table 3). For the populations of the Northern group, traces of restricted gene flow with isolation by distance were noted, which is concordant with the previously mentioned hypothesis that these populations survived long periods of the Pleistocene in small pools. On the other hand, the southern populations inhabit localities that were on the sea shore even during the glacial maxima. The importance of the Pleistocene epoch for the demographic history of the investigated species is highlighted by dating the onset of presumable population expansion of the Narta, Nin and Pag populations to that epoch.

Implications for conservation

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.

ACKNOWLEDGEMENTSTop

We gratefully acknowledge the contribution of our colleagues in collecting samples, especially Pero Tutman, Armin Pallaoro, Meta Povž, Aleksandar Joksimović, Saimir Bequiraj, Luka Šupraha, Filip Bukša and Draško Holcer. We are also grateful to Aljoša Duplić for his help in obtaining financial support. This study was financed by the Croatian Ministry of Science, Education and Sport (Project nos. 119-1782739-1233 and 119-0000000-3184) and the State Institute for Nature Protection.

REFERENCESTop

Alfirević S. 1965. Geology of the Adriatic. Matica Hrvatska, Split. (in Croatian)

Astraldi M., Balopoulos S., Candela J., et al. 1999. The role of straits and channels in understanding the characteristics of Mediterranean circulation. Prog. Oceanogr. 44: 65-108.
http://dx.doi.org/10.1016/S0079-6611(99)00021-X

Bakran-Petricioli T., Jakl Z. 2007. Marine habitats. Handbook for inventarization and monitoring. Ministry of Culture, State Department for Nature Protection, Croatia, Zagreb. (in Croatian)

Bardakci F., Tatar N., Hrbek T. 2004. Genetic relationships between Anatolian species and subspecies of Aphanius based on RAPD markers. Biologia (Bratislava) 59: 559-566.

Clement M., Posada D., Crandall K.A. 2000. TCS: a computer program to estimate gene genealogies. Mol. Ecol. 9: 1657-1660.
http://dx.doi.org/10.1046/j.1365-294x.2000.01020.x

Cognetti G., Maltagliati F. 2000. Biodiversity and adaptive mechanisms in brackish water fauna. Mar. Pollut. Bull. 40: 7-14.
http://dx.doi.org/10.1016/S0025-326X(99)00173-3

Domingues V.S., Bucciarelli G., Almada V.C., et al. 2005. Historical colonization and demography of the Mediterranean damselfish, Chromis chromis. Mol. Ecol. 14: 4051-4063.
http://dx.doi.org/10.1111/j.1365-294X.2005.02723.x

Drummond A.J., Suchard M.A., Xie D., et al. 2012. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29: 751-761.
http://dx.doi.org/10.1093/molbev/mss075

Dulčić J., Soldo A., Jardas I. 2005. Adriatic fish biodiversity and review of bibliography related to Croatian small-scale coastal fisheries. AdriaMed Technical Documents nº15: 103-125.

Ferrito V., Pappalardo A.M., Canapa A., et al. 2013. Mitochondrial phylogeography of the killifish Aphanius fasciatus (Teleostei, cyprinodontidae) reveals highly divergent Mediterranean populations. Mar. Biol. 160: 3193-3208.
http://dx.doi.org/10.1007/s00227-013-2307-4

Gysels E.S., Helemans B., Patarnello T., et al. 2004. Current and historic gene flow of the sand goby Pomatoschistus minutus on the European Continental Shelf and in the Mediteranean Sea. Biol. J. Linn. Soc. 83: 561-576.
http://dx.doi.org/10.1111/j.1095-8312.2004.00411.x

Hall T.A. 1999. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucl. Acids. Symp. Ser. 41: 95-98.

Hrbek T., Meyer A. 2003. Closing of the Tethys Sea and the phylogeny of Eurasian killifishes (Cyprinodontiformes: Cyprinodontidae). J. Evol. Biol. 16: 17-36.
http://dx.doi.org/10.1046/j.1420-9101.2003.00475.x

Kryštufek B., Buzan E.V., Hutchinson W.F., et al. 2007. Phylogeography of the rare Balkan endemic Martino’s vole, Dinaromys bogdanovi, reveals strong differentiation within the western Balkan Peninsula. Mol. Ecol. 16: 1221-1232.
http://dx.doi.org/10.1111/j.1365-294X.2007.03235.x

Lee W.-J., Conroy J., Howell W.H., et al. 1995. Structure and evolution of teleost mitochondrial control regions. J. Mol. Evol. 41: 54-66.
http://dx.doi.org/10.1007/BF00174041

Librado P., Rozas J. 2009. DNAsp v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451-1452.
http://dx.doi.org/10.1093/bioinformatics/btp187

Lipej L., Dulčić J. 2010. Checklist of the Adriatic Sea Fishes. Zootaxa 2589: 1-92.

Ljubisavljević K., Arribas O., Džukić G., et al. 2007. Genetic and morphological differentiation of mosor rock lizards, Dinarolacerta mosorensis (Kolombatovic, 1886), with the description of new species from the Prokletije Mountain Massif (Montenegro) (Squamata: Lacertidae). Zootaxa 1613: 1-22.

Maggio T., Lo Brutto S., Garoia F., et al. 2009. Microsatellite analysis of red mullet Mullus barbatus (Perciformes, Mullidae) reveals the isolation of the Adriatic Basin in the Mediterranean Sea. ICES J. Mar. Sci. 66: 1883-1891.
http://dx.doi.org/10.1093/icesjms/fsp160

Matić-Skoko S., Peharda M., Pallaoro A., et al. 2005. Species composition, seasonal fluctuations, and residency of inshore fish assemblages in the Pantan estuary of the eastern middle Adriatic. Acta Adriat. 46 (2): 201-212.

Meyer A. 1993. Evolution of mitochondrial DNA in fishes. In: Hochachka P.W., Mommsen T.P. (eds), Biochemistry and Molecular Biology of Fishes, vol 2. Elsevier Science Publisher, Amsterdam, pp 1-38.

Panchal M. 2007. The automation of nested clade phylogeographic analysis. Bioinformatics 23: 509-510.
http://dx.doi.org/10.1093/bioinformatics/btl614

Parenti L.R., Tigano C. 1993. Polymorphic skeletal characters in Aphanius fasciatus Teleostei: Cyprinodontiformes. Copeia 4: 1132-1137.
http://dx.doi.org/10.2307/1447094

Perdices A., Carmona J.A., Fernández-Delgado C., et al. 2001. Nuclear and mitochondrial data reveal high genetic divergence among Atlantic and Mediterranean populations of the Iberian killifish Aphanius iberus (Teleostei: Cyprinodontidae). Heredity 87: 314-324.
http://dx.doi.org/10.1046/j.1365-2540.2001.00888.x

Pirazzoli P.A. 2000. Sea-Level Changes - the Last 20,000 Years. Wiley and Sons Ltd. Chichester, pp 1-211.

Podnar M., Mayer W., Tvrtković N. 2004. Mitochondrial phylogeography of the Dalmatian wall lizard, Podarcis melisellensis (Lacertidae). Organ. Div. Evol. 4: 307-317.
http://dx.doi.org/10.1016/j.ode.2004.04.004

Posada D., Crandall K.A., Templeton A.R. 2000. GeoDis: a program for the cladistic nested analysis of the geographical distribution of genetic haplotypes. Mol. Ecol. 9: 487-488.
http://dx.doi.org/10.1046/j.1365-294x.2000.00887.x

Ramos-Onsins S.E., Rozas J. 2002. Statistical properties of new neutrality tests against population growth. Mol. Biol. Evol. 19: 2092-2100.
http://dx.doi.org/10.1093/oxfordjournals.molbev.a004034

Rocco L., Ferrito V., Costagliola D., et al. 2007. Genetic divergence among and within four Italian populations of Aphanius fasciatus (Teleostei, Cyprinodontiformes). Ital. J. Zool. 74: 371-379.
http://dx.doi.org/10.1080/11250000701451225

Rogers A.R., Harpending H. 1992. Population growth makes waves in the distribution of pairwise genetic differences. Mol. Biol. Evol. 9: 552-569.

Ronquist F., Huelsenbeck J.P. 2003. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19: 1572-1574.
http://dx.doi.org/10.1093/bioinformatics/btg180

Schneider-Jacoby M., Schwarz U., Sackl P., et al. 2006. Rapid assessment of the Ecological Value of the Bojana-Buna Delta (Albania/Montenegro). Euronatur, Radolfzell.
http://www.euronatur.org/Publications.411.0.html

Schmidt T.R., Gold, J.R. 1993. Complete sequence of the mitochondrial cytochrome b gene in the cherryfin shiner, Lytrurus roseipinnis (Teleostei: Cyprinidae). Copeia 3: 880-883.
http://dx.doi.org/10.2307/1447258

Swofford D.L. 2002. PAUP*: Phylogenetic Analysis Using Parsimony (*and Other Methods), Version 4. Sinauer Associates.

Templeton A.R. 1998. Nested clade analyses of phylogeographic data: testing hypotheses about gene flow and population history. Mol. Ecol. 7: 381-397.
http://dx.doi.org/10.1046/j.1365-294x.1998.00308.x

Templeton A.R. 2005. Inference Key for the Nested Haplotype Tree Analysis of geographical Distances. Available from URL:
http://darwin.uvigo.es/

Templeton A.R., Crandall K.A., Sing C.F. 1992. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics 132: 619-633.

Triantafyllidis A., Leonardos I., Bista I., et al. 2007. Phylogeography and genetic structure of the Mediterranean killifish Aphanius fasciatus (Cyprinodontidae). Mar. Biol. 152: 1159-1167.
http://dx.doi.org/10.1007/s00227-007-0760-7

SUPPLEMENTARY MATERIAL

The following material is available through the online version of this article and at the following link:
http://www.icm.csic.es/scimar/supplm/sm04198esm.pdf

Table S1. – The list of localities of the origin of sequences retrieved from GenBank (N, number of sequences). Country codes: SLO, Slovenia; CRO, Croatia; MNE, Montenegro; ALB, Albania; ITA, Italy; MLT, Malta; TUN, Tunisia; GRE, Greece; ESP, Spain; ALG, Algeria; IRN, Iran.

Table S2. – Haplotype frequencies for cyt B and mt DNA CR from the Adriatic localities. NG and SG explain populations belonging to the Northern and Southern groups, respectively.