Craniofacial disorders represent some of the most common birth defects affecting ¼ of the human population worldwide. Many causal mutations implicated in craniofacial malformations are also responsible for generating the natural craniofacial trait variation observed between species. The advancement of next generation sequencing technologies has allowed for rapid acquisition and identification of genome-wide variants mediating adaptive craniofacial differences between species. Rock dwelling and sand dwelling cichlid fishes of Lake Malawi (LM) exhibit extensive variations in craniofacial morphologies, specifically in jaw size and length making them an outstanding model for the study of craniofacial morphological evolution. Here, we performed whole-genome re-sequencing of 16 rock and sand dwelling LM cichlids and identified regions of genomic differentiation between the groups. We show evidence that (1) regions of genomic differentiation are enriched for genes involved in pathways and processes relating to cranial neural crest cell differentiation and development (2) alternately fixed SNPs between rock and sand species are within and around crucial genes of the cranial neural crest cell gene regulatory network and (3) alternately fixed SNPs between rock and sand species are also within enhancer elements regulating craniofacial and neural crest development. Taken together, our comparative genomic analysis provides a global picture of the key genomic changes likely modulating species-specific craniofacial differences between rock dwelling and sand dwelling LM cichlids.