Detection of footprints of historical natural selection on quantitative traits in cross-sectional data sets is challenging, especially when the number of populations to be compared is small and the populations are subject to strong random genetic drift. We extend a recent Bayesian multivariate approach to differentiate between selective and neutral causes of population differentiation by the inclusion of habitat information. The extended framework allows one to test for signals of selection in two ways: by comparing the patterns of population differentiation in quantitative traits and in neutral loci, and by comparing the similarity of habitats and phenotypes. We illustrate the framework using data on variation of eight morphological and behavioral traits among four populations of nine-spined sticklebacks (Pungitius pungitius). In spite of the strong signal of genetic drift in the study system (average FST = 0.35 in neutral markers), strong footprints of adaptive population differentiation were uncovered both in morphological and behavioral traits. The results give quantitative support for earlier qualitative assessments, which have attributed the observed differentiation to adaptive divergence in response to differing ecological conditions in pond and marine habitats.