Phenotypic covariation among suites of traits may constrain or promote diversification both within and between species, yet few studies have empirically tested this relationship. In this study, we investigate whether phenotypic covariation of craniofacial traits is associated with diversification in an adaptive radiation of pupfishes found only on San Salvador Island, Bahamas (SSI). The radiation includes generalist, durophagous, and lepidophagous species. We compared phenotypic variation and covariation (i.e., the P matrix) between (1) allopatric populations of generalist pupfish from neighboring islands and estuaries in the Caribbean, (2) SSI pupfish allopatric lake populations with only generalist pupfish, and (3) SSI lake populations containing the full radiation in sympatry. Additionally, we examine patterns observed in the P matrices of two independent lab-reared F2 hybrid crosses of the two most morphologically distinct members of the radiation to make inferences about the underlying mechanisms contributing to the variation in craniofacial traits in SSI pupfishes. We found that the P matrix of SSI allopatric generalist populations exhibited higher levels of mean trait correlation, constraints, and integration with simultaneously lower levels of flexibility compared to allopatric generalist populations on other Caribbean islands and sympatric populations of all three species on SSI. We also document that while many craniofacial traits appear to result from additive genetic effects, variation in key traits such as head depth, maxilla length, and lower jaw length may be produced via non-additive genetic mechanisms. Ultimately, this study suggests that differences in phenotypic covariation significantly contribute to producing and maintaining organismal diversity.