In evolutionary quantitative genetics, the genetic variance–covariance matrix,
G
, and the vector of directional selection gradients,
β
, are key parameters for predicting multivariate selection responses and genetic constraints. Historically, investigations of
G
and
β
have not overlapped with those dissecting the genetic basis of quantitative traits. Thus, it remains unknown whether these parameters reflect pleiotropic effects at individual loci. Here, we integrate multivariate genome-wide association study (GWAS) with
G
and
β
estimation in a well-studied system of multivariate constraint: sexual selection on male cuticular hydrocarbons (CHCs) in
Drosophila serrata
. In a panel of wild-derived re-sequenced lines, we augment genome-based restricted maximum likelihood to estimate
G
alongside multivariate single nucleotide polymorphism (SNP) effects, detecting 532 significant associations from 1 652 276 SNPs. Constraint was evident, with
β
lying in a direction of
G
with low evolvability. Interestingly, minor frequency alleles typically increased male CHC-attractiveness suggesting opposing natural selection on
β
. SNP effects were significantly misaligned with the major eigenvector of
G
,
g
max
, but well aligned to the second and third eigenvectors
g
2
and
g
3
. We discuss potential factors leading to these varied results including multivariate stabilizing selection and mutational bias. Our framework may be useful as researchers increasingly access genomic methods to study multivariate selection responses in wild populations.
Evolution of a predator-induced, nonlinear reaction norm