Amino acid changes in mitochondrial (mt) oxidative phosphorylation (OXPHOS) genes have been suggested as a key adaptation to environmental variation. Here, we analyzed 416 sequences of ATPase synthase 6 (MT-ATP6) and NADH dehydrogenase 2 (MT-ND2) in 22 different hare (Lepus) species from across a wide range of habitats and climates. We used site- and branch-based methods to test for positive selection on specific codons and lineages. We found four codons in MT-ATP6 and five in MT-ND2 under positive selection, affecting several species lineages. We investigated the association of protein variants at each locus with climate zone, using multinomial generalized linear models (glm), including species, regions, historical introgression events, and the co-occurring protein variant at the other locus as additional explanatory variables. A significant climate effect as based on the Koppen climate classification was observed for MT-ND2 protein variants as translated from our nucleotide sequences. Moreover, MT-ND2 protein variants were significantly affected by the co-occurring MT-ATP6 protein variant in the same mtDNA molecule. Contrary to the expectation for non-recombining mitochondrial DNA molecules, the presence of an evolutionarily relatively ancestral protein variant at one locus was associated with a relatively derived protein at the other locus in the same mitochondrial molecule, respectively. The relative evolutionary status of a protein variant was evaluated according to its positions relative to the respective out-group protein variant in a network analysis of nucleotide sequences. All our results suggest a complex effect of various climatic parameters acting on multiple mtOXPHOS genes in a co-adaptive way, favoring combinations of ancestral and derived variants.