Reproductive isolation is often achieved when genes that are neutral or beneficial in their genomic background become functionally incompatible in a foreign genome, causing inviability, sterility or low fitness in hybrids. Recent studies suggest that mitonuclear interactions are among the initial incompatibilities to evolve at early stages of population divergence across taxa. Yet, the genomic architecture of mitonuclear incompatibilities has rarely been elucidated. We employ an experimental evolution approach starting with low-fitness F2 interpopulation hybrids of the copepod Tigriopus californicus, in which frequencies of compatible and incompatible nuclear alleles change in response to a alternate mitochondrial background. After about nine generations, we observe a generalized increase in population size and in survivorship, suggesting efficiency of selection against maladaptive phenotypes. Whole genome sequencing of evolved populations showed some consistent allele frequency changes across three replicates of each reciprocal cross, but markedly different patterns between mitochondrial backgrounds. In only a few regions (~6.5% of the genome), the same parental allele was overrepresented irrespective of the mitochondrial background. About 33% of the genome showed allele frequency changes consistent with divergent selection, with the location of these genomic regions strongly differing between mitochondrial backgrounds. The dominant allele matched the mitochondrial background in 87 and 89% of these genomic regions, consistent with mitonuclear coadaptation. These results suggest that mitonuclear incompatibilities have a complex polygenic architecture that differs between populations, potentially generating genome wide barriers to gene flow between closely related taxa.

Sampling

To understand the genetic architecture of mitonuclear incompatibilities we focused on two well-studied populations of Tigriopus californicus: San Diego (SD) and Santa Cruz (SC). These populations were suitable for a competition experiment between divergent nuclear alleles under a fixed mitochondrial background because previous studies have shown that: i) marked mitochondrial divergence (22.17%; Pereira et al. 2016) has driven adaptive evolution targeting genes functionally interacting with the mitochondria (Barreto & Burton 2013a), ii) 91% of the nuclear polymorphisms are fixed among parental populations (Pereira et al. 2016) allowing an accurate estimation of allelic frequencies in hybrids (Lima et al. 2019), iii) F2 hybrids show breakdown at all physiological and life history traits known to be associated with mitonuclear incompatibilities (Ellison & Burton 2006; a; b; Barreto & Burton 2013b; Barreto et al. 2015), and iv) strong selection during F2 development results in significant deviations from the expected Mendelian inheritance (Foley et al. 2013; Healy & Burton 2020).

Stock populations were collected in San Diego (32◦ 44′ 41.17′′ N, 117◦ 15′ 19.43′′ W) and Santa Cruz (36◦ 56′ 58.32′′ N, 122◦ 2′ 48.98′′ W). Multiple cultures from each site were maintained in 400 mL beakers at common garden at 20°C with a 12:12 L:D photoperiod, in filtered seawater mixed with ground “Algae Wafers” (Kyorin Co., Himeji, Japan). Beakers were monthly mixed to promote panmixia within populations and medium was renewed.

Experimental evolution

T. californicus is sexually dimorphic. An adult male clasps an immature female until her terminal molt, when she is then inseminated and released. Females can mate only once and use stored sperm to fertilize sequential clutches of eggs that can add to several hundred progeny (Vittor 1971). This species lacks heteromorphic sex chromosomes and recombination occurs only in males (Burton et al. 1981). We obtained virgin females by separating clasped pairs, and produced F1 with similar mixed nuclear genome and fixed mitochondrial backgrounds; i.e. SD mitochondrial background on SD♀ x SC♂ cross, and SC mitochondrial background on SC♀ x SD♂. For each reciprocal cross, we outcrossed F1s to produce recombinant F2 hybrids that were allowed to mate randomly.

Each experimental line started with 100 outbred F2 gravid females. Lines evolved under these conditions for nine months with overlapping generations, replenishing the growing medium monthly. Since on average females reach adulthood in 2-3 weeks and produce multiple egg clutches until they are 4-6 weeks old, this experimental design corresponds to approximately nine generations of experimental evolution up to F11. This procedure was replicated 10 and 7 times for the SC and SD mitochondrial backgrounds, respectively. We followed the same procedure to generate one control line with fully matched nuclear and mitochondrial genomes for each parental population (i.e. SC♀ x SC♂ and SD♀ x SD♂).

Fitness recovery

Relatively small experimental populations may lead to strong genetic drift, and conversely to limited response to selection imposed by the fixed mitochondrial backgrounds. If selection is strong relative to drift in evolved lines we expect an increase in productivity and associated recovery in one or multiple fitness traits associated with mitonuclear incompatibilities (Ellison & Burton 2008b). To test these hypotheses, at the end of the experimental evolution, we measured: census size (as the number of adults after), fecundity (as the number of nauplii larvae hatching from the first clutch of a female), and survivorship (as the fraction of nauplii surviving to 14 days). Fecundity measurements were replicated between 4 to 12 times, depending on the number of available virgin females, and survivorship was replicated between 10 to 28 times, depending on the number of available gravid females. To monitor how average fitness varied along the course of the experiment, we have also measured survivorship 3 to 8 additional times, using 4 to 12 replicates. Additionally, we measured these two fitness traits for the initial reciprocal F2 hybrids and for the pure parental populations, as a reference for fitness breakdown and recovery respectively. We estimated mean plus ±1SE. We tested for significant hybrid breakdown by comparing fitness of the F2 hybrid with its maternal population, using a Mann–Whitney U-test and an alpha of 0.05 in R 2.15.1 (R Development Core Team). We tested for significant recovery in lines for which the mean reached or passed the reference parental fitness, adjusting the P-value when multiple comparisons occur at the same time.

Allelic frequency change

To provide insights into the genetic architecture of mitonuclear incompatibilities we examined allelic frequency change across the whole genome in lines evolving under each mitochondrial background, relative to the initial F2 females. We used the pool-seq approach, which has recently been established for evolve-and-resequencing studies (Schlotterer et al. 2015), including in T. californicus (Lima & Willett 2018; Lima et al. 2019; Healy et al. 2020; Griffiths et al. 2020).

For lines evolving under each mitochondrial background, we selected three replicates that showed larger census size, no further decrease in fecundity, and some recovery in survivorship. We pooled 200 adult individuals, or the maximum available, and extracted genomic DNA using the phenol chloroform protocol and digestion with RNAse (Sambrook & Russell 2001). For the initial F2 females, we also extracted DNA from pools of 100 adult outcross F2 females that gave rise to the experimental populations, after removing their egg sacs. All samples were sequenced in the Illumina HiSeq 2500 platform with 100 bp paired-end libraries. Reads were trimmed for quality using PoPoolation2 (Kofler et al. 2011), discarding bases with Phred quality scores lower than 25 and keeping reads of at least 50 bp after trimming.

We only considered SNPs that were fixed between parental populations and therefore that can be used to determine the ancestry of the nuclear alleles favored in either mitochondrial background. We used the bioinformatics pipeline established by (Lima & Willett 2018; Lima et al. 2019) (see references for details), and the syntenic reference genomes of SD and SC, where more than 90% of the genome is anchored to the 12 chromosomes (Barreto et al. 2018). In short, first, we made the two reference genomes equivalent in length and accuracy by adding N’s to any position where either parental reference had an “N”. Second, we established a list of fixed SNPs by i) performing reciprocal mapping of reads from one population (Barreto et al. 2018) against the reference of the other, ii) considering only “fixed” positions where all mapped reads showed an alternative nucleotide to the reference (coverage ≥ 15), and iii) comparing reciprocal mapping to keep only SNPs that were “fixed” in both mappings. Third, we mapped our reads for each hybrid dataset to both parental references using BWA-MEM with default parameters (Li & Durbin 2009) and keeping reads that mapped with MAPQ score > 20. When mapping reads from hybrids to one of the parental references, there is a known bias towards an over-representation of reads of the reference population (Lima & Willett 2018). By averaging the read counts for each SNP, from the alignments to each parental population, this reference-population bias is overcome. Allele counts for the fixed SD and SC alleles and respective frequencies were determined using PoPoolation2 (Kofler et al. 2011) (minimum coverage of minor allele ≥ 4). Finally, data for all cross datasets were compared and only SNP that passed all filters up to this point, for both reciprocal crosses, were kept for statistical analysis of allele frequency change. These SNPs are expected to show high levels of variance, due to the variation on sequencing coverage and sampling alleles from the pool of individuals (Lima & Willett 2018). We averaged allele counts (and frequencies) in non-overlapping windows of 500 consecutive SNPs. We refer to these as “genomic windows”.

In each mitochondrial background separately, we determined which of these genomic windows show a significant allelic frequency change by comparing the initial F2 hybrids to the evolved hybrids, using the Cochran-Mantel-Haenszel (CMH) test as implemented in the mantelhaen.test function in the statistical R package (R Development Core Team). This test operates on 2x2 contingency tables (times the number of replicates), comparing the counts for SD and SC alleles at the beginning and at end of the experiment, and thus it is highly sensitive to variation in sequencing coverage. To account for differences in coverage between the initial and evolved hybrids, we estimated initial counts of SD and SC alleles by multiplying the allelic frequency estimated in the initial F2 hybrid pool by the total coverage observed at each evolved hybrid separately (i.e. counts for SD + SC). This approach results in a P-value per genomic window, which reflects the probability of rejecting the null hypothesis of similar allelic frequency during experimental evolution in a given mitochondrial background.

We used a z-score of 2 to identify genomic windows showing higher P-values relative to the average values observed genome wide, consistent with stronger selection driving allelic frequency change. We have estimated three z-score threshold based on three distributions of P-values: i) from both crosses taken together (corresponding to a -log10 (P-value) of 3.24), which allows identifying differences in the strength of selection between the two mitochondrial backgrounds; ii) only from crosses evolving under the SC mitochondria (corresponding to a -log10 (P-value) of 4.17), which allows identifying targets of selection within the SC mitochondrial background;  and iii) only from crosses evolving under the SD mitochondria (corresponding to a -log10 (P-value) of 1.73), which allows identifying targets of selection within the SD mitochondrial background. Since adjacent genomic windows are physically linked, we looked for consecutive genome windows that are under selection. We refer to these as genomic regions. Genomic regions were classified as under selection through the following criteria: i) at least 10 genomic windows with -log10 (P-value) ³ 1; ii) with at least one window above a z-score of 2; and iii) ending with one window -log10 (P-value) < 1. To represent which parental allele dominates during the experimental evolution, we classified genomic regions as blue for overrepresentation of the SC nuclear allele and red for overrepresentation of SD allele. By comparing allelic frequency change in reciprocal crosses, we identified genomic regions fitting three categories: “uniform selection”, where the same allele dominates irrespective of the mitochondrial background; “divergent selection” where the dominant nuclear allele matches the mitochondria; and “antagonistic selection” where the dominant allele does not match the mitochondria. While we chose to use z-scores as threshold for selection, the results were qualitatively similar when considering a threshold based off the top 5% of the distribution of P-values.