Cumaná guppy inhabiting a stream within the Cariaco drainage has been shown to be distinct in terms of male sexual colouration, and partially reproductively isolated form a common guppy population inhabiting a stream in the San Juan drainage . The two drainages in north-eastern Venezuela are adjacent to each other, but are separated by a mountain range, Cordillera de la Costa, which may have formed a barrier to gene flow. It has been hypothesised that this barrier may have facilitated the phenotypic differentiation between guppy populations inhabiting separated drainages . To test this hypothesis, we studied guppy populations inhabiting these two drainages using both molecular (mtDNA and microsatellites) and morphometric approaches.
In accordance with the hypothesis that the Cordillera de la Costa constitutes a barrier to gene flow, the grouping of populations by drainage explained a significant proportion of genetic variation in both mtDNA and microsatellites. STRUCTURE results generally supported the distinction between the Cariaco and San Juan drainage guppies. In the analysis with two genetic clusters imposed, the populations clustered mostly by drainage (see Figure 2B), the exception being the populations CAS, GUAR and SF, located near the western and eastern borders of the Cariaco drainage. DAPC yielded similar results, with GUAR and CAS showing the greatest admixture of genes from the San Juan cluster, and all other populations clustering in accordance with drainage membership (see Figure 3A). Multidimensional scaling of the matrix of pairwise F
ST showed a similar pattern of separation between drainages, except one Cariaco drainage population (CAS) which is separated from San Juan drainage cluster along the first dimension (see Additional file 2: Figure S1).
Microsatellite differentiation between the Cariaco and San Juan drainages (7.8% for all populations sampled, and 5.93% for populations sampled in 2011) was larger than that reported by Suk and Neff  for differentiation between guppies in the Caroni and Oropuche drainages in Trinidad (4.6%), a level of differentiation proposed by Schories et al. to represent separate Poecilia species. For comparison, in the smooth (Lissotriton vulgaris) and Carpathian (L. montandoni) newts which constitute well defined species, 7.79% microsatellite variation was attributed to interspecific differentiation , a value comparable to the one reported here. However, different races of the genus Heliconius, which are considered examples of incipient speciation, were found to differ by over 30% at microsatellite loci . These examples imply considerable variation in the degree of differentiation at microsatellite loci among incipient/young species. However, microsatellites generally have not supported clustering together of populations according to the region (see Additional file 3: Figure S2), providing no evidence of species status to populations from the Cariaco region.
In the mtDNA phylogeny, two main lineages were distinguished; one composed exclusively of haplotypes from the Cariaco drainage, the other comprising mostly haplotypes found in guppies from the San Juan drainage, but containing also a well-supported Cariaco sub-lineage B. Both the Cariaco drainage lineage C and sub-lineage B received a high bootstrap support, indicating its distinctiveness and substantial divergence from the San Juan lineage (Figure 4), as reported by Alexander et al. for the mtDNA control region. Assuming mutation rate typical for many vertebrates (divergence of 2% per million years; ), the net divergence between the two main lineages of haplotypes (0.035) would be consistent with the 600,000 years of separation between the two drainages starting with the uplifting of Cordillera de la Costa, although this is only a rough estimate and should be treated with caution.
While mtDNA and microsatellite divergence estimates reported above support the hypothesis of considerable isolation and differentiation between guppies inhabiting streams on the opposite sides of the Cordillera de la Costa, STRUCTURE results do not corroborate such a clear distinction. The Evanno et al. method indicated eight, not two, as the most probable number of genetic clusters. While this result should be treated with caution, as STRUCTURE is not well suited for the analysis of data exhibiting hierarchical population structure, DAPC yielded the same optimal number of clusters. Furthermore, a greater number of clusters than two is consistent with the results of AMOVA, which showed that most genetic variation is distributed between populations within drainages. The fact that the highest degree of admixture was observed in populations near the borders of the drainages (CAS, GUAR, SF, see Additional file 5: Figure S5) strongly indicates that this pattern is due to gene flow rather than homoplasy. Gene flow is further corroborated by mtDNA introgression, which shows similar geographic pattern.
Neighbouring populations and populations sampled from within the same stream mostly fall into the same cluster (e.g. CAS and GUAR, both sampled from Rio Casanay, and VL and ACH or PA and RSJ, population pairs sampled from adjacent streams). Also populations CC, CCAR and WC, which are located in the Cumaná flood plain, were assigned to a common cluster. This pattern is reflected in significant, although weak, isolation by distance.
Interestingly, the population CUM, (ca. 1 km distance from CC which was sampled in 2011), clusters with a more distant population PER, which was sampled during the same season as CUM, 10 years after the first three populations were sampled. This genetic similarity between geographically more distant populations, but temporally closer sampling events indicates that temporal changes in genetic composition of populations may be rapid and have a considerable effect on population structure of guppies in north-eastern Venezuela. Recent simulations  showed that large temporal allele frequency fluctuations can arise as a consequence of relatively small sizes of guppy populations . Despite this, drainage explained a similar percentage of variation in microsatellites when samples from 2001 were excluded, and for mtDNA the structure was even more pronounced in this limited dataset. Therefore, the overall picture of significant differentiation between drainages was not biased due to a temporal gap in sampling events from the Cariaco drainage.
Surprisingly, SF located at the western border of the region did not fall into any group composed of Cariaco drainage populations in the STRUCTURE analysis, and instead shared the cluster with the geographically distant San Juan drainage populations. A possible explanation could be that introgression from a western common guppy population occurred; indeed, individuals with the haplotype from haplogroup A, typical of San Juan drainage guppies, were found in populations located west of SF (MH and JR, unpublished data).
Our results suggest a recent increase in introgression of San Juan genes into Cariaco drainage guppies (please note, that we are using the term introgression to indicate gene flow between genetically differentiated groups within a species; we do not suggest that we are dealing with two separate species). Two pairs of populations in our data, sampled from very close locations, YCUAL/SF and GUAR/CAS (in each pair sampled in 2001 and 2011, respectively), corroborate the idea of rapid invasion happening in few years. In the STRUCTURE analysis with two main genetic clusters, YCUAL showed more genetic similarity to guppies in the Cariaco drainage, while SF sampled 10 years later in close proximity showed a high proportion of the San Juan genes. In the second case, population GUAR, sampled in 2001, clustered to a great extent with San Juan populations but also showed considerable admixture with Cariaco genes, whereas CAS, sampled in 2011 from the same stream, grouped with the San Juan drainage populations. A similar situation has been reported for Trinidadian guppies, where substantial admixture was found between populations in adjacent streams, but belonging to different drainages (; based on 7 microsatellite loci). Furthermore, in populations at the borders of the Cariaco and San Juan drainages, CAS and GUAR, only mitochondrial haplotypes typical for San Juan drainage were present, and one was also found in SF (see Additional file 1: Table S7, Additional file 5: Figure S5). This clear geographical pattern strongly suggests introgression, rather than incomplete lineage sorting, as an underlying mechanism. A recent SNP-based study also indicated introgression from Poza Azufre (SanJuan drainage) to Cumaná (see Figure 3 in ).
Our data also provide evidence that some introgression may have occurred in a more distant past. In a phylogram of mitochondrial haplotypes within the San Juan lineage (Figure 4), a Cariaco sub-lineage was clearly distinguished, supported by a high bootstrap value. The presence of this lineage may be explained by a past introgression of a San Juan haplotype into Cariaco populations, which then evolved, in isolation from San Juan guppies, within these populations.
The reasons for the San Juan drainage guppy’s invasion into the Cariaco drainage are not clear. The pattern of introgression, most pronounced near the borders of the Cariaco drainage, suggests the role of natural processes and not of human-induced introductions. We are not aware of any abrupt environmental changes in this area, but occasional floods occurring in the region may favour occasional migration events . AMOVA results based on microsatellite loci indicate the proportion of variation explained by drainage was nearly twice as high as that estimated in the analysis in which the highest level of hierarchy was defined by mtDNA haplotypes. This suggests that when the geographic barrier is crossed, reproductive barriers do not prevent introgression of mtDNA.
The analysis of body shape and colour patterns that distinguish the Cumaná guppy from other guppy populations  showed that only the population CC represented all the characteristics which phenotypically define the Cumaná male morphotype (presence of double sword, presence of black crescents, larger tail orange area and narrower caudal peduncle, relative to typical P. reticulata). Some other populations from Cariaco drainage exhibited individual Cumaná characteristics, but not their complete combination (see Figures 3B and 5). Individuals from the WC and CUM populations, located near the CC population, had an increased number of black crescents, while the WC population revealed an increase in the mean orange tail area and, together with CAS, a relatively high occurrence of the double sword, compared to other guppy populations measured. The only trait which significantly differed between the Cariaco and San Juan drainages was the presence of black crescents. Thus overall, there is little evidence that most guppies inhabiting Cariaco region can be equated with Cumaná morphotype, as suggested by Poeser . Indeed, most individuals from CC population were assigned to cluster 1 in the DAPC based on phenotypic traits, whereas vast majority of individuals from other populations from both San Juan and Cariaco drainages, were assigned to cluster 2 (see Figure 3B). Consequently, there was no significant difference between drainages in probability of assignment to a drainage. This result clearly shows that morphological differences are not aligned with a geographic barrier between two drainages. Rather, the CC population represents the typical Cumaná morphotype, with other traits, except for black crescents, occurring only in some populations in the region.
That we found no evidence for original designation of P. wingei according to drainage  does not exclude the possibility that Cumaná population could be genetically and morphologically distinct enough to be considered a different species. However, our results do not provide evidence for genetic differentiation, as CC population, representing Cumaná morphotype, clustered together in the STRUCTURE with CCAR which represents the common guppy morph (see Figures 2A and 5). Furthermore, phenotypic analyses showed that populations in the San Juan drainage also contained individuals with some ‘Cumaná traits’, for example, double swords were observed in the Poza Azufre and Rio San Juan populations, and high mean orange area in Rio Juan Sanchez population. This pattern is visible on the graphical representation of DAPC, where individuals belonging to the cluster 1, dominated by the ‘Cumaná-type’ phenotypic characteristics, occur sporadically throughout all populations (see Figure 3B). Likewise, CC population contained individuals assigned to 'common guppy' cluster 2 (Figure 3B). Overall, our data do not support existence of two species which can be easily separated based on male morphology. The lack of support for the existence of two clusters in DAPC analysis of morphological traits further strengthens this conclusion.
There are two possible explanations for the patterns we observed. First, the Cumaná city area constituted the source of Cumaná morphotype, and migrants from this area spread widely throughout the Cariaco drainage. The second is that the whole Cariaco drainage was the original range of the Cumaná type guppy, but the populations located closer to the border with the San Juan guppies have been invaded by immigrants from the San Juan drainage.
As most genetic variance was distributed among populations rather than between drainages, it seems that limited gene flow between different streams, to a greater degree than between drainages, may have facilitated the evolution of divergent female preferences of females form Cumaná . Apart from limited gene flow, genetic divergence between populations may be due to drift, as effective sizes of many guppy populations are <1000 . A recent theoretical model showed that under such conditions, speciation by sexual selection is particularly likely .
Once divergent mating preferences evolve, they may constitute a powerful mechanism maintaining divergence in sexual traits. Genetic clustering of populations from the Cumaná flood plain (Figure 2A), representing both Cumaná (CC population) and common guppy morphs (CCAR) suggests that female choice has imposed differential rates of introgression among genes that do or do not encode for traits related to biological divergence. Divergent female preferences have also been shown to facilitate reproductive isolation among guppy populations differing in the strength of selection from predators . Thus, although sexual selection can under some scenarios also constrain speciation , our results are consistent with a traditional view that sexual selection may facilitate the emergence of reproductive isolation . However, we have also found that phenotypic traits typical for Cumaná morphotype are present in other populations and vice versa, which indicates that the barrier to the flow of genes associated with traits subject to mating preferences is incomplete. More intensive sampling effort in Cumaná region will be necessary to determine if morphological homogenisation with respect to male sexual traits has been occurring in consequence.