Compatibility between the nuclear (nDNA) and mitochondrial (mtDNA) genomes is important for organismal health. However, its significance for major evolutionary processes such as speciation is unclear, especially in vertebrates. We previously identified a sharp mtDNA-specific sequence divergence between morphologically indistinguishable chameleon populations (Chamaeleo chamaeleon recticrista) across an ancient Israeli marine barrier (Jezreel Valley). Since mtDNA introgression and gender-based dispersal were ruled out, we hypothesized that mtDNA spatial division was maintained by mito-nuclear functional compensation. Here, we studied RNA-seq generated from each of ten chameleons representing the north and south populations and identified candidate non-synonymous substitutions (NSS) matching the mtDNA spatial distribution. The most prominent NSS occurred in 14 nuclear DNA-encoded mitochondrial proteins. Increased chameleon sample size (N=70) confirmed the geographic differentiation in POLRMT, NDUFA5, ACO1, LYRM4, MARS2 and ACAD9. Structural and functionality evaluation of these NSS revealed high functionality. Mathematical modeling suggested that this mito-nuclear spatial divergence is consistent with hybrid breakdown. We conclude that our presented evidence and mathematical model underline mito-nuclear interactions as a likely role-player in incipient speciation in vertebrates.