Introduction:
Speciation occurs when populations accumulate genetic differences that cause reproductive incompatibilities in hybrids, leading to maladaptive phenotypes such as inviability, sterility, or to reduced fitness relative to parentals. Bateson (Bateson 1909), Dobzhansky (Dobzhansky 1937), and Muller (Muller 1942) have proposed that such incompatibilities would involve two or more loci, so that substitutions that are adaptive or neutral in their own genomic background can be functionally incompatible with alleles that are present in a foreign genomic background; the so-called Bateson-Dobzhansky-Muller incompatibilities or BDMIs. Although this model has guided research on the genomic basis of speciation for more than 70 years, only within the last decade have genomic regions (and in few cases genes) causing full sterility or inviability been identified in taxa as diverse as yeast, flies, mice and plants (Presgraves 2010). Such genes are found to evolve rapidly, as expected for genes caught up in open-ended molecular evolutionary arms race, and are often involved in co-evolution between host and pathogen, or between selfish genes and suppressors (Presgraves 2010). Particularly in studies with highly divergent species pairs, BDMIs seem to have a complex genomic architecture that involves tens of genes spread throughout the genome (Tang & Presgraves 2009; Schumeret al. 2018). Yet, full sterility or inviability evolve later during the speciation process, and thus it is unclear whether those genes reflect initial barriers that arise between diverging populations or whether they reflect other evolutionary forces that arose after reproductive isolation was established. Therefore, understanding the genetic architecture of incompatibilities causing partial reproductive isolation, such as reduced fertility, remains an important task in evolutionary biology (Corbett-Detig et al. 2013; Turner & Harr 2014; Rafati et al. 2018).
Recent studies in taxa as diverse as arthropods, nematodes, vertebrates, yeast, and angiosperms (reviewed in Burton et al. 2013; Sloanet al. 2017; Hill et al. 2018) have shown that the co-evolution between the mitochondrial and the nuclear genome often results in partially reduced fitness in hybrids between closely related taxa (i.e. species, subspecies, or even populations). This observation has led several authors to suggest that mitonuclear BDMIs play a disproportional role at early stages of reproductive isolation relative to the more widely studied nuclear-nuclear BDMIs (Burton & Barreto 2012; Hill 2016). To understand why such a general pattern would arise across species, one needs to consider the origin and evolution of the mitochondrion in eukaryotes. Some 1.5 billion years ago, when a proteobacterial heterotroph (proto-mitochondrion; Martijn et al.2018) became the obligatory endosymbiont of an archaebacterial methanogen (proto-nucleus), this symbiosis resulted in a strong functional partition (Rivera et al. 1998), where the proto-mitochondrion became responsible for metabolic functions and the proto-nucleus for transcription and translation. Because intracellular clonal replication favors smaller genomes (Taylor et al. 2002; Clark et al. 2012), a race for replication of the proto-mitochondrion lead to the reduction in the endosymbiont genome size, both through the loss of redundant genes and through gene transfer to the host genome (Selosse et al. 2001; Timmis et al.2004). Such selective pressure is still observed today in the genome evolution of intracellular endosymbionts in insects (Wernegreen 2002), in the recent transposition of mitochondrial genes into the nuclear genome (Hazkani-Covo et al. 2010), or even in the deletion of essential mitochondrial genes with significant fitness costs for the organism (Taylor et al. 2002). Currently, in most metazoans, over 1,000 proteins are housed in the mitochondrion (Calvo & Mootha 2010), but only 13 of those remain coded in the mitochondrial genome. Hence, mitochondrial function is dependent upon nuclear-encoded proteins, many of which interact closely with mtDNA encoded proteins, RNAs or DNA-binding sites. This mitonuclear cooperation sets the stage for a tight evolutionary arms race between the nuclear and mitochondrial genomes to preserve functionality within eukaryotic evolutionary lineages in essential cell functions, such as respiration and mitochondrial protein synthesis, that have large impacts on organismal fitness. Such mitonuclear coadaptation becomes exposed in F2 hybrids, where independently evolving mitochondrial and nuclear alleles are forced to interact, often being functionally incompatible, disrupting organelle function and, consequently, causing hybrid breakdown (Burtonet al. 2013).
Contrary to nuclear-nuclear BDMIs, mitonuclear incompatibilities are expected to be highly asymmetric for a variety reasons, as described in the model of compensatory coevolution (Rand et al. 2004). First, because of the mode of replication of its circular genome, mutation rate (µ ) is higher in the mitochondrial relative to the nuclear genome. Second, due to its lack of recombination and matrilineal inheritance, mitochondrial genes have an effective population size (Ne ) that is 4 times smaller than that of nuclear genes, resulting in higher rates of fixation via genetic drift and conversely a reduced efficiency of selection. These two processes result in the consistent observation of faster evolution rates of mitochondrial relative to nuclear genes, which are 2-fold in drosophilids, 20-fold in ungulates and up to 40-fold in primates (Osada & Akashi 2012). Even though most de novo mutations are partially deleterious, the relatively higher levels of genetic drift allow the fixation of these new mutations. This asymmetry of µ and Ne leads to a higher fixation rate for weakly deleterious mutations in the mitochondrial genomes relative to the nuclear genome; a process known as “Muller’s ratchet” (Lynch & Blanchard 1998). This accumulation of deleterious mutations in the mitochondrial genome elicits compensatory mutations in the interacting nuclear genes (Osada & Akashi 2012; Barreto & Burton 2013a; Sloan et al. 2014) but not in the rest of the nuclear genome, maintaining the stability and function of mitonuclear protein complexes. Although such mitochondrial deleterious mutations are effectively silenced within a population, they become exposed in interpopulation crosses where coadapted mitonuclear complexes become mismatched in hybrids, contributing for the establishment of genetic barriers between recently diverged taxa (Sloan et al. 2017). The magnitude of such early barriers to gene flow directly depends on the genetic architecture of mitonuclear BDMIs, which remains unknown across species.
Although mitonuclear incompatibilities are now recognized to play a general role in establishing reproductive isolation between emerging species (Reinhardt et al. 2013), their effect is most visible in taxa presenting exceptionally high mitochondrial evolution rates, such as the copepod Tigriopus californicus (Willett 2012). Allopatric divergence between populations of this species resulted in parallel patterns of genomic divergence that are consistent with mitonuclear coevolution happening within independent populations. Protein coding genes from the mitochondria evolve 2 to 14 times faster than those from the nuclear genome (Pereira et al. 2016), suggesting that mitochondrial genes drive intragenomic coevolution. Nuclear encoded proteins that functionally interact with the mitochondria evolve more rapidly than non-interacting nuclear encoded proteins (Barreto et al. 2018), suggesting that selection favoring compensatory mutations targets specific nuclear genes. Finally, experimental interpopulation F2 hybrids show that fitness breakdown in multiple life history traits (such as fecundity, survivorship and developmental time) scales with mitochondrial divergence (Burton 1990a; Edmands 1999). Notably, F2 fitness breakdown is rescued in maternal backcrosses (Ellison & Burton 2008b), where mitochondrial and nuclear coevolving units are rematched, demonstrating that fitness breakdown is caused by interactions between the mitochondria and unknown nuclear loci. Despite such an evident parallelism in mitonuclear coevolution within populations, evolutionary processes conditioned by local habitat are non-parallel. For example, ecological trade-offs along latitudinal gradients have resulted in differential adaptation to temperature (Willett 2010; Hong & Shurin 2015; Pereira et al. 2017) and to salinity (Leong et al.2017). Moreover, smaller Ne at southern populations has resulted in stronger genetic drift relative to northern populations (Pereiraet al. 2016), and potentially in a faster accumulation of deleterious mutations genome wide.
Recent studies with interpopulation hybrids of T. californicushave established that mitonuclear incompatibilities result in fitness breakdown at various organizational levels. While heterozygous F1s are vigorous, mismatched mitonuclear complexes in F2 and inbred lines results in: reduced mitochondrial function (Ellison & Burton 2008a), reduced ATP production (Ellison & Burton 2006), elevated oxidative damage to DNA (Barreto & Burton 2013b), upregulation of pathways involved in physiologic stress (Barreto et al. 2015), and breakdown at multiple life history traits (fecundity, survivorship and developmental time; Burton 1990a; Edmands 1999). It is still unclear whether such generalized fitness breakdown cascade from mitonuclear BDMIs involving few nuclear genes with generalize effect, or involving multiple regions spread throughout the nuclear genome. It is worth noting that this species lacks sex chromosomes, and instead has polygenic sex determination (Voordouw & Anholt 2002; Alexander et al. 2015). Therefore, factors leading to the disproportionate role of sex chromosomes in the evolution of BDMIs reported in many species (Presgraves 2018), does not apply in this case.
Previous efforts to understand the genetic architecture of mitonuclear BDMIs in T. californicus have focused on deviations from the expected Mendelian inheritance in F2 hybrids surviving to adulthood (Foley et al. 2013; Lima & Willett 2018; Lima et al.2019), in hybrid swarms with recovered fitness (Pritchard & Edmands 2012) or among high- and low-fitness hybrids within an F2 cohort (Healy & Burton 2020). Although those generally confirm a prominent role of mitonuclear incompatibilities relative to nuclear-nuclear BDMIs (but see Pritchard et al. 2011), they strongly differ in the estimated number, location and effect size of nuclear loci interacting with the mitochondria (ranging from few loci to all 12 chromosomes). Understanding the genomic architecture of mitonuclear BDMIs requires a whole-genome approach, over multiple generations of hybridization and selection, and under divergent selection associated to alternative mitochondrial backgrounds. To address this, we employed an experimental evolution approach initiated with replicate populations of low-fitness F2 interpopulation hybrids. Over the course of the experiment, parental nuclear alleles competed in alternate mitochondrial backgrounds (i.e. treatment) established by reciprocal crosses, for about nine generations (Fig. 1). To assess the efficiency of natural selection in small experimental populations we tested for increase in population size and for fitness recovery in female fecundity and nauplii survivorship, which are known to be associated with mitonuclear incompatibilities. Next, to find regions likely responding to selection (uniform or divergent between treatments), we identified genomic regions with significant allelic frequency change during the experiment. Finally, to tease apart nuclear genomic regions involved in mitonuclear incompatibilities we identified regions that differ among alternative mitochondrial backgrounds. We found genomic regions likely involved in mitonuclear BDMIs on multiple chromosomes and that these regions differ between mitochondrial backgrounds, suggesting that mitonuclear incompatibilities have a complex and asymmetric genetic architecture.