The identification of the genes involved in morphological variation in nature is still a major challenge. Here we explore a new approach: we combine 178 samples from a natural hybrid zone between two subspecies of the house mouse (Mus musculus domesticus and Mus musculus musculus), and high coverage of the genome (~145K SNPs) to identify loci underlying craniofacial shape variation. Due to the long history of recombination in the hybrid zone, high mapping resolution is anticipated. The combination of genomes from subspecies allows the mapping of both, variation within subspecies and intersubspecific differences, thereby increasing the overall amount of causal genetic variation than can be detected. Skull and mandible shape were measured using 3D landmarks and geometric morphometrics. Using principle component axes as phenotypes, and a linear mixed model accounting for genetic relatedness in the mapping populations, we identified 9 genomic regions associated with skull and 10 with mandible shape. High mapping resolution (median size of significant regions = 148 kb) enabled identification of single or few candidate genes in most cases. Some of the genes act as regulators or modifiers of signaling pathways relevant for morphological development and bone formation, including several with known craniofacial phenotypes in mice and humans. The significant associations combined explain 13% and 7% of the skull and mandible shape variation. In addition, a positive correlation was found between chromosomal length and proportion of variation explained. Our results suggest a complex genetic architecture for shape traits, and support a polygenic model.This article is protected by copyright. All rights reserved.