A total of 3.06 million to 15.98 million clean reads passed the data filtering for each sample with an average value of 6.66 million. Among them, 95.17–98.38% were successfully mapped to the reference genome. The Stacks analysis for all samples generated 300,683 RAD loci, the number of which is almost twice that of estimated recognition sites (160523) of EcoRI in the reference genome (two RAD loci should be associated with each restriction site). One admixed individual was identified and removed from the American eel population. A total of 4,311,758 SNPs from 328,631 RAD loci representing 26,229,220 bp were identified in European eel, and 1,829,324 SNPs from 318,424 RAD loci representing 25,414,634 bp were obtained from American eel. After removing SNPs with extreme coverage or deviation from HWE, 4,214,493 and 1,817,481 SNPs remained with an average depth per site of 15.9 and 17.5 for European and American eel, respectively. The values of genome-wide nucleotide diversity (π) were 0.0060 in European eel and 0.0068 in American eel, comparable to estimates reported by Nikolic et al. (2020), also showing slightly higher π in American eel.
Demographic history reconstruction
The folded SFS of European eel and American eel are shown in Fig. 1, both of which exhibit clear L-shapes. Demographic histories of the two species, inferred using Stairway Plot 2 are shown in Fig. 2. Assuming a generation length of 13 years for European eel, the demographic history coinciding with and prior to the presumptive time of speciation was aligned for the two species when a generation length of 10.7 years was assumed for American eel. Conversely, assuming the same generation length for both species led to shifted timing of demographic histories (see Supplementary Information Fig. S1).
Compared to PSMC analysis (Li & Durbin 2011) based on individual genomes (Nikolic et al. 2020), Stairway Plot 2 provides higher resolution and more detailed histories. The time range of histories for both species encompasses ca. 5.00 million years ago (Mya) to the present. The historical Ne values since species divergence are estimated to range from 39,769 to 3,984,743 for European eel, and from 49,793 to 4,999,485 for American eel, with the maximum being about 100 times higher than the minimum for both species. With one notable exception the trajectories for both species are quite similar. Hence, Ne appears stable until ca. 170,000 generations ago, when a decline occurs that represents the minimum Ne for both species. This is followed by an increase of Ne starting ca. 160,000 generations ago. The time of this event corresponds well with the divergence time and hence time of speciation (150,000-160,000 generations ago) inferred using δaδi (Gutenkunst et al. 2009) in Nikolic et al. (2020). The second major event involves strong declines starting ca. 90,000 generations. ago This corresponds to the time (87,000–92,000 generations) of onset of secondary contact between the species, inferred using δaδi (Nikolic et al. 2020). Subsequently, Ne of American eel increases strongly and remains high and stable towards the present. European eel, in contrast, shows a much lower increase, remains stable from ca.70,000 to 40,000 generations ago, possibly with a decline towards the end of this interval. From ca. 40,000 generations ago towards the present, Ne has remained stable after a drastic expansion.
Environmental factors influencing demographic history
Nikolic et al. (2020) provided evidence for speciation of Atlantic eels involving vicariance, followed by isolation and later secondary contact. Our analyses show that both divergence and secondary contact coincide with periods of population decline. Declines at the time of divergence could simply reflect a previously common gene pool splitting into two, whereas the declines at the onset of secondary contact are more surprising; intuitively it could be expected that strong population growth rather than declines would lead to more overlap of spawning regions and hence secondary contact. The life cycle of Atlantic eels is intimately linked to ocean currents, particularly the Gulf Stream (Tesch 2003). It has been suggested that initial divergence of the species was due to ocean current changes leading to isolation of the spawning regions (Avise et al. 1990; Nikolic et al. 2020), possibly mediated by Heinrich Events; massive discharges of icebergs during glaciations, leading to decreased salinity and weakening of the Gulf Stream (Heinrich 1988; Hodell et al. 2008). The second decline could again result from Heinrich Events which caused ocean current changes bringing the spawning regions back into contact. In contrast, the subsequent time period with highly different demographic histories of the species is more likely to reflect factors influencing one species and not the other, hence not involving ocean current changes. Geomagnetism is considered highly important for navigation of eels during their larval spawning migrations (Durif et al. 2021), and hence changes to the Earth’s magnetic fields could have negative effects. The Brunhes-Matuyama magnetic reversal occurred ca. 780,000 years ago (Bassinot et al. 1994), which roughly coincides with the time of divergent demographic histories of the species (Fig. 2). If this has negatively affected Atlantic eels, the effect would be expected to be stronger for European eel (with spawning migration > 5,000 km) compared to American eel (spawning migration > 1,500 km). We note, however, that this event does not entirely overlap with the time span during which the species show different demographic histories, and it can at most be one of several factors involved. Other, as yet unidentified factors affecting the European but not North American continent are also possible.
Based on analysis of mitogenome sequences, it has been suggested that speciation of Atlantic eels occurred ca. 3.38 Mya, coinciding with the closure of the Panama Gateway (Jacobsen et al. 2014b). However, there are no imprints in the demographic histories coinciding with this time (Figure. 1). Strikingly, there are also no discernible demographic effects of cycles of glacial and interglacial periods and even relatively recent drastic events like the Last Glacial Maximum (Fig. 2), as otherwise evident in demographic histories inferred in a range of organisms including freshwater fishes, birds and mammals using PSMC (Moura et al. 2014; Liu et al. 2016; Nadachowska-Brzyska et al. 2016; Pujolar et al. 2017; Mays et al. 2018). The result may seem counterintuitive, given that subfossil data suggests that European eel was absent from the northern part of its present distributional range during the Last Glacial Maximum (Kettle et al. 2008), and similar patterns are expected for American eel. An explanation could be that the Gulf Stream was displaced further South during glaciations, and juvenile eels could have been swept towards more southern regions of North Africa and South America (Kettle et al. 2008), leading to different distributional ranges but without population declines.
Current declines in the context of historical fluctuations
Our results suggest that the current declines of European and American eels do not coincide with general historical trends of declines. On the contrary, both species are suggested to have been relatively stable for at least the last hundreds of thousands of years, and the current declines seem unprecedented in the species’ history. For instance, the 90–99% decrease in European eel recruitment only find parallels when considering events that occurred > 2 Mya. However, such conclusions should of course be treated with caution. First, our results depend on assumptions of generation length and mutation rate, but it would require unrealistic changes of these to lead to qualitatively different conclusions regarding timing of events, and relative magnitudes of Ne through time would remain the same. Second, current declines are estimated as census population sizes (N) encompassing juvenile recruitment or spawning escapement, whereas estimated demographic histories are based on Ne. Comparison of trends at N and Ne assumes a constant ratio between the two, but there is evidence for increased Ne/N ratios in some cases of declining populations (Palstra & Ruzzante 2008), possibly reflecting density-dependent effects. However, there is virtually no knowledge of Ne/N ratios and their constancy or variation in eels. Third and finally, methods like Stairway Plot 2 capture major trends of demographic history and would not be likely to detect fluctuations occurring at the scale of decades or centuries. Nevertheless, the population dynamics and especially long generation lengths of Atlantic eels mean that recovery from low population sizes will be slow. Hence, it has been shown that even if all fishing for European eel was abandoned it would still take 80 years to achieve recovery comparable to pre-1980 levels (Astrom & Dekker 2007). This increases the possibilities that drastic historical declines, if they had occurred, would last sufficiently long to leave detectable imprints in genomes.
In total, our results suggest that the current drastic declines of European and American eel occur against a backdrop of long-term stability. Hence, the species may now be entering unchartered territories, and concerns about low densities of spawners in the huge spawning region resulting in Allee effects must be considered very real.