The fish used in the phylogenetic study (Table 1, Fig. 1) were collected from different locations in the sebkha: holes 71 (23.273778N, 15.917255W), 121 (23.272678N, 15.921623W), 129 (23.273244, 15.920921W), 130 (23.273296N, 15.920892W), and 148 (23.273961N, 15.920653W); the fish used in the microsatellite study were collected from four localities: water holes 35 (23.276394N, 15.915238W), 37 (23.276436N, 15.915553W), and 121, and the swamp (23.282684N, 15,909980W). The authorization to collect was provided by the “Haut Commissariat aux eaux et forêts et à la lutte contre la désertification” (HCEFLCD-Maroc). All fishing and tissue processing was carried out under existing collaborative arrangements between the Scientific Institute (Mohammed V University) in Rabat and the NGO Nature Initiative in Dakhla for the purposes of biodiversity research. As there is no Institutional Animal Care and Use Committee or equivalent animal ethics committee in Morocco, no formal approval was obtained for this work. Hundreds of fish live in colonized water holes (84 of the 161 existing water holes host either fish, crustaceans or molluscs), meaning that this species is not sensitive to light sampling of a few individuals. Nevertheless, to minimize the impact on these populations, we only took specimen fin clips and preserved them in a solution of 90% ethanol. To minimize the suffering of the individuals studied, live specimens were captured using seine nets and immediately anaesthetized with MS-222 (Tricaine-S, Western Chemicals Inc) before fin clip extraction. After sampling, the fish were kept for one hour in a bucket protected from light to allow them to recover, before re-introducing them into their natural habitat. The fish sampled in Senegal (Hann Bay off the coast of Dakar) or in Côte d’Ivoire (Layo in the Ebrié Lagoon; Grand Lahou; San Pedro and Fresco on the Sassandra River; the Marahoué River) were collected from the catches of local fishermen.

Approximately 50 mg of the sample tissue (fin clip) were sheared into fine pieces before being digested at 55 °C overnight using 10 mL proteinase K (10 mM/mL) in 190 mL of an extraction buffer solution (1 M Tris, 0.5 M NaCl₂, 1% SDS). DNA was then extracted from each sample following the protocol described for genomic DNA extraction for PCR-based techniques [8]. The DNA extracted was suspended in sterile double distilled water and stored at –20 °C until PCR amplification.

A 1147-bp fragment containing the complete ND2 gene was amplified in each sample using two primers: ND2F 5′-CATACCCCAAACATGTTGGT-3′ (forward) and ND2R 5′-GGAGATTTTCACTCCCGCTTA-3′ (reverse). Amplifications were performed in a final volume of 50 mL containing 0.25 mM of MgCl₂, 0.2 mM of each dNTP, 1 mM of each primer, 5 mL of 10 × buffer and 10 units of Taq polymerase (Promega). The PCR reaction conditions were: 3 min of pre-heating at 94 °C followed by 35 cycles of 30 s at 94 °C, 30 s at 48 °C and 1 min at 72 °C, with a final elongation phase at 72 °C for 5 min.

Additional sequences obtained from GenBank were also included in the analysis (Table 1). Sequences from Coptodon zillii, a closely related species [9], sequences from Sarotherodon melanotheron and S. galilaeus, two other “tilapia” species, and sequences from Chromidotilapa guentheri and Tylochromis leonensis Cichlid species were added as an outgroup to root the trees according to [10]. The sequences were aligned manually using BioEdit 5.09 [11]. The aligned sequences were analyzed using the Maximum Likelihood (ML) Minimum Evolution (ME) and Distance Method (DM, the Kimura two-parameter model [12], followed by the Neighbor Joining methods [13]) using MEGA (Molecular Evolutionary Genetics Analysis) software version 5.1 [14]. Before starting the analysis, an evolutionary model for ML and ME was selected by MEGA 5.1 [15], using the Bayesian information criterion (BIC) [16]. The model with the lowest BIC scores was considered as best describing the substitution pattern. Supports for inferred clades were obtained through the non-parametric bootstrap [17] with 2000 replicates.

One microgram of total genomic DNA was isolated from the individuals and sent to GenoScreen, France (www.genoscreen.com) to develop microsatellite libraries through 454 GS-FLX Titanium pyrosequencing of enriched DNA libraries, as described in [18]. Briefly, total DNA was enriched for AG, AC, AAC, AAG, AGG, ACG, ACAT, and ATCT repeat motifs and subsequently amplified. The PCR products were purified and quantified, and GsFLX libraries were then constructed in accordance with the manufacturer's protocols (Roche Diagnostics) and sequenced on a GsFLX-PTP. The bioinformatics programme QDD [19] was used to filter for redundancy, resulting in a final set of sequences from which it was possible to design primers. Finally, of the 5347 sequences including a microsatellite motif, 128 primer sets were designed and a sub-group of 47 primer pairs was tested for amplification on 8 DNA samples. Primer sets were discarded if they failed to amplify or led to multiple fragments. Then, 24 microsatellite loci from validated primer sets were selected for a polymorphism study on 15 DNA samples. PCR amplifications were performed in 25 μL reactions containing 20 ng of template DNA, 1 × reaction buffer, 37.5 pmol of MgCl₂, 6 pmol of dNTP, 10 pmol of fluorescent primer, 10 pmol of primer, and 1 U Taq of polymerase (FastStart–Roche Diagnostics). The PCR cycling consisted of an initial denaturation at 95 °C for 10 min, followed by 40 cycles: denaturation at 95 °C for 30 s, annealing at 55 °C for 30 s, and extension at 72 °C for 1 min and final extension at 72 °C for 10 min. Each microsatellite amplification was diluted with H₂O (1:50), mixed with Hi-Di Formamide and GeneScan 500 LIZ Size Standard (Applied Biosystems). The fragments were separated using an Applied Biosystems 3730XL DNA Analyser. The alleles were scored using GeneMapper v5.0 software (Applied Biosystems). Finally, 16 markers were selected for population genotyping (Table 2).

Microsatellite data were checked for scoring errors due to stuttering or large allele dropout and the presence of null alleles, using MICRO-CHECKER software [20]. The intra-population genetic variability was measured by estimating the observed heterozygosity (HO), the expected heterozygosity (HE), and the inbreeding coefficient (FIS) using GENEPOP software version 3.4 [21]. FST and associated probabilities were estimated using a Markov chain method in accordance with [22]. The length of the Markov chain involved a burn-in period of 1000 iterations and 100 batches of 1000 iterations thereafter. A Factorial Correspondence Analysis (FCA) was carried out with GENETIX software version 4.05 [23] in order to investigate the relationships between individuals. This type of analysis can explain a maximum amount of genetic variation using a minimum number of factors and can also provide the means for visualizing the genetic relationships between populations. The GENETIX software version 4.05 [23] was also used to determine whether or not there was a correlation between geographical distances and genetic distances between in the different samples.

The clustering approach implemented in STRUCTURE v2.3.3 [24] was used to infer the number of potential population units within the sebkha. We chose the admixture model with correlated allele frequencies and, because our sampling scheme involved collecting many individuals from a few discrete distant locations [25], we used the sampling location as prior information [26]. We used default values for all the other parameters in the software. Each run consisted of a burn-in period of 10⁵ Markov chain Monte Carlo (MCMC) iterations, followed by 10⁶ MCMC iterations. We carried out 10 replicate runs for each value of the number (K) of clusters, set between 2 and 4 (i.e. the number of samples).

Eighty-nine new ND2 sequences (1025 bp) were obtained during this study (Table 1), representing only 13 different haplotypes, A to M (Fig. 3). Two were present in the Imlili population (A and B), two in the Hann population (C and D) and eleven in the Côte d’Ivoire populations (E to M). All these sequences were aligned with GenBank sequences from C. zillii, Sarotherodon melanotheron, S. galilaeus, Chromidotilapa guentheri, and Tylochromis leonensis, and formed a matrix composed of 96 different sequences. This sequence matrix contained 372 variable sites, from which 206 were parsimony informative.

According to the different values of the corrected Akaike Information Criterion (AIC) obtained with MEGA 5.1, the optimum model for sequence evolution was the HKY + G model [27]. This model was used in the ML and ME analyses. Phylogenetic relationships among all the haplotypes observed based on the ML, ME or DM methods were congruent.

Fig. 3 shows the ME consensus tree obtained, on which the bootstrap values observed with the three methods are indicated. When rooted using Chromidotilapia guentheri and Tylochromis leonensis the two Sarotherodon species (S. melanotheron and S. galilaeus) appeared to be the sister group to the Coptodon species. Within the Coptodon clade, Coptodon zillii from Cairo and C. zillii of unknown origin were grouped together, forming the sister group to all the C. guineensis haplotypes. Within the C. guineensis clade, haplotypes from Kouilou (Congo) and another of unknown origin appeared basal, while the different haplotypes from Côte d’Ivoire, on the one hand, and the Imlili and Hann haplotypes, on the other one, appeared more derived (forming a monophyletic group that presented more recent genetic divergence). All these groups were supported by high bootstrap values varying from 79% to 100%.

The numbers of alleles observed in the 16 loci studied in the different samples are summarized in Table 3. Three loci were monomorphic M11, M14, and M17. A total of 54 different alleles was detected at the 13 polymorphic loci. The largest number of alleles was detected for locus M-41 (eight alleles) and the lowest for the loci M-07, M-20, M-23 and M-38 (two alleles). The mean number of alleles per locus varied between 2.9 in the swamp sample and 3.5 in the water hole 37 sample.

Using MICRO-CHECKER software, no scoring error caused by stuttering, no large allele dropout and no null alleles were detected, except at locus M27 in the population from the swamp. In this population, an excess of homozygotes was observed, suggesting that one null allele may have been present in this locus. No other genotype frequencies departure from Hardy Weinberg expectations was observed.

Inter-population differentiation as measured by F_ST [23], and the geographical distances between them are presented in Table 4. All pairwise F_ST values were between 0.031 and 0.167 and statistically significant (p < 1%), indicating significant differentiation among populations. Moreover, and even if the small number of samples did not make it possible to achieve a statistical test, it was worth noting (Fig. 4) that there was a differentiation by distance pattern. The more distant the populations were, the more genetically differentiated they were, while on the contrary, the closest populations were the least genetically differentiated. A mantel test carried out using GENETIX [23] produced Z = 1495.69. The 23 possible permutations in the matrix (4! − 1) produced only one value ≥ to the observed Z value, indicating the existence of a differentiation signal for distance relationships. In terms of the mean number of migrants exchanged per population at each generation, using the formula of Wright (1969) formula (N_m = (1 – F_ST)/4·F_ST), these numbers vary between 0.122 and 7.89.

With the programme GENETIX [23], we carried out a three-dimensional FCA (Fig. 5) that helped visualize these differentiations. Individuals from the swamp sample appeared differentiated from those one of the three other samples on the first axis. Individuals from the water hole 121 were mainly differentiated from individuals from the two other water holes, 35 and 37, on the second axis, while individuals from these two water holes were partially separated on the third axis.

STRUCTURE analyses [24] of all the individuals sampled within the sebkha provided consistent results over the 10 runs tested for each K value (Fig. 6). The natural logarithm for the likelihood of the data ln(P(X|K)) increased from K = 2 (–2350.7) to K = 4 (–2248.0), which was the maximal. For K = 2, all the individuals from the swamp and water hole 121 were clearly distributed in two different populations (Fig. 6), while the individuals from water holes 35 and 37 were not assigned to either of them. For K = 3, the individuals from the swamp and water hole 121 were still assigned to two different populations as well as the individuals from water hole 35, which constituted a third population. The individuals from water hole 37 were still not assigned to any of these three populations. When K was set to 4, all the individuals from the swamp and water holes 121, 35 and 37 were respectively assigned to different populations.