Eleven haplotypes were recognized among the 117 T. zillii D-loop sequences obtained (Table 1). The most widespread haplotype (one) was found in all samples except for the Kishon River (Kishon system), Ein Feshkha (Dead Sea sytem), and Nitzanim (coastal introduced) samples. The two fish farm specimens also possessed this haplotype. The Ein Feshkha and Ne'ot HaKikar samples (both from the Dead Sea system) were monomorphic, but for different haplotypes (Table 1; Figure 2). The Kinneret sample (Jordan River system) showed greater haplotype diversity than any other (h = 0.69, sd = 0.12). As a consequence of the high haplotype diversity, the effective size estimation of this population is also higher than all the others (θ_s = 2.11, sd = 1.05), with the Kinneret being the largest most continuous body of water examined. Among the artificial ponds' samples, that from Ramat Gan was the most diversified (h = 0.56, sd = 0.09). With the exception of the monomorphic Dead Sea samples, the Kishon River sample provided the smallest θ_s estimation (θ_s = 0.35, sd = 0.35) (Table 1).

When including both natural and artificial habitat samples, AMOVA revealed that within-population variation explained about half (48.8 %) of the total genetic variation. The other approximate half of the variation (48.3 %) was explained by variation among populations within the various aquatic systems (Φ_sc = 0.50, P < 10^-5). Variation among aquatic systems was thus small (2.9 %) and non significant (Φ_ct = 0.03, P =0.35). Pairwise Φ_st indices as well as Nei’s corrected average numbers of pairwise difference (Table 2) indicated that samples from artificial ponds (Nitzanim and Ramat-Gan) significantly diverged from all other samples and from one another, including those from the coastal system (i.e., the natural system within which the artificial ponds were located). The sample from Ne'ot HaKikar (Dead Sea system) appeared to differ from the second Dead Sea sample (Ein Feshkha) but not from some samples from the Jordan River and coastal water systems. Taken together with its recent discovery, the population of Ne'ot HaKikar is suspected of being a recent introduction. These findings are reflected in the haplotype network (Figure 2).

When the three samples that are presumed to have been artificially introduced (Ne'ot HaKikar, Nitzanim and Ramat-Gan) were excluded from the AMOVA, a significant difference among some of the systems emerged (Φ_ct = 0.58, P = 0.013). Analysis of Nei’s corrected average numbers of pairwise differences and pairwise Φ_st indices revealed that the Dead Sea sample of Ein Feshkha differed from the others (2.08 < D < 2.61, P < 10^-5). This was also the case for the Kishon River (1.46 < D < 2.08, P <10^-5). However, the divergence between the Jordan River and coastal systems was not significant (D = 0.02, P = 0.25).

In the nested MANOVA of 22 meristic characters (Table 3) the factors "population within aquatic systems" and "aquatic system" significantly partitioned the variance of the meristic counts (F₁₀₅, ₅₃₈ = 4.19, P < 10^-5). MANOVA results were not affected by the inclusion or removal of samples from the introduced populations. A scatterplot of the first and second discriminant functions (accounting for 36% and 25.3% of the variation, respectively) clearly separated the Ein Feshkha sample from all others. In this plot (Figure 3), individuals of the coastal sample from Ein Afek and the (probably introduced) sample from Ne'ot HaKikar partly overlap. The highest correct classification rates, according to meristic counts, were of the Ein Feshkha sample (93.8%) and the Kishon sample (90%), while the percentage of correct classification of the other samples ranged from 56 to 83. Two characters were found to influence the first and second discriminant functions, respectively: the position of the longest pectoral ray; and the number of scales between the lower lateral line and the anal fin. The mean values and standard deviation of these characters in each sample are presented in Table 4..

The data sets supporting the results of this article are available in GenBank under accession numbers EU163705-EU163723 and FJ613474-FJ613479 and in the Dryad Repository: http://dx.doi.org/10.5061/dryad.29g768sc.

The research was approved by the Israel Nature and Parks Authority (permits 2006/25631, 2007/28858, 2008/31863). Fish were sampled using a seine net. Seven habitats were sampled: a) two coastal system habitats (Ein Afek on the Na'aman River, one sampling site; and the Taninim River, two sampling sites), b) the Kishon River (one sampling site), c) two habitats of the Jordan River system (Lake Kinneret, two sampling sites; and the Beit She'an Valley springs, three sampling sites), d) two Dead Sea system habitats (the Ein Feshkha Nature Reserve, two sampling sites; and the Ne'ot HaKikar Salt Flat, one sampling site). The sampling sites of each water body were at most five kilometers apart. Due to the small sample sizes, data from the same water body were pooled (i.e., they were treated as a single population) in the statistical analyses (AMOVA - Analysis of Molecular Variance and MANOVA - Multivariate Analysis of Variance). Consequently, our design did not allow us to analyze population structure within water bodies, although we cannot exclude the possibility that such a structure exists. We test here only the existence of structure among water basins.

Two artificial recreational ponds were also sampled (the Ramat-Gan park pond and the Nitzanim Nature Reserve pond, Figure 1). These artificial ponds are in the coastal basin. In addition, we obtained two T. zillii individuals from the HaMa'apil fish breeding farm. Unfortunately, other breeding farms did not allow us to examine their fish stock. For the morphological analysis, we also included specimens from the Steinhardt National Collection of Natural History, Zoological Museum at Tel-Aviv University (hereafter NCNH). Sample sizes and mean standard body lengths are provided in Table 4.

A liver sample was removed from each specimen collected (Table 4) and fixed in 70% ethanol for DNA extraction. The fish bodies were also preserved in 70% ethanol for meristic counts and submitted to the NCNH. While the individuals sampled were used for both morphological and genetic analysis, we extended our morphological sample size with the addition of fish from the NCNH. Unfortunately, the museum specimens were found to be unsuitable for DNA extraction. Before including these data in our morphological analysis, we verified that the date of capture and methods of conservation were not significant factors, using the non-parametric MANOVA test [¹⁶] as implemented in the software PAST [¹⁷]. Sample sizes are provided in Table 4. Further sampling would have potentially put the small studied populations at risk.

For each individual, 20–40 mg of liver tissue preserved in ethanol was homogenized with 500 μl TNES-Urea buffer [¹⁸] and 12.5 μl of proteinase K for 3 hours at 37°C. Following homogenization, DNA was extracted using a standard phenol-chloroform protocol followed by ethanol–sodium acetate precipitation. We selected the mitochondrial control region (D – loop) as our molecular marker, as in other studies of cichlid populations [^19-21]. The control region was amplified using the primers Ormt 449 up [²²] and a newly designed primer mit-tRNA-phe (5'-AGGGYCYATCTTAACATCTTCAGTG-3'). Amplified fragments were directly sequenced (to avoid cloning artifact) on an ABI PRISM 3100 (Applied Biosystems) genetic analyzer. The fragments were sequenced on both strands. Additionally, each polymorphic site and indel was confirmed manually by looking at the corresponding region on the chromatogram. Haplotypes were submitted to GenBank under accession numbers EU163705-EU163723 and FJ613474-FJ613479. Sequences were aligned using ClustalX [²³].The resulting alignment was then refined manually. The final alignment included 862 positions with two indels [Additional file 1.

The software Arlequin 3.0 [²⁴] was used to compute several diversity indices (Table 1): the haplotype diversity h[²⁵], the nucleotide diversity π[²⁵,²⁶] and two estimators of θθ_s, and θ_π. θ is the effective size of the population multiplied by the mutation rate. The θ_s estimate of θ is based on the relationship between the number of segregating sites and the sample size [²⁷], while θ_π is the mean number of pairwise differences in the sample [²⁶]. Both estimators assume the infinite site model. Since the mutation rate for the T. zillii control region sequence in unknown, the θ estimators can only be used to infer the relative effective sizes of the populations and consequently the relative effect of drift. Arlequin 3.0 was also used to conduct an AMOVA and to compute fixation index estimates (Φ estimates) [²⁸]. The total molecular variance was partitioned into four "aquatic systems" (i.e., coastal, Kishon, Jordan River and Dead Sea), "Populations within aquatic system" and "within population" components. The sample sizes did not allow the examination of population structure within water bodies. Two analyses were conducted. In the first, samples from artificial habitats (i.e., Ramat-Gan and Nitzanim) were considered as coastal, since these habitats are situated within the coastal basin, and the Ne'ot HaKikar sample, which has been suspected of having an aquacultural origin, was considered to belong to the Dead Sea system. In the second, samples of populations, which were probably introduced (i.e., Ramat-Gan, Nitzanim and Ne'ot HaKikar) were excluded. The HaMa'apil sample was excluded from the AMOVA analysis due to its small size (n = 2). Pairwise Φ_st values were tested with an exact contingency test as implemented in Arlequin 3.0. When aquatic systems or populations were found to be significantly different, Nei’s corrected average numbers of differences D[²⁹] were computed between systems or populations with the program Arlequin, using Jukes-Cantor distances. Finally, a median-joining network [³⁰] was reconstructed with the program Network 4.2.0.1 (available at http://www.fluxus-engineering.com). In all analyses, gaps were treated as a fifth character.

Twenty-two meristic characters were counted manually (Table 3). The data were logarithmically transformed to eliminate size-related differences [³¹]. The meristic counts were analysed using a nested MANOVA to examine whether the variation of morphological traits was significantly partitioned by geographical factors (i.e., among populations and among aquatic systems) discriminant analysis was used to visualize the morphological variation among populations [^32-34]. These analyses were conducted with the software Statistica (StatSoft) and PAST. Morphological data were deposited in the Dryad Repository: http://dx.doi.org/10.5061/dryad.29g768sc.

From left to right: haplotype number (Hap.), the polymorphic sites in the Control region sequence (polymorphic sites), the number of occurrences of each haplotype in each sample. Diversity parameters : n = sample size, h = haplotype diversity, π = nucleotide diversity, Ө_s and Ө_π -estimates of population effective size multiplied by the mutation rate. *Introduced population.

Above diagonal: Nei’s corrected average numbers of pairwise differences. Below diagonal: pairwise Φ_st indices. *P < 0.05, **P < 0.01, ***P < 0.001 and ^ns non-significant. Introduced populations are indicated by an asterisk.

* Introduced population, SL = mean standard body length, sd = standard deviation of standard body length. HaMa'apil sample was not measured due to small sample size. TR3 = Number of scales between the lower lateral line and the anal fin. LPR = The position of the longest pectoral ray