Comparing the patterns of population divergence using both neutral genetic and phenotypic traits provides an opportunity to examine the relative importance of evolutionary mechanisms in shaping population differences. We used microsatellite markers to examine population genetic structure in the Mexican sailfin mollies Poecilia velifera and P. petenensis. We compared patterns of genetic structure and divergence to that in two types of phenotypic traits: morphological characters and mating behaviours. Populations within each species were genetically distinct, and conformed to a model of isolation by distance, with populations within different geographical regions being more genetically similar to one another than were populations from different regions. Bayesian clustering and barrier analyses provided additional support for population separation, especially between geographical regions. In contrast, none of the phenotypic traits showed any type of geographical pattern, and population divergence in these traits was uncorrelated with that found in neutral markers. There was also a weaker pattern of regional differences among geographical regions compared to neutral genetic divergence. These results suggest that while divergence in neutral traits is likely a product of population history and genetic drift, phenotypic divergence is governed by different mechanisms, such as natural and sexual selection, and arises at spatial scales independent from those of neutral markers.