Divergent selection at ecologically important traits is thought to be a major factor driving phenotypic differentiation between populations. To elucidate the role of different evolutionary processes shaping the variation in gill-raker number of European whitefish (Coregonus lavaretus sensu lato) in the Baltic Sea basin, we assessed the relationships between genetic and phenotypic variation among and within three whitefish ecotypes (sea-, river and lake-spawners). To generate expected neutral distribution of FST and to evaluate whether highly variable microsatellite loci resulted in deflated FST estimates compared to less variable markers we performed population genetic simulations under finite island and hierarchical island models. The genetic divergence observed among (FCT=0.010) and within (FST=0.014–0.041) ecotypes was rather low. The divergence in gill-raker number, however, was substantially higher between sea- and river-spawners compared to observed microsatellite data and simulated neutral baseline (PCT>FCT). This suggests that the differences in gill-raker number between sea- and river-spawners are likely driven by divergent natural selection. We also found strong support for divergent selection on gill-raker number among different populations of sea-spawners (PST>FST), most likely caused by highly variable habitat use and diverse diet. The putative role of divergent selection within lake-spawners initially inferred from empirical microsatellite data, was not supported by simulated FST distributions. This work provides a first formal test of divergent selection on gill-raker number in Baltic whitefish, and demonstrates the usefulness of population genetic simulations to generate informative neutral baselines for PST–FST analyses helping to disentangle the effects of stochastic evolutionary processes from natural selection.