Inter-population differences in otolith morphology are genetically encoded in the killifish Aphanius fasciatus ( Cyprinodontiformes )

Inter-population differences in otolith shape, morphology and chemistry have been used effectively as indicators for stock assessment or for recognizing environmental adaptation in fishes. However, the precise parameters that affect otolith morphology remain incompletely understood. Here we provide the first direct support for the hypothesis that inter-population differences in otolith morphology are genetically encoded. The study is based on otolith morphology and two mitochondrial markers (D-loop, 16S rRNA) of three natural populations of Aphanius fasciatus (Teleostei: Cyprinodontidae) from Southeast Tunisia. Otolith and genetic data yielded congruent tree topologies. Divergence of populations likely results from isolation events in the course of the Pleistocene sea level drops. We propose that otolith morphology is a valuable tool for resolving genetic diversity also within other teleost species, which may be important for ecosystem management and conservation of genetic diversity. As reconstructions of ancient teleost fish faunas are often solely based on fossil otoliths, our discoveries may also lead to a new approach to research in palaeontology.


INTRODUCTION
Otoliths are aragonitic mineralizations that are arranged in three pairs in the inner ear of teleost fishes, where they play an important role in the senses of hearing and balance (Popper et al. 2005).The saccular otoliths (termed otoliths in the following) provide morphological characters that can be used in family, genus, and species discrimination (e.g.Smale et al. 1995, Volpedo and Echeverría 2000, Tuset et al. 2008).In addition, otoliths can be found abundantly as fossils and are used to reconstruct ancient teleost fish diversity, zoogeography and evolution (e.g.Nolf 1995, Girone and Nolf 2009, Schwarzhans and Bratishko 2011).
Inter-population variability is known to occur widely in otoliths, especially with regard to size and contour, daily and annual growth rings, trace elements, and isotopic compositions.Fourier analysis and landmarks have been used to quantify otolith variation (size and contour) between species, populations, and even stocks (e.g.Castonguay et al. 1991, Torres et al. 2000, Monteiro et al. 2005, Mérigot et al. 2007).Differences in otolith chemistry have been used for stock discrimination, analyses of population structure, reconstruction of the environmental history, and ecosystem monitoring (e.g.Campana 1999, Volpedo and Cirelli 2006, Woods et al. 2010).However, only a few studies have addressed the question whether inter-population differences in otolith traits are genetically encoded or the result of differences in certain habitat parameters (Torres et al. 2000, Stransky et al. 2008, Lombarte et al. 2010, Vignon and Morat 2010).
Mitochondrial DNA is a useful tool in population studies for various reasons, including maternal transmission, limited length and elevated evolution rate (Brown et al. 1979, Wilson et al. 1985).The mitochondrial control region (D-loop, Brown 1983) has been found to be an excellent population marker and has thus been used frequently in the analysis of phylogeography in marine fishes (e.g.Bernardi 2000, Stepien et al. 2001, Astolfi et al. 2005).In addition, many studies have used ribosomal RNA genes because of their highly repetitive nature, ease of manipulation and biological importance in order to elucidate the evolutionary and demographic history of populations (e.g.Sollner-Webb and Mougey 1991, Ghigliotti et al. 2008, Torres-Machorro et al. 2010).
Here we study whether inter-population differences in otolith morphology of a tooth-carp species (Teleostei: Cyprinodontidae) are consistent with data obtained from mitochondrial markers.Aphanius fasciatus (Valenciennes, in Humboldt and Valenciennes 1821) is used as a model because tooth-carps and killifishes are particularly suitable for the study of microevolutionary processes (Villwock 1994, Fuller 2008, Martin and Wainwright 2011).A. fasciatus is widely distributed in coastal and brackish-water habitats of the central and eastern Mediterranean Sea (Wildekamp 1993).A specific attribute of A. fasciatus is that it has no economic use, so the populations are not manipulated but rather reflect the natural distribution of the genotype (Maltagliati 1999).A. fasciatus therefore serves as model organism in many studies on genetic structures and levels of differentiation based on different methods, i.e. osteological characters (Tigano et al. 1999, 2001, Ferrito et al. 2003, 2007), genetic data (Hrbek and Meyer 2003, Triantafyllidis et al. 2007, Pappalardo et al. 2008), both osteological and genetic data (Maltagliati et al. 2003, Tigano et al. 2004, 2006), and otoliths (Reichenbacher et al. 2007).These studies have revealed that A. fasciatus from Italy, Greece and Northern Tunisia is genetically structured by geographic isolation and that a clear link exists between genetic structures and the Pleistocene history of the Mediterranean Sea.Based on this, we hypothesize that genetic structuring is present and recognizable in mitochondrial markers in A. fasciatus from SE Tunisia (of which genetic differentiation has not previously been studied) and, given that otolith morphology is genetically encoded, is also recognizable in differences in otolith morphology.We tested this hypothesis based on three natural A. fasciatus populations using otolith morphometry, D-loop and 16S rRNA genes to assess levels of divergence and compare the degree of genetic variability.

Sampling and sites
Adult individuals, identifiable by the typical pigmentation, with standard lengths (SL) ranging from 25 to 43 mm were collected by hand nets from the coastal sites Sfax and Luza in the Gulf of Gabes and the brackish site Hamdoun on the Dkhila coast (fed by the freshwater source Oued Hamdoun) (Fig. 1, Table 1).Female and male individuals were manually sorted by their distinctive colour pattern and pigmentation and preserved in 99.9% ethanol.The relatively high number of specimens for the otolith study is due to the fact that we had to check whether sex dimorphism is reflected in otolith morphology as well, and that artefacts due to restricted sample sizes had to be avoided.Dissected otoliths are deposited in the Bavarian State Collection for Palaeontology and Geology in Munich, Germany (accession number BSPG 2009 X).

Otolith preparation and analyses
Otolith dissection and preparation followed standard procedures.Study of otolith morphology was based on SEM images (LEO 1430 VP at the Zoological State Collection Munich); otolith terminology follows Nolf (1985) and Tuset et al. (2008) (Fig. 2A).For morphometric analyses, digital images were captured using a Leica DFC 295 camera and the Leica Image Access Software (IMAGIC 1000, Imagic Bildverarbeitung AG, Glattbrugg, Switzerland).Three angles and eight linear distances were measured (Fig. 2B) and ten otolith variables were calculated from the measurements; they describe specific segments of the otolith (Table 2A; for details see Reichenbacher et al. 2007).

Molecular analysis
The molecular study included mitochondrial DNA D-loop (D-loop) and 16S ribosomal RNA (16S rRNA) analyses.Total DNA was extracted from muscle tissue preserved in ethanol, using a Wizard Genomic DNA extraction kit following the manufacturer's protocol (Promega).The concentration of extracted DNA was spectrophotometrically estimated.
For the D-loop analysis, a 378 bp fragment of the mitochondrial control region was amplified from total DNA extracts using newly designed primers: AFDF: 5'-CCCCGCCGCCC ATCAATAAT-3' and AFDR: 5'-CCAGGAATAATTCACTAAGTGC-3'.D -loop sequence of Aphanius fasciatus (GenBank accession AM884570) was used to make the alignment that allowed the definition of these primers.
For the 16S rRNA gene analysis, a fragment of 154 bp was amplified from total DNA extracts using also newly designed primers: AFrF 5'-GTA ATC CAG GTC AGT TTC TAT CT-3' and AFrR 5'-GTA ATC CAG GTC AGT TTC TAT CT-3'.The primers were designed on the basis of the 16S ribosomal RNA sequence of Aphanius fasciatus (GenBank accession EF640854).
Polymerase chain reaction (PCR) amplification was performed on total genomic DNA in a total volume of 50 μl containing: 20ng DNA, 2.5mM dNTPs, Taq Buffer (1X), 50mM MgCl 2 , 25 pmol from each primer and 0.5 U of Taq DNA polymerase (Madison, Promega).The final volume (50 μl) was adjusted by adding sterile distilled water.Negative controls were performed for all reactions.
The cycling profile was 94°C for 3 min, followed by 35 cycles at 94°C for 1 min, 56°C for 1 min, 72°C for 1 min, and a final 10-min extension step at 72°C.
The amplified DNA segments encoding D-loop and 16S rRNA genes were purified using the "Wizard PCR preps DNA purification kit" according to the manufacturer's instructions (Promega, Madison, WI) and then sequenced.Cycle sequencing was performed by Macrogen Inc. using Automated Applied Biosystems (AB) sequencing and the Taq Dye Deoxy Terminator cycle sequencing kit.Chromatograms and alignments were visually checked and verified.
JModelTest (Posada 2008) was run to determine the most suitable model of DNA evolution to consider for our set of sequences through the five model selection strategies available in the program.neighbour-joining and maximum likelihood analyses were performed using SEAVIEW (Gouy et al. 2010), and PHYML on line at the ATGC Montpellier bioinformatics platform (v3.0,Guindon and Gascuel 2003), respectively.The evolutionary relationships among sequences were also estimated with the Bayesian Markov chain Monte Carlo (MCMC) phylogenetic analyses (BI) as included in MrBayes v3.1 using the default priors (Huelsenbeck andRonquist 2001, Ronquist andHuelsenbeck 2003).Three heated chains and a single cold chain were employed in all MCMC analyses, and runs were initiated with random trees.Two independent MCMC runs were conducted with 2 million generations per run, and trees and parameters were sampled every 100 generations.Stationarity was assessed by examining the average standard deviation of split frequencies and the potential scale reduction factor (Ronquist et al. 2005).For each run, the first 25% of sampled trees were discarded as burn-in.Bayesian posterior probabilities were used to assess branch support of the MCMC tree.
The genetic variation within groups was then estimated using basic statistics.Haplotype (h) and nucleotide (π in percentage) diversities and their standard deviations (±SD) were estimated using DNASP (v4.10.9,Rozas et al. 2003).The MEGA software version 3.1 (Kumar et al. 2004) was used to estimate genetic distances (Tamura and Nei 1993).
Finally, F ST values were calculated and we used analysis of molecular variance (AMOVA), as implemented in ARLEQUIN (version 3, Excoffier et al. 2005) to quantify the genetic variation among groups, among populations within groups, and within populations.The variance components of the different hierarchical levels were tested statistically by nonparametric randomization tests using 10 000 permutations.
Monophyly of trees was tested by creating a monophyletic backbone tree for each of the clades.The backbone tree was then used as a constraint tree for a maximum likelihood search using the parameters from the unconstrained maximum likelihood tree.The monophyletic signal was tested with the Shimodaira-Hasegawa test (Shimodaira and Hasegawa, 1999) using PAUP* v4.0b10 (Swofford 2003), with RELL (Kishino et al. 1990) (resampling estimated log-likelihood) optimization and 10 000 bootstrap replicates.

Results of otolith analysis
All studied A. fasciatus otoliths are characterized by a triangular to slightly rounded contour (Fig. 3).The rostrum is present and has a rounded or slightly pointed tip.The antirostrum is shorter than the rostrum and rounded, the excisura ostii is wide with a deep, acute notch.The sulcus is located in a median position and subdivided into a small, slightly deep-ened, funnel-like shaped ostium and a long, straight or posteriorly slightly bent cauda covered with colliculi.The inner face of all otoliths is planar, the outer face slightly convex.

Within-population differences
Based on qualitative comparison of otolith images (Fig. 3) and evaluation of standard deviations of otolith variables (Table 1A), the variation of the otolith morphology within a population is low for both the Sfax (Fig. 3a-l) and Corsica (Fig. 3k1-p1) fishes.Within the Luza population (Fig. 3m-x), a few otoliths possess a rounded rather than triangular contour (Fig. 3s) or an unusually short rostrum (Fig. 3x), so the variation within the Luza population is considered as low to moderate.In contrast, distinctive intra-population differences are visible in the otoliths from the Hamdoun fishes (Fig. 3y-j1).These otoliths include a "normal" otolith type (Fig. 3y-z, e1, g1-i1) similar to that seen in Sfax, Luza and Corsica individuals, and a second, more elongate, narrow-triangular type (Fig. 3b1-d1, j1); intermediate morphologies are also present (Fig. 3a1, f1).The relatively high standard deviations of some otoliths variables, i.e.L/H index, relative rostrum height, excisura angle, and relative dorsal length (Table 2A), corroborate these observations for the Hamdoun fishes.

Inter-population differences
Based on qualitative comparison of otolith images (Fig. 3) and uni-and multivariate statistics (Table 2), no clear differences in otolith morphology and otolith variables can be recognized between the populations from Sfax and Luza.In contrast, three otolith variables are significantly different between Sfax/Luza and Hamdoun, and four otoliths variables are significantly different between Sfax/Luza and Corsica.Moreover, a single otolith variable is different between the individuals from Hamdoun and Corsica (Table 2B).The great similarity between the Sfax and Luza otoliths, but their clear differentiation from the Hamdoun and Corsica otoliths, is additionally shown by a dendrogram using the Euclidean distance as a measure of dissimilarity and the "between groups linkage method" as the clustering algorithm (Fig. 4).Moreover, the canonical discriminant analysis reveals a high classification success for each group (Sfax/Luza, Hamdoun, Corsica) and thus clearly supports the importance of otolith variables for deciphering inter-population differences (Table 3).

Results of molecular analysis
A total of 16 sequences of 378 bp was obtained for the D-loop of A. fasciatus (Genbank accession numbers: JX406312 to JX406327).Among them, eight different haplotypes were identified.Fourteen sites were variable and 11 were parsimony-informative.The TPM1uf model with a gamma distribution shape parameter equal to 0.016 and the proportion of invariable sites equal to 0.169 was the best evolutionary model.The nucleotides frequencies were 37.01%, 17.38%, 15.79% and 29.82% for A, C, G and T, respectively.The different phylogenetic analyses provided very similar topologies and for the sake of simplicity we only present the phylogenetic tree corresponding to the maximum likelihood analysis (Fig. 5A).The A. fasciatus sequences were distributed among two clades that are separated by a moderate value of Tamura and Nei genetic divergence (2.7±0.08%).Clade I consists of the Sfax and Luza specimens, while clade II comprises the Hamdoun specimens.Genetic diversity in clade II (Nucleotide diversity Pi=0.00±0.000;Haplotype diversity Hd=0.99±0.000) is slightly higher than in clade I (Pi=0.0005±0.000;Hd=0.94±0.070).The mean F ST value between clades I and II is 0.88 (p=0.001).Most of the variation is explained between the two clades (89.7%), whereas variation within clades is smaller (11.75%) (Fig. 5A).
A total of 16 sequences of 154 bp was obtained for the 16S rRNA of A. fasciatus (pending submission to treebase).Among them, four different haplotypes were identified.Fifteen sites were variable and two were parsimony informative.The F81 model with a gamma distribution shape parameter equal to 0.016 and the proportion of invariable sites equal to 0.00 was the best evolutionary model.The nucleotides frequencies were 31.9%,27.47%, 19.7% and 20.93% for A, C, G and T, respectively.The different phylogenetic analyses resulted in very similar topologies and for the sake of simplicity we only represent the phylogenetic tree corresponding to the maximum likelihood analysis (Fig. 5B).The A. fasciatus sequences were distributed among two clades, of which clade I includes the Sfax and Luza specimens, and clade II comprises the Hamdoun specimens.The clades are separated by a moderate value of Tamura and Nei genetic divergence (0.7±0.007%), and genetic diversity in clade I (Pi=0.0029±0.001;Hd=0.417±0.191) is slightly higher than in clade II (Pi=0.002±0.000;Hd=0.286±0.019).The mean F ST value between clades I and II is 0.79 (p=0.001).Most of the variation is explained among the two clades (80.07%), whereas variation within the clades is smaller (20.72%).
For both mitochondrial markers (D-loop and 16Sr-RNA), Shimodaira-Hasegawa tests based on the maximum likelihood tree indicate that clade II (Hamdoun) represents a monophyletic clade (P<0.05).

DISCUSSION
Fish populations inhabiting similar environments and interconnected through migration and gene flow usually display largely similar phenotypic and genetic traits (e.g.Carvalho 1993).Aphanius fasciatus is widely distributed along the southeastern coast of Tunisia.Hence, gene flow, resulting in mixing between populations and stabilization of the basic genome, can be expected regardless of whether the species is known for having a rather sedentary life history, with large demersal eggs and without larval dispersal stages (Maltagliati 1999).However, both otolith and mitochondrial marker analyses (D-loop, 16S rRNA) are indicative of restricted gene flow between the A. fasciatus individuals from Sfax/ Luza and Hamdoun.Moreover, the dendrogram and trees reveal congruent topologies (Figs 4-5).In the 16S rRNA-based tree, the monophyly of clade II (Hamdoun) is slightly less well supported, which might be due to the reduced substitution rate in 16S in comparison with the mitochondrial control region (D-loop).Another line of evidence in support of a genetic basis for population differences in otolith morphology is the absence of otolith differences between the A. fasciatus individuals from the heavily polluted site Sfax (Messaoudi et al. 2009) and the non-polluted site Luza.
These results raise the question as to what caused the genetic divergence between the A. fasciatus populations from SE Tunisia?The occurrence of genetic divergence within a species is affected by many factors, including population size, time since isolation, and porosity of the isolating barrier or mechanism (Frankham 1995).Genetic drift is faster if populations are small or the isolating barrier is very effective (Leis et al. 2011).A factor that significantly influenced the physical connection between coastal sites in the Mediterranean Sea was the changing global sea level between glacial and interglacial climate conditions during the Pleistocene (1,8 million to 11000 years ago) (Lambeck et al. 2002).Genetic differentiation among populations that may be linked with Pleistocene sea level falls are known for several teleost species from the Mediterranean Sea, including the silverside Atherina boyeri Risso, 1810 (Milana et al. 2012 and references therein) and the goby Pomatoschistus tortonesei (Mejri et al. 2009).According to these studies, distinctive oceanic currents in the Mediterranean Sea represent hydrographic barriers for coastal species and maintained their differentiation for long periods of time.In the case of A. fasciatus, the strong sea level drop (up to -120 m) during the last Pleistocene glaciation (20000 years ago) was suggested to be responsible for the genetic divergence between populations from Italy (Rocco et al. 2007).Moreover, the sea level fall during the Early Pleistocene (ca.1.7 million years ago) was favoured as an explanation for the genetic divergence between A. fasciatus populations from Sicily, northern Tunisia and Malta (Tigano et al. 2006).
We therefore hypothesize that the sea level fall during the Early or Late Pleistocene and the current oceanographic conditions have affected the genetic structure of A. fasciatus from SE Tunisia.Notably, the allopatric divergence that developed during times of isolation was strong enough to prevent gene flow during the subsequent sea level rise when the climate became warm and wet again.As mentioned above, the Hamdoun site (inhabited by clade II) is fed by freshwater, whereas the sites of Sfax and Luza (inhabited by clade I) are strictly coastal.We therefore suggest that the specific salinity and hydrological parameters of the Hamdoun site have triggered the genetic divergence in this population.This assumption concurs with results from previous studies, in which specific salinity parameters have been detected as the causative agents of strong selective pressure on organisms that eventually led to reproductive isolation in spite of possible intermingling during larval and adult life stages (Bekkevold et al. 2005, Fuller 2008, Williams et al. 2008).
If our assumption that inter-population differences in otolith morphology are genetically encoded is correct, then the results of our study also indicate genetic divergence between A. fasciatus from Corsica (Fanjo Delta) and SE Tunisia (Sfax/Luza).A likely explanation for this pattern is that the respective populations are separated by the oceanic currents of the Sicily Strait (Astraldi et al. 1999, 2002, Béranger et al. 2005), which is known to represent a breakpoint to gene flow in the Mediterranean basin (Mejri et al. 2009).On the other hand, the otolith similarities between A. fasciatus from Corsica and the northernmost of the here studied SE-Tunisian sites, i.e.Hamdoun, are difficult to explain without conducting additional molecular studies as the Strait of Sicily probably prevented gene flow also between these sites.Also the development of new markers (e.g.nuclear genome) would be useful in order to elucidate in detail the micro-evolutionary processes acting in Aphanius fasciatus populations.
The results of this study shed new light on previous work on the differences in otolith morphology between populations of Aphanius iberus (Reichenbacher and Sienknecht 2001) and A. dispar (Reichenbacher et al. 2009a, 2009b, Teimori et al. 2012a, 2012b).In these studies, otolith differences have been interpreted as indicating allopatric genetic divergence, but a test of this interpretation by molecular data analyses has not been conducted.We assume that a genetic basis in otolith morphology is indicative for population differences also in A. iberus and A. dispar and probably also in several other species of Aphanius.Our study therefore may also serve as a starting point for the development of a new utility of otolith morphology for analysing genetic structures between populations of Aphanius species and also between other species of teleost fishes.Such analyses are important for ecosystem management and conservation of genetic diversity (e.g.Clarke et al. 2011, Leis et al. 2011).In addition, this new importance of otoliths can provide an innovative approach for studying evolution and phylogeography of fossil teleost fish faunas, which are often preserved solely based on otoliths.
fig. 1. -Geographic overview of the study sites.
fig. 5. -Maximum likelihood phylogenetic relationships for 16 specimens of Aphanius fasciatus from Sfax, Luza and Hamdoun (SE Tunisia) based on D-loop sequences (A) and 16S RNA sequences (B).Numbers beside nodes indicate bootstrap values (% > 50) obtained by the neighbour-joining, maximum likelihood and posterior probabilities for Bayesian analysis, respectively.

tabLe 1 .
-Details of the sites and samples.SL, standard length.

tabLe 2 .
-Values of otolith variables (mean values and standard deviations) in the studied Aphanius fasciatus populations (A) and list of otolith variables that differ significantly between the studied populations (one-way ANOVA with Duncan post hoc test, p<0.001) (B); n, number of specimens.Luza vs. Hamdoun Rel.rostrum length, rel.antirostrum length, rel.medial length Sfax/Luza vs. Corsica Length-height-index, rel.antirostrum height, rel.rostrum length, posterior angle Hamdoun vs. Corsica Rel.medial length fig. 4. -Relationships between Aphanius fasciatus populations as suggested by cluster analysis based on otolith variables (Euclidean distances of mean values ± standard deviations).tabLe 3. -Classification matrix of the CDA (jackknifed) based on the otolith variables of the studied Aphanius fasciatus populations.The otoliths from Sfax and Luza were merged.