Wandering behaviour prevents inter and intra oceanic speciation in a coastal pelagic fish

Wandering behaviour prevents inter and intra oceanic speciation in a coastal pelagic fish www.nature.com/scientificreports OPEN Wandering behaviour prevents inter and intra oceanic speciation in a coastal pelagic fish Received: 31 January 2017 1,2 1 1 1 Gonçalo Silva , Regina L. Cunha , Ana Ramos & Rita Castilho Accepted: 26 April 2017 Small pelagic fishes have the ability to disperse over long distances and may present complex Published: xx xx xxxx evolutionary histories. Here, Old World Anchovies (OWA) were used as a model system to understand genetic patterns and connectivity of fish between the Atlantic and Pacific basins. We surveyed 16 locations worldwide using mtDNA and 8 microsatellite loci for genetic parameters, and mtDNA (cyt b; 16S) and nuclear (RAG1; RAG2) regions for dating major lineage-splitting events within Engraulidae family. The OWA genetic divergences (0–0.4%) are compatible with intra-specific divergence, showing evidence of both ancient and contemporary admixture between the Pacific and Atlantic populations, enhanced by high asymmetrical migration from the Pacific to the Atlantic. The estimated divergence between Atlantic and Pacific anchovies (0.67 [0.53–0.80] Ma) matches a severe drop of sea temperature during the Günz glacial stage of the Pleistocene. Our results support an alternative evolutionary scenario for the OWA, suggesting a coastal migration along south Asia, Middle East and eastern Africa continental platforms, followed by the colonization of the Atlantic via the Cape of the Good Hope. Periodic climatic events affect the evolution of the species, shaping their biogeographic and macroecological pattern. S s peciation in the marine realm is mostly related to geological and climatic events that have occurred 2, 3 throughout different periods. Pleistocene climatic events promoted fluctuations in sea surface temperatures, sea level, ice sheet coverage and changes on the global oceanographic circulation patterns, having impact on living species distribution, diversity, population structure and speciation ). ( C e. o gn . r te em f. 4 porary changes on oceanographic features such as in frontal systems, upwelling events and environmental transitions may also have implications on species distribution patterns and life history traits (ref. 4 and references therein). Although cycli- cal periodicity of range shifts may enhance secondary contacts and prevent spe , t cih a e i tio so n lation experienced by some peripheral populations may also promote differentiation. Species that diverged during the Pleistocene often exhibit shallow divergence due to their recent isolation or incipient sp . The e a cia pt pio ar n ent absence of physical barriers in the marine environment creates additional difficulties to explain vicariant all -opatric or peri patric speciation reported during this per . N io od netheless, organisms oen s ft how limited distributional ranges imposed by intrinsic physiological constraints, dispersal ability or the incapacity to adapt to new environments. Throughout the Pleistocene climatic oscillations, organisms were able to cross temporarily interrupted barriers during specific periods of time. Transitions over soft barriers include long-distance migrations through warm Equatorial waters (e.g. ref. 8) or interoceanic migrations (e.g. ref. 9). Here, we used the Old World Anchovies (OWA) species comp a lex s a case study of coastal pelagic fishes with high dispersal ability to analyse evolutionary relationships and the level of connectivity between the Atlantic and Pacific basins. The OWA are coastal fish species distributed along offshore areas above the continental platforms of the Atlantic and Pacific oceans, restricted sea basins in the Mediterranean Sea, the Baltic and the Black seas as well as inshore environments such as estuaries, inlets an . Thi d b s sa pys ecies complex comprises five nominal species: the Japanese anchovy Engraulis japonicus and the Australian anchovy E. australis in the western Pacific; the Cape anchovE. c y apensis in the southeastern Atlantic Ocean; the silver a E. eu ncho rv ysy tole in the western Atlantic Ocean; the European anchovy E. encrasicolus in the eastern Atlantic, Baltic Sea, Mediterranean Sea and the Black Sea. More recently, the white anchovy E. albidus was added to the OWA group, based o - n ecologi cal, morphological and genetic differences between inshore and pelagic populations in the Mediterra nean Sea (Fig. 1), but its specific status remains uncerta . Thin e original description of OWA nominal species was mainly 1 2 Centre of Marine Sciences, CCMAR, University of Algarve, Gambelas, 8005-139, Faro, Portugal. Present address: MARE – Marine and Environmental Sciences Centre; ISPA - Instituto Universitário, Rua Jardim do Tabaco 34, 1149-041, Lisboa, Portugal. Correspondence and requests for materials should be addressed to G.S. (email: gsilva@ispa.pt) Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 1 www.nature.com/scientificreports/ Figure 1. Present-day distribution of nominal Old World Anchovies s; b pel cies ack dots represent sample locations; raw map was downloaded from https://freevectormaps.com/ and edited in Adobe Illustrator CS5.1 (Adobe Systems Inc., CA, USA). based on traditional taxonomy and geography, with taxa longitudinally and latitudinally disjunct and hence clas- sified as separated species. Despite a large body of literature focusing on this group, the phylogenetic relation- ship between the OWA species remain poorly understood. This group is thought to have diverged recently from the remaining Engraulidae at about one million yea . rs ago e O Th WA exhibit shallow genetic divergences between putative species and several conflicting aspects in the taxonomy and molecular classification of the group exist. Based on morphological characters, Whitehead and his colleagues proposed that the whole group should be considered a single species and molecular studies revealed that some of the described species have no genetic support. No significant genetic differences were found between E. encrasicolus and E. eurystole and the existence of shared mtDNA haplotypes between E. encrasicolus and E. 15 16 capensis or between E. australis and E. japonicus, seriously compromise the current taxonomy of this group. The lack of genetic divergence among disjunct distributions contrasts with partial reproductive isolation and 12, 17 parallel genetic differentiation among sympatric ecological morpho o typ r w esith the complex population 14, 18–22 structure observed at regional spatial scales in the European anch . A ov ls yo, E. encrasicolus mtDNA is 13–15 divided in two lineages (clade A and B), but the paraphyly of these mtDNA lineages generated further de . bate Pleistocene climatic swings influenced anchovies range shifts, lead to trans-equatorial dispersals during epi- 15, 16 16 sodes of global cooling or through deep cold wat a er nd promoted inter-oceanic migratio . A nn s chovies are thought to have colonised the Atlantic Ocean through the southern India . n An O cce ho a v n ies in the Atlantic Ocean likely experienced several extinction-colonization cycles at the extremes of the distribution driven by 14, 16 Pleistocene climate shifts . Moreover, the western Atlantic was colonized after the last glacial maximum (LGM) from anchovies dispersing from western African populations of the European a. R nce hcen ovy tly, a mitochondrial-based analysis of E. encrasic (mi olustochondrial clade B) indicated that some of these individ- uals were found under selection, possibly related to adaptation to cold tem . Ipn t era ht e P ur acific O es cean, the Japanese and Australian anchovies diverged between 105 ka (thousand years ago) an . The g d 420 ka enetic signature of Pacific anchovies revealed persistence on separated hemispheres over several glacial cycles, although more recent dispersals were identi. fie u Th ds far, studies involving OWA phylogenetic assessments were based on 15, 16 allozymes and on a small fragment of the mitochondrial cytochrome b gene (cyt b; .521 bp) In this study, we used a large portion (1044 bp) o f the cyt a b nd 8 nuclear microsatellites (as described in Silva et al. ) to further analyse evolutionary relationships within OWA and provide a novel perspective on the level of connectivity between the Atlantic and Pacific anchovies. We dated main lineage splitting events within Engraulidae and propose a new biogeographic scenario for the OWA. Results Population structure, differentiation and connectivity. Multilocus genotypes from the 16 locations from the Atlantic Ocean and Pacific Ocean were obtained for 462 anchovies (Supplement S1 a). Th ry Fig e num .  ber of alleles per locus varied from 19 (locus Ee2–91b) to 85 (locus Ee10) over all locations (Supplem S1 en ).t ary Table  Mean allelic richness, standardized for comparison across a minimum common sample size of nine individuals, ranged from 6.9 (Senegal) to 9.3 (USA) in the Atlantic Ocean and from 7.4 to 9.4 in the Pacific Ocea1 n (T ). able  POWSIM analysis indicated at least a 95% probability (based on the proportion of significant chi-squared tests) of detecting an F value of as low as 0.001 using the microsatellite data set (Supplemen S2 t). E ary Fig xpec.  ted ST heterozygosity (H ) varied between 0.797 (Senegal) and 0.894 (USA) in the Atlantic Ocean and between 0.849 and 0.899 in the Pacific Ocean, while the observed heterozygosit) va y (H ried between 0.653 (Guinea-Bissau) and 0.825 (Portugal north) in the Atlantic Ocean and between 0.699 and 0.742 in the Pacific Ocea1 n (T ). able  The DAPC results show that E. encrasicolus, E. eurystole and E. capensis are closer to each other than to the E. australis and E. japonicus clusters (Fig 2a .  ). The probability of assignment of individuals to their nominal species (Fig.  2b) shows that E. eurystole and E. capensis individuals were mostly assigned to the E. encr asicolus clusters with probabilities close to 0.8, and E. au a sn trd a lE is. japonicus were mostly assigned to their own clusters Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 2 www.nature.com/scientificreports/ Mitochondrial Cytochrome b Microsatellites Location Code Nominal species Long Lat Year N n h π N Aavg A = 9 A = 18 Effnum H H h r r O E Norway NO E. encrasicolus 10.6 59.0 2007 24 17 0.953 0.009 40 13.13 7.38 10.25 6.973 0.714 0.850 Poland PL E. encrasicolus 16.5 54.6 2008 9 7 0.917 0.014 9 8.38 8.38 — 6.081 0.736 0.864 English Channel EC E. encrasicolus 0.1 50.8 2007 27 18 0.963 0.010 45 13.00 7.25 9.75 6.596 0.729 0.837 Bay of Biscay BB E. encrasicolus −2.9 43.5 2007 23 19 0.980 0.015 45 16.00 7.88 10.38 6.926 0.762 0.837 Portugal - North PN E. encrasicolus −8.8 40.7 1998 25 24 0.997 0.013 45 17.25 8.25 11.50 8.383 0.825 0.882 Portugal - South PS E. encrasicolus −8.4 37.1 2007 29 27 0.995 0.011 43 17.88 7.63 12.63 8.588 0.777 0.880 Canary Islands CA E. encrasicolus −15.0 28.3 1999 24 23 0.996 0.007 42 17.75 8.38 12.13 9.929 0.792 0.888 Senegal SN E. encrasicolus −17.6 14.8 1999 25 25 1.000 0.006 37 12.50 6.88 9.38 4.973 0.693 0.797 Guinea-Bissau GU E. encrasicolus −14.2 9.7 2006 20 20 1.000 0.006 19 13.13 8.38 12.88 6.751 0.653 0.868 Ghana GH E. encrasicolus 0.0 5.6 2008 25 25 1.000 0.006 27 15.50 8.88 12.88 8.664 0.766 0.883 Namibia NM E. capensis 11.7 −17.2 2007 24 17 0.920 0.010 32 14.38 8.88 11.63 8.375 0.769 0.875 South Africa SA E. capensis 21.0 −34.7 2007 13 10 0.923 0.019 21 12.63 9.13 11.50 8.148 0.768 0.887 USA US E. eurystole −66.1 41.5 2006 12 9 0.909 0.004 18 13.50 9.25 13.50 7.978 0.793 0.894 Japan JP E. japonicus 139.9 35.6 2006 24 24 1.000 0.011 30 17.13 7.75 13.00 9.014 0.699 0.869 Australia AU E. australis 151.0 −35.0 2008 34 33 0.998 0.007 44 21.38 9.38 13.63 10.530 0.742 0.899 New Zealand NZ E. australis 175.0 −36.7 2005 35 15 0.709 0.001 34 14.38 7.38 11.38 6.771 0.717 0.849 Total 373 269 0.993 0.023 531 47.63 32.00 39.63 6.997 0.746 0.866 Table 1. Sample locations, sample abbreviations, collection dates, sample sizes and summary statistics for a 1044b p sequence fragment of the mtDNA cytochrome b and eight nuclear microsatellites of Old World Anchovies (OWA). Long, longitude; Lat, latitude; N, number of individ, n ua um ls; b ner of haplotypes; h, haplotype diverstity; π, nucleotide diversit, a y; v Aa era vg ge number of alleles, , aAr llelic richness; Effnum , Effective number of alleles; , o H bserved mean heterozigozity, exp ; H ected mean heterozygozity. 0 E ab 0.8 0.6 0.4 0.2 NO EC BB PN PS CA SN GU GH NM SA US JP AU NZ 0.8 0.6 0.4 E. encrasicolus 0.2 E. capensis E. eurystole individuals Cluster 1 E. japonicus Cluster 2 E. australis Figure 2. (a) Discriminant analysis of principal components (DAPC) of multi-locus Old World Anchovies genotypes; individual genotypes appear as circles; ellipses represent the centre of dispersion of each putative species (see Fig.  1 and Table  1). Horizontal and vertical axes are the first two principal components, respectively; (b) scatter plot of a pri -defin ori ed nominal species; (c) individuals assigned to their genetic cluster without forcing them into pre-determined groups. also with high probabilities. The a posterior D i APC analysis (Fig 2c .  ) shows two clusters, one comprising 167 E. encrasicolus individuals, and all E. eurystole, E. capensis, E. australis and E. japonicus, and another with 157 E. encrasicolus. Estimates with microsatellite data of the mutation scaled population size parameter were higher in the Pacific Ocean than in the Atlantic Ocean (Θ = 11.6 and Θ = 98.4). Our comparison of candidate models of gene Atlantic Pacific flow between populations, clearly rejects panmixia, and showed that migration from the Pacific to the Atlantic (model 2) fitted our microsatellite data best (T 2), b ab u le  t with extremely low gene flow ( = M 1.5). A total of 373 individuals from five OWA putative species and 16 locations were analysed for mitochondrial cytb gene, yielding 269 haplotypes. Haplotype diver hs )i w ty a ( s generally high, ranging from 0.917 to 1.000 in the eastern Atlantic Ocean (from Norway to South Africa), 0.909 in the western Atlantic, and from 0.709 to 1.000 in the Pacific Ocean (Table  1). Nucleotide diversity (π) was low, ranging from 0.4% (USA) to 1.9% (South Africa) in Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 3 membership probability PL www.nature.com/scientificreports/ Model Marker Models parameters Bézier dBézier Probability ATL <−> Model 1 **** −5339.10 0.00 1 PAC Model 2 ATL <− PAC *0** −5396.31 −57.20 1.4342E-25 mtDNA Model 3 ATL −> PAC **0* −5375.88 −36.77 1.071E-16 Model 4 (ATL + PAC) *0*0 −5524.26 −185.15 3.87457E-81 ATL <−> Model 1 **** −329230.64 −232312.46 0 PAC Model 2 ATL <− PAC *0** −96918.17 0.00 1 Microsatellites Model 3 ATL −> PAC **0* −205619.62 −108701.45 0 Model 4 (ATL + PAC) *0*0 −384554.60 −287636.43 0 Table 2. Bayes factors model comparison of migration models for Old World Anchovies between the Atlantic (ATL) and the Pacific (PAC) oceans. Model parameters code as follows: asterisk (*) indicates that a particular migration rate was estimated by the model and 0 indicates that no migration was allowed. The first sign indicates theta for Atlantic, the second sign migration to Atlantic, the third sign theta for Pacific and the forth sign migration to Pacific. E. encrasicolus E. capensis E. eurystole E. japonicus E. australis E. encrasicolus 0.000 0.001 0.000 0.003 E. capensis 0.000 0.001 0.000 0.003 E. eurystole 0.001 0.001 0.001 0.003 E. japonicus 0.000 0.000 0.001 0.002 E. australis 0.003 0.003 0.004 0.002 Table 3. Estimates of net evolutionary divergence between putative species of Old World Anchovies (below diagonal) and standard error values (above the diagonal). the Atlantic Ocean and from 0.1% to 1.1% in the Pacific Ocean (T 1a ). E ble  volutionary net divergence between OWA putative species varied between 0–0.4% (T3 a). ble  e h Th aplotype network revealed four main clades, two in Atlantic Ocean and two in the Pacific Oce 3). an (Fig .  Clades from the two oceanic basins were separated by a minimum of 30–32 mutations, with two individuals from South Africa clustering within the northern Pacific clade. Within the Atlantic Ocean, putaE. en tive s crp as ecies i- colus, E. capensis and E. eurystole were not reciprocally monophyletic, but alternatively fit on European anchovy clades A and B. The southern Pacific clade includes haplotypes mostly from E. australis, but also two haplotypes fromE. ja ponicus that are separated nine mutations from the most frequent haplotype of this clade. es Th e likely represent intermediate haplotypes between the two clades. e Th four clades exhibited different haplotype patterns: E. encrasicolus clade A and E. australis were characterized by multiple star-like radiations with relatively shallow genetic divergences; E. encrasicolus clade B and E. japonicus lacked distinct star patterns and exhibited possi- bly unsampled or extinct haplotypes (Fig 3). P.  airwise-taxa and pairwise-locatio , n GF and D indicated ST ST EST that genetic differentiation among putative species was low, with the Atlantic putative species (E. encrasicolus, E. capensis and E. eurystole) and the japanese anchovies (E. japonicus) showing the lowest genetic differentiation between sea basins (Supplementary Fig S3)..  Estimates with mtDNA data of the mutation scaled population size parameter were the same between ocean basins (Θ = 0.1). Our comparison of candidate models of gene flow between populations, clearly rejects panmixia and showed that the full gene flow model (model 1) fitted our mtDNA data best (T 2), w ab it le  h highly asymmetri- cal immigration between ocean basins, with the Pacific population providing five times as many immigrants into the Atlantic Ocean than vice-vers= a ( 14.5 M vs. 2.9, respectively). Phylogenetic analyses and Divergence Time Estimates. e Th resulting topologies inferred from the BI and ML analyses based on the two data sets (with and without the mutation in the codon 368 of the analysed portion of cyt b) were identical (Fig 4 an . d Supplementary Figs  S4 and S5). Bayesian (−ln L = 21,088.60) and ML (−ln L = 20,986.77) analyses arrived at similar topologie 4 s a (F nid g S .  upplementary FiS g4 . , respectively) and therefore BI is shown in Supplementary Information (Supplementa S4 r). P y Fig ot .  ential Scale Reduction Factors in the BI analysis were about 1.00 and the average ESS value for all parameters (1385.01) was significantly larger than 200, which indicates convergence of the runs. In ML analysi4 s (Fig ), four m .  ain monophyletic groups were retrieved within the OWA Engrau sl p ip s . complex corresponding to E. japonicu E. a s, ustralis, a third clade includ- ing specimens assigned to E. encrasicolus clade A and a fourth clade including specimens assigned to E. encra- sicolus clade B. Engraulis eurystole specimens grouped within E. encrasicolus clade A, and E. capensis specimens grouped both within clade A and B of E. encrasicoE lu n sg . raulis encrasicolus clades A and B were well supported in both analyses (clade A: Bayesian posterior probabilit= ies, B 87, a Pn P d ML bootstrap proportions, B = 59; P clade B: BPP = 100, BP = 89). Phylogenetic relationships within Engraulidae, either with ML or BI approaches, were Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 4 www.nature.com/scientificreports/ Figure 3. Minimum-spanning tree for the Old World Anchovies, based on the mitochondrial cytochrome b (1044 bp, 373 individuals). The colour and the size of the circles represent the geographic source (according to Fig. 1, with the exception of New Zealand haplotypes which are represented here in light brown) and frequency of each haplotype, respectively. The smallest colored circles represent a singleton haplotype. Black circles represent either extant unsampled sequences or extinct ancestral sequences. The length of the lines connecting haplotypes is proportional to the number of mutational differences separating the haplotypes, except when otherwise indicated. Figure 4. Phylogenetic relationships of Old World Anchovies Engrau . m lisaxim spp um likelihood analyses inferred from a fragment of 1044 bp o f cyt b. Maximum likelihood bootstrap values for major supported clades larger 50% are shown above branches. Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 5 www.nature.com/scientificreports/ Figure 5. (a) Bayesian dating analysis of the Engraulidae family based on the concatenation of mitochondrial (cytochrome b: 1131 bp and 16S rRNA: 800b p) and nuclear (RAG1: 1480 bp a nd RAG2: 1221 bp) gene fragments; (b) inset of the Old World Anchovies (OWA); numbers (1) and (3) are fossil calibrations, while number (2) represents a geologic event calibration; (A) (B) (C) (D) are relevant nodes ages; all branches show BPP > 75% except the node signed with * (53%). generally unresolved and it was not possible to confidently establish the sister group of t sp hp e . co Eng m rau plex lis (Fig. 4 and Supplementary Fig S4 . ). Beast dating analysis based on the concatenated data set of mtD a N n A (c d 16S) a yt b nd nuclear introns (RAG1 and RAG2) fragments of the family Engraulidae estimated the age of the most recent common ancestor (MRCA) of the OWA at 0.67 M a ago [0.53–0.80] (Fig.  5). The age of the MRCA of the Atlantic anchovies was estimated at 0.40 [0.27–0.52] Ma, while the split of the Pacific anchovies was estimated at about 0.51 [0.39–0.63] Ma. Discussion Genetic divergences detected between Atlantic and Pacific anchovies are more compatible with intraspecific divergence (net divergen = ce 0.2%; Table  3) and indicate ongoing gene flow between oceanic basins. These results suggest that reproductive isolation between OWA putative species is not complete, and revealed the existence of regional variants or incipient species. The reconstructed genetic patterns within the OWA species complex found in the present study are in agreement with the findings of Whitehead and his colle ba a sgues ed on morphological characters that consider the OWA as a single species. Regardless of the increased sequence length and the addi- tion of more taxa, the precise origin of the OWA remains uncertain. The colonisation of the Atlantic Ocean likely occurred via South Africa, with anchovies dispersing across the northern Indian Ocean along the continental platforms of South Asia, Middle East and eastern Africa. The Atlantic and Pacific OWA revealed independent demographic histories but with contemporary gene flow. Before addressing the main interpretations and con- clusions of these results, two main caveats must be addressed. First, data from the Pacific Ocean rely on only three sampling sites. Second, despite the rare records in the eastern African coast, no intermediate locations were sampled between the Atlantic and the Pacific oceans, thus we may have a restricted representation of the genetic diversity from the Indo-Pacific region. Nevertheless, the results presented in this study provide compelling evi- dence to set a biogeographical scenario for the OWA. Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 6 www.nature.com/scientificreports/ Previous phylogenetic analyses based on mitochondrial da bt : 521 a (c b yt p) recovered two lineages (clades 13, 15, 16 A and B) within E. encrasicolus that did not group together . The assumption of paraphyly for this species was uniquely based on a neighbour joining analysis that grouped the clade B w japoi n t ic h uE s w . ith relatively low statistical support (BP 53%). A more recent study showed that some individuals belonging to the E. en crasicolus mitochondrial clade B were under selec a tio nd t n herefore it is pivotal to understand if selection plays a - n impor tant role in shaping the evolutionary and demographic history of this complex. All phylogenetic analyses (ML and BI) presented here that included more sequences (55 Engraulidae species) and a larger fragb m p) en tt (1 han 044 previous studies, returned very similar results, regardless the effect of selection over the analysed portion of the cyt b. Moreover, unlike previous results, our ML analyses retrieved E. encrasicolus as monoph4 y), w lethi ic (Fig le .  the BI results do not reject it (SupplementaS4 ry Fig ). .  Our Bayesian dating analysis estimated that the divergence between the Pacific and the Atlantic OWA 13, 16 occurred during the Pleistocene at 0.67 Ma (Fig . 5). This divergence is more recent than previously estimate. d Grant and his colleagues used a fixed mutation rate of 1.9%/Ma based on the separation between Cetengraulis edentulus and C. mysticetus determined by the closure of the Panama seaway. We also used the closure of the Panama Isthmus as a calibration event but added two fossil calibrations that were modelled with a lognormal distribution, which may explain the differences between age estimates from both studies. Our results indicated that the OWA coalesced during the late Pleisto 5cen ). De (Fig espite t .  he few observed exceptions, such as those observed e.g. in the genus Gobiodo th n at are compatible with a Pleistocene divergence time frame [0.01–2.60 M a], most contemporary reef fishes show an earlier origin back to the Pliocene [2.60– 5.30 Ma] or Miocene [5.30–23 M a] (reviewed in ref. ). W 5 ithin the Pacific Ocean, E. japonicus and E. australi s diverged at 0.51 [0.39–0.63 M a] (Fig. 5), earlier than previously estimated [0.11–0.42 Ma] . To the best of our knowledge, values of genetic divergence between OWA putative species based on cyt b (0.0– 0.4%; Table  3) constitute the lowest percentage of sequence divergence reported for intrageneric fis . h species Microsatellite data also showed a high degree of admixture between putative species divided into two clusters (Fig. 2c) that largely coincide with the northern part of the Atlantic Ocean (E. encrasicolus and E. eurystole) and the southeastern Atlantic and Pacific Oceans (E. encrasic E. c olua s, pensis, E. japonicus and E. australis). The DAPC (Fig. 2) also shows the existence of contemporary gene flow between anchovies from the Atlantic and Pacific basins. Low intraspecific genetic differentiation between individuals from distant locations was also detected on 27–29 other coastal fish species . The climate oscillation on relatively short time scales from decades to hundreds of thousands of years promotes shifts in distribution ranges, abundance 4 a(r nd efr . eferences therein) and cyclical population extinctions and recolonizat t io hn as t prevent the formation of deep lineages in small pelagics such as sardines and anchovies . Our results do not unequivocally support the origin of the OWA in the Pacific Ocean as no clear sister group of the OWA emerges from the ML and BI phylogenetic analyses of th r e c egyt io n b (Fig.  4 and Supplementary Fig. S4). Also, we cannot conclude if the Atlantic OWA resulted from a single colonization event or if they are the product of more than one invasion on different timesc g ailes ven the low statistical support of the node corresponding to the Atlantic MRCA in all phylogenetic trees retrieved with different met 4 a ho nd d s (Fig.  Supplementary Fig S4 .  ). Regardless the number of colonisations, the Pacific anchovies may have reached the Atlantic Ocean by three possible colonization pathways: (1) via Cape Horn in South America, (2) across the Bering strait in the Arctic Ocean, and/or (3) coastal dispersal throughout the Indian Ocean via the Cape of Good Hope in South Africa. Dispersal through the Cape Horn was unlikely because the OWA do not occur on the coastline of South America and there are no paleontological records in the area of species belonging to this group. Colonization of the Atlantic Ocean through the Bering Strait would only have been possible during an interglacial period when the ice sheets melted and temperatures raised up between 0 a °C. H nd 2 ow ever, the estimated origin of their MRCA occurred at 0.67 Ma (Fig . 5), which coincided with a glacial period. Moreover, OWA from the north Atlantic and north Pacific are not genetically clos2 e (Figs  and 3; Supplementary Fig S3 . ). The third hypothetical dispersal route from the Pacific to the Atlantic is via the Cape of the Good Hope in South Africa. There is evidence supporting this alternative route around the South Africa gateway: (1) morphological similarities between Pacific and Atlantic Ocean OWA and (2) the predominant direction of OWA migration is from the Pacific to the Atlantic Ocean (Table  2). This dispersal most likely occurred throughout an interglacial period. The estimated divergence that occurred at 0.67 Ma (O WA MRCA; Fig. 5) between anchovies from the Indian and the Atlantic waters was probably promoted by the Günz glacial period that prevented dispersal between the two ocean basins. Studies on the evolution of the OWA performed thus far already suggested this route as the most likely for the Atlantic 15, 16, 31 15, 16 colonisation . Grant and his colleagues postulated that anchovies used this migration path through open-ocean, via temperate Indian Ocean and South Africa. These authors considered that this route would be the only possible colonisation pathway for temperate-water species as occurs with other small pelagic fish (e.g. sardines) , given the high sea surface temperatures (SST’s) in the coastlines of eastern Africa and southern Asia. Until now, few records of OWA were reported in warmer ha,b b it ua tt r secent studies indicated that higher SST’s are not a limiting factor for the OWA, as they can be found in tropical waters both in the Atlantic and Pacific 10, 14, 23 16 oceans . The colonisation pathway proposed by Grant and Bowen would imply more than 8000 km o f open-ocean migration between the western Australia and South Africa without any stepping-stones. Anchovies are coastal pelagic fishes that live in average up to three years and it would be extremely unlikely that anchovies survived to a trans-oceanic migration of this magnitude in a single generation with no s . Ot po en-o povce era sn areas usually exhibit low productivity and scarce food res , w our hic ces h would create additional difficulties to large-scale migrations. Moreover, anchovies from Australia and southeastern Atlantic are not closely related in 15, 16 any of the phylogenetic analyses performed thus fa2 r (Figs  and 3) . Alternatively, we propose a dispersal route of the OWA from the Pacific to the Atlantic across the continental platforms of the northern Indian Ocean, eastern African coast and South Africa, impelled by the North Equatorial Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 7 www.nature.com/scientificreports/ and Agulhas Currents (Supplementary Fig S6). Th.  e existence of shared microsatellite alleles (Fig 2) an . d haplo- types (Fig.  3) between Atlantic and Pacific OWA revealed by our data, points to the existence of recurrent migra- tion events and contemporary gene flow between ocean basins. We propose that the Atlantic colonisers were likely seeded from an Indo-Pacific pool that could have included both north and south Pacific anchovy ances- tors, departing from Philippines and Indonesia, using Somalia, Mauritius and Seychelles as stepp. ing-stones Pleistocene colonisations from the Pacific to the Atlantic through the Cape of the Good Hope were inferred for other coastal fish, mostly for tropical species (e.g. ref. 32). Trans-oceanic dispersals of tropical fish between the Indian and the Atlantic oceans are generally conditioned by the cold Benguela upwelling in the western South Africa, but fluctuations on the intensity of this current along the Pleistocene allowed punctuated episodes of dispersal. Apparently, trans-oceanic dispersal of OWA from the Indian to the Atlantic Oceans is not limited by the Agulhas and Angola warm currents, nor the Benguela cold current due to the apparent wide temperature range tolerance of anchovy species. Climatic fluctuations of mid/late Quaternary are likely to have contributed to anchovies demographic expansions/contractions as seen in other species 4 an(r d r ef ef . erences therein) lead- ing to extinction-colonisation cycles at the extreme of the distribution areas and to population connectivity/ differentiation. Conclusions Comprehensive information on coastal fish genetic diversity and dispersal provides a wide-scale perspective of marine species connectivity and evolution. The results of our survey of the OWA complex Engr.au he lir se s pp highlight that (i) coastal fish may disperse along large distances and frequently cross biogeographic barriers while maintaining high levels of genetic diversity and low genetic differentiation among populations; (ii) the dispersal route from the Pacific to the Atlantic Ocean is shared by other coastal fish species (e.g. ref. 32), increasing the sup- port for a common model of intra-specific long distance dispersal; (iii) the permeability of biogeographical barri- ers depends on the relaxation of environmental conditions of ocean fronts; (iv) differentiation patterns depend on an intricate relationship between common geographic distances among populations, main circulation patterns, punctuated effectiveness of the biogeographic soft barriers promoted by climate oscillations and on ecological and biological traits at intraspecific level. Material and Methods Ethics statement. No specific permits were required for the field studies described here, since fish wer- e pur chased at fish markets or were collected on scientific cruises with fishing procedures. We confirm that the study locations were not privately owned or protected, and the field sampling activities did not involve endangered or protected species beyond the focal species. Sample collection, DNA extraction and PCR amplification. Samples of OWA likely representing five putative species were collected from 16 locations from both Atlantic and Pacific oce 1 a an ns d (Fig Tab.  le  1). e Th identification of the sampled specimens was based on their geographical origin, given the lack of morpholog- ical differentiation and diagnosing characters between putative s . Fi pes cies h were purchased at small coastal fish markets, as artisanal fishermen do not venture far, or were collected on scientific cruises (see acknowledge- ments). A small portion of white muscle or fin was preserved in 96% ethanol and st20 or °C. T ed at o− tal DNA extraction, polymerase chain reaction (PCR), purification of the PCR product, sequencing for a fragment of the mitochondrial cyt b (1044 bp) and microsatellite genotyping of eight loci were performed as describe d in Silva et al. . Sequences of the nuclear recombination-activating genes RAG1 (1480 bp) and R AG2 (1221b p) and of the mitochondrial 16S rRNA (800 bp) f rom OWA putative species used in the Bayesian dating analysis were obtained as described in Bloom & Lovejo.y Sequence alignment, population structure, differentiation and connectivity. To determine con- temporary genetic structuring and individual assignments based on the autosomal microsatellite data set, we used 36 37 discriminate analysis of principal components (DAPC) implemented in the adegenet p oac f R 2.15.3 kage and following Silva et. al . DAPC was chosen over Bayesian clustering methods because this method is model free does not assuming Hardy–Weinberg equilibrium or linkage disequilibrium being more appropriate for situations where such assumptions are not met, as is oen t ft he case with ancho . vies e Th program POWSIM 4.0 was used to evaluate statistical power for detecting pairwise genetic differentia- tion at F levels ranging from 0.00 to 0.10. We simulated the divergence of three subpopulations, corresponding ST to (1) European anchovy (E. encrasicolu E. c s, apensis and E. eurystole), (2) E. japonicus and (3) E. australis, from a single ancestral population through genetic drift to a given ov va era ll ue defin l F ed by controlling effective ST population size (e) N and number of generations (t). To best reflect the assumingly Nl e ao rf ga e nchovy, we let Ne = 10 000 and varied ftrom 0 to 2078 for simulating different levels of differentiation. Aer t ft he simulation, each subpopulation was sampled a = t n 381 and divergence from genetic homogeneity was tested w χ-exac ith t test. This procedure was repeated 100 times and the proportion of significant outcomes was used to estimate statistical power for detecting pairwise genetic differentiation. Cytb sequences were aligned usiC ng l ustalX 2.0.3 with default settings, implemente G de in ne ious 5.4 , checked and trimmed manually. Sequences were reduced to haplotypes using Collaps . N eum 1.2ber of individuals (N ), number of haplotypes () a n nd haplotype ( ) a h nd nucleotide diversities ( ) wπ ere calculated in Arlequin 3.5.1.2 using the cyt b data set. Summary statistics, number of individua ), a ls ( vN erage number of alleles (A ), observed heterozygosity ( ) H avg O and expected heterozygosity ( ) w H ere calculated for each location and for each locus with Gen . N oet di ve evolutionary divergence between putative species of OWA was calculated on M us E iG ng A t 5 he Maximum Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 8 www.nature.com/scientificreports/ Composite Likelihood model. The rate variation among sites was modelled with a gamma distribution (shape parameter = 1.48). To examine the relationship between mitochondrial haplotypes, a minimum spanning network was con- 42 45 structed witAr h lequin 3.5.1.2 and visualized wit Haps h tar . Pairwise genetic differentiation was estimated 46 47 48 with G and Jost’s D value , both within and between putative species, following Pennin.gs feo t al r st_est est mtDNA and using the R package Diversity for microsatellites. 50, 51 Estimation of migration rates. We used the coalescent-based program Migrate-N to compare different biogeographic hypotheses for the past and present migration of OWA between the Atlantic and the Pacific oceans. We conducted the analyses with two sets of data, mitochondrial DNA and microsatellites, struc- tured into two groups according to geographical regions: southern Atlantic Ocean (pooled southern locations Senegal, Guinea-Bissau, Ghana, Namibia, Angola and South Africa) and Pacific Ocean (Japan, Australia and New Zealand). We tested four variations of the two-population (Atlantic-Pacific) migration model: bidirectional migration (full model, four parameters), strict Atlantic to Pacific migration (three parameters), strict Pacific to Atlantic migration (three parameters) and panmictic model that assumes the Atlantic and Pacific are part of a panmictic population (one parameter). Testing the directionality of gene flow is justified because the dominant ocean current between the ocean basins, the Agulhas flow, runs westerly from the Indian to the Atlantic Ocean 52, 53 and is thought to play a limiting role in marine dispersal in the opposite dir . Ie ni ct tio ian l values were calcu- lated using F. Mutation rates were set to be constant among loci. The Migrate-N run parameters were calibrated ST on the full model for convergence of the Markov chain Monte Carlo sampling method. The prior distributions were uniform for mutation-scaled population size parameter), t s th het at a a (r θe four times the product of the effective population size and the mutation rate, and mutation-scaled migration rates M, that is, immigration rate scaled by the mutation rate, over the ran = g e o 0.05–0.5 a f θ nd prior migration rat= e M 0–100 f or mtDNA, and θ = 10–100 and M= 0–50 for microsatellites. These settings resulted in converged posterior distributions with a clear maximum for each estimate. The Bayesian run for mtDNA consisted of one long chain with a total of 15 million states visited and 50,000 states recorded for the generation of posterior distribution histograms for each locus aer di ft scarding the first 10,000 genealogies as burn-in; for all loci, a total of 48 million states were visited and 160,000 samples were recorded. For all the analyses we used an adaptive heating scheme with four simulta- neous chains using different acceptance ratios (temperature settings were 1.0; 1.5; × 10 3.0; ); t 1h e analyses were run on a cluster computer using 4 compute nodes per run. The Bayesian run for microsatellites consisted of one long chain with a total of 6 million states were visited and 20,000 states were recorded for the gen - eration of pos terior distribution histograms for each locus aer di ft scarding the first 10,000 genealogies as burn-in; for all loci, a total of 48 million states were visited and 160,000 samples were recorded. For all the analyses we used an adaptive heating scheme with four simultaneous chains using different acceptance ratios (temperature settings were 1.0; 1.5; 3.0; 1,000,000.0); the analyses were run on a cluster computer using four compute nodes per run. Overall loci information was combined into a single estimate by Bézier approximation of the thermodynamic scores as described by Beerli & Palczewsk . W i e averaged the Bézier score over three different runs and used as input to evaluate multiple models using Bayes fac . tors Phylogenetic and dating analyses. Phylogenetic relationships within Engraulidae were based on a frag- ment of the mitochondrial c g yt en be (121 taxa, corresponding to 55 Engraulidae species; 1044 bp). A t least one representative species of Engraulidae per genus (except Papueng ) a ra nu dli n s ine randomly chosen specimens from each of the five currently recognized OWA species were included in the phylogenetic analyses, with the exception of E. capensis from which only four specimens were available (accession numbers in Supplementary Table  S2). According to Lavoué et al . , we selected the following outgroup species: Chirocentrus dC or lu ab p, ea harengus, Denticeps clupeoides, Ilisha africana, Sardina pilchardus, Sundasalanx mekongensis. The Akaike 57 58 Information Criterion (AIC i)mplemented in Modeltest 3,. s 7 elected the GT+ RI+Γ as the evolutionary model that best-fitted the data set. The inferred parameters were used in maximum likelihood (ML) and Bayesian Inference (BI) analyses. BI analyses were conducted with MrBayes 3.2.1 . Metropolis-coupled Markov chain Monte Carlo (MCMC) analyses were ran for 20,000,000 generations with sample frequency of 2000. Final trees were calculated aer a b ft urnin of 1,000 generations. PhyML 3.0 was used to estimate the ML tree and to test by non-parametric bootstrapping the robustness of the inferred trees using 1,000 pseudoreplicates. 15, 16 Previous work recovered E. encrasicolus as paraphyletic (specimens assigned to the species grouped into two clades that did not group together). To test if natural selection could interfere with phylogenetic inference, we performed all phylogenetic analyses using the above data set that only included individuals that do not show any evidence of being under selection and repeated the procedures using another data set (116 t bp) t axa; 1044 hat included individuals presenting a mutation in codon 368 of the cyt b as identified in S.il.va et al To estimate the OWA origin and date lineage-splitting events within Engraulidae we used a Bayesian relaxed molecular-clock approach as implemented in Beas t ba 2.1.3 sed on a concatenated dataset of four p - ar tial fragments of mtDNA (cyt : 1131 b bp; 16S: 800 bp) and nuclear (RAG1: 1480 bp; R AG2: 1221 bp) genes. We included sequence data of 49 Engraulidae lineages/t , li axa kely representing 14 genera (out of 16: the monospe- cific Lycothrissa and Papuengraulis genera are not represented), from which 5 are OWA (accession numbers in Supplementary TabS3 le  ). According to our ML and BI analyses, E. capensis an E. eu d rystole are conspecific with E. encrasicolus, and two clades (hereaer c ft lade A and clade B) were recovered within the latter. Hence, to perform the dating analysis we only selected a single representative of E. encrasicolus from clade A and two from clade B, both without the mutation at codon 368. We used the BirthDeath model for the tree prior that assumes that at any point in time, every lineage undergoes speciation at rate λ or goes extinc, a t a n t ra d thr te eμ e calibration points. One refers to the earliest record of Engraulidae [6–12] million years (Ma) from the Miocene - Lower Pliocene of Cyprus . The second calibration corresponds to age estimated for the divergence between Anchovia clu peoides Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 9 www.nature.com/scientificreports/ 13, 64 and A. macrolepidota [2.8–3.1] Ma due to the closure of the Panama seawa . Thy e third calibration corre- sponds to E. japonicus [2–0] Ma from Kokubu group, Japa . C n alibrations using the two fossils were modelled with a lognormal distribution, where 95% of the prior weight fell within the geological interval in which each fossil was discovered. For the Engraulidae [12–6] Ma, the parameters of the lognormal calibration prior were: 95% interval: mean in real space: 1.4, offset: 6.0 and log stdev: 1.0. For E. japonicus [2–0] Ma, the parameters of the lognormal calibration prior were: 95% interval: mean in real space: 0.465, offset: 0 and log stdev: 1.0. For the divergence between Anchovia clupeoide an s d A. macrolepidota we used a calibration according to Lessios . e t al where the closure of the isthmus of Panama occurred between 3.1–2.8 Ma. Log normal calibration was set to: 95% interval: mean in real space: 0.071, offset: 2.8 and log stdev: 1.0. MCMC analyses were run for 20,000,000 genera- tions with a sample frequency of 20,000, following a discarded burn-in of 2,000,000 steps. The convergence to the stationary distributions was confirmed by inspection of the MCMC samples using Tracer 1.6 . ML, BI and dating analyses were performed on the R2C2 research group cluster facility provided by the IT Department of the University of Algarve. Data Accessibility. Sequences of the mitochondrial c (1044 yt b bp) were deposited in GenBank (acces- sion numbers: JQ716609–JQ716731, JQ716748–JQ716756, JX683020–JX683113, KF601435–KF601478 and KJ007642–KJ007734). Sequences of the nuclear recombination-activating genes RAG1 and RAG2 and of the mitochondrial 16S rRNA used in the Bayesian dating analysis were deposited in GenBank (accession num- bers: KX824115–KX824124). The accession numbers of the sequences retrieved from GenBank correspond- ing to the remaining Engraulidae specimens used in all phylogenetic analyses are indicatS2 ed in T and a S3 b les  (SI). Genotypes of individuals obtained from microsatellites are deposited in Dryad repository (http://dx.doi. org/10.5061/dryad.r4f8v). References 1. Jansson, R. & Dynesius, M. The fate of clades in a world of recurrent climatic change: Milankovitch oscillations and evolution. Annual Review of Ecology and Systematics 33, 741–777, doi:10.1146/annurev.ecolsys.33.010802.150520 (2002). 2. Briggs, J. C. Antitropical distribution and evolution in the Indo-West Pacific Ocean. Systemati 36 c B, 237–247, iology doi:10.2307/2413064 (1987). 3. Briggs, J. C. Antitropicality and vicariance. Systematic Bi36 olog , 206–207, do y i:10.2307/2413269 (1987). 4. Henriques, R., Potts, W. M., Santos, C. V., Sauer, W. H. H. & Shaw, P. W. Population connectivity and phylogeography of a coastal fish, Atractoscion aequidens (Sciaenidae), across the Benguela current region: evidence of an ancient vicariant even 9t. , Plos One e87907, doi:10.1371/journal.pone.0087907 (2014). 5. Rocha, L. A. & Bowen, B. W. Speciation in coral-reef fishes. Journal of Fish Biol72 og, 1101–1121, do y i:10.1111/j.1095- 8649.2007.01770.x (2008). 6. Palumbi, S. R. Marine speciation on a small planet. Trends in Ecology and Evo 7l, 114–118, do ution i:10.1016/0169-5347(92)90144-Z (1992). 7. Mirams, A. G. K., Treml, E. A., Shields, J. L., Liggins, L. & Riginos, C. Vicariance and dispersal across an intermittent barrier: population genetic structure of marine animals across the Torres Strait land bridg30 e. C , 937–949, do oral Reefs i:10.1007/s00338- 011-0767-x (2011). 8. Burridge, C. P. Antitropicality of Pacific fishes: molecular insights. Environmental Biology o f65 Fi, 151–164, shes doi:10.1023/A:1020040515980 (2002). 9. Bowen, B. W., Muss, A., Rocha, L. A. & Grant, W. S. Shallow mtDNA coalescence in Atlantic pygmy angelfishes (genus Ce ) ntropyge indicates a recent invasion from the Indian Ocean. Journal of Her 97 ed , 1–12, do ity i:10.1093/jhered/esj006 (2006). 10. Whitehead, P. J., Nelson, W. S. & Wongratana, T. Clupeoid fishes of the World (Suborder Clupeoidei). Vol. 7 (FAO, 1988). 11. Borsa, P., Collet, A. & Durand, J. D. Nuclear-DNA markers confirm the presence of two anchovy species in the Mediterranean. Comptes Rendus Biologies 327, 1113–1123, doi:10.1016/j.crvi.2004.09.003 (2004). 12. Montes, I. et al. Transcriptome analysis deciphers evolutionary mechanisms underlying genetic differentiation between coastal and offshore anchovy populations in the Bay of Biscay. Marine Bio 163 logy, 205, doi:10.1007/s00227-016-2979-7 (2016). 13. Grant, W. S., Lecomte, F. & Bowen, B. W. Biogeographical contingency and the evolution of tropical anchovies (genus ) Cetengraulis from temperate anchovies (genus Engrau ). J lo isurnal of Biogeography 37, 1352–1362, doi:10.1111/j.1365-2699.2010.02291.x (2010). 14. Silva, G., Horne, J. B. & Castilho, R. Anchovies go north and west without losing diversity: post-glacial range expansions in a small pelagic fish. Journal of Biogeography 41, 1171–1182, doi:10.1111/jbi.12275 (2014). 15. Grant, W. S., Leslie, R. W. & Bowen, B. W. Molecular genetic assessment of bipolarity in the anchovy genu . J s oEn urn ga rau l ol fi F s ish Biology 67, 1242–1265, doi:10.1111/j.1095-8649.2005.00820.x (2005). 16. Grant, W. S. & Bowen, B. W. Living in a tilted world: climate change and geography limit speciation in Old World anchovies (Engraulis; Engraulidae). Biological Journal of the Linnean Society 88 , 673–689 (2006). 17. Le Moan, A., Gagnaire, P. A. & Bonhomme, F. Parallel genetic divergence among coastal–marine ecotype pairs of European anchovy explained by differential introgression aer s ft econdary contact. Molecular Ec 25 ol , 3187–3202, do ogy i:10.1111/mec.13627 (2016). 18. Jemaa, S., Bacha, M., Khalaf, G. & Amara, R. Evidence for population complexity of the European anchovy (Engraulis encrasicolus) along its distributional range. Fisheries Res 168 earc, 109–116, do h i:10.1016/j.fishres.2015.04.004 (2015). 19. Kristoer ff sen, J. B. & Magoulas, A. Population structure of anchovy Engraulis encrasicolus L. in the Mediterranean Sea inferred from multiple methods. Fisheries Resear 91 ch , 187–195, doi:10.1016/j.fishres.2007.11.024 (2008). 20. Magoulas, A., Castilho, R., Caetano, S., Marcato, S. & Patarnello, T. Mitochondrial DNA reveals a mosaic pattern of phylogeographical structure in Atlantic and Mediterranean populations o En f a gr n ac u h li o s e vy n ( crasicolus). Molecular Phylogenetics and Evolution 3, 734–746, doi:10.1016/j.ympev.2006.01.016 (2006). 21. Sanz, N., Garcia-Marín, J. L., Viñas, J., Roldán, M. & Pla, C. Spawning groups of European anchovy: population structure and management implications. ICES Journal of Marine Scien 65 ce, 1635–1644, do i:10.1093/icesjms/fsn128 (2008). 22. Viñas, J. et al. Genetic population structure of European anchovy in the Mediterranean Sea and the northeast Atlantic Ocean using sequence analysis of the mitochondrial DNA control region. ICES Journal of Marine S 7, 391–397, do cience i:10.1093/icesjms/fst132 (2013). 23. Silva, G., Lima, F. P., Martel, P. & Castilho, R. Thermal adaptation and clinal mitochondrial DNA variation of European anchovy. Proceedings of the Royal Society of London B: Biological Sciences 281, 20141093, doi:10.1098/rspb.2014.1093 (2014). 24. Liu, J. X. et al. Late Pleistocene divergence and subsequent population expansion of two closely related fish species, Japanese anchovy (Engraulis japonicus) and Australian anchovy (Engraulis austr ). al M is olecular Phylogenetics and Evolution 40, 712–723, doi:10.1016/J. Ympev.2006.04.019 (2006). Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 10 www.nature.com/scientificreports/ 25. Duchene, D., Klanten, S. O., Munday, P. L., Herler, J. & van Herwerden, L. Phylogenetic evidence for recent diversic fi ation of obligate coral-dwelling gobies compared with their host corals. Molecular Phylogenetics and Ev69 olu, 123–132, do tion i:10.1016/j. ympev.2013.04.033 (2013). 26. Cárdenas, L. et al. Origin, diversification, and historical biogeography of the genus Tra (P chu er ru cif s ormes: Carangidae). Molecular Phylogenetics and Evolution 35, 496–507, doi:10.1016/j.ympev.2005.01.011 (2005). 27. DiBattista, eJ t. al. Twisted sister species of pygmy angelfishes: discordance between taxonomy, coloration, and phylogenetics. Coral Reefs 31, 839–851, doi:10.1007/s00338-012-0907-y (2012). 28. Horne, J. B., van Herwerden, L., Choat, J. H. & Robertson, D. R. High population connectivity across the Indo-Pacific: congruent lack of phylogeographic structure in three reef fish congeners. Molecular Phylogenetics and Evo 49 lut , 629–638, do ion i:10.1016/j. ympev.2008.08.023 (2008). 29. Reece, J. S., Bowen, B. W., Smith, D. G. & Larson, A. Comparative phylogeography of four Indo-Pacific moray eel species (Muraenidae) reveals comparable ocean-wide genetic connectivity despite v fi e-fold differences in available adult habitat. Marine Ecology Progress Series 437, 269–277, doi:10.3354/meps09248 (2011). 30. Schwartzlose, R. A. et al . Worldwide large-scale fluctuations of sardine and anchovy populations. South African Journal of Marine Science 21, 289–347, doi:10.2989/025776199784125962 (1999). 31. Grant, W. S. & Bowen, B. W. Shallow population histories in deep evolutionary lineages of marine fishes: insights from sardines and anchovies and lessons for conservation. e J Th ournal of Hered89 ity, 415–426,  doi:10.1093/jhered/89.5.415 (1998). 32. Lessios, H. A. & Robertson, D. R. Speciation on a round planet: phylogeography of the goatfish genus MulloJidic ourn hal th of ys. Biogeography 40, 2373–2384, doi:10.1111/jbi.12176 (2013). 33. Sigman, D. M. & Hain, M. P. The biological productivity of the ocean. Nature Education Knowle3 d, 21 (2012). ge 34. Marlow, J. R., Lange, C. B., Wefer, G. & Rosell-Melé, A. Upwelling intensification as part of the Pliocene-Pleistocene climate transition. Science 290, 2288–2291, doi:10.1126/science.290.5500.2288 (2000). 35. Bloom, D. D. & Lovejoy, N. R. The evolutionary origins of diadromy inferred from a time-calibrated phylogeny for Clupeiformes (herring and allies). Proceedings of the Royal Society B: Biological Scienc 281 es , doi:10.1098/rspb.2013.2081 (2014). 36. Jombart, T. Adegenet: a R package for the multivariate analysis of genetic Bm ioi a nfor rker ms. atics 24, 1403–1405, doi:10.1093/ bioinformatics/btn129 (2008). 37. R Development Core Team. R: A language and environment for statistical computing. (2013). 38. Zarraonaindia, I., Pardo, M. A., Iriondo, M., Manzano, C. & Estonba, A. Microsatellite variability in European anchovy (Engraulis encrasicolus) calls for further investigation of its genetic structure and biogeography. ICES Journal of Ma66 rin , 2176–2182, e Science doi:10.1093/icesjms/fsp187 (2009). 39. Ryman, N. & Palm, S. POWSIM: a computer program for assessing statistical power when testing for genetic differentiation. Molecular Ecology Notes 6, 600–602, doi:10.1111/j.1471-8286.2006.01378.x (2006). 40. Geneious v. 5.1. (Biomatters Ltd, available at: http://www.geneious.com, Auckland, New Zealand), (2011). 41. Posada, D. Collapse 1.2: Describing haplotypes from sequence alignments http://darwin.uvigo.es/software/collapse.html (2004). 42. Excoe ffi r, L. & Lischer, H. E. L. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources10 , 564–567, doi:10.1111/j.1755-0998.2010.02847.x (2010). 43. Meirmans, P. G. & Van Tienderen, P. H. GENOTYPE and GENODIVE: two programs for the analysis of genetic diversity of asexual organisms. Molecular Ecology Notes 4, 792–794 (2004). 44. Tamura, K.e t al. MEGA5: Molecular Evolutionary Genetics Analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Molecular Biology and Evolutio 28 n , 2731–2739, doi:10.1093/Molbev/Msr121 (2011). 45. Teacher, A. G. F. & Griffiths, D. J. HapStar: automated haplotype network layout and visualization. Molecular Ecology R 11, esources 151–153, doi:10.1111/j.1755-0998.2010.02890.x (2011). 46. Hedrick, P. W. & Goodnight, C. A standardized genetic differentiation measure. Ev ol 59 u, 1633–1638, do tion i:10.1554/05-076.1 (2005). 47. Jost, L. G(ST) and its relatives do not measure differentiation. Molecular Ec 17, 4015–4026, do ology i:10.1111/J.1365- 294x.2008.03887.X (2008). 48. Pennings, P. S., Achenbach, A. & Foitzik, S. Similar evolutionary potentials in an obligate ant parasite and its two host species. Journal of Evolutionary Biology 24, 871–886, doi:10.1111/j.1420-9101.2010.02223.x (2011). 49. Keenan, K., McGinnity, P., Cross, T. F., Crozier, W. W. & Prodöhl, P. A. diveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors. Methods in Ecology and E 4, 782–788, do volution i:10.1111/2041- 210x.12067 (2013). 50. Beerli, P. & Felsenstein, J. Maximum likelihood estimation of migration rates and effective population numbers in two populations using a coalescent approach. Genet 152 ics , 763–773 (1999). 51. Beerli, P. & Felsenstein, J. Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach. Proceedings of the National Academy of Scien 9 c8 e, 4563–4568, do s i:10.1073/pnas.081068098 (2001). 52. Ivanova, E. V. The global thermohaline paleocirculation. (Springer, 2009). 53. Lutjeharms, J. R. E. In e s Th ea: the global coastal ocean Vol. 14B e s Th ea - ideas and observations on the progress in the study of the seas (eds Allan, R. Robinson & Kenneth, Brink) 783–834 (Harvard University Press, 2006). 54. Beerli, P. & Palczewski, M. Unified framework to evaluate panmixia and migration direction among multiple sampling locations. Genetics 185, 313–326, doi:10.1534/genetics.109.112532 (2010). 55. Bloomquist, E. W., Lemey, P. & Suchard, M. A. Three roads diverged? Routes to phylogeographic infer Tr en ence ds . in Ecology & Evolution 25, 626–632, 10.1016/j.tree.2010.08.010 (2010). 56. Lavoué, S., Miya, M. & Nishida, M. Mitochondrial phylogenomics of anchovies (family Engraulidae) and recurrent origins of pronounced miniaturization in the order Clupeiformes. Molecular Phylogenetics and Evo 56 lu, 480–485, do tion i:10.1016/j. ympev.2009.11.022 (2010). 57. Akaike, H. A new look at the statistical model identifications. IEEE Transactions on automatic c 19 o, 716–723 (1974). ntrol 58. Posada, D. & Crandall, K. A. Modeltest: testing the model of DNA substituition. Bioi nfor 14, 817–818 (1998). matics 59. Ronquist, Fe. t al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology 61, 539–542, doi:10.1093/sysbio/sys029 (2012). 60. Guindon, S. et al. New Algorithms and Methods to Estimate Maximum-Likelihood Phylogenies: Assessing the Performance of PhyML 3.0. Systematic Biology 59, 307–321, doi:10.1093/sysbio/syq010 (2010). 61. Bouckaert, R e. t al. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Com 10 pu, e1003537, do t Biol i:10.1371/ journal.pcbi.1003537 (2014). 62. Gernhard, T. The conditioned reconstructed process. J Theor Biol 253, 769–778, doi:10.1016/j.jtbi.2008.04.005 (2008). 63. Grande, L. & Nelson, G. Interrelationships of fossil and recent anchovies (Teleostei: Engrauloidea) and description of a new species from the Miocene of Cyprus. American Museum Novitat2826 es  , 1–16 (1985). 64. Lessios, H. A., Kessing, B. D., Robertson, D. R. & Paulay, G. Phylogeography of the pantropical sea urchin in r Euce ida lart is ion to land barriers and ocean currents. Evol u 53 ti, 806–817 (1999). on 65. Yabumoto, Y. Pleistocene clupeid and engraulidid fishes from Kokubu group in Kagoshima Prefecture, Japan. Bulletin of Kitakyushu Museum Natural History 8, 55–74 (1988). 66. Rambaut, A., Suchard, M. A., Xie, D. & Drummond, A. J. Tracer v1.6. Available f h ro tm tp://beast.bio.ed.ac.uk/Tracer (2014). Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 11 www.nature.com/scientificreports/ Acknowledgements We would like to thank to A. Gonçalves and M. Ruano (Portugal), K. Rowling (Australia) and R. Tesker (New Zealand), for collecting and shipping samples (see also refs 14 and 23 acknowledgments). We thank N. Coelho, C. Patrão and M. Valente for technical help. GS doctoral grant (SFRH/BD/36600/2007) was supported by the Portuguese Foundation for Science & Technology (FCT) and presently he is supported through the FCT strategic project UID/MAR/04292/2013 granted to MARE. RLC was supported by a FCT post-doctoral fellowship (SFRH/ BPD/65830/2009). The work was funded by FCT strategic plan UID/Multi/04326/2013 granted to CCMAR. Author Contributions G.S. and R.C. planned the study. G.S. organized the sampling. G.S. and A.R. carried out the molecular work. G.S., R.L.C. and R.C. performed data analyses. G.S. and R.C. made the figures. G.S., R.L.C. and R.C. wrote the manuscript. All authors read and approved the final version of the manuscript. Additional Information Supplementary information accompanies this paper at doi:10.1038/s41598-017-02945-0 Competing Interests: The authors declare that they have no competing interests. Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Cre- ative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not per- mitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. © The Author(s) 2017 Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 12 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Scientific Reports Springer Journals

Wandering behaviour prevents inter and intra oceanic speciation in a coastal pelagic fish

Free
12 pages
Loading next page...
 
/lp/springer_journal/wandering-behaviour-prevents-inter-and-intra-oceanic-speciation-in-a-7BypooRkbb
Publisher
Nature Publishing Group UK
Copyright
Copyright © 2017 by The Author(s)
Subject
Science, Humanities and Social Sciences, multidisciplinary; Science, Humanities and Social Sciences, multidisciplinary; Science, multidisciplinary
eISSN
2045-2322
D.O.I.
10.1038/s41598-017-02945-0
Publisher site
See Article on Publisher Site

Abstract

www.nature.com/scientificreports OPEN Wandering behaviour prevents inter and intra oceanic speciation in a coastal pelagic fish Received: 31 January 2017 1,2 1 1 1 Gonçalo Silva , Regina L. Cunha , Ana Ramos & Rita Castilho Accepted: 26 April 2017 Small pelagic fishes have the ability to disperse over long distances and may present complex Published: xx xx xxxx evolutionary histories. Here, Old World Anchovies (OWA) were used as a model system to understand genetic patterns and connectivity of fish between the Atlantic and Pacific basins. We surveyed 16 locations worldwide using mtDNA and 8 microsatellite loci for genetic parameters, and mtDNA (cyt b; 16S) and nuclear (RAG1; RAG2) regions for dating major lineage-splitting events within Engraulidae family. The OWA genetic divergences (0–0.4%) are compatible with intra-specific divergence, showing evidence of both ancient and contemporary admixture between the Pacific and Atlantic populations, enhanced by high asymmetrical migration from the Pacific to the Atlantic. The estimated divergence between Atlantic and Pacific anchovies (0.67 [0.53–0.80] Ma) matches a severe drop of sea temperature during the Günz glacial stage of the Pleistocene. Our results support an alternative evolutionary scenario for the OWA, suggesting a coastal migration along south Asia, Middle East and eastern Africa continental platforms, followed by the colonization of the Atlantic via the Cape of the Good Hope. Periodic climatic events affect the evolution of the species, shaping their biogeographic and macroecological pattern. S s peciation in the marine realm is mostly related to geological and climatic events that have occurred 2, 3 throughout different periods. Pleistocene climatic events promoted fluctuations in sea surface temperatures, sea level, ice sheet coverage and changes on the global oceanographic circulation patterns, having impact on living species distribution, diversity, population structure and speciation ). ( C e. o gn . r te em f. 4 porary changes on oceanographic features such as in frontal systems, upwelling events and environmental transitions may also have implications on species distribution patterns and life history traits (ref. 4 and references therein). Although cycli- cal periodicity of range shifts may enhance secondary contacts and prevent spe , t cih a e i tio so n lation experienced by some peripheral populations may also promote differentiation. Species that diverged during the Pleistocene often exhibit shallow divergence due to their recent isolation or incipient sp . The e a cia pt pio ar n ent absence of physical barriers in the marine environment creates additional difficulties to explain vicariant all -opatric or peri patric speciation reported during this per . N io od netheless, organisms oen s ft how limited distributional ranges imposed by intrinsic physiological constraints, dispersal ability or the incapacity to adapt to new environments. Throughout the Pleistocene climatic oscillations, organisms were able to cross temporarily interrupted barriers during specific periods of time. Transitions over soft barriers include long-distance migrations through warm Equatorial waters (e.g. ref. 8) or interoceanic migrations (e.g. ref. 9). Here, we used the Old World Anchovies (OWA) species comp a lex s a case study of coastal pelagic fishes with high dispersal ability to analyse evolutionary relationships and the level of connectivity between the Atlantic and Pacific basins. The OWA are coastal fish species distributed along offshore areas above the continental platforms of the Atlantic and Pacific oceans, restricted sea basins in the Mediterranean Sea, the Baltic and the Black seas as well as inshore environments such as estuaries, inlets an . Thi d b s sa pys ecies complex comprises five nominal species: the Japanese anchovy Engraulis japonicus and the Australian anchovy E. australis in the western Pacific; the Cape anchovE. c y apensis in the southeastern Atlantic Ocean; the silver a E. eu ncho rv ysy tole in the western Atlantic Ocean; the European anchovy E. encrasicolus in the eastern Atlantic, Baltic Sea, Mediterranean Sea and the Black Sea. More recently, the white anchovy E. albidus was added to the OWA group, based o - n ecologi cal, morphological and genetic differences between inshore and pelagic populations in the Mediterra nean Sea (Fig. 1), but its specific status remains uncerta . Thin e original description of OWA nominal species was mainly 1 2 Centre of Marine Sciences, CCMAR, University of Algarve, Gambelas, 8005-139, Faro, Portugal. Present address: MARE – Marine and Environmental Sciences Centre; ISPA - Instituto Universitário, Rua Jardim do Tabaco 34, 1149-041, Lisboa, Portugal. Correspondence and requests for materials should be addressed to G.S. (email: gsilva@ispa.pt) Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 1 www.nature.com/scientificreports/ Figure 1. Present-day distribution of nominal Old World Anchovies s; b pel cies ack dots represent sample locations; raw map was downloaded from https://freevectormaps.com/ and edited in Adobe Illustrator CS5.1 (Adobe Systems Inc., CA, USA). based on traditional taxonomy and geography, with taxa longitudinally and latitudinally disjunct and hence clas- sified as separated species. Despite a large body of literature focusing on this group, the phylogenetic relation- ship between the OWA species remain poorly understood. This group is thought to have diverged recently from the remaining Engraulidae at about one million yea . rs ago e O Th WA exhibit shallow genetic divergences between putative species and several conflicting aspects in the taxonomy and molecular classification of the group exist. Based on morphological characters, Whitehead and his colleagues proposed that the whole group should be considered a single species and molecular studies revealed that some of the described species have no genetic support. No significant genetic differences were found between E. encrasicolus and E. eurystole and the existence of shared mtDNA haplotypes between E. encrasicolus and E. 15 16 capensis or between E. australis and E. japonicus, seriously compromise the current taxonomy of this group. The lack of genetic divergence among disjunct distributions contrasts with partial reproductive isolation and 12, 17 parallel genetic differentiation among sympatric ecological morpho o typ r w esith the complex population 14, 18–22 structure observed at regional spatial scales in the European anch . A ov ls yo, E. encrasicolus mtDNA is 13–15 divided in two lineages (clade A and B), but the paraphyly of these mtDNA lineages generated further de . bate Pleistocene climatic swings influenced anchovies range shifts, lead to trans-equatorial dispersals during epi- 15, 16 16 sodes of global cooling or through deep cold wat a er nd promoted inter-oceanic migratio . A nn s chovies are thought to have colonised the Atlantic Ocean through the southern India . n An O cce ho a v n ies in the Atlantic Ocean likely experienced several extinction-colonization cycles at the extremes of the distribution driven by 14, 16 Pleistocene climate shifts . Moreover, the western Atlantic was colonized after the last glacial maximum (LGM) from anchovies dispersing from western African populations of the European a. R nce hcen ovy tly, a mitochondrial-based analysis of E. encrasic (mi olustochondrial clade B) indicated that some of these individ- uals were found under selection, possibly related to adaptation to cold tem . Ipn t era ht e P ur acific O es cean, the Japanese and Australian anchovies diverged between 105 ka (thousand years ago) an . The g d 420 ka enetic signature of Pacific anchovies revealed persistence on separated hemispheres over several glacial cycles, although more recent dispersals were identi. fie u Th ds far, studies involving OWA phylogenetic assessments were based on 15, 16 allozymes and on a small fragment of the mitochondrial cytochrome b gene (cyt b; .521 bp) In this study, we used a large portion (1044 bp) o f the cyt a b nd 8 nuclear microsatellites (as described in Silva et al. ) to further analyse evolutionary relationships within OWA and provide a novel perspective on the level of connectivity between the Atlantic and Pacific anchovies. We dated main lineage splitting events within Engraulidae and propose a new biogeographic scenario for the OWA. Results Population structure, differentiation and connectivity. Multilocus genotypes from the 16 locations from the Atlantic Ocean and Pacific Ocean were obtained for 462 anchovies (Supplement S1 a). Th ry Fig e num .  ber of alleles per locus varied from 19 (locus Ee2–91b) to 85 (locus Ee10) over all locations (Supplem S1 en ).t ary Table  Mean allelic richness, standardized for comparison across a minimum common sample size of nine individuals, ranged from 6.9 (Senegal) to 9.3 (USA) in the Atlantic Ocean and from 7.4 to 9.4 in the Pacific Ocea1 n (T ). able  POWSIM analysis indicated at least a 95% probability (based on the proportion of significant chi-squared tests) of detecting an F value of as low as 0.001 using the microsatellite data set (Supplemen S2 t). E ary Fig xpec.  ted ST heterozygosity (H ) varied between 0.797 (Senegal) and 0.894 (USA) in the Atlantic Ocean and between 0.849 and 0.899 in the Pacific Ocean, while the observed heterozygosit) va y (H ried between 0.653 (Guinea-Bissau) and 0.825 (Portugal north) in the Atlantic Ocean and between 0.699 and 0.742 in the Pacific Ocea1 n (T ). able  The DAPC results show that E. encrasicolus, E. eurystole and E. capensis are closer to each other than to the E. australis and E. japonicus clusters (Fig 2a .  ). The probability of assignment of individuals to their nominal species (Fig.  2b) shows that E. eurystole and E. capensis individuals were mostly assigned to the E. encr asicolus clusters with probabilities close to 0.8, and E. au a sn trd a lE is. japonicus were mostly assigned to their own clusters Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 2 www.nature.com/scientificreports/ Mitochondrial Cytochrome b Microsatellites Location Code Nominal species Long Lat Year N n h π N Aavg A = 9 A = 18 Effnum H H h r r O E Norway NO E. encrasicolus 10.6 59.0 2007 24 17 0.953 0.009 40 13.13 7.38 10.25 6.973 0.714 0.850 Poland PL E. encrasicolus 16.5 54.6 2008 9 7 0.917 0.014 9 8.38 8.38 — 6.081 0.736 0.864 English Channel EC E. encrasicolus 0.1 50.8 2007 27 18 0.963 0.010 45 13.00 7.25 9.75 6.596 0.729 0.837 Bay of Biscay BB E. encrasicolus −2.9 43.5 2007 23 19 0.980 0.015 45 16.00 7.88 10.38 6.926 0.762 0.837 Portugal - North PN E. encrasicolus −8.8 40.7 1998 25 24 0.997 0.013 45 17.25 8.25 11.50 8.383 0.825 0.882 Portugal - South PS E. encrasicolus −8.4 37.1 2007 29 27 0.995 0.011 43 17.88 7.63 12.63 8.588 0.777 0.880 Canary Islands CA E. encrasicolus −15.0 28.3 1999 24 23 0.996 0.007 42 17.75 8.38 12.13 9.929 0.792 0.888 Senegal SN E. encrasicolus −17.6 14.8 1999 25 25 1.000 0.006 37 12.50 6.88 9.38 4.973 0.693 0.797 Guinea-Bissau GU E. encrasicolus −14.2 9.7 2006 20 20 1.000 0.006 19 13.13 8.38 12.88 6.751 0.653 0.868 Ghana GH E. encrasicolus 0.0 5.6 2008 25 25 1.000 0.006 27 15.50 8.88 12.88 8.664 0.766 0.883 Namibia NM E. capensis 11.7 −17.2 2007 24 17 0.920 0.010 32 14.38 8.88 11.63 8.375 0.769 0.875 South Africa SA E. capensis 21.0 −34.7 2007 13 10 0.923 0.019 21 12.63 9.13 11.50 8.148 0.768 0.887 USA US E. eurystole −66.1 41.5 2006 12 9 0.909 0.004 18 13.50 9.25 13.50 7.978 0.793 0.894 Japan JP E. japonicus 139.9 35.6 2006 24 24 1.000 0.011 30 17.13 7.75 13.00 9.014 0.699 0.869 Australia AU E. australis 151.0 −35.0 2008 34 33 0.998 0.007 44 21.38 9.38 13.63 10.530 0.742 0.899 New Zealand NZ E. australis 175.0 −36.7 2005 35 15 0.709 0.001 34 14.38 7.38 11.38 6.771 0.717 0.849 Total 373 269 0.993 0.023 531 47.63 32.00 39.63 6.997 0.746 0.866 Table 1. Sample locations, sample abbreviations, collection dates, sample sizes and summary statistics for a 1044b p sequence fragment of the mtDNA cytochrome b and eight nuclear microsatellites of Old World Anchovies (OWA). Long, longitude; Lat, latitude; N, number of individ, n ua um ls; b ner of haplotypes; h, haplotype diverstity; π, nucleotide diversit, a y; v Aa era vg ge number of alleles, , aAr llelic richness; Effnum , Effective number of alleles; , o H bserved mean heterozigozity, exp ; H ected mean heterozygozity. 0 E ab 0.8 0.6 0.4 0.2 NO EC BB PN PS CA SN GU GH NM SA US JP AU NZ 0.8 0.6 0.4 E. encrasicolus 0.2 E. capensis E. eurystole individuals Cluster 1 E. japonicus Cluster 2 E. australis Figure 2. (a) Discriminant analysis of principal components (DAPC) of multi-locus Old World Anchovies genotypes; individual genotypes appear as circles; ellipses represent the centre of dispersion of each putative species (see Fig.  1 and Table  1). Horizontal and vertical axes are the first two principal components, respectively; (b) scatter plot of a pri -defin ori ed nominal species; (c) individuals assigned to their genetic cluster without forcing them into pre-determined groups. also with high probabilities. The a posterior D i APC analysis (Fig 2c .  ) shows two clusters, one comprising 167 E. encrasicolus individuals, and all E. eurystole, E. capensis, E. australis and E. japonicus, and another with 157 E. encrasicolus. Estimates with microsatellite data of the mutation scaled population size parameter were higher in the Pacific Ocean than in the Atlantic Ocean (Θ = 11.6 and Θ = 98.4). Our comparison of candidate models of gene Atlantic Pacific flow between populations, clearly rejects panmixia, and showed that migration from the Pacific to the Atlantic (model 2) fitted our microsatellite data best (T 2), b ab u le  t with extremely low gene flow ( = M 1.5). A total of 373 individuals from five OWA putative species and 16 locations were analysed for mitochondrial cytb gene, yielding 269 haplotypes. Haplotype diver hs )i w ty a ( s generally high, ranging from 0.917 to 1.000 in the eastern Atlantic Ocean (from Norway to South Africa), 0.909 in the western Atlantic, and from 0.709 to 1.000 in the Pacific Ocean (Table  1). Nucleotide diversity (π) was low, ranging from 0.4% (USA) to 1.9% (South Africa) in Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 3 membership probability PL www.nature.com/scientificreports/ Model Marker Models parameters Bézier dBézier Probability ATL <−> Model 1 **** −5339.10 0.00 1 PAC Model 2 ATL <− PAC *0** −5396.31 −57.20 1.4342E-25 mtDNA Model 3 ATL −> PAC **0* −5375.88 −36.77 1.071E-16 Model 4 (ATL + PAC) *0*0 −5524.26 −185.15 3.87457E-81 ATL <−> Model 1 **** −329230.64 −232312.46 0 PAC Model 2 ATL <− PAC *0** −96918.17 0.00 1 Microsatellites Model 3 ATL −> PAC **0* −205619.62 −108701.45 0 Model 4 (ATL + PAC) *0*0 −384554.60 −287636.43 0 Table 2. Bayes factors model comparison of migration models for Old World Anchovies between the Atlantic (ATL) and the Pacific (PAC) oceans. Model parameters code as follows: asterisk (*) indicates that a particular migration rate was estimated by the model and 0 indicates that no migration was allowed. The first sign indicates theta for Atlantic, the second sign migration to Atlantic, the third sign theta for Pacific and the forth sign migration to Pacific. E. encrasicolus E. capensis E. eurystole E. japonicus E. australis E. encrasicolus 0.000 0.001 0.000 0.003 E. capensis 0.000 0.001 0.000 0.003 E. eurystole 0.001 0.001 0.001 0.003 E. japonicus 0.000 0.000 0.001 0.002 E. australis 0.003 0.003 0.004 0.002 Table 3. Estimates of net evolutionary divergence between putative species of Old World Anchovies (below diagonal) and standard error values (above the diagonal). the Atlantic Ocean and from 0.1% to 1.1% in the Pacific Ocean (T 1a ). E ble  volutionary net divergence between OWA putative species varied between 0–0.4% (T3 a). ble  e h Th aplotype network revealed four main clades, two in Atlantic Ocean and two in the Pacific Oce 3). an (Fig .  Clades from the two oceanic basins were separated by a minimum of 30–32 mutations, with two individuals from South Africa clustering within the northern Pacific clade. Within the Atlantic Ocean, putaE. en tive s crp as ecies i- colus, E. capensis and E. eurystole were not reciprocally monophyletic, but alternatively fit on European anchovy clades A and B. The southern Pacific clade includes haplotypes mostly from E. australis, but also two haplotypes fromE. ja ponicus that are separated nine mutations from the most frequent haplotype of this clade. es Th e likely represent intermediate haplotypes between the two clades. e Th four clades exhibited different haplotype patterns: E. encrasicolus clade A and E. australis were characterized by multiple star-like radiations with relatively shallow genetic divergences; E. encrasicolus clade B and E. japonicus lacked distinct star patterns and exhibited possi- bly unsampled or extinct haplotypes (Fig 3). P.  airwise-taxa and pairwise-locatio , n GF and D indicated ST ST EST that genetic differentiation among putative species was low, with the Atlantic putative species (E. encrasicolus, E. capensis and E. eurystole) and the japanese anchovies (E. japonicus) showing the lowest genetic differentiation between sea basins (Supplementary Fig S3)..  Estimates with mtDNA data of the mutation scaled population size parameter were the same between ocean basins (Θ = 0.1). Our comparison of candidate models of gene flow between populations, clearly rejects panmixia and showed that the full gene flow model (model 1) fitted our mtDNA data best (T 2), w ab it le  h highly asymmetri- cal immigration between ocean basins, with the Pacific population providing five times as many immigrants into the Atlantic Ocean than vice-vers= a ( 14.5 M vs. 2.9, respectively). Phylogenetic analyses and Divergence Time Estimates. e Th resulting topologies inferred from the BI and ML analyses based on the two data sets (with and without the mutation in the codon 368 of the analysed portion of cyt b) were identical (Fig 4 an . d Supplementary Figs  S4 and S5). Bayesian (−ln L = 21,088.60) and ML (−ln L = 20,986.77) analyses arrived at similar topologie 4 s a (F nid g S .  upplementary FiS g4 . , respectively) and therefore BI is shown in Supplementary Information (Supplementa S4 r). P y Fig ot .  ential Scale Reduction Factors in the BI analysis were about 1.00 and the average ESS value for all parameters (1385.01) was significantly larger than 200, which indicates convergence of the runs. In ML analysi4 s (Fig ), four m .  ain monophyletic groups were retrieved within the OWA Engrau sl p ip s . complex corresponding to E. japonicu E. a s, ustralis, a third clade includ- ing specimens assigned to E. encrasicolus clade A and a fourth clade including specimens assigned to E. encra- sicolus clade B. Engraulis eurystole specimens grouped within E. encrasicolus clade A, and E. capensis specimens grouped both within clade A and B of E. encrasicoE lu n sg . raulis encrasicolus clades A and B were well supported in both analyses (clade A: Bayesian posterior probabilit= ies, B 87, a Pn P d ML bootstrap proportions, B = 59; P clade B: BPP = 100, BP = 89). Phylogenetic relationships within Engraulidae, either with ML or BI approaches, were Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 4 www.nature.com/scientificreports/ Figure 3. Minimum-spanning tree for the Old World Anchovies, based on the mitochondrial cytochrome b (1044 bp, 373 individuals). The colour and the size of the circles represent the geographic source (according to Fig. 1, with the exception of New Zealand haplotypes which are represented here in light brown) and frequency of each haplotype, respectively. The smallest colored circles represent a singleton haplotype. Black circles represent either extant unsampled sequences or extinct ancestral sequences. The length of the lines connecting haplotypes is proportional to the number of mutational differences separating the haplotypes, except when otherwise indicated. Figure 4. Phylogenetic relationships of Old World Anchovies Engrau . m lisaxim spp um likelihood analyses inferred from a fragment of 1044 bp o f cyt b. Maximum likelihood bootstrap values for major supported clades larger 50% are shown above branches. Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 5 www.nature.com/scientificreports/ Figure 5. (a) Bayesian dating analysis of the Engraulidae family based on the concatenation of mitochondrial (cytochrome b: 1131 bp and 16S rRNA: 800b p) and nuclear (RAG1: 1480 bp a nd RAG2: 1221 bp) gene fragments; (b) inset of the Old World Anchovies (OWA); numbers (1) and (3) are fossil calibrations, while number (2) represents a geologic event calibration; (A) (B) (C) (D) are relevant nodes ages; all branches show BPP > 75% except the node signed with * (53%). generally unresolved and it was not possible to confidently establish the sister group of t sp hp e . co Eng m rau plex lis (Fig. 4 and Supplementary Fig S4 . ). Beast dating analysis based on the concatenated data set of mtD a N n A (c d 16S) a yt b nd nuclear introns (RAG1 and RAG2) fragments of the family Engraulidae estimated the age of the most recent common ancestor (MRCA) of the OWA at 0.67 M a ago [0.53–0.80] (Fig.  5). The age of the MRCA of the Atlantic anchovies was estimated at 0.40 [0.27–0.52] Ma, while the split of the Pacific anchovies was estimated at about 0.51 [0.39–0.63] Ma. Discussion Genetic divergences detected between Atlantic and Pacific anchovies are more compatible with intraspecific divergence (net divergen = ce 0.2%; Table  3) and indicate ongoing gene flow between oceanic basins. These results suggest that reproductive isolation between OWA putative species is not complete, and revealed the existence of regional variants or incipient species. The reconstructed genetic patterns within the OWA species complex found in the present study are in agreement with the findings of Whitehead and his colle ba a sgues ed on morphological characters that consider the OWA as a single species. Regardless of the increased sequence length and the addi- tion of more taxa, the precise origin of the OWA remains uncertain. The colonisation of the Atlantic Ocean likely occurred via South Africa, with anchovies dispersing across the northern Indian Ocean along the continental platforms of South Asia, Middle East and eastern Africa. The Atlantic and Pacific OWA revealed independent demographic histories but with contemporary gene flow. Before addressing the main interpretations and con- clusions of these results, two main caveats must be addressed. First, data from the Pacific Ocean rely on only three sampling sites. Second, despite the rare records in the eastern African coast, no intermediate locations were sampled between the Atlantic and the Pacific oceans, thus we may have a restricted representation of the genetic diversity from the Indo-Pacific region. Nevertheless, the results presented in this study provide compelling evi- dence to set a biogeographical scenario for the OWA. Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 6 www.nature.com/scientificreports/ Previous phylogenetic analyses based on mitochondrial da bt : 521 a (c b yt p) recovered two lineages (clades 13, 15, 16 A and B) within E. encrasicolus that did not group together . The assumption of paraphyly for this species was uniquely based on a neighbour joining analysis that grouped the clade B w japoi n t ic h uE s w . ith relatively low statistical support (BP 53%). A more recent study showed that some individuals belonging to the E. en crasicolus mitochondrial clade B were under selec a tio nd t n herefore it is pivotal to understand if selection plays a - n impor tant role in shaping the evolutionary and demographic history of this complex. All phylogenetic analyses (ML and BI) presented here that included more sequences (55 Engraulidae species) and a larger fragb m p) en tt (1 han 044 previous studies, returned very similar results, regardless the effect of selection over the analysed portion of the cyt b. Moreover, unlike previous results, our ML analyses retrieved E. encrasicolus as monoph4 y), w lethi ic (Fig le .  the BI results do not reject it (SupplementaS4 ry Fig ). .  Our Bayesian dating analysis estimated that the divergence between the Pacific and the Atlantic OWA 13, 16 occurred during the Pleistocene at 0.67 Ma (Fig . 5). This divergence is more recent than previously estimate. d Grant and his colleagues used a fixed mutation rate of 1.9%/Ma based on the separation between Cetengraulis edentulus and C. mysticetus determined by the closure of the Panama seaway. We also used the closure of the Panama Isthmus as a calibration event but added two fossil calibrations that were modelled with a lognormal distribution, which may explain the differences between age estimates from both studies. Our results indicated that the OWA coalesced during the late Pleisto 5cen ). De (Fig espite t .  he few observed exceptions, such as those observed e.g. in the genus Gobiodo th n at are compatible with a Pleistocene divergence time frame [0.01–2.60 M a], most contemporary reef fishes show an earlier origin back to the Pliocene [2.60– 5.30 Ma] or Miocene [5.30–23 M a] (reviewed in ref. ). W 5 ithin the Pacific Ocean, E. japonicus and E. australi s diverged at 0.51 [0.39–0.63 M a] (Fig. 5), earlier than previously estimated [0.11–0.42 Ma] . To the best of our knowledge, values of genetic divergence between OWA putative species based on cyt b (0.0– 0.4%; Table  3) constitute the lowest percentage of sequence divergence reported for intrageneric fis . h species Microsatellite data also showed a high degree of admixture between putative species divided into two clusters (Fig. 2c) that largely coincide with the northern part of the Atlantic Ocean (E. encrasicolus and E. eurystole) and the southeastern Atlantic and Pacific Oceans (E. encrasic E. c olua s, pensis, E. japonicus and E. australis). The DAPC (Fig. 2) also shows the existence of contemporary gene flow between anchovies from the Atlantic and Pacific basins. Low intraspecific genetic differentiation between individuals from distant locations was also detected on 27–29 other coastal fish species . The climate oscillation on relatively short time scales from decades to hundreds of thousands of years promotes shifts in distribution ranges, abundance 4 a(r nd efr . eferences therein) and cyclical population extinctions and recolonizat t io hn as t prevent the formation of deep lineages in small pelagics such as sardines and anchovies . Our results do not unequivocally support the origin of the OWA in the Pacific Ocean as no clear sister group of the OWA emerges from the ML and BI phylogenetic analyses of th r e c egyt io n b (Fig.  4 and Supplementary Fig. S4). Also, we cannot conclude if the Atlantic OWA resulted from a single colonization event or if they are the product of more than one invasion on different timesc g ailes ven the low statistical support of the node corresponding to the Atlantic MRCA in all phylogenetic trees retrieved with different met 4 a ho nd d s (Fig.  Supplementary Fig S4 .  ). Regardless the number of colonisations, the Pacific anchovies may have reached the Atlantic Ocean by three possible colonization pathways: (1) via Cape Horn in South America, (2) across the Bering strait in the Arctic Ocean, and/or (3) coastal dispersal throughout the Indian Ocean via the Cape of Good Hope in South Africa. Dispersal through the Cape Horn was unlikely because the OWA do not occur on the coastline of South America and there are no paleontological records in the area of species belonging to this group. Colonization of the Atlantic Ocean through the Bering Strait would only have been possible during an interglacial period when the ice sheets melted and temperatures raised up between 0 a °C. H nd 2 ow ever, the estimated origin of their MRCA occurred at 0.67 Ma (Fig . 5), which coincided with a glacial period. Moreover, OWA from the north Atlantic and north Pacific are not genetically clos2 e (Figs  and 3; Supplementary Fig S3 . ). The third hypothetical dispersal route from the Pacific to the Atlantic is via the Cape of the Good Hope in South Africa. There is evidence supporting this alternative route around the South Africa gateway: (1) morphological similarities between Pacific and Atlantic Ocean OWA and (2) the predominant direction of OWA migration is from the Pacific to the Atlantic Ocean (Table  2). This dispersal most likely occurred throughout an interglacial period. The estimated divergence that occurred at 0.67 Ma (O WA MRCA; Fig. 5) between anchovies from the Indian and the Atlantic waters was probably promoted by the Günz glacial period that prevented dispersal between the two ocean basins. Studies on the evolution of the OWA performed thus far already suggested this route as the most likely for the Atlantic 15, 16, 31 15, 16 colonisation . Grant and his colleagues postulated that anchovies used this migration path through open-ocean, via temperate Indian Ocean and South Africa. These authors considered that this route would be the only possible colonisation pathway for temperate-water species as occurs with other small pelagic fish (e.g. sardines) , given the high sea surface temperatures (SST’s) in the coastlines of eastern Africa and southern Asia. Until now, few records of OWA were reported in warmer ha,b b it ua tt r secent studies indicated that higher SST’s are not a limiting factor for the OWA, as they can be found in tropical waters both in the Atlantic and Pacific 10, 14, 23 16 oceans . The colonisation pathway proposed by Grant and Bowen would imply more than 8000 km o f open-ocean migration between the western Australia and South Africa without any stepping-stones. Anchovies are coastal pelagic fishes that live in average up to three years and it would be extremely unlikely that anchovies survived to a trans-oceanic migration of this magnitude in a single generation with no s . Ot po en-o povce era sn areas usually exhibit low productivity and scarce food res , w our hic ces h would create additional difficulties to large-scale migrations. Moreover, anchovies from Australia and southeastern Atlantic are not closely related in 15, 16 any of the phylogenetic analyses performed thus fa2 r (Figs  and 3) . Alternatively, we propose a dispersal route of the OWA from the Pacific to the Atlantic across the continental platforms of the northern Indian Ocean, eastern African coast and South Africa, impelled by the North Equatorial Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 7 www.nature.com/scientificreports/ and Agulhas Currents (Supplementary Fig S6). Th.  e existence of shared microsatellite alleles (Fig 2) an . d haplo- types (Fig.  3) between Atlantic and Pacific OWA revealed by our data, points to the existence of recurrent migra- tion events and contemporary gene flow between ocean basins. We propose that the Atlantic colonisers were likely seeded from an Indo-Pacific pool that could have included both north and south Pacific anchovy ances- tors, departing from Philippines and Indonesia, using Somalia, Mauritius and Seychelles as stepp. ing-stones Pleistocene colonisations from the Pacific to the Atlantic through the Cape of the Good Hope were inferred for other coastal fish, mostly for tropical species (e.g. ref. 32). Trans-oceanic dispersals of tropical fish between the Indian and the Atlantic oceans are generally conditioned by the cold Benguela upwelling in the western South Africa, but fluctuations on the intensity of this current along the Pleistocene allowed punctuated episodes of dispersal. Apparently, trans-oceanic dispersal of OWA from the Indian to the Atlantic Oceans is not limited by the Agulhas and Angola warm currents, nor the Benguela cold current due to the apparent wide temperature range tolerance of anchovy species. Climatic fluctuations of mid/late Quaternary are likely to have contributed to anchovies demographic expansions/contractions as seen in other species 4 an(r d r ef ef . erences therein) lead- ing to extinction-colonisation cycles at the extreme of the distribution areas and to population connectivity/ differentiation. Conclusions Comprehensive information on coastal fish genetic diversity and dispersal provides a wide-scale perspective of marine species connectivity and evolution. The results of our survey of the OWA complex Engr.au he lir se s pp highlight that (i) coastal fish may disperse along large distances and frequently cross biogeographic barriers while maintaining high levels of genetic diversity and low genetic differentiation among populations; (ii) the dispersal route from the Pacific to the Atlantic Ocean is shared by other coastal fish species (e.g. ref. 32), increasing the sup- port for a common model of intra-specific long distance dispersal; (iii) the permeability of biogeographical barri- ers depends on the relaxation of environmental conditions of ocean fronts; (iv) differentiation patterns depend on an intricate relationship between common geographic distances among populations, main circulation patterns, punctuated effectiveness of the biogeographic soft barriers promoted by climate oscillations and on ecological and biological traits at intraspecific level. Material and Methods Ethics statement. No specific permits were required for the field studies described here, since fish wer- e pur chased at fish markets or were collected on scientific cruises with fishing procedures. We confirm that the study locations were not privately owned or protected, and the field sampling activities did not involve endangered or protected species beyond the focal species. Sample collection, DNA extraction and PCR amplification. Samples of OWA likely representing five putative species were collected from 16 locations from both Atlantic and Pacific oce 1 a an ns d (Fig Tab.  le  1). e Th identification of the sampled specimens was based on their geographical origin, given the lack of morpholog- ical differentiation and diagnosing characters between putative s . Fi pes cies h were purchased at small coastal fish markets, as artisanal fishermen do not venture far, or were collected on scientific cruises (see acknowledge- ments). A small portion of white muscle or fin was preserved in 96% ethanol and st20 or °C. T ed at o− tal DNA extraction, polymerase chain reaction (PCR), purification of the PCR product, sequencing for a fragment of the mitochondrial cyt b (1044 bp) and microsatellite genotyping of eight loci were performed as describe d in Silva et al. . Sequences of the nuclear recombination-activating genes RAG1 (1480 bp) and R AG2 (1221b p) and of the mitochondrial 16S rRNA (800 bp) f rom OWA putative species used in the Bayesian dating analysis were obtained as described in Bloom & Lovejo.y Sequence alignment, population structure, differentiation and connectivity. To determine con- temporary genetic structuring and individual assignments based on the autosomal microsatellite data set, we used 36 37 discriminate analysis of principal components (DAPC) implemented in the adegenet p oac f R 2.15.3 kage and following Silva et. al . DAPC was chosen over Bayesian clustering methods because this method is model free does not assuming Hardy–Weinberg equilibrium or linkage disequilibrium being more appropriate for situations where such assumptions are not met, as is oen t ft he case with ancho . vies e Th program POWSIM 4.0 was used to evaluate statistical power for detecting pairwise genetic differentia- tion at F levels ranging from 0.00 to 0.10. We simulated the divergence of three subpopulations, corresponding ST to (1) European anchovy (E. encrasicolu E. c s, apensis and E. eurystole), (2) E. japonicus and (3) E. australis, from a single ancestral population through genetic drift to a given ov va era ll ue defin l F ed by controlling effective ST population size (e) N and number of generations (t). To best reflect the assumingly Nl e ao rf ga e nchovy, we let Ne = 10 000 and varied ftrom 0 to 2078 for simulating different levels of differentiation. Aer t ft he simulation, each subpopulation was sampled a = t n 381 and divergence from genetic homogeneity was tested w χ-exac ith t test. This procedure was repeated 100 times and the proportion of significant outcomes was used to estimate statistical power for detecting pairwise genetic differentiation. Cytb sequences were aligned usiC ng l ustalX 2.0.3 with default settings, implemente G de in ne ious 5.4 , checked and trimmed manually. Sequences were reduced to haplotypes using Collaps . N eum 1.2ber of individuals (N ), number of haplotypes () a n nd haplotype ( ) a h nd nucleotide diversities ( ) wπ ere calculated in Arlequin 3.5.1.2 using the cyt b data set. Summary statistics, number of individua ), a ls ( vN erage number of alleles (A ), observed heterozygosity ( ) H avg O and expected heterozygosity ( ) w H ere calculated for each location and for each locus with Gen . N oet di ve evolutionary divergence between putative species of OWA was calculated on M us E iG ng A t 5 he Maximum Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 8 www.nature.com/scientificreports/ Composite Likelihood model. The rate variation among sites was modelled with a gamma distribution (shape parameter = 1.48). To examine the relationship between mitochondrial haplotypes, a minimum spanning network was con- 42 45 structed witAr h lequin 3.5.1.2 and visualized wit Haps h tar . Pairwise genetic differentiation was estimated 46 47 48 with G and Jost’s D value , both within and between putative species, following Pennin.gs feo t al r st_est est mtDNA and using the R package Diversity for microsatellites. 50, 51 Estimation of migration rates. We used the coalescent-based program Migrate-N to compare different biogeographic hypotheses for the past and present migration of OWA between the Atlantic and the Pacific oceans. We conducted the analyses with two sets of data, mitochondrial DNA and microsatellites, struc- tured into two groups according to geographical regions: southern Atlantic Ocean (pooled southern locations Senegal, Guinea-Bissau, Ghana, Namibia, Angola and South Africa) and Pacific Ocean (Japan, Australia and New Zealand). We tested four variations of the two-population (Atlantic-Pacific) migration model: bidirectional migration (full model, four parameters), strict Atlantic to Pacific migration (three parameters), strict Pacific to Atlantic migration (three parameters) and panmictic model that assumes the Atlantic and Pacific are part of a panmictic population (one parameter). Testing the directionality of gene flow is justified because the dominant ocean current between the ocean basins, the Agulhas flow, runs westerly from the Indian to the Atlantic Ocean 52, 53 and is thought to play a limiting role in marine dispersal in the opposite dir . Ie ni ct tio ian l values were calcu- lated using F. Mutation rates were set to be constant among loci. The Migrate-N run parameters were calibrated ST on the full model for convergence of the Markov chain Monte Carlo sampling method. The prior distributions were uniform for mutation-scaled population size parameter), t s th het at a a (r θe four times the product of the effective population size and the mutation rate, and mutation-scaled migration rates M, that is, immigration rate scaled by the mutation rate, over the ran = g e o 0.05–0.5 a f θ nd prior migration rat= e M 0–100 f or mtDNA, and θ = 10–100 and M= 0–50 for microsatellites. These settings resulted in converged posterior distributions with a clear maximum for each estimate. The Bayesian run for mtDNA consisted of one long chain with a total of 15 million states visited and 50,000 states recorded for the generation of posterior distribution histograms for each locus aer di ft scarding the first 10,000 genealogies as burn-in; for all loci, a total of 48 million states were visited and 160,000 samples were recorded. For all the analyses we used an adaptive heating scheme with four simulta- neous chains using different acceptance ratios (temperature settings were 1.0; 1.5; × 10 3.0; ); t 1h e analyses were run on a cluster computer using 4 compute nodes per run. The Bayesian run for microsatellites consisted of one long chain with a total of 6 million states were visited and 20,000 states were recorded for the gen - eration of pos terior distribution histograms for each locus aer di ft scarding the first 10,000 genealogies as burn-in; for all loci, a total of 48 million states were visited and 160,000 samples were recorded. For all the analyses we used an adaptive heating scheme with four simultaneous chains using different acceptance ratios (temperature settings were 1.0; 1.5; 3.0; 1,000,000.0); the analyses were run on a cluster computer using four compute nodes per run. Overall loci information was combined into a single estimate by Bézier approximation of the thermodynamic scores as described by Beerli & Palczewsk . W i e averaged the Bézier score over three different runs and used as input to evaluate multiple models using Bayes fac . tors Phylogenetic and dating analyses. Phylogenetic relationships within Engraulidae were based on a frag- ment of the mitochondrial c g yt en be (121 taxa, corresponding to 55 Engraulidae species; 1044 bp). A t least one representative species of Engraulidae per genus (except Papueng ) a ra nu dli n s ine randomly chosen specimens from each of the five currently recognized OWA species were included in the phylogenetic analyses, with the exception of E. capensis from which only four specimens were available (accession numbers in Supplementary Table  S2). According to Lavoué et al . , we selected the following outgroup species: Chirocentrus dC or lu ab p, ea harengus, Denticeps clupeoides, Ilisha africana, Sardina pilchardus, Sundasalanx mekongensis. The Akaike 57 58 Information Criterion (AIC i)mplemented in Modeltest 3,. s 7 elected the GT+ RI+Γ as the evolutionary model that best-fitted the data set. The inferred parameters were used in maximum likelihood (ML) and Bayesian Inference (BI) analyses. BI analyses were conducted with MrBayes 3.2.1 . Metropolis-coupled Markov chain Monte Carlo (MCMC) analyses were ran for 20,000,000 generations with sample frequency of 2000. Final trees were calculated aer a b ft urnin of 1,000 generations. PhyML 3.0 was used to estimate the ML tree and to test by non-parametric bootstrapping the robustness of the inferred trees using 1,000 pseudoreplicates. 15, 16 Previous work recovered E. encrasicolus as paraphyletic (specimens assigned to the species grouped into two clades that did not group together). To test if natural selection could interfere with phylogenetic inference, we performed all phylogenetic analyses using the above data set that only included individuals that do not show any evidence of being under selection and repeated the procedures using another data set (116 t bp) t axa; 1044 hat included individuals presenting a mutation in codon 368 of the cyt b as identified in S.il.va et al To estimate the OWA origin and date lineage-splitting events within Engraulidae we used a Bayesian relaxed molecular-clock approach as implemented in Beas t ba 2.1.3 sed on a concatenated dataset of four p - ar tial fragments of mtDNA (cyt : 1131 b bp; 16S: 800 bp) and nuclear (RAG1: 1480 bp; R AG2: 1221 bp) genes. We included sequence data of 49 Engraulidae lineages/t , li axa kely representing 14 genera (out of 16: the monospe- cific Lycothrissa and Papuengraulis genera are not represented), from which 5 are OWA (accession numbers in Supplementary TabS3 le  ). According to our ML and BI analyses, E. capensis an E. eu d rystole are conspecific with E. encrasicolus, and two clades (hereaer c ft lade A and clade B) were recovered within the latter. Hence, to perform the dating analysis we only selected a single representative of E. encrasicolus from clade A and two from clade B, both without the mutation at codon 368. We used the BirthDeath model for the tree prior that assumes that at any point in time, every lineage undergoes speciation at rate λ or goes extinc, a t a n t ra d thr te eμ e calibration points. One refers to the earliest record of Engraulidae [6–12] million years (Ma) from the Miocene - Lower Pliocene of Cyprus . The second calibration corresponds to age estimated for the divergence between Anchovia clu peoides Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 9 www.nature.com/scientificreports/ 13, 64 and A. macrolepidota [2.8–3.1] Ma due to the closure of the Panama seawa . Thy e third calibration corre- sponds to E. japonicus [2–0] Ma from Kokubu group, Japa . C n alibrations using the two fossils were modelled with a lognormal distribution, where 95% of the prior weight fell within the geological interval in which each fossil was discovered. For the Engraulidae [12–6] Ma, the parameters of the lognormal calibration prior were: 95% interval: mean in real space: 1.4, offset: 6.0 and log stdev: 1.0. For E. japonicus [2–0] Ma, the parameters of the lognormal calibration prior were: 95% interval: mean in real space: 0.465, offset: 0 and log stdev: 1.0. For the divergence between Anchovia clupeoide an s d A. macrolepidota we used a calibration according to Lessios . e t al where the closure of the isthmus of Panama occurred between 3.1–2.8 Ma. Log normal calibration was set to: 95% interval: mean in real space: 0.071, offset: 2.8 and log stdev: 1.0. MCMC analyses were run for 20,000,000 genera- tions with a sample frequency of 20,000, following a discarded burn-in of 2,000,000 steps. The convergence to the stationary distributions was confirmed by inspection of the MCMC samples using Tracer 1.6 . ML, BI and dating analyses were performed on the R2C2 research group cluster facility provided by the IT Department of the University of Algarve. Data Accessibility. Sequences of the mitochondrial c (1044 yt b bp) were deposited in GenBank (acces- sion numbers: JQ716609–JQ716731, JQ716748–JQ716756, JX683020–JX683113, KF601435–KF601478 and KJ007642–KJ007734). Sequences of the nuclear recombination-activating genes RAG1 and RAG2 and of the mitochondrial 16S rRNA used in the Bayesian dating analysis were deposited in GenBank (accession num- bers: KX824115–KX824124). The accession numbers of the sequences retrieved from GenBank correspond- ing to the remaining Engraulidae specimens used in all phylogenetic analyses are indicatS2 ed in T and a S3 b les  (SI). Genotypes of individuals obtained from microsatellites are deposited in Dryad repository (http://dx.doi. org/10.5061/dryad.r4f8v). References 1. Jansson, R. & Dynesius, M. The fate of clades in a world of recurrent climatic change: Milankovitch oscillations and evolution. Annual Review of Ecology and Systematics 33, 741–777, doi:10.1146/annurev.ecolsys.33.010802.150520 (2002). 2. Briggs, J. C. Antitropical distribution and evolution in the Indo-West Pacific Ocean. Systemati 36 c B, 237–247, iology doi:10.2307/2413064 (1987). 3. Briggs, J. C. Antitropicality and vicariance. Systematic Bi36 olog , 206–207, do y i:10.2307/2413269 (1987). 4. Henriques, R., Potts, W. M., Santos, C. V., Sauer, W. H. H. & Shaw, P. W. Population connectivity and phylogeography of a coastal fish, Atractoscion aequidens (Sciaenidae), across the Benguela current region: evidence of an ancient vicariant even 9t. , Plos One e87907, doi:10.1371/journal.pone.0087907 (2014). 5. Rocha, L. A. & Bowen, B. W. Speciation in coral-reef fishes. Journal of Fish Biol72 og, 1101–1121, do y i:10.1111/j.1095- 8649.2007.01770.x (2008). 6. Palumbi, S. R. Marine speciation on a small planet. Trends in Ecology and Evo 7l, 114–118, do ution i:10.1016/0169-5347(92)90144-Z (1992). 7. Mirams, A. G. K., Treml, E. A., Shields, J. L., Liggins, L. & Riginos, C. Vicariance and dispersal across an intermittent barrier: population genetic structure of marine animals across the Torres Strait land bridg30 e. C , 937–949, do oral Reefs i:10.1007/s00338- 011-0767-x (2011). 8. Burridge, C. P. Antitropicality of Pacific fishes: molecular insights. Environmental Biology o f65 Fi, 151–164, shes doi:10.1023/A:1020040515980 (2002). 9. Bowen, B. W., Muss, A., Rocha, L. A. & Grant, W. S. Shallow mtDNA coalescence in Atlantic pygmy angelfishes (genus Ce ) ntropyge indicates a recent invasion from the Indian Ocean. Journal of Her 97 ed , 1–12, do ity i:10.1093/jhered/esj006 (2006). 10. Whitehead, P. J., Nelson, W. S. & Wongratana, T. Clupeoid fishes of the World (Suborder Clupeoidei). Vol. 7 (FAO, 1988). 11. Borsa, P., Collet, A. & Durand, J. D. Nuclear-DNA markers confirm the presence of two anchovy species in the Mediterranean. Comptes Rendus Biologies 327, 1113–1123, doi:10.1016/j.crvi.2004.09.003 (2004). 12. Montes, I. et al. Transcriptome analysis deciphers evolutionary mechanisms underlying genetic differentiation between coastal and offshore anchovy populations in the Bay of Biscay. Marine Bio 163 logy, 205, doi:10.1007/s00227-016-2979-7 (2016). 13. Grant, W. S., Lecomte, F. & Bowen, B. W. Biogeographical contingency and the evolution of tropical anchovies (genus ) Cetengraulis from temperate anchovies (genus Engrau ). J lo isurnal of Biogeography 37, 1352–1362, doi:10.1111/j.1365-2699.2010.02291.x (2010). 14. Silva, G., Horne, J. B. & Castilho, R. Anchovies go north and west without losing diversity: post-glacial range expansions in a small pelagic fish. Journal of Biogeography 41, 1171–1182, doi:10.1111/jbi.12275 (2014). 15. Grant, W. S., Leslie, R. W. & Bowen, B. W. Molecular genetic assessment of bipolarity in the anchovy genu . J s oEn urn ga rau l ol fi F s ish Biology 67, 1242–1265, doi:10.1111/j.1095-8649.2005.00820.x (2005). 16. Grant, W. S. & Bowen, B. W. Living in a tilted world: climate change and geography limit speciation in Old World anchovies (Engraulis; Engraulidae). Biological Journal of the Linnean Society 88 , 673–689 (2006). 17. Le Moan, A., Gagnaire, P. A. & Bonhomme, F. Parallel genetic divergence among coastal–marine ecotype pairs of European anchovy explained by differential introgression aer s ft econdary contact. Molecular Ec 25 ol , 3187–3202, do ogy i:10.1111/mec.13627 (2016). 18. Jemaa, S., Bacha, M., Khalaf, G. & Amara, R. Evidence for population complexity of the European anchovy (Engraulis encrasicolus) along its distributional range. Fisheries Res 168 earc, 109–116, do h i:10.1016/j.fishres.2015.04.004 (2015). 19. Kristoer ff sen, J. B. & Magoulas, A. Population structure of anchovy Engraulis encrasicolus L. in the Mediterranean Sea inferred from multiple methods. Fisheries Resear 91 ch , 187–195, doi:10.1016/j.fishres.2007.11.024 (2008). 20. Magoulas, A., Castilho, R., Caetano, S., Marcato, S. & Patarnello, T. Mitochondrial DNA reveals a mosaic pattern of phylogeographical structure in Atlantic and Mediterranean populations o En f a gr n ac u h li o s e vy n ( crasicolus). Molecular Phylogenetics and Evolution 3, 734–746, doi:10.1016/j.ympev.2006.01.016 (2006). 21. Sanz, N., Garcia-Marín, J. L., Viñas, J., Roldán, M. & Pla, C. Spawning groups of European anchovy: population structure and management implications. ICES Journal of Marine Scien 65 ce, 1635–1644, do i:10.1093/icesjms/fsn128 (2008). 22. Viñas, J. et al. Genetic population structure of European anchovy in the Mediterranean Sea and the northeast Atlantic Ocean using sequence analysis of the mitochondrial DNA control region. ICES Journal of Marine S 7, 391–397, do cience i:10.1093/icesjms/fst132 (2013). 23. Silva, G., Lima, F. P., Martel, P. & Castilho, R. Thermal adaptation and clinal mitochondrial DNA variation of European anchovy. Proceedings of the Royal Society of London B: Biological Sciences 281, 20141093, doi:10.1098/rspb.2014.1093 (2014). 24. Liu, J. X. et al. Late Pleistocene divergence and subsequent population expansion of two closely related fish species, Japanese anchovy (Engraulis japonicus) and Australian anchovy (Engraulis austr ). al M is olecular Phylogenetics and Evolution 40, 712–723, doi:10.1016/J. Ympev.2006.04.019 (2006). Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 10 www.nature.com/scientificreports/ 25. Duchene, D., Klanten, S. O., Munday, P. L., Herler, J. & van Herwerden, L. Phylogenetic evidence for recent diversic fi ation of obligate coral-dwelling gobies compared with their host corals. Molecular Phylogenetics and Ev69 olu, 123–132, do tion i:10.1016/j. ympev.2013.04.033 (2013). 26. Cárdenas, L. et al. Origin, diversification, and historical biogeography of the genus Tra (P chu er ru cif s ormes: Carangidae). Molecular Phylogenetics and Evolution 35, 496–507, doi:10.1016/j.ympev.2005.01.011 (2005). 27. DiBattista, eJ t. al. Twisted sister species of pygmy angelfishes: discordance between taxonomy, coloration, and phylogenetics. Coral Reefs 31, 839–851, doi:10.1007/s00338-012-0907-y (2012). 28. Horne, J. B., van Herwerden, L., Choat, J. H. & Robertson, D. R. High population connectivity across the Indo-Pacific: congruent lack of phylogeographic structure in three reef fish congeners. Molecular Phylogenetics and Evo 49 lut , 629–638, do ion i:10.1016/j. ympev.2008.08.023 (2008). 29. Reece, J. S., Bowen, B. W., Smith, D. G. & Larson, A. Comparative phylogeography of four Indo-Pacific moray eel species (Muraenidae) reveals comparable ocean-wide genetic connectivity despite v fi e-fold differences in available adult habitat. Marine Ecology Progress Series 437, 269–277, doi:10.3354/meps09248 (2011). 30. Schwartzlose, R. A. et al . Worldwide large-scale fluctuations of sardine and anchovy populations. South African Journal of Marine Science 21, 289–347, doi:10.2989/025776199784125962 (1999). 31. Grant, W. S. & Bowen, B. W. Shallow population histories in deep evolutionary lineages of marine fishes: insights from sardines and anchovies and lessons for conservation. e J Th ournal of Hered89 ity, 415–426,  doi:10.1093/jhered/89.5.415 (1998). 32. Lessios, H. A. & Robertson, D. R. Speciation on a round planet: phylogeography of the goatfish genus MulloJidic ourn hal th of ys. Biogeography 40, 2373–2384, doi:10.1111/jbi.12176 (2013). 33. Sigman, D. M. & Hain, M. P. The biological productivity of the ocean. Nature Education Knowle3 d, 21 (2012). ge 34. Marlow, J. R., Lange, C. B., Wefer, G. & Rosell-Melé, A. Upwelling intensification as part of the Pliocene-Pleistocene climate transition. Science 290, 2288–2291, doi:10.1126/science.290.5500.2288 (2000). 35. Bloom, D. D. & Lovejoy, N. R. The evolutionary origins of diadromy inferred from a time-calibrated phylogeny for Clupeiformes (herring and allies). Proceedings of the Royal Society B: Biological Scienc 281 es , doi:10.1098/rspb.2013.2081 (2014). 36. Jombart, T. Adegenet: a R package for the multivariate analysis of genetic Bm ioi a nfor rker ms. atics 24, 1403–1405, doi:10.1093/ bioinformatics/btn129 (2008). 37. R Development Core Team. R: A language and environment for statistical computing. (2013). 38. Zarraonaindia, I., Pardo, M. A., Iriondo, M., Manzano, C. & Estonba, A. Microsatellite variability in European anchovy (Engraulis encrasicolus) calls for further investigation of its genetic structure and biogeography. ICES Journal of Ma66 rin , 2176–2182, e Science doi:10.1093/icesjms/fsp187 (2009). 39. Ryman, N. & Palm, S. POWSIM: a computer program for assessing statistical power when testing for genetic differentiation. Molecular Ecology Notes 6, 600–602, doi:10.1111/j.1471-8286.2006.01378.x (2006). 40. Geneious v. 5.1. (Biomatters Ltd, available at: http://www.geneious.com, Auckland, New Zealand), (2011). 41. Posada, D. Collapse 1.2: Describing haplotypes from sequence alignments http://darwin.uvigo.es/software/collapse.html (2004). 42. Excoe ffi r, L. & Lischer, H. E. L. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources10 , 564–567, doi:10.1111/j.1755-0998.2010.02847.x (2010). 43. Meirmans, P. G. & Van Tienderen, P. H. GENOTYPE and GENODIVE: two programs for the analysis of genetic diversity of asexual organisms. Molecular Ecology Notes 4, 792–794 (2004). 44. Tamura, K.e t al. MEGA5: Molecular Evolutionary Genetics Analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Molecular Biology and Evolutio 28 n , 2731–2739, doi:10.1093/Molbev/Msr121 (2011). 45. Teacher, A. G. F. & Griffiths, D. J. HapStar: automated haplotype network layout and visualization. Molecular Ecology R 11, esources 151–153, doi:10.1111/j.1755-0998.2010.02890.x (2011). 46. Hedrick, P. W. & Goodnight, C. A standardized genetic differentiation measure. Ev ol 59 u, 1633–1638, do tion i:10.1554/05-076.1 (2005). 47. Jost, L. G(ST) and its relatives do not measure differentiation. Molecular Ec 17, 4015–4026, do ology i:10.1111/J.1365- 294x.2008.03887.X (2008). 48. Pennings, P. S., Achenbach, A. & Foitzik, S. Similar evolutionary potentials in an obligate ant parasite and its two host species. Journal of Evolutionary Biology 24, 871–886, doi:10.1111/j.1420-9101.2010.02223.x (2011). 49. Keenan, K., McGinnity, P., Cross, T. F., Crozier, W. W. & Prodöhl, P. A. diveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors. Methods in Ecology and E 4, 782–788, do volution i:10.1111/2041- 210x.12067 (2013). 50. Beerli, P. & Felsenstein, J. Maximum likelihood estimation of migration rates and effective population numbers in two populations using a coalescent approach. Genet 152 ics , 763–773 (1999). 51. Beerli, P. & Felsenstein, J. Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach. Proceedings of the National Academy of Scien 9 c8 e, 4563–4568, do s i:10.1073/pnas.081068098 (2001). 52. Ivanova, E. V. The global thermohaline paleocirculation. (Springer, 2009). 53. Lutjeharms, J. R. E. In e s Th ea: the global coastal ocean Vol. 14B e s Th ea - ideas and observations on the progress in the study of the seas (eds Allan, R. Robinson & Kenneth, Brink) 783–834 (Harvard University Press, 2006). 54. Beerli, P. & Palczewski, M. Unified framework to evaluate panmixia and migration direction among multiple sampling locations. Genetics 185, 313–326, doi:10.1534/genetics.109.112532 (2010). 55. Bloomquist, E. W., Lemey, P. & Suchard, M. A. Three roads diverged? Routes to phylogeographic infer Tr en ence ds . in Ecology & Evolution 25, 626–632, 10.1016/j.tree.2010.08.010 (2010). 56. Lavoué, S., Miya, M. & Nishida, M. Mitochondrial phylogenomics of anchovies (family Engraulidae) and recurrent origins of pronounced miniaturization in the order Clupeiformes. Molecular Phylogenetics and Evo 56 lu, 480–485, do tion i:10.1016/j. ympev.2009.11.022 (2010). 57. Akaike, H. A new look at the statistical model identifications. IEEE Transactions on automatic c 19 o, 716–723 (1974). ntrol 58. Posada, D. & Crandall, K. A. Modeltest: testing the model of DNA substituition. Bioi nfor 14, 817–818 (1998). matics 59. Ronquist, Fe. t al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology 61, 539–542, doi:10.1093/sysbio/sys029 (2012). 60. Guindon, S. et al. New Algorithms and Methods to Estimate Maximum-Likelihood Phylogenies: Assessing the Performance of PhyML 3.0. Systematic Biology 59, 307–321, doi:10.1093/sysbio/syq010 (2010). 61. Bouckaert, R e. t al. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Com 10 pu, e1003537, do t Biol i:10.1371/ journal.pcbi.1003537 (2014). 62. Gernhard, T. The conditioned reconstructed process. J Theor Biol 253, 769–778, doi:10.1016/j.jtbi.2008.04.005 (2008). 63. Grande, L. & Nelson, G. Interrelationships of fossil and recent anchovies (Teleostei: Engrauloidea) and description of a new species from the Miocene of Cyprus. American Museum Novitat2826 es  , 1–16 (1985). 64. Lessios, H. A., Kessing, B. D., Robertson, D. R. & Paulay, G. Phylogeography of the pantropical sea urchin in r Euce ida lart is ion to land barriers and ocean currents. Evol u 53 ti, 806–817 (1999). on 65. Yabumoto, Y. Pleistocene clupeid and engraulidid fishes from Kokubu group in Kagoshima Prefecture, Japan. Bulletin of Kitakyushu Museum Natural History 8, 55–74 (1988). 66. Rambaut, A., Suchard, M. A., Xie, D. & Drummond, A. J. Tracer v1.6. Available f h ro tm tp://beast.bio.ed.ac.uk/Tracer (2014). Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 11 www.nature.com/scientificreports/ Acknowledgements We would like to thank to A. Gonçalves and M. Ruano (Portugal), K. Rowling (Australia) and R. Tesker (New Zealand), for collecting and shipping samples (see also refs 14 and 23 acknowledgments). We thank N. Coelho, C. Patrão and M. Valente for technical help. GS doctoral grant (SFRH/BD/36600/2007) was supported by the Portuguese Foundation for Science & Technology (FCT) and presently he is supported through the FCT strategic project UID/MAR/04292/2013 granted to MARE. RLC was supported by a FCT post-doctoral fellowship (SFRH/ BPD/65830/2009). The work was funded by FCT strategic plan UID/Multi/04326/2013 granted to CCMAR. Author Contributions G.S. and R.C. planned the study. G.S. organized the sampling. G.S. and A.R. carried out the molecular work. G.S., R.L.C. and R.C. performed data analyses. G.S. and R.C. made the figures. G.S., R.L.C. and R.C. wrote the manuscript. All authors read and approved the final version of the manuscript. Additional Information Supplementary information accompanies this paper at doi:10.1038/s41598-017-02945-0 Competing Interests: The authors declare that they have no competing interests. Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Cre- ative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not per- mitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. © The Author(s) 2017 Scientific Repo R ts | 7: 2893 | DOI:10.1038/s41598-017-02945-0 12

Journal

Scientific ReportsSpringer Journals

Published: Jun 6, 2017

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off