Meta-analysis of northeast Atlantic marine taxa shows contrasting phylogeographic patterns following post-LGM expansions

Literature search

To compare the phylogeography of benthic and demersal organisms across the northeast Atlantic, we undertook a meta-analysis of molecular phylogeographic studies. A literature search was conducted by TLJ and JRS using Web of Science (Thomson Reuters) in February 2015. Search terms were chosen to maximise the inclusion of articles reporting on population structure and phylogeography from the northeast Atlantic. The following sets of Boolean search terms were submitted to the Advanced Search Tool: (1) gene flow OR population structure OR genetic diversity OR phylogeograph*; (2) marine OR intertidal OR subtidal OR estuar*; and (3) Atlantic. Titles and abstracts were screened by TLJ and JRS to identify suitable articles within the scope of this study and only articles that matched the following criteria were retained: (a) organisms were fully marine or estuarine throughout their life history (diadromous species were excluded); (b) studies of temporal changes, hybridisation or introgression from closely related species were omitted; (c) the study included at least three sampling sites from within the northeast Atlantic (Fig. 1—sites outside of this area were not considered); (d) datasets contained a minimum of five individuals per site and a total sample size of at least 50; and (e) the study included latitude and longitude of the sampling sites or a detailed description or map which provided sufficient detail to determine the geographical location of sample origins. The studies were reviewed independently by TLJ and JRS and there were no disputes regarding inclusion or rejection that needed adjudication. Given the proven utility of mtDNA in comparative phylogeography (e.g., Patarnello, Volckaert & Castilho, 2007; Ni et al., 2014) and the availability of these data in the public domain, a meta-analysis was conducted using published mtDNA gene sequences.

Data reanalysis

A standardised methodology was implemented to ensure that all mtDNA datasets were reanalysed in a consistent and directly comparable manner. Data analyses in the original studies were far from consistent, particularly with respect to the analysis of haplotype networks and historical demography. The majority of studies reported information about population structure, however, in several instances the studies included additional samples outside of the northeast Atlantic in their analysis. Therefore, standardised tests of population structure were undertaken de novo for each species. Sites that were genetically homogeneous (as described by the original authors) and which were spatially close or situated in the same geographical region were combined in some datasets. This ensured that phylogeography within and across seas was examined in this meta-analysis. Population differentiation was examined using global values of Jost’s D (Jost, 2008) and F_ST (Weir & Cockerham, 1984) using the fastDivPart function from the R package diveRsity (Keenan et al., 2013; R Core Team, 2016) and significance was assessed using 10,000 permutation replicates.

To examine the genealogical relationships within species, haplotype networks were constructed using the haploNet function from the R package pegas (Paradis, 2010). Tajima’s D (Tajima, 1989), Fu’s F_S (Fu, 1997) and Ramos-Onsins’ R ₂ (Ramos-Onsins & Rozas, 2002) neutrality tests were performed in DnaSP v5.10 (Librado & Rozas, 2009) to determine whether each species carried a signal that deviated from neutrality (significance was assessed using 10,000 bootstrap replicates). Mismatch analyses (frequency of pairwise nucleotide-site differences between sequences) were carried out using the population growth-decline model in DnaSP to further examine the demographic history, and Harpending’s raggedness index (r) (Harpending, 1994) was used to evaluate the fit of the observed distribution to the growth-decline model (10,000 bootstrap replicates). A non-significant index suggests that the observed data have a relatively good fit to the growth-decline model. In contrast, a significant index is indicative of a stable population which is typically thought to show a ‘ragged’, multi-modal mismatch (Harpending, 1994).

The equation t = τ/(2µk) was used to estimate the timing of a population expansion (t), where τ is the date of the expansion measured in units of mutational time (Tau –estimated using DnaSP), µis the mutation rate per site per year and k is the sequence length. In addition, Bayesian Skyline Plots (BSPs) were run using BEAST2 v2.5.0 (Drummond et al., 2005; Bouckaert et al., 2014). BEAST2 uses a Markov chain Monte Carlo (MCMC) sampling procedure to estimate effective population size (N_e) through time based on the temporal distribution of coalescences in gene genealogies. For each dataset, the substitution model was selected using bModelTest (Barido-Sottani et al., 2018), which uses reversible jump MCMC that allows the Markov chain to jump between states representing different possible substitution models. A strict clock and a coalescent Bayesian Skyline prior was implemented. Each run consisted of 100 million steps with a burn-in of one million and parameters were sampled every 10,000 steps. Chain convergence and BSPs were analysed with Tracer v1.7.1 (Rambaut et al., 2018).

Recent studies have shown that the use of mutation rates derived from ancient calibration dates or from phylogenetic analyses may not be appropriate for studies at the population level (Ho et al., 2008; Ho et al., 2011). In this study, therefore, mutation rates were chosen based on the most recent calibration date available for the closest taxonomic relative (Table S1). In published studies where a mutation rate was not specified, the genetic distance provided by the study was divided by the date of the calibration event (in Myr) to obtain a % mutation rate per Myr. For cases where only calibration dates older than 5 Myr were available for the species and gene of interest, a three-fold correction in mutation rate was applied to the original rate to control for the potential time-dependency of molecular rates. This adjustment was implemented because rates have been found to vary by three to six-fold for several marine species when calibration dates younger than 5 Myr vs. older dates have been tested (Crandall et al., 2012; Laakkonen, Strelkov & Väinölä, 2015). A range of mutation rates based on the rates reported by previous studies were used to calculate a minimum, maximum and average time estimate since a population expansion.

Literature search

The initial search using Boolean terms identified 1,120 articles, which was reduced to 56 articles after the titles and abstracts were examined and the search criteria were applied (Fig.  S1). The final database for the meta-analysis consisted of mtDNA gene sequence data from 21 studies (Table 1); some studies from the previous step were not included due to the use of RFLPs in mtDNA or because some mtDNA datasets were not publicly available. The final database spanned several taxonomic groups, with fishes, molluscs and crustaceans accounting for the majority of species (81%). The most common mitochondrial gene across all studies was cytochrome oxidase I (COI), followed by cytochrome b (Cyt b), the control region (CR) and the intergenic spacer region (IGS). COI was the most commonly used gene for invertebrate studies, IGS for macroalgae, and studies of fish used either the CR or the Cyt b gene.

Genetic structure

Sixteen species showed significant global Jost’s D and F_ST values, indicative of population differentiation (Table 2), while the remaining five species showed little evidence of population differentiation. Across the 21 datasets, four different types of haplotype network were putatively identified based on the structure of the networks (Fig. 2) (all haplotype networks are presented in Fig. S2):

(i) A ‘Star’ network (Fig. 2A), in which a single, widespread haplotype is typically positioned at the centre of the network and is thought to be the ancestral haplotype. Additional haplotypes are linked to this dominant haplotype by a single (or a few) mutational step(s), suggesting these haplotypes are the product of recent mutation events. Eight species showed this type of relationship (Celleporella hyalina, Conger conger, Nassarius nitidus, Nassarius reticulatus, Palinurus elephas, Pelvetia canaliculata, Pomatoschistus microps and Raja clavata). In one case, the dominant haplotype had far fewer connections than a low-frequency haplotype in the network, making it difficult to distinguish the centre of the network with confidence (Pomatoschistus microps);

(ii) A ‘Complex star’ network (Fig. 2B), in which there are multiple high-frequency haplotypes and connections. Six species showed this type of relationship (Carcinus maenus, Cerastoderma edule, Maja brachydactyla, Pomatoschistus minutus, Solea solea, Symphodus melops);

(iii) A ‘Reciprocally monophyletic’ network (Fig. 2C), in which more than one lineage is apparent and each lineage is linked by a long branch associated with numerous mutations. Four species showed this type of relationship (Macoma balthica, Modiolus modiolus, Owenia fusiformis and Pectinaria koreni);

(iv) A ‘Complex mutational’ network (Fig. 2D), in which some branches were separated by a very large number of mutations, while other branches had contrarily one or two mutations. Three species showed this type of relationship (Dicentrarchus labrax, Labrus bergylta and Neomysis integer). In most cases, a dominant haplotype was present and was presumed to be the ancestral form. However, Neomysis integer presented an unusual network in which a distinct ancestral haplotype was not apparent and the centre of the haplotype network was not readily distinguishable.

Historical demography

Historical demography was inferred for each species based on the observed mismatch distribution, neutrality tests and the raggedness index (Table 2). Four main types of mismatch distributions were observed: unimodal, skewed unimodal, multimodal and bimodal (Fig. 3) (all mismatch distributions are presented in Fig. S3). Unimodal is associated with a sudden population expansion (e.g., Maja brachydactyla; Fig. 3A), and skewed unimodal is generally associated with a recent expansion or bottleneck (e.g., Nassarius reticulatus; Fig. 3B). Multimodal (e.g., Labrus bergylta; Fig. 3C) and bimodal (e.g., Macoma balthica; Fig. 3D) are usually associated with constant population size. However, previous research has suggested that bimodal peaks may indicate the presence of two distinct lineages (e.g., Alvarado-Bremer et al., 2005), which would potentially violate the assumptions of coalescent theory if analysed as one ‘genetic’ population. In this case, the first peak would represent intra-clade pairwise differences, whereas the second peak would likely represent more ancient inter-clade pairwise differences (Fig. 3D). For each instance of bimodality, the haplotype network was inspected for evidence of two or more lineages. The networks indicated that more than one distinct lineage was evident for all bimodal mismatches (Macoma balthica, Modiolus modiolus, Owenia fusiformis and Pectinaria koreni) and, therefore, mismatch analysis and neutrality tests were carried out on each lineage separately. These analyses were not conducted for lineage 2 of Modiolus modiolus due to the small number of individuals (N = 3) comprising this lineage.

Neutrality statistics for testing the drift–mutation equilibrium (Tajima’s D, F_S and R ₂) were found to be contrasting between species (Table 2). These tests tended to be significant for species that showed a star-shaped network and for which the mismatch graph was unimodal or skewed unimodal. This supported evidence that a signal of rapid population expansion was detected; however, a selective sweep can also produce the same genetic signal. Harpending’s r suggested that two datasets departed from a model of demographic expansion (Table 2), but inspection of the mismatch graphs and neutrality tests indicated there was strong evidence to support a rapid population expansion (or selective sweep) in both datasets. No signatures of rapid population expansion were detected in five species (Celleporella hyalina, Macoma balthica, Nassarius nitidus, Neomysis integer and Raja clavata), suggesting a stable constant population size.

For the remaining 19 datasets (16 species, 19 including lineages), a historic population expansion was assumed and the timing of the expansion was estimated (Fig. 4). All expansions were found to take place during the Pleistocene or the Holocene epoch. Estimated timings for 17 datasets were after or overlapped the earliest estimate for the LGM (∼26.5 Ka). Expansion estimates for one fish (Labrus bergylta) and one lineage of the polychaete Owenia fusiformis pre-dated the LGM but were still positioned during the last glacial period. Bayesian Skyline Plots (Fig. 5) were generally consistent with the results from the mismatch analyses. Among the 17 datasets for which from the mismatch analyses expansion times were estimated to have occurred after the LGM, a rise in N_e post-LGM was apparent in 15 of these datasets, but the strength of the increase varied across datasets. In comparison to the mismatch analysis, the BSP for L. bergylta (Fig. 5F) and O. fusiformis lineage 2 (Fig. 5P) indicated a population expansion after the earliest estimate for the LGM as opposed to pre-dating the LGM. In addition, although the mismatch analyses inferred a post-LGM expansion for M. modiolus lineage 1 (Fig. 5M) and O. fusiformis lineage 3 (Fig. 5Q), BSPs generally suggested N_e was constant after the LGM.

Demographic history

Demographic history was variable across species in the northeast Atlantic, as evidenced by both the diverse structuring of the haplotype networks and the observed mismatch distributions within species. The presence of one or more lineages and the complexity of mutational patterns in several networks suggested some species have undergone pronounced changes in their demography and genealogy. Connections with large mutation steps separating some haplotypes are indicative of deep phylogenetic splits in the genealogies and suggests the persistence of old populations in these species. Accumulating new mutations is a relatively slow process and, therefore, sufficient time since coalescence must have elapsed to facilitate these large sequence divergences (Avise, 2009).

In the northeast Atlantic, the LGM has often been viewed as a possible explanation for discrepancies in genealogies and for rapid population expansions via recolonisation as glaciers started to retreat from their maximum positions (Hewitt, 2004). In this study, we detected rapid expansions in many different taxa, of which the majority were estimated to occur after the LGM. This supports evidence for post-LGM expansions, possibly from periglacial refugia (Maggs et al., 2008) or via recolonisation of areas previously affected by the Northern Hemisphere ice sheets. These results are in contrast to the northeast Pacific where regional persistence during the LGM appeared to be common in rocky-shore organisms (Marko et al., 2010). The conclusions of several previous studies reanalysed in this meta-analysis also detected rapid expansions (e.g., Jolly et al., 2006; Sotelo et al., 2008; Larmuseau et al., 2009); however, the authors of these studies estimated the dates of these expansions to have occurred pre-LGM. This discrepancy could be due to the differences in mutation rates, whereby the original authors typically used rates derived from ancient calibrations, while in this study we attempted to use more recent calibration dates to correct for the potential time-dependency of molecular rates (Ho et al., 2011).

Of course, we acknowledge that the signal of deviation from neutrality we detected may, in some cases, be the result of a selective sweep and not a rapid expansion. This signal could be distinguished by incorporating multi-locus data; nevertheless, given that a variety of species in this study showed similar genealogical patterns consistent with demographic expansion, it seems likely that most of them did indeed experience demographic changes associated with the end of the LGM, rather than selective sweeps. Moreover, distinctive haplotypes were found in several species networks (Pelvetia canaliculata, Pomatoschistus minutus, Owenia fusiformis and Pectinaria koreni) to the south of where the Eurasian ice sheet is proposed to have extended during the LGM (Fig. 1). This finding suggests populations of these species may have survived in southern glacial refugia; though, as pointed out by some of the original authors, deep sequence divergences in some species (e.g., O. fusiformis and P. korena) and the lack of a species-specific molecular clock calibration makes inferences about refugia challenging (Jolly et al., 2005; Jolly et al., 2006).

It is difficult to suggest an explanation for the two expansions estimated to have pre-dated the LGM (using mismatch analysis), but which fall within the last glacial period. This pattern of pre-LGM expansion has also been reported in a number of previous studies for a variety of marine taxa (e.g., Hoarau et al., 2007; Marko et al., 2010; Ni et al., 2014; Almada et al., 2017). One potential explanation for this pattern is that sea level during the last glacial cycle did not decrease uniformly towards the level observed at the LGM, but oscillated rapidly over a period of 60 Ka to 30 Ka (see Fig. 3A in Lambeck, Esat & Potter, 2002). Therefore, it may be possible that we are detecting the signature of a population expansion during one of these sudden increases in sea level during the last glacial period. Alternatively, as the BSP analysis inferred a post-LGM expansion for these two datasets, this could be a limitation associated with the mismatch analysis approach, which does not consider genealogy, and may, therefore, produce a less precise estimation. In addition, the sample of genetic diversity for this species may not be representative (Karl et al., 2012) or the genetic signal we detected may have been the result of a selective sweep and not a rapid expansion.

The use of single marker mtDNA genealogies and coalescence theory can introduce challenges associated with the interpretation of data and these limitations should be acknowledged (Karl et al., 2012). For example, the populations under study may have experienced multiple episodes of growth and decline; however, only the most recent expansion event can be detected using coalescence analysis and, in some cases, these events may not be sufficiently severe to be detected (Karl et al., 2012). In addition, coalescent histories can differ amongst loci because they can experience mutation and drift independently. Therefore, analysis of a single gene only gives insight into the coalescent history of that locus, which may not always be representative of population history. Analysis of multiple loci and genomics would help to alleviate these concerns, and would likely provide enhanced resolution for exploring the phylogeography of northeast Atlantic marine fauna.

Although population expansions were detected in a number of species in this study and also in the wider literature, populations of other marine species, including five from this study, have been found to remain stable throughout the LGM. As previously reported, not all coastal marine taxa appear prone to demographic changes during or after ice ages (Janko et al., 2007; Marko et al., 2010; Olsen et al., 2010). It is also important to acknowledge that earlier events in the Pleistocene and more ancient events that pre-date the Pleistocene may have helped shape the contemporary patterns of genealogical structure observed in this study.