Oxidative phosphorylation depends on the interaction of protein complexes in the inner mitochondrial membrane. In most animals, the four core complexes mediating OXPHOS are encoded by the 13 protein-encoding mitochondrial genes and an estimated 80 nuclear-encoded genes (Nicholls and Ferguson, 2013). The respirasome (comprising complexes I, III, and IV) uses the energy liberated during electron transfer from NADH to O₂ to drive proton-translocation to the intermembrane space and thus establish a proton gradient across the inner mitochondrial membrane (Figure 1) (Gu et al., 2016; Letts et al., 2016; Wu et al., 2016). The ATP synthase (complex V) uses the proton-motive force thus generated to chemiosmotically drive ATP synthesis (Allegretti et al., 2015; Hahn et al., 2016). Nuclear-encoded proteins also serve as electron carriers (e.g., cytochrome c), alternative electron inputs (e.g., complex II), and assembly factors (e.g., NDUFC1, SURF1) throughout the OXPHOS system (Mashkevich et al., 1997; Nicholls and Ferguson, 2013; Stroud et al., 2016).

Intimate associations between mitochondrial- and nuclear-encoded subunits are required for the electron transport chain of the respirasome and ATP synthase activity for efficient mitochondrial ‘coupling,’ i.e., the ratio of ATP synthesis per unit substrate and O₂ consumed (Lowell and Spiegelman, 2000). The intimacy of such interactions is particularly reflected by OXPHOS complex I. This enzyme couples electron transfer through its nuclear-encoded hydrophilic arm to proton translocation through its four mitochondrially encoded proton pumps. This depends on long-range conformational changes mediated through protein–protein interactions of the mitochondrially encoded subunits with core and supernumerary nuclear-encoded subunits (Fiedorczuk et al., 2016; Zhu et al., 2016). Genetic studies have demonstrated that assembly of this complex depends strictly on 39 of its 45 subunits (Stroud et al., 2016), and even single amino acid substitutions can alter coupling efficiency (Mimaki et al., 2012; Gershoni et al., 2014). In addition to serving as the primary electron input in the respiratory chain (Nicholls and Ferguson, 2013), OXPHOS complex I is also the main site of cellular reactive oxygen species (ROS) production (Murphy, 2009).

The central physiological importance of OXPHOS means that mitonuclear compatibilities are required for optimal fitness. Even minor biochemical inefficiencies can have major fitness consequences for an organism by modulating their energetic efficiency and oxidative stress levels. There are therefore strong selective pressures to maintain optimal mitonuclear interactions in the OXPHOS system (Rand et al., 2004; Dowling et al., 2008; Burton et al., 2013).

There are multiple lines of evidence that OXPHOS is important for local adaptation. Experimental approaches with model organisms have allowed researchers to test an impressive array of mitonuclear combinations and assess their functional effects under a wide set of conditions (Dowling et al., 2007; Arnqvist et al., 2010; Paliwal et al., 2014; Ma et al., 2016; Mossman et al., 2016). These efforts have revealed that mismatched mitonuclear interactions (gene–gene interactions) can have profound consequences, such as reduced metabolic performance, fecundity, and lifespan. When mitonuclear combinations are assessed in multiple environments (e.g., diet, temperature or hypoxia), interaction effects are commonly context dependent (gene–gene-environment interactions; Koevoets et al., 2012; Hoekstra et al., 2013).

In the wild, OXPHOS traits have been correlated with a wide range of environmental pressures, including heat stress (Morales et al., 2016a,b), cold stress (Cheviron et al., 2014; Stier et al., 2014a,b), nutrient limitation (da Fonseca et al., 2008), and hypoxia (da Fonseca et al., 2008; Scott et al., 2011). Consistently, there is evidence for positive selection and climate-linked differences in the sequences and expression of OXPHOS genes in a range of animal species with wide biogeographic ranges (Mishmar et al., 2003; Ruiz-Pesini et al., 2004; Toews and Brelsford, 2012; Garvin et al., 2015; Morales et al., 2016a). Accordingly, genes encoding the OXPHOS machinery are frequently candidates for positive selection. This likely reflects the high levels of mitochondrial DNA variation within and among populations combined with the selection pressures for optimally adapted phenotypes (Gershoni et al., 2009; Bar-Yaacov et al., 2012).

As mechanistic understanding of OXPHOS activity continues to improve, it should be increasingly possible to make more specific predictions of what kinds of protein-level changes should be adaptive. For example, it is hypothesized that the coupling efficiencies of OXPHOS complexes are closely linked to adaptive thermal biology. OXPHOS generates chemical energy and heat in proportions that depend on the coupling efficiency of the respirasome and ATP synthase (Lowell and Spiegelman, 2000). It is proposed that the heat produced from less-coupled mitochondria may be particularly beneficial for adaptation of endothermic organisms to colder environments. In contrast, heat production may be deleterious in warm environments, necessitate higher nutritional intake, and is associated with high oxidative stress due to increased ROS production (Figure 1) (Pörtner et al., 1998; Brand, 2000; Somero, 2002; Fangue et al., 2009; Stier et al., 2014a,b). It is important to account for variation in such predictions among organisms and environments. For example, contrary to the expectations for endotherms, cold adaptation in fishes is linked to higher mitochondrial densities in muscle (White et al., 2011), and so the associated high energy demands for synthesis and maintenance of mitochondria may favor genotypes with high coupling efficiency (Pavlova et al., 2017).

Mitonuclear co-evolution should be enforced under strong selection given the complex interactions and essential functions mediated by OXPHOS (Burton et al., 2013). Challenges to positive co-evolution include the fast mutation rate of the mitochondrial genome (due to its proximity to ROS production, high rate of replication and lack of efficient repair mechanisms), typically maternal inheritance, and lack of recombination, which generate a mutation load that the nuclear genome must counter by compensatory mutation (Rand et al., 2004; Lynch et al., 2006; Osada and Akashi, 2012; Havird et al., 2015; Havird and Sloan, 2016). In addition, mitonuclear co-evolution can be disrupted by mechanisms that generate genetic variation or promote gene flow. In particular, substantial gene flow and recombination in nuclear DNA will tend to break up optimally functioning allele combinations of co-adapted genes in each sexual generation (Rand et al., 2004; Burton and Barreto, 2012; Burton et al., 2013). Accordingly, we propose that nuclear genomic architecture should tend to evolve to suppress recombination and prevent the segregation of genome regions that mediate epistatic functions of nuclear-encoded mitochondrial genes.

With improvements in techniques to explore genome structure, examples are building of how genomic architecture can drive evolutionary adaptation, for example the ‘supergene’ of 125 genes associated with differences in male mating strategies in birds (Küpper et al., 2016). Natural selection can locate co-adapted loci in genome areas of low recombination (Schwander et al., 2014; Thompson and Jiggins, 2014) or promote genomic clustering of synergistically adaptive alleles (Yeaman, 2013) so that they can be co-inherited and/or co-regulated. The three-dimensional organization of the genome can dictate how and which loci should be subject to genome changes that will favor their co-location (Lanctôt et al., 2007; Wijchers and de Laat, 2011). Reduced recombination among co-adapted genes (increasing their co-inheritance) can occur through the evolution of recombination modifiers or chromosomal re-arrangements, such as transposition of a gene to a location close to a co-adapted gene, movement of co-adapted genes toward a centromere or into a region within an inversion between diverged lineages (Figure 1) (Rieseberg, 2001; Butlin, 2005; Kirkpatrick and Barton, 2006; Yeaman and Whitlock, 2011; Yeaman, 2013; Ortiz-Barrientos et al., 2016).

To date, mitonuclear genomic architecture (encompassing mitochondrial-nuclear and nuclear-nuclear) remains relatively underexplored, except in the context of biased co-transmission of mitochondrial and nuclear-encoded mitochondrial genes on sex chromosomes. Because mtDNA, due to its maternal inheritance, accumulates mutations that are deleterious in males (mother’s curse), selection to restore fitness in males that drives compensatory evolution of nuclear-encoded mitochondrial genes could be expected to prevent concentrations of such genes on female-linked chromosomes (Havird and Sloan, 2016). Results for different taxa variously support overrepresentation, underrepresentation, or unbiased distribution of nuclear-encoded mitochondrial genes on X and Z chromosomes with respect to autosomes, supporting multiple theories of the distribution of nuclear-encoded mitochondrial genes: co-adaptation, sexual conflict and sexual selection (Drown et al., 2012; Hill and Johnson, 2013; Dean et al., 2014, 2015; Hill, 2014; Rogell et al., 2014). More recently, we uncovered an autosomal genomic island of divergence associated with mitonuclear interactions in a passerine (Figure 1) (Morales et al., 2016a). This genomic island of divergence is implicated in maintaining deep mitochondrial divergence between two parapatric lineages connected by nuclear gene flow. Observations of the fluidity of positioning of nuclear-encoded mitochondrial genes in some systems (but see Dean et al., 2015) raise questions about the role of genome organization in rates of co-evolution between mitochondria and nuclear genes (Hill, 2014). We contend that exploration of genomic architecture may be crucial for understanding mitonuclear co-evolution, and also vice versa given the crucial roles and considerable number of genes concerned with mitochondrial function (Pagliarini et al., 2008).

Detection of natural selection is one of the most contentious and active fields in evolutionary biology. Here we highlight, in the context of mitonuclear co-evolution, the more general issues that are explored in depth elsewhere (Nielsen, 2005; Haasl and Payseur, 2016; Manel et al., 2016; Stephan, 2016).

The power to detect candidate loci that evolve under natural selection rests on the molecular tools available for a given system: reduced representation genomic scans (SNPs), sequence-based genomic scans (candidate genes, exome-sequencing, or RNA-sequencing), whole genome re-sequencing, and/or physical linkage maps (Manel et al., 2016; Stephan, 2016). The key limitation to detecting natural selection in the wild is that several ecological and evolutionary processes can leave a similar signature to selection and lead to a high rate of false positives. Confounding factors include demographic processes (e.g., population size change and structure), background selection, and heterogeneous mutation and recombination rates. Given that mitochondrial and nuclear genomes can have largely independent evolutionary histories (e.g., different introgression patterns and mutation load), knowledge of the demographic history of the study system is especially useful to interpret patterns of mitonuclear co-evolution (Bar-Yaacov et al., 2015; Morales et al., 2016a,b; Pereira et al., 2016; Sloan et al., 2016). Given that mitonuclear co-evolution is likely to respond to environmental variation (Burton et al., 2013), approaches to detecting selection that rely on gene-environment associations could be particularly useful to identify candidate loci under selection (Rellstab et al., 2015; Forester et al., 2016).

A common starting point in the search for signatures of selection in mitonuclear co-evolution is sequencing full mitochondrial genomes. A family of methods proven to be especially useful in the context of mitogenome evolution are codon-based approaches, which rely on the estimation of the non-synonymous to synonymous ratio (ω = d_N/d_S) [HyPhy and Datamonkey (Pond and Frost, 2005); PAML (Yang, 2007)]. There are multiple examples in the literature of how complementary codon-based approaches have been combined to discriminate positive and relaxed purifying selection in mitogenome-encoded OXPHOS components (da Fonseca et al., 2008; Garvin et al., 2011; Morales et al., 2015; Wertheim et al., 2015; Pavlova et al., 2017). However, these types of methods have important limitations: they require data across multiple species, or sequences that are reasonably diverged, and only allow selection inference within coding regions (Gonçalves da Silva, 2017). It is important to consider this last limitation since mitonuclear incompatibilities have been mapped to non-coding regulatory genes, non-coding sequences such as transfer RNAs and the mitochondrial control region (Montooth et al., 2009; Meiklejohn et al., 2013; Rollins et al., 2016; Jhuang et al., 2017).

A natural follow-up is to look for evidence of natural selection in nuclear-encoded mitochondrial genes and signals of mitonuclear co-evolutionary adaptation (Mishmar et al., 2006; Gagnaire et al., 2012; Bar-Yaacov et al., 2015; Pereira et al., 2016). For example, such approaches have identified supernumerary and assembly factors of OXPHOS complex I implicated in local adaptation (Garvin et al., 2016; Morales et al., 2016a).

Recent advances in understanding structure-function relationships in oxidative phosphorylation enable better prediction of how genetic substitutions affect mitochondrial function. Largely as a result of recent advances in cryo-electron microscopy, complete atomic-resolution structures of all components in the mammalian electron transport chain are now available, including the mitochondrially co-encoded complex I (Figure 1) (Fiedorczuk et al., 2016; Zhu et al., 2016), complex III (Iwata et al., 1998), complex IV (Tsukihara et al., 1996), and the respirasome supercomplex (Gu et al., 2016; Letts et al., 2016; Wu et al., 2016). In addition, near-complete structures of yeast ATP synthase are also available (Allegretti et al., 2015; Hahn et al., 2016).

With these newly available protein structures, it is now possible to map the locations of subunits and amino acids predicted to be under selection using protein visualization software (Pettersen et al., 2004) and to develop homology models using public servers (Källberg et al., 2012; Kelley et al., 2015). Such approaches have been used to predict the mechanistic effects of amino acid substitutions observed in OXPHOS subunits across diverse species (Scott et al., 2011; Finch et al., 2014; Caballero et al., 2015; Zhang and Broughton, 2015; Campana et al., 2016; Morales et al., 2016a). This approach was first highlighted by Garvin et al. (2011) who detected an amino acid under positive selection in the Pacific salmon in an unusual region of OXPHOS complex I: the piston-like horizontal helix (helix HL) of ND5. Meta-analysis has since suggested that this helix is the most common region of positive selection in the mitogenomes of diverse animal species (Garvin et al., 2014). While the function of the helix HL remains unresolved, it is hypothesized to influence coupling by propagating conformational changes from proximal to distal proton pumps; hence fine-tuning its properties may have adaptive consequences for heat and energy production (Torres-Bacete et al., 2011; Sazanov, 2014).

In animal systems, it is challenging to validate experimentally that certain amino acid substitutions affect mitochondrial function. Due to their membrane localisation, multi-subunit cofactor-bound composition, and complex assembly pathways, OXPHOS complexes are incompatible with recombinant protein expression and can rarely be purified natively (Nicholls and Ferguson, 2013). However, a suite of physiological techniques enable us to measure the activities, kinetics, and efficiencies of OXPHOS complexes. For example, classical respirometry techniques enable measurement of the rates of substrate oxidation or oxygen consumption by whole cells or purified mitochondria; it is possible to calculate mitochondrial coupling efficiencies and to probe the activities of specific protein complexes by systematically comparing basal respiration rates with those in the presence of specific agonists, inhibitors, and uncouplers (Figure 1) (Nicholls and Ferguson, 2013; Stier et al., 2013, 2016; Toews et al., 2014). Well-established assays using lysed cells also enable measurement of potentially relevant parameters such as the expression levels, protein content, and kinetic parameters of individual OXPHOS complexes (Scott et al., 2011).

Another development that enhances the ability to measure biologically relevant mitochondrial function is the recent discovery that non-mammalian animals harbor functional mitochondria in their erythrocytes. This presents options for non-destructive sampling of wild populations (Stier et al., 2013, 2015, 2016). These cellular measurements of mitochondrial respiration can be complemented with whole-organism measurements of basal and maximal metabolic rates; while rarely adopted in mitonuclear ecology, such approaches may have value for understanding relationships between nutritional intake and energy expenditure (White et al., 2007, 2011; Halsey and White, 2010).

Taking an approach amenable to experimental manipulations, laboratory-based crossings have also been used to assess the effects of intraspecific and interspecific mitonuclear compatibilities using individuals sampled from wild populations from different environments. As elaborated in a case study below, there are several examples of how crosses have been combined with measurements of enzymatic, cellular, or organismal performance to consolidate genotype-to-phenotype links (Edmands and Burton, 1999; Willett and Burton, 2001; Arnqvist et al., 2010; Chang et al., 2015; Dordevic et al., 2016). In the few study systems where cell lines can be established, experimental cellular approaches can provide valuable functional insights in study systems (Blier et al., 2006), for example through the construction of mitonuclear hybrid cell systems (cybrids) that allow testing of the effects of mitogenome variation on fitness in a constant nuclear background (e.g., Barrientos et al., 1998; Moreno-Loshuertos et al., 2006; Dingley et al., 2014). Within these kinds of manipulative approaches, as well as more broadly, rapid advances in molecular genomics, including the increasing tractability of RNA sequencing, are facilitating investigations of the roles, mechanisms and evolutionary genomics of genes of interest in adaptation and divergence of wild species (Harrisson et al., 2014; Havird and Sloan, 2016; Latorre-Pellicer et al., 2016).

The level of resolution that can be reached in examining genomic architecture depends on the genomic resources available. Full resolution requires assembled genomes and physical linkage maps, rarely available for wild organisms. However, powerful approximations can be made by mapping genetic variants of interest (e.g., candidate genes or loci under selection) on to a reference annotated genome of the same species or a close relative with known conserved synteny. Mapped variation provides the presumed order, position and identity of loci of interest (e.g., nuclear-encoded mitochondrial genes).

Linkage (gametic) disequilibrium (LD) is the non-random association of alleles at different loci within individuals. These correlations can arise through genes being near each other on a chromosome, via population subdivision, and driven by epistatic selection. Accordingly, LD can arise between markers on different nuclear chromosomes, and between the mitochondrial and nuclear genomes (Sloan et al., 2015). Linkage disequilibrium has proven a powerful tool for studying the genomic architecture of population divergence, local adaptation, and reproductive isolation (Nosil et al., 2009; Servedio, 2009; Smadja and Butlin, 2011). Natural selection can favor the evolution of high LD when multiple loci that influence a trait experience the same divergent selection (Nielsen, 2005). Strong LD between nuclear-encoded mitochondrial alleles could signal co-adapted genes responding to the same selective drivers, which may or may not be maintained by genomic architecture favoring reduced levels of recombination. Significant clustering of nuclear genes encoding mitochondrial or chloroplast proteins in Arabidopsis has been demonstrated (Alexeyenko et al., 2006). Similar analysis for animal taxa has rarely been performed, not least because fully assembled genomes are unavailable for many of the organisms for which mitonuclear co-evolution might be relevant, but some significant mitochondrial-nuclear LD is present in humans (Sloan et al., 2015). As genomic resources for non-model organisms expand, we expect to see more studies of mitonuclear genomic architecture (e.g., genomic re-arrangements of nuclear-encoded mitochondrial genes), as is increasingly the case for other traits linked to reproductive isolation (Lowry and Willis, 2010; Jones et al., 2012; Egan et al., 2015).

Even without a reference genome, LD can be estimated by analyzing LD clustering with LDna (Kemppainen et al., 2015). This tool finds clusters of loci with similarly high levels of LD independently of their position in the genome. This is a valuable first step to studying the genomic architecture of organisms without genomic resources. It is important to note that high LD can also emerge through processes unrelated to selection, such as population history and structure, which should be accounted for (Mangin et al., 2012; Goicoechea et al., 2015). Demographic factors, however, should impact many loci across the genome, so significant excesses of high LD among nuclear genes with mitochondrial functions is indicative of non-neutral processes (Morales et al., 2016a). Comparative genomic approaches are also recommended to investigate whether re-arranged or ancestral genomic architectures in closely related taxa are more or less prone to evolution of mitonuclear interactions, and whether nuclear-encoded mitochondrial gene re-arrangements can be favored to reduce recombination between locally co-adapted alleles and can promote adaptive genetic divergence (Faria et al., 2011; Yeaman, 2013).

One of the richest examples of the role of mitochondrial evolution in local adaptation comes from studies of the bar-headed goose Anser indicus. During its well-documented migrations over the Himalayas, this bird sustains high metabolic rates as a result of a multitude of physiological adaptations (Bishop et al., 2015). Comparative physiological studies show that adaptations at multiple levels of organization, from protein activity to organ morphology, enable this species to enhance O₂ supply and modulate O₂ usage compared to low-altitude geese species (Scott et al., 2009, 2011).

Among the potentially adaptive differences observed is a difference in the substrate kinetics of OXPHOS complex IV in bar-headed geese compared to other species. Cardiac muscle measurements show that this enzyme has a comparatively high affinity but low activity for its nuclear-encoded substrate cytochrome c. Scott et al. (2011) proposed that this may be adaptively relevant by enabling the mitochondrion to maintain redox balance in response to limitations and fluctuations in their O₂ supply during their extreme flights. The authors inferred the genotypic basis of these changes by comparing the sequences of the mtDNA and nuclear-encoded mitochondrial genes from complex IV subunits between low- and high-altitude species. This revealed several non-synonymous substitutions in bar-headed geese, including a unique W116R substitution in COX3, as well as subtle differences in the expression levels of the complex (Scott et al., 2011).

Despite these strong phenotype-to-genotype links, it remains to be determined how these substitutions affect the biochemistry of the complex. Homology modeling suggests that the W116R mutation disrupts inter-subunit interactions in complex IV, but it is unclear if and how this affects the binding of cytochrome c (Scott et al., 2011). These gaps reflect the major challenges associated with purifying this enzyme complex for kinetic or structural characterisation. A striking contrast is provided by studies on why the O₂-affinity of hemoglobin is so high in bar-headed geese; the relative ease of purifying this protein from erythrocytes has facilitated structural studies showing that a single amino acid substitution markedly shifts the cooperative behavior of the hemoglobin tetramer and in turn modulates O₂ affinity (Zhang et al., 1996; Liang et al., 2001).

These studies of the bar-headed goose offer a strong example of integrated research leading to an understanding of the mechanistic links between genes affecting mitochondrial function and ecophysiological phenotype. While there are multiple classes of gene that might be expected to contribute to the adaptive phenotype (Scott et al., 2011), there is no particular expectation that strong genome architecture is required to promote the co-evolution of these: the species does not interbreed with another, and there are no apparent differently adapted lineages within the bar-headed goose, so gene flow should not disrupt co-adapted combinations.

Through extensive studies, Burton and colleagues have demonstrated that mitonuclear discordance has contributed to allopatric population divergence of Tigriopus californicus copepods. There is extraordinary genetic differentiation between populations of this intertidal copepod across geographic barriers in the Californian coast, with mtDNA divergence exceeding 18% between Santa Cruz and San Diego populations (Burton and Lee, 1994; Burton, 1998). Elegantly designed aquarium experiments revealed that F₁ hybrids from the T. californicus populations are viable, but subsequent generations exhibit a range of fitness defects and reduced ATP production rates (Burton and Lee, 1994; Ellison and Burton, 2006, 2008). Consistent with a mitochondrial origin, maternal but not paternal backcrossing can restore fitness of progeny (Ellison and Burton, 2008).

Targeted sequencing revealed that there are high levels of divergence in genes encoding key determinants of the mitochondrial electron transport chain. While most of these substitutions appear to be neutral, ω-based approaches provided strong evidence for positive selection for substitutions in mitochondrially encoded complex IV and its nuclear-encoded substrate cytochrome c (Rawson et al., 2000; Willett and Burton, 2004). Interpopulation crossing experiments substantiated this by showing that mitonuclear compatibility was required for optimal complex IV activity (Edmands and Burton, 1999) and that the mitotype modulated segregation ratios of cytochrome c (Willett and Burton, 2001). The authors went further by validating these predictions using biochemical approaches. Intrapopulation pairs of complex IV and cytochrome c consistently showed up to fourfold higher activity than did interpopulation pairs (Rawson and Burton, 2002). Moreover, recombinantly produced cytochrome c variants (reflecting different nuclear backgrounds) interacted differentially with complex IV in tissue homogenates (reflecting different mitochondrial backgrounds) (Harrison and Burton, 2006). This work proved that single substitutions are sufficient to cause mitonuclear incompatibilities in wild populations.

Despite these accomplishments, the evolutionary processes and pressures that result in allopatric population divergence remain under investigation. Pereira et al. (2016) recently approximated the contribution of genetic drift and natural selection in T. californicus divergence by comparing whole-transcriptome sequences of allopatric populations at different stages of divergence (Pereira et al., 2016). They found that the pattern of shared polymorphism could be partially explained by genetic drift, as lower effective population sizes led to less shared polymorphism between populations, and higher mutation load. However, natural selection possibly drives accelerated evolution of some genes, including nuclear-encoded mitochondrial ones. The authors predict that genomic architecture should regulate the efficiency of selection and the impact of drift, for example by modulating recombination rates, but this prediction was not tested.

Our studies on the population structure of eastern yellow robin Eopsaltria australis have emphasized the importance of studying genomic architecture (Morales et al., 2016a). This songbird is one of multiple animals that maintains functional mitonuclear interactions despite discordance between its mitochondrial and nuclear genomes (Toews and Brelsford, 2012). Population genetic data have shown that the major axis of nuclear DNA differentiation runs north-south through the species range in Eastern Australia, whereas mitochondrial DNA has diverged into two mitolineages in the perpendicular coastal-inland direction (Pavlova et al., 2013; Morales et al., 2016a). Coalescent analyses suggest that the two genomes initially differentiated together in a north–south direction during the early Pleistocene, but their evolutionary history became separated due to two mitochondrial introgression events in the mid-to-late Pleistocene (Pavlova et al., 2013; Morales et al., 2016b). The two mitolineages show sharp climate-correlated differences in their distributions. This suggests that the mitochondrial introgression and resulting divergence were driven by natural selection (Morales et al., 2016b).

There is evidence of positive selection for non-conservative amino acid differences in the proton pumps ND4 and ND4L of OXPHOS complex I between the mitolineages of the eastern yellow robin. These polymorphisms are predicted to cause differences in electrostatic subunit-subunit interactions and in turn influence coupling efficiencies of the complex, though this remains to be validated experimentally (Morales et al., 2015). Comparison of fixation indexes in the nuclear genomes between eastern yellow robin populations across their biogeographic range revealed the existence of two genomic islands of divergence against a background of low differentiation (Morales et al., 2016a). One of these islands, located on autosome 1A (Figure 1), is statistically overrepresented with nuclear-encoded genes with predicted mitochondrial functions; among them are three complex I supernumerary subunits proposed to be functionally linked to the mitochondrially encoded ND4 and ND4L genes. Moreover, markers within the genomic island of divergence exhibit very strongly elevated LD, suggesting genome architecture that promotes reduced recombination between nuclear-encoded genes with mitochondrial functions. Further research will disentangle whether mitonuclear co-evolution promoted the evolution of this particular genomic architecture or pre-existing genomic architecture enabled this tight mitonuclear co-evolution.