Introduction

Elucidating the different types of pre- or post-zygotic reproductive isolating barriers is crucial for understanding speciation (Ramsey et al., 2003; Coyne and Orr 2004; Nosil et al., 2005; Rieseberg and Willis, 2007; Lowry et al., 2008). Prezygotic mechanisms impede mating through, for example, habitat, temporal or sexual isolation, or mechanical isolation because of genitalic incompatibility (Nosil et al., 2005; Lessios, 2011). On the other hand, postzygotic selection acts after the zygote has been formed because of reduced fitness of hybrids that may eventually lead to reinforcement and speciation. Here a distinction is made between extrinsic and intrinsic postzygotic barriers (Mayr, 1963; Nosil et al., 2005; Seehausen et al., 2014). In extrinsic postzygotic isolation, hybrid phenotypes are maladapted to the environment relative to parental phenotypes (Nosil et al., 2005; Schluter, 2009). Intrinsic postzygotic barriers involve genetic incompatibility when suboptimal allelic combinations are brought together in hybrids. Deleterious combinations will lead to reduced fitness or infertility, as demonstrated both theoretically (Gavrilets, 2003, 2004; Coyne and Orr, 2004) and empirically (Brideau et al., 2006; Presgraves, 2007; Maheshwari and Barbash, 2011). A special case, the so-called cytonuclear incompatibility, involves mitochondrial and nuclear genes (Burton et al., 2013). Genes from both genomes interact within protein complexes of the oxidative phosphorylation pathway that is responsible for energy (adenosine triphosphate (ATP)) production (Saraste 1999; Ballard and Whitlock, 2004). Thus, coadaptation within populations or species may result in mismatches and protein dysfunction in hybrids leading to reduced fitness (Burton et al., 2013). Evidence for such mismatches has been shown in several organisms, including flies (Meiklejohn et al., 2013), copepods (Ellison and Burton, 2006) and monkey flowers (Fishman and Willis 2006). Nonetheless, the absence of coevolution is not always deleterious. For example, a recent study of the redbelly dace (Chrosomus eos) showed that introgression of a foreign mitogenome, of up to 10 Myr of independent evolution, is not necessarily deleterious but may rather have beneficial effects (Deremiens et al., 2015). In total, few good cases of cytonuclear incompatibility exist outside the laboratory and more empirical evidence is needed to understand the significance of this mechanism in speciation (Burton et al., 2013).

Here we investigate pre- and post-zygotic barriers between the two sister species of North Atlantic eels: European (Anguilla anguilla) and American eel (A. rostrata) (Tesch, 2003). Both species are panmictic (Als et al., 2011; Côté et al., 2013), show high effective population sizes (Côté et al., 2013; Pujolar et al., 2013; Jacobsen et al., 2014b) and spawn in the thermal fronts of the southern Sargasso Sea (McCleave et al., 1987; Tesch, 2003). Their spawning shows extensive overlap in space (Schmidt, 1923; Tesch, 2003; Munk et al., 2010) and time, with peak spawning time in February–April for American and March–June for European eel (McCleave et al., 1987; Tesch, 2003; Miller et al., 2015). After spawning, the larvae (leptocephali) are transported by ocean currents back to their respective feeding grounds (Schmidt, 1923), where they enter fresh or brackish water as glass eels. After a period of 4–20 years as yellow eels, they metamorphose into silver eels that migrate back to the Sargasso Sea to spawn and subsequently die (Tesch, 2003). The spawning and larval migrations are highly differentiated in Atlantic eels with a spawning migration of ∼5000 km for European eel and 2500 km for American eel (Aoyama, 2009) and a subsequent duration of larval phases of ∼2 years and 7–9 months (Tesch, 2003), respectively, although this is still debated (Tesch, 2003; Zenimoto et al., 2011).

Mitogenome sequencing suggests that species divergence of Atlantic eels was initiated ∼3.3 million years ago (Myr) (Jacobsen et al., 2014b). However, gene flow still occurs, as hybrids have been reported in several studies (Avise et al., 1990; Albert et al., 2006; Gagnaire et al., 2009; Pujolar et al., 2014a, 2014b). Thus, although temporal and spatial separation at the spawning grounds likely constitutes prezygotic barriers, they are not absolute. Nonetheless, the frequency of hybridization is unknown and, to date, no F1 hybrids have been found in the Sargasso Sea, where merely two second-generation backcrosses previously have been detected based on 86 species-diagnostic (FST ⩾0.95) single-nucleotide polymorphisms (SNPs) in samples from 2007 (Pujolar et al., 2014a). The finding of only backcross individuals but no F1 hybrids may be explained by temporally varying hybridization (Pujolar et al., 2014a). This could result from differences in temporal overlap of spawning time between the species across years, potentially mediated by annual differences in the location of the thermal fronts that are highly dynamic across years (Tesch, 2003; Ullman et al., 2007).

F1 hybrids have only been observed in Iceland (Avise et al., 1990; Albert et al., 2006; Pujolar et al., 2014a). This is possibly a consequence of an intermediate larval migration behavior, where F1 hybrids either actively or passively get transported to Iceland, approximately located longitudinally intermediate between North America and Europe (Avise et al., 1990; Pujolar et al., 2014a). In continental Europe and North America, only few later-generation backcrosses have been observed (Albert et al., 2006; Pujolar et al., 2014b), suggesting strong postzygotic barriers. One such barrier has been proposed to be cytonuclear incompatibility (Gagnaire et al., 2012), based on the finding of significant inter-species nonsynonymous differentiation in genes involved in the ATP synthase complex (Gagnaire et al., 2012). This protein complex is part of the oxidative phosphorylation chain and catalyzes the synthesis of ATP from adenosine diphosphate (Saraste, 1999). The two genes showing highest nonsynonymous differentiation in Atlantic eels are the mitochondrial ATP6 (ATP synthase F0 subunit 6) and the nuclear ATP5c1 (ATP synthase F1 subunit gamma) genes. Thus, coevolution likely has occurred between these genes and mismatches in hybrids are expected to be associated with postzygotic selection. Such selection may be linked to the possibility of completing the migratory loop; either the spawning migration or the subsequent larval phase. Both life-history traits may be associated with energetics and are significantly correlated with amino acid changes in ATP6 in freshwater eels (Anguillidae) (Jacobsen et al., 2015). Although European eel can be reared in vitro through the preleptocephali stage (Sørensen et al., 2016), the same is not yet possible for the American eel and none of them can presently be reared artificially. Thus, hybrids can only be studied from wild-caught specimens (Tesch, 2003). Here recent studies based on restriction site associated DNA (RAD) sequencing (Pujolar et al., 2014b) and analysis of a subset of species-diagnostic SNPs (Pujolar et al., 2014a) have detected hybrids from beyond the F1 generation backcross level. Screening for unusual ATP6–ATP5c1 combinations in such hybrids can be used for investigating the possibility of cytonuclear incompatibility in Atlantic eels (Blier et al., 2001, Gagnaire et al., 2012).

The use of species-diagnostic SNPs to detect hybrids may, however, involve potential problems. Genome-wide differentiation between American and European eel is highly heterogeneous with an average FST of 0.041. Nevertheless, thousands of SNPs show fixed differences and are candidates for being under divergent selection between the species (Jacobsen et al., 2014a). Hence, it is likely that selection acts on these loci in hybrids. As such, whereas F1 hybrids should be heterozygous for all diagnostic loci, selection and coadaptation may distort genotype frequencies in post-F1 hybrids, leading to increased inaccuracy in estimating numbers of generations since initial hybridization.

Here we addressed the following questions. (1) Is the degree of hybridization temporally varying in the Sargasso Sea? (2) Does selection affect diagnostic SNPs and thereby categorization of post-F1 hybrids? (3) Are ATP6–ATP5c1 combinations in post-F1 hybrids in accordance with the hypothesis of cytonuclear incompatibility partly underlying postzygotic barriers between the species? Using a panel of 96 near-diagnostic SNPs (Pujolar et al., 2014a) we analyzed eel larvae sampled in the Sargasso Sea in 2014 (N=460), in addition to new samples from Morocco (N=21) and the Faroe Islands (N=30). Data were used for assessing temporal variation in hybridization and selection in combination with previously published data, in particular from larvae collected in the Sargasso Sea in 2007 (N=92) and Icelandic eels collected in the early 2000s (N=159). We expected that F1 hybrids should be overrepresented compared with backcrosses at the spawning grounds because of postzygotic barriers, in case of continuous gene flow between the species. Furthermore, we sequenced the ATP6 and ATP5c1 genes from 19 hybrids in order to test for species-specific ATP6–ATP5c1 matches and thus indirectly for cytonuclear incompatibility. These hybrids were identified in the present and previous studies (Jacobsen et al., 2014a; Pujolar et al., 2014a, 2014b) by analyzing species-diagnostic SNPs using STRUCTURE (Pritchard et al., 2000; Falush et al., 2003, 2007) and NEWHYBRIDS (Anderson and Thompson, 2002). Overall, our study provides new insights into the pre- and post-zygotic barriers involved in speciation in the case of North Atlantic eels in particular, and speciation in the oceans in general.

Materials and methods

Identifying species and hybrids among eel larvae from the Sargasso Sea

A total of 472 candidate Anguilla sp. larvae (leptocephali), sampled in the Sargasso Sea from 18 March to 15 April 2014, were analyzed in this study. They represent 47 different localities along 7 transects at longitude: 68.5°W, 65.5°W, 62.7°W, 59.5°W, 57.0°W, 53.5°W and 50°W (Figures 1 and 2a) and cover the main part of the spawning ground shared by the two species that historically has been ranging from ca. 70 to 58°W (Miller et al., 2015). A single sampling station at longitude 51.8°W was also included. Sampling was conducted using a ring net (diameter: 3.5 m, length: 25 m and mesh size of 560 μm). Larvae were stored in separate eppendorf tubes in either ethanol (96%) or RNAlater. In addition, new samples also included glass eels collected in Morocco (Oved Sebou; N=21) and yellow eels collected in the Faroe Islands (Streymnes; N=15 and Miðvágur; N=15) in 2011 (Figure 1 and Table 1).

Figure 1
figure 1

Map showing all localities of new and previously published data analyzed in this study. Black circles denote localities with new data; gray circles denote localities with both new and already published data and white squares show localities with only data already published. See Table 1 for information on exact sampling location and sample size.

Full size image
Figure 2
figure 2

Maps showing (a) sampling localities of the leptocephali larvae caught during the 2014 expedition and genotyped in this study, and (b) the frequencies of the two species at each of the localities sampled during the 2014 expedition. The two red circles denote the localities from which hybrids (bAR × AA) were detected in 2014 and red asterisks show the positions of two hybrids found among larvae sampled in 2007. The numbers in brackets denote sample sizes of analyzed larvae for each locality.

Full size image
Table 1 Overview over sampling localities, sample sizes and sequencing approaches of new and reanalyzed data
Full size table

DNA was extracted using E.Z.N.A. purifications columns (OMEGA Bio-Tek, Norcross, GA, USA). Species identification was based on PCR of the mitochondrial cytochrome b gene (CytB) according to Trautner (2006). A total of 96 species-diagnostic SNPs (FST ⩾0.95) developed by Pujolar et al. (2014a) based on RAD sequencing were then genotyped on 96.96 Dynamic Arrays (Fluidigm Corporation, San Francisco, CA, USA) using the Fluidigm EP1 instrumentation according to the manufacturer’s recommendations. In addition, genotype data from Pujolar et al. (2014a), who analyzed 280 European eels for the same 96 SNPs, were reanalyzed in the study. These data consisted of eel larvae collected in the Sargasso Sea (N=92) in 2007 (caught between 70 and 64°W), yellow and glass eels collected in Iceland (N=159) in 2000–2003 and yellow eels collected in the Faroe Islands (N=29) in 2011 (for details on exact sampling localities see Pujolar et al., 2014a). In total, the data set encompassed 803 individuals (Figure 1and Table 1).

All individuals were analyzed using NEWHYBRIDS (Anderson and Thompson, 2002). No prior information on the origin of individuals was used and an assignment threshold of >0.95 was applied. The software assigns individuals to different hybrid classes by calculating the probability of the individual belonging to a specific hybrid category. In this study we used the same hybrid categories as in Pujolar et al. (2014a): pure European eel (AA), pure American eel (AR), first-generation hybrids AA × AR (F1), F1 × F1 (F2), AA × F1 (bAA), AR × F1 (bAR), bAA × AA, bAA × AR, bAR × AR, bAR × AA, bAA × F1 and bAR × F1. This means that, for example, bAR × AA is an individual resulting from a pure European eel (AA) mating with a hybrid backcross (bAR) that again is the result of a pure American eel (AR) mating with an F1 hybrid.

Seven hybrids were assigned to the ‘unusual’ bAR × AA class and showed seven loci that were homozygous in all individuals for the European eel allele. As this is unlikely to happen by chance, and may indicate selection, we examined these loci further. First, we calculated the probability of observing the homozygote genotype in all seven individuals. Assuming Hardy–Weinberg proportions we first estimated genotype frequencies in American and European eel based on the allele frequencies observed in each of the seven biallelic loci. Subsequently, we used these frequencies to estimate genotype frequencies in F1 hybrids, backcrosses to American eel (bAR) and bAR × AA hybrids and finally the individual probabilities (see Supplementary Notes S1 and S2). Second, the positions of all SNPs were matched to the predicted complementary DNA annotation file of the European eel draft genome (Henkel et al., 2012) (http://www.zfgenomics.com/sub/eel). This allowed us to investigate: (1) whether SNPs were located within genes (introns+exons) or in noncoding regions; and (2) the possible function of the genes in which those SNPs were located. We were particularly interested in testing whether these genes matched the Gene Ontology (GO) terms ‘development’ and ‘phosphorylation’ as the recent study of Jacobsen et al. (2014a) showed these two GO term groups to be overrepresented in FST outlier SNPs between European and American eel.

Cytonuclear incompatibility in hybrids

A total of 19 hybrid individuals were analyzed for the mitochondrial ATP6 and interacting nuclear ATP5c1 genes (Table 2). These individuals corresponded to all previously analyzed hybrid individuals from which DNA or tissue was still available. The individuals belonged to samples genotyped based on either the Fluidigm 96.96 Dynamic Arrays (N=12) (Pujolar et al., 2014a; this study) or RAD sequencing (N=7) (Pujolar et al., 2014b). The seven RAD sequenced individuals from Pujolar et al. (2014b) were previously analyzed using STRUCTURE and the ‘gensback’ (G) option. This option allows testing whether an individual has an immigrant ancestor in the last G generations (Pritchard et al., 2000; Falush et al., 2003, 2007). However, the study by Pujolar et al. (2014b) only used a G of five (corresponding to the fourth-generation backcross category), and thus some individuals potentially having an older immigrant ancestor may have been wrongly assigned. As erroneous hybrid assignment may bias the assessment of possible cytonuclear selection, we reanalyzed the data using a G option of 9 (corresponding to the eight-generation backcross category) and furthermore evaluated assignment power in STRUCTURE. This is described in detail in Supplementary Note S3.

Table 2 Results of ATP6 and ATP5c1 sequencing from the 19 hybrid samples
Full size table

Both the mitochondrial ATP6 gene that contains five species-diagnostic nonsynonymous changes and the part of the nuclear interactor ATP5c1 that covers two quasi-diagnostic nonsynonymous SNPs positioned in the seventh exon were sequenced. Primers developed by Gagnaire et al. (2012) were used (Supplementary Table S1). PCR was performed in 20 μl reactions containing 10 μl RedTaq mix (SIGMA-ALDRICH, Brøndby, Denmark), 0.25 mM of each primer and 50–100 ng of DNA. The amplification parameters were: 94 °C for 5 min, followed by 38 cycles (94 °C for 30 s, 60 °C or 30 s, 72 °C for 45 s) and finally 5 min at 72 °C. Sequencing was conducted using the commercial service provided by Macrogen (Amsterdam, Holland). The forward primers Ang-ATP6-L and ATP-ATP5c1-6L were used for sequencing. Sequence identity was assessed using BLAST (Altschul et al., 1997) (http://blast.ncbi.nlm.nih.gov) for the ATP6 gene and by evaluating the genotypes at the two quasi-diagnostic SNPs for ATP5c1 gene.

For all analyzed hybrids, the probability of obtaining the observed ATP5c1 genotypes was calculated based on the assigned hybrid classes, and the observed species-specific allele frequencies of ATP5c1 from Gagnaire et al. (2012) (American=1/120 and European eel=76/0 (the numbers correspond to the frequency of European and American alleles observed in each species)). Hardy–Weinberg proportions were assumed. The subsequent P-values were used to infer deviations from neutral expectations and in combination with the mitogenome data used to assess possible cytonuclear incompatibility.

Results

Species distribution in the Sargasso Sea and identification of hybrids

Species identification of the eel larvae collected in 2014, based on the mitochondrial cytochrome b gene, revealed a higher proportion of American than European eel in the two westernmost transects (Figures 2a and b). The frequency of European eel increased in the eastward direction and reached 100% for the four easternmost transects (Figure 2b). At transects 1–4, larvae from both species were found, indicating some spatial overlap of spawning grounds. Among 460 larvae, 358 showed the European eel CYTB gene, whereas 102 showed the American eel CYTB gene. The difference in the number of European and American eel caught is to some extent related to sampling, as the easternmost localities sampled are expected to be dominated by European eel. A qualitatively similar distribution of species in the Sargasso Sea was observed in 2007 (Munk et al., 2010), although the hybrid larvae were caugth slightly more westwards in 2007 (Figure 2b).

Of the 96 SNPs genotyped on the Fluidigm array, 11 showed either no variation or had high levels of missing data across individuals and were omitted from subsequent analyses. The final SNP data set encompassed 85 species-diagnostic SNPs. A total of 12 larvae showed a high proportion of missing data and were omitted from all further analyses. Five of these represented larvae of other deep sea (non-Anguilla) eel species spawning in the Sargasso Sea, as further validated by sequencing the mitochondrial 16S rRNA gene using the L1854 and H3059 primers (Aoyama et al., 2000) Supplementary Table S2). A total of 460 eel larvae from 2014 were analyzed.

Among the 460 larvae sampled in 2014, a total of 2 individuals were identified as hybrids. Both showed the European eel mitochondrial genotype (Figure 2b). Using NEWHYBRIDS, the two hybrids were classified as second-generation backcrosses (bAR × AA) between a first-generation backcross (American eel × F1 hybrid) male and a pure European eel female. This hybrid class was the same, as for the two hybrid larvae sampled in 2007 (sample L125 and L211, Figure 2b) (Pujolar et al., 2014a), that were also reanalyzed in this study. No hybrids were found among the new samples analyzed from Morocco or the Faroe Islands. The percentage of bAR × AA hybrids in the Sargasso Sea was lower in the 2014 sample compared with 2007: 0.43% (2/460) vs 2.17% (2/92). This was also the case when only using the more abundant European eel samples: 0.56% (2/358) vs 4.08% (2/49). However, in both cases the difference was not significant (Fisher’s exact test; P=0.077 and P=0.136).

Among the reanalyzed Icelandic hybrids, 13 were classified as F1 hybrids, 1 as bAA and 3 as bAR × AA, as previously found by Pujolar et al. (2014a) analyzing the same individuals. Further analyses of the 7 bAR × AA hybrids (3 from Iceland and 4 from the Sargasso Sea) showed that at 7 of the analyzed 85 loci, all individuals were homozygous for the European eel-specific allele (see Supplementary Table S5). The probability that all 7 individuals are homozygous at the same locus is between 1.09 × 10−4 and 4.15 × 10−5 depending on the locus (Table 3, see Supplementary Notes S1 and S2 for specific details about the calculations) and therefore highly unlikely to occur by chance. Six out of seven of these SNPs were located inside annotated genes. Among these genes, four included GO terms related to development, and one gene related to phosphorylation (Supplementary Table S3). The frequency associated with different types of gene categories was higher among the seven SNPs but did not differ significantly compared with the rest of the analyzed SNPs (Supplementary Table S4).

Table 3 Allele frequencies of the seven candidate loci for selection and estimated genotype frequencies in bAR × AA hybrids
Full size table

Assessment of cytonuclear selection and incompatibility

Out of the 19 hybrids analyzed, 7 were assigned to the F1 hybrid category, 5 to the second-generation backcross category (bAR × AA) and 7 to the third-generation backcross category or later (Table 2). Estimation of hybrid classes for the RAD sequenced individuals is described in detail in Supplementary Note S3. Out of all hybrids, 3 had American and 16 had European mitochondrial DNA ATP6 haplotypes. Observed genotypes at the nuclear ATP5c1 gene overall followed neutral expectations given the estimated hybrid classes (Figure 3). However, one individual inferred to be a sixth generation or later American backcross (in the following denoted ⩾6th) (RBG10 Table 2 and Figure 3) showed an introgressed European mitogenome and was heterozygote at ATP5c1 (see Supplementary Note S4 and Supplementary Tables S11 and S12 for the calculations of the probabilities). Calculation of the probability for one individual of being an ATP5c1 heterozygote in the sixth or later generation American backcross level is low, with an estimated probability of ⩽0.032 (Figure 3). Overall, all individuals showed a match between mitochondrial and nuclear genotype with one exception: one individual sampled in Europe (LEG29; Table 2) showed a complete mismatch between the ATP6 and ATP5c1 genotypes. This individual was inferred to be a ⩾6th-generation European backcross and was homozygous for the European ATP5c1 genotype but showed an American ATP6 haplotype (Figure 3).

Figure 3
figure 3

Expected and observed genotype frequencies at ATP5c1 of the examined hybrids, given the estimated hybrid classes. The estimated expected frequencies are based on assumed Hardy–Weinberg proportions and previously reported allele frequencies of the ATP5c1 gene in American eel (Gagnaire et al., 2012). The probabilities of the observed genotype combinations are shown to the right and the individual showing ATP6–ATP5c1 mismatch is denoted by an asterisk. Sample sizes (N) are shown in the respective histograms. AA, European eel (A. Anguilla); AR, American eel (A. rostrata).

Full size image

Discussion

Temporally varying hybridization in North Atlantic eels

One interesting finding in our study is the lack of first-generation (F1) hybrids found in the Sargasso Sea in 2014 despite the large number of eel larvae analyzed (N=460). This is in accordance with the study of Pujolar et al. (2014a) and Als et al. (2011), who also reported no F1 hybrids in the Sargasso Sea in 2007, although relying on smaller sample sizes (N=92 and 388). Surprisingly, only bAR × AA hybrids were observed in the Sargasso Sea in both 2007 and 2014. Although the frequency of these hybrids decreased from 2007 to 2014, the differences were not significant and the result does not in itself favor temporally varying hybridization. Nonetheless, although F1 hybrids of Atlantic eel may show an increase in fitness, the so-called hybrid vigor (see, for example, Templeton, 1986; Johnson et al., 2010), post-F1 hybrids are likely to be selected against because of extrinsic selection, or intrinsic selection because of, for example, genetic incompatibilities. Indeed, the presence of strong postzygotic barriers seems supported in Atlantic eels that, despite low average genetic differentiation, show many diagnostic and near-diagnostic SNPs distributed across the genome (Jacobsen et al., 2014a). The results in this study also support a postzygotic component of selection in hybrids, as 7 of the analyzed 85 loci showed strong evidence for selection for the European eel-specific allele. More direct support for postzygotic selection can be found in Iceland. Here both amplified fragment length polymorphism- and SNP-based studies of glass and yellow eels report more than twofold higher frequencies of F1 hybrids compared with combined first and second backcrosses (Albert et al., 2006; Pujolar et al., 2014a). Besides being the only place where F1 hybrids have been found, Iceland is also the only locality where first- (bAA) and second-generation backcrosses (bAR × AA) have been detected outside the Sargasso Sea. Although this may be because of differences in hybrid frequencies between Iceland and Continental Europe and North America, there is currently no evidence for this. In fact, out of several studies investigating hybridization in Atlantic eels from the mainland (excl. Iceland) using microsatellite (Als et al., 2011: N=1010), amplified fragment length polymorphism (Albert et al., 2006: N=379) or SNP markers (Pujolar et al., 2014a, 2014b: N=310), none have observed any such individuals. Given these observations, we find it likely that F1 hybrids should exist in higher frequency compared with later-generation backcrosses if hybridization is continuous over time. As such, we find that the apparent lack of F1 hybrids in the Sargasso Sea in 2007 and 2014 is surprising and best explained by temporally varying hybridization.

Temporally varying hybridization may be mediated by differences in temporal or spatial overlap of spawning between the species over time. This could be a result of annual variation in the location of the thermal fronts (Pujolar et al., 2014a) that show considerable differences in intensity and geographical position across years (Tesch, 2003; Ullman et al., 2007; Munk et al., 2010). Indeed, temperature has increased steadily over the past decades in the Sargasso Sea (Bonhommeau et al., 2008; Huffard et al., 2014), where the 22.5 °C isotherm moreover has moved northwards since ∼1970 (Friedland et al., 2007). As the 22.5 °C isotherm is generally considered to be near the northern limit of spawning of Atlantic eels (Kleckner and McCleave 1988; Tesch and Wegner, 1990), these changes are likely to be associated with a more northern displacement of the shared spawning grounds. As such, the distances that American and European eel need to cover to reach the Sargasso Sea may have changed recently, potentially affecting time of arrival and hybridization, resulting in fewer hybrids being produced. Temporally varying hybridization has also been suggested by Albert et al. (2006) to explain the decrease in proportion of hybrids observed in Iceland from 2000 to 2003 when comparing glass eel samples, a pattern that was further supported when comparing different cohorts comprising glass eels and yellow eels. Year-to-year variation in population density could also explain the lack of F1 hybrids in our study. In this sense, both European and American eel have experienced drastic declines in the past ca. 30 years (Åström and Dekker, 2007), possibly linked to inland pollution, dams and fisheries (Busch and Braun, 2014; Laporte et al., 2016). This decline suggests possible reduced densities of spawning eels that could lead to a reduced spatial overlap and hybridization (Albert et al., 2006; Jacobsen et al., 2014b). As environmental change within the Sargasso Sea itself may affect eel population sizes (Friedland et al., 2007; Bonhommeau et al., 2008), these two scenarios may not be mutually exclusive and the lack of hybrids could be a consequence of different environmental effects.

Sampling bias in either space or time may also explain the lack of F1 hybrids in our study. However, although a distinct hybrid zone in theory could go unnoticed, there is little evidence for this. If it did indeed exist, it would be expected to be present in the core area of the two overlapping spawning areas that historically has been ranging from 70 to 58°W (Miller et al., 2015), and that was extensively sampled during expeditions in both 2007 and 2014 (Als et al., 2011; Pujolar et al., 2014a; this study: Figure 2). Another possibility is that hybridization occurs as a last option, when individuals are unable to locate conspecifics. If so, F1 hybrids should be the product of either European eels arriving too early or American eel arriving too late at the spawning area, or when both species arrive outside the primary spawning season. A recent analysis of all published records of larvae (<10 mm) collected in the Sargasso Sea suggests that the spawning season extends from 13 February to 27 April in American eel and 27 February to 21 July in European eel, with rare spawning events possibly occurring outside these periods (Miller et al., 2015). Thus, given that sampling took place from March to April during the 2007 and 2014 expeditions to the Sargasso Sea, a temporal sampling bias cannot be ruled out. A way of investigating this possibility would be to analyze more samples from Iceland that so far is the only place where F1 hybrids have been observed (Avise et al., 1990; Albert et al., 2006; Gagnaire et al., 2009; Pujolar et al., 2014a). If the apparent lack of F1 hybrids in 2007 and 2014 is true, then this pattern should also be reflected for the same cohorts in Iceland.

Selection at loci in natural hybrids

Identification of hybrids was based on a set of SNPs that showed almost fixed differences between the two species and are likely to mark chromosomal regions under divergent selection in Atlantic eels (Jacobsen et al., 2014a; Pujolar et al., 2014a). The use of genetic markers under selection enhance the ability to determine the population of origin of individuals within species (Nielsen et al., 2012) and detect post-F1 hybrids (Pujolar et al., 2014a). However, it may also involve potential problems when used for hybrid assignment by affecting the genotypic proportions that hybrid identification is based on. Even if near fixation of different alleles in the two species reflects selection, it does not infer whether selection is weak or strong. In fact, even weak directional selection may lead to fixation given sufficient time and low genetic drift and gene flow. Specifically in our study, all putative second-generation bAR × AA backcrosses were homozygous for the European eel allele at the same 7 SNPs (out of 85 SNP loci). This observation is highly unlikely to happen by chance (Table 3) and strongly suggests a role of selection. We find the most parsimonious explanation to be that the seven SNPs are linked, or part of, to functional alleles that are codominant and subject to strong selection between the species. This could result in all genotypes being homozygous, whereas at the same time the genotypes at other loci are more in accordance with neutral expectations, given the specific crosses underlying the hybrids. However, although this may explain the specific genotypes at these loci, omission of the loci would not cause the hybrids to be assigned to a less complex hybrid class than bAR × AA, such as, for example, first-generation AA backcross (F1 × AA). Nonetheless, it demonstrates the potential pitfalls when using such markers in hybrid assessment and suggests an increasing bias when attempting to identify hybrids several generations back in time.

Further investigation of the genomic positions of the 7 SNPs supports that selection is playing a role as (1) 6 out of the 7 SNPs were located in genes (as opposed to being randomly located in noncoding regions) and (2) 4 genes matched GO terms related to development and one gene matched GO terms related to phosphorylation/energy production. This is in concordance with the recent study of Jacobsen et al. (2014a) investigating speciation in North Atlantic eels, in which those two GO term categories were overrepresented among genes associated with candidate SNPs for differential selection showing FST=1 when comparing European and American eel samples.

Cytonuclear incompatibility in North Atlantic eels

The cytonuclear incompatibility hypothesis states that there needs to be a match between coadapted interacting nuclear and mitochondrial genes that might otherwise lead to protein malfunction and consequently reduced fitness (Blier et al., 2001). Results from our study are overall in accordance with this hypothesis as hybrid individuals generally showed a match between the putatively coadapted species-specific nuclear and mitochondrial alleles. However, in most cases the observed ATP5c1 genotypes, found within each hybrid category, followed neutral expectations and do thus neither support nor contradict the cytonuclear hypothesis. Interestingly, however, one American glass eel backcross (RBG10; Table 2 and Figure 3) sampled in Riviere Blanche was heterozygous for the nuclear ATP5c1 gene, but showed an introgressed European mitochondrial DNA ATP6 gene. Considering that this individual was estimated to belong to a ⩾6th-generation backcross, the probability of being heterozygous at ATP5c1 is low (P⩽0.032, Figure 3). This points toward selection for a match between corresponding ATP6 haplotypes and ATP5c1 alleles and thus supports the hypothesis of cytonuclear incompatibility.

On the other hand, one glass eel collected in Ireland and identified as a ⩾6th-generation European backcross (LEG29; Table 2 and Figure 3) represented a complete mismatch, as it showed the American mitochondrial haplotype but was homozygous for the European nuclear allele. This demonstrates that the combination of an American ATP6 haplotype and two European ATP5c1 alleles is not immediately lethal, and this is in opposition to the expectations of cytonuclear incompatibility. This can potentially be explained if incompatibility is incomplete. In such a case, individuals showing mismatch could be viable, although fitness is reduced. Nonetheless, with only one individual in clear favor and only one that contradicts the expected pattern of mitonuclear incompatibility, this possibility cannot be verified and the hypothesis of cytonuclear incompatibility cannot be fully accepted or rejected based on the results in this study.

Nonetheless, the finding of an ATP6–ATP5c1 mismatch is surprising, as the finding of significant interspecies nonsynonymous differentiation in genes involved in the ATP synthase complex (Gagnaire et al., 2012), in combination with studies showing strong evidence for positive selection within ATP6 (Jacobsen et al., 2014b, 2015), is difficult to explain without cytonuclear incompatibility. However, three explanations may account for the finding of a complete ATP6–ATP5c1 mismatch even in case of complete cytonuclear incompatibility: (1) cytonuclear selection might first act in a later life stage beyond the glass eel stage, (2) compensatory interactions involving other nuclear subunits of the ATP synthase protein complex might occur or (3) unidirectional incompatibility. Given the function of the ATP synthase complex, ATP6–ATP5c1 mismatch is expected to lead to reduced energy production (Saraste, 1999). However, it is possible that compensatory mechanisms exist early in life. For example, Drosophila mutants show metabolic compensation through increased glycolytic flux, ketogenesis and Krebs’ cycle activity despite ATP6 dysfunction, although severely deleterious phenotypes are already observed after few days and mutants die prematurely (Celotto et al., 2011). If a similar mechanism acts in Atlantic eels, ATP6–ATP5c1 mismatch may be observed in juvenile individuals, whereas selection may act against this genetic combination later in life. Compensatory mechanisms may also exist, involving other nuclear subunits. In addition to ATP5c1, 13 other nuclear genes are involved in the ATP synthase protein complex (Ballard and Whitlock, 2004). Although transcriptome analysis has shown that ATP5c1 is the most divergent subunit in terms of diagnostic nonsynonymous substitutions between European and American eels, diagnostic changes also exist in at least four other subunits (Gagnaire et al., 2012). Thus, ATP synthase function and hence postzygotic selection may depend on the genotypes of these subunits in combination with ATP5c1. A last possibility is that cytonuclear incompatibility is unidirectional in Atlantic eels, with ATP6–ATP5c1 mismatch being fatal in American eel, but not in European eel. This could potentially explain the observed asymmetrical patterns of introgression, from American to European eel, reported in some studies (Gagnaire et al., 2009; Wielgoss et al., 2014). On the other hand, if unidirectional cytonuclear incompatibility occurred we would expect introgression of the American mitogenome in European eel. However, little evidence of this exists. In fact, the individual sequenced in this study (LEG29) showing ATP6–ATP5c1 mismatch is, to our knowledge, the only European individual to date that show the American mitogenome (see, for example, Avise et al., 1990; Daemen et al., 2001; Jacobsen et al., 2014b) besides F1 hybrids from Iceland (Avise et al., 1990; Albert et al., 2006; this study). As such, we find the possibility of unidirectional cytonuclear incompatibility unlikely.

Conclusion

In conclusion, our results provide new insights into the pre- and post-zygotic barriers acting between the two species of Atlantic eels. The lack of F1 hybrids, but presence of a few later-generation hybrids in the Sargasso Sea in both 2007 and 2014, suggests that interbreeding between the two species is restricted to some years, whereas in other years only backcrosses occur that are the result of interbreeding in past generations. We find the most likely explanation to be interannual variation in the location of the thermal fronts in which spawning takes place that may subsequently lead to increased or decreased overlap of spawning time of the two species. Moreover, overrepresentation of homozygotes for some SNPs used for hybrid assessment was observed that we interpret as divergent selection between species being particularly strong at some loci. On one side, this provides evidence for postzygotic barriers, but on the other side complicates accurate identification of post-F1 hybrid classes. Finally, we find some evidence supporting the suggestion of cytonuclear incompatibility involving the mitochondrial ATP6 gene and the nuclear ATP5c1 interactor gene (Gagnaire et al., 2012), but mismatches between the two are not necessarily lethal and the result should be interpreted with precautions until more data are available. The results obtained in this study are not only important for understanding speciation in Atlantic eels, but also for understanding the variety of isolating barriers acting in nature (Coyne and Orr, 2004; Nosil et al., 2005) and especially in the marine environment (Puebla, 2009). Here, organisms like many fishes and invertebrates show several traits in common with Atlantic eels: they are difficult to rear in captivity, show large population sizes, have no obvious geographic barriers to gene flow and show long larval phases (Palumbi, 1994). All these features promote low genetic drift, efficiency of even low selection and increasing possibilities for speciation with gene flow (Palumbi, 1994; Riginos and Victor, 2001; Jacobsen et al., 2014a). Thus, both the difficulties in studying these organisms and the overall evolutionary forces involved in their speciation may be quite similar to the situation in Atlantic eels.

Data archiving

The sequence data are deposited in Genbank accessions numbers KX818059–KX818096. The FLUIDIGM SNP and the RAD SNP data sets are deposited in DRYAD: DOI:10.5061/dryad.v5m24.