Abstract

Larval development has strong impacts on dispersal potential and gene flow among populations of marine invertebrates. However, Pleistocene glaciations have also played an important role in shaping population structure in benthic taxa in the Northern Hemisphere, even those with planktotrophic larvae. Each glacial advance tended to fragment species distributions, often separating populations for long periods and setting the stage for their differentiation. This study examines patterns of sequence divergence of the mitochondrial cytochrome c oxidase subunit I gene in North American populations of two bivalve species complexes, Hiatella arctica s. l. and Macoma balthica s. l., with complementary data from the nuclear internal transcribed spacer-2 (ITS2) gene for the latter. Our results confirm the presence of two known species from the M. balthica complex in Canada, but also provide evidence for a third clade in Atlantic Canada. Our study confirms that the H. arctica complex in Canada contains at least four species, with support for a novel clade (Hiatella N) in the northeastern Pacific. Our results extend the range of a previously identified Hiatella clade (K) to include the northwestern Atlantic. Both M. balthica s. l. and H. arctica s. l. have broad Holarctic distributions and planktotrophic larvae, but this work reveals differences in phylogeographic structure and genetic diversity.

INTRODUCTION

Life-history attributes, vicariance events and environmental tolerances all play a role in determining where species occur (Reid, 1990; Hewitt, 2000). However, contemporary patterns of population structure in the Northern Hemisphere can only be understood by considering dispersal capacity and past population subdivision in response to glacial events and how these factors impacted population isolation and differentiation (e.g. Wares & Cunningham, 2001). Patterns of population structure in marine invertebrates with a sessile adult stage are interesting because dispersal potential in these taxa is largely dictated by the presence or absence of a planktonic larval stage (Meehan, 1985; Lee & Boulding, 2009). Some larvae spend weeks in the plankton and passive dispersal by oceanic currents can cause enough gene flow to ensure near panmixis across the species’ range (Jablonski, 1986; Arndt & Smith, 1998; Lee & Boulding, 2009). Although planktotrophic species may demonstrate less regional genetic divergence than those lacking a planktonic stage, substantial population structure occurs in some species with planktonic larvae. For instance, Macoma balthica shows clear genetic divergence between populations in the northwestern and northeastern Atlantic (Väinölä, 2003; Nikula, Strelkov & Väinölä, 2007).

Glacial cycles had a dramatic impact on Canadian marine environments. During glacial maxima, sea levels declined by up to 170 m and the Arctic Ocean was covered by ice. The northwestern Atlantic was heavily impacted by the last glacial maximum (LGM), but the Cordilleran ice sheet on the west coast produced lesser impacts (Warner, Mathewes & Clague, 1982; Bernatchez & Wilson, 1998; Rohling et al., 1998; Hewitt, 2000; Mandryk et al., 2001; Marko, 2004). In fact, Marko et al. (2010) suggested that many rocky-shore species in the northeastern Pacific were able to persist during the LGM. In addition, the recurrent opening and closing of the Bering Strait linked to glacial cycles promoted the intermittent isolation and exchange of species between the Pacific and Atlantic Oceans (Briggs, 1970; Vermeij, 1991; Taylor & Dodson, 1994; Dodson et al., 2007). During each glacial advance, species retreated to refugia and then expanded their range during the subsequent interglacial (Hewitt, 2000). For example, studies of mtDNA diversity in two species with planktonic larvae—the capelin, Mallotus villosus, and the green sea urchin, Strongylocentrotus droebachiensis—demonstrated their persistence in glacial refugia (Addison & Hart, 2005; Dodson et al., 2007). Because each glacial advance tended to fragment species distributions in a consistent way, populations were repeatedly separated for prolonged periods, setting the stage for their differentiation (Hewitt, 1996; Maggs et al., 2008; Dapporto, 2009). As species expanded their range during interglacials, divergent lineages re-established contact at sites along Canada's coasts (Harper & Hart, 2007). Species with similar life history attributes, but differing ecological preferences, have undoubtedly responded differently to glacial cycles (Bernatchez & Wilson, 1998).

Holarctic species are ideal candidates for determining how glacial cycles have shaped contemporary patterns of genetic variation. This study targets two wide-ranging bivalves, Hiatella arctica and Macoma balthica, both with a planktonic larval phase. There is evidence that each of these taxa is best viewed as a species complex. Coan (1971) recognized M. petalum as a sibling species of M. balthica from the northwestern Atlantic, and genetic analyses later revealed divergent lineages of M. balthica in the eastern and western North Atlantic (Meehan, 1985; Meehan, Carlton & Wenne, 1989). Based on DNA sequence analysis, Väinölä (2003) confirmed the status of M. petalum as a separate species and uncovered four lineages in the M. balthica complex. Nikula et al. (2007) corroborated these results using additional molecular data, and further discovered that M. b. balthica likely consists of two lineages (C and D). Both studies highlighted the role of Pleistocene glacial cycles in shaping genetic structure in this complex. Lastly, a recent study by Layton, Martel & Hebert (2014) explored DNA barcode variation in marine molluscs and highlighted three M. balthica clades in Canada.

Hiatella arctica also has a complex taxonomic history with two species names (H. pholadis and H. striata) synonymized with it (Marshall & Gofas, 2015). Taxonomic decisions are complicated because adults of Hiatella display considerable phenotypic plasticity, shell shape varying with the substrate and microhabitat type (Alison & Marincovich, 1982). Recent work examining variation of the mitochondrial cytochrome c oxidase subunit I (COI) gene has highlighted the occurrence of multiple Hiatella clades across North America (Layton et al., 2014) with at least 11 taxa in the Northern Hemisphere as a whole (Laakkonen, Strelkov & Väinölä, 2015). As in the case of M. balthica, Layton et al. (2014) highlighted cryptic diversity in the H. arctica complex, but did not examine the patterning of genetic diversity in these groups. The present study addresses this gap by using DNA barcodes (COI sequences) to compare patterns of phylogeographic structure in Canadian populations of H. arctica s. l. and M. balthica s. l., with a secondary goal of advancing understanding of the taxonomic status of their component lineages. We expand on prior work by both incorporating additional localities in Canada and sequences of the nuclear internal transcribed spacer-2 (ITS2) gene for the M. balthica complex.

METHODS

Sample collection, vouchers and data deposition

From 4 to 60 specimens of the Macoma balthica and Hiatella arctica complexes were collected at 33 sites in Alaska, British Columbia, Labrador, Manitoba, New Brunswick, Newfoundland, Nova Scotia and Prince Edward Island between 2007 and 2014 (Fig. 1). Specimen details, sequences and trace files are available from the Dataset dx.doi.org/10.5883/DS-COIPOP on BOLD (Ratnasingham & Hebert, 2007) and specimens are held at the Biodiversity Institute of Ontario. Sequences have also been deposited in GenBank (accession numbers: KP977574-KP977969, KP977970-KP978014). Specimens were obtained from rock crevices, algal mats and holdfasts in the intertidal zone and also from subtidal habitats using otter trawls, dredges and SCUBA diving. A lateral incision was made along each shell to ensure proper preservation of internal tissues when transferred to 90% ethanol. Specimens were then stored at −20 °C until ready for tissue sampling. All specimens were examined at the Canadian Museum of Nature to ensure that they matched the morphological description of H. arctica and M. balthica. Subspecies designations for M. balthica follow Väinölä (2003) and Nikula et al. (2007), while Hiatella lineage assignments follow Laakkonen et al. (2015). For simplicity, we employ species names in the broad sense throughout the rest of this section.

Figure 1.

Collection sites for Hiatella arctica s. l. and Macoma balthica s. l. Sample sizes for each species at each site are shown in pie diagrams.

Figure 1.

Collection sites for Hiatella arctica s. l. and Macoma balthica s. l. Sample sizes for each species at each site are shown in pie diagrams.

DNA amplification and sequencing

Tissue from the adductor muscle was removed from each specimen and placed in cetyltrimethylammonium bromide (CTAB) lysis buffer solution with proteinase K. The samples were then incubated for 12 h at 56 °C before DNA extraction was carried out using a manual glass-fibre plate method (Ivanova, Fazekas & Hebert, 2008). Following incubation, the DNA was eluted with 40 µl of ddH2O. After re-suspension, 2 µl of each DNA extract was placed into a well in a separate plate with 18 µl ddH2O to ensure the dilution of salts or mucopolysaccharides that could inhibit PCR. Both COI and ITS2 were amplified from M. balthica, while only COI was amplified from H. arctica. Species-specific primer sets were used for COI: HiaF1/HiaR1: AAGTTGTAATCATCGAGATATTGG and TAGACTTCTGGGTGCCCGAAAAACCA for H. arctica and MMacF1/LepR1: CTTTTATTAGCTGCACCTGATAT and TAAACTTCTGGATGTCCAAAAAATCA for M. balthica (S. Prosser, unpublished). Universal ITS2 animal primers were also used for M. balthica: CAS5p8sFc: TGAACATCGACATTTYGAACGCACAT and CAS28sB1d: TTCTTTTCCTCCSCTTAYTRATATGCTTAA (Ji, Zhang & He, 2003). Each well was filled with 2 µl of diluted DNA and the following reagents to generate a 12.5 µl PCR reaction mix: 6.25 µl 10% trehalose, 2 µl ddH2O, 1.25 µl 10× PCR buffer, 0.625 µl MgCl2 (50 mM), 0.125 µl of each forward and reverse primer (10 µM), 0.0625 µl dNTP (10 mM) and 0.06 µl Platinum Taq polymerase. The thermocycling regime consisted of one cycle of 1 min at 94 °C, 40 cycles of 40 s at 94 °C, 40 s at 52 °C and 1 min at 72 °C, and finally 5 min at 72 °C. An E-GelH 96 (Invitrogen) was used to check 3 µl of each PCR product and reactions that generated an amplicon were bidirectionally sequenced using PCR primers and BigDye v. 3.1 on an ABI 3730xl DNA Analyzer (Applied Biosystems). Sequences were manually edited using CodonCode (CodonCode Corporation) and an amino acid alignment for COI was generated by eye in MEGA v. 6 (Tamura et al., 2013). All ITS2 sequences were aligned with Kalign2 (Lassmann, Frings & Sonnhammer, 2009).

Data analysis

Neighbor-joining (NJ) trees for COI (both species) and ITS2 (M. balthica) were constructed in MEGA v. 6 using a Tamura-3-parameter substitution model and 1,000 bootstrap replicates (Tamura, 1992; Tamura et al., 2013). A model test in MEGA v. 6 returned the lowest Bayesian Information Criterion value for the Tamura-3-parameter substitution model for both the Macoma and Hiatella datasets. The M. balthica and H. arctica NJ trees were rooted with outgroups (Macoma inquinata and Panopea generosa, respectively) and each COI cluster showing more than 2% divergence was treated as a distinct clade. These NJ trees also included 13 and 11 COI sequences from previously identified lineages from Väinölä (2003) and Laakkonen et al. (2015), respectively. The COI sequences generated by Laakkonen et al. (2015) only partially overlap (105 bp) with the barcode region. Clustering patterns in each NJ tree were compared with those in the corresponding median-joining haplotype network that was constructed in Network v. 4.6.1 (fluxus-engineering.com; Bandelt, Forster & Röhl, 1999). Haplotype networks were recreated in TCS v. 1.21 to identify ancestral haplotypes (Clement, Posada & Crandall, 2000).

Arlequin v. 3.5 (Excoffier & Lischer, 2010) was employed to examine patterns of genetic structure in each species using COI sequences. Haplotype and nucleotide diversity for each population was calculated with pairwise differences, and Tajima's D test of neutrality with 1,000 simulated samples was used to test for evidence of nonneutral evolution and population size fluctuation (Nei, 1987; Tajima, 1989). The proportion of unique haplotypes was also quantified to ascertain the geographic location(s) with the highest genetic exclusivity. An analysis of molecular variance (AMOVA) was conducted with pairwise differences and significance was tested with 1,000 permutations (Kimura, 1980). The AMOVA results were used to determine whether the majority of genetic variation in each species existed within or between populations, as a measure of regional differentiation (Excoffier & Lischer, 2010). Fixation indices (FST) were estimated with 1,000 permutations to determine the partitioning of variance (Weir & Cockerham, 1984). Slatkin's linearized FST (Slatkin, 1993) was plotted against geographic distance to test for isolation-by-distance (IBD) and ultimately to infer whether genetic variation among populations reflects long-term historical divergence or geographic distance (Slatkin, 1993; Kyle & Boulding, 2000; Marko, 2004; Keever et al., 2009). In order to test the significance of IBD, a Mantel test with 1,000 permutations was conducted in Arlequin v. 3.5 (P-values <0.05 were treated as significant) after geographical distance between sites was calculated using Google Earth.

RESULTS

Sequence and haplotype diversity

The 197 COI and 45 ITS2 sequences from the Macoma balthica complex ranged in length from 377 to 655 bp (mean = 429 bp) and 210 to 349 bp (mean = 326 bp), respectively. The 199 COI sequences of Hiatella arctica ranged in length from 550 to 655 bp (mean = 648 bp). All sequences of variable length were included in the construction of NJ trees and haplotype networks, but M. balthica s. l. and H. arctica s. l. sequences were trimmed to 378 and 588 bp for the Arlequin analysis. All sequences contained less than 1% ambiguous bases and lacked stop codons and double peaks. The considerable length variation in sequences for the M. balthica complex reflected the need to employ internal primers to recover sequences from some specimens with degraded DNA. Maximum and mean intraspecific divergences in the H. arctica complex were 23.5 and 6.9%, respectively, while maximum and mean intraspecific divergences in the M. balthica complex were 15.4 and 8.8%, respectively (Fig. 2).

Figure 2.

Intraspecific COI sequence divergence (% K2P) for Hiatella arctica s. l. (A) and Macoma balthica s. l. (B).

Figure 2.

Intraspecific COI sequence divergence (% K2P) for Hiatella arctica s. l. (A) and Macoma balthica s. l. (B).

A NJ tree indicated that four lineages were present within H. arctica s. l. (Fig. 3). Three corresponded to haplotypes I, K and L of Laakkonen et al. (2015), while the fourth from British Columbia and southeastern Alaska was new and was designated as Hiatella lineage N. Hiatella K was previously reported from multiple locations in Europe (Laakkonen et al., 2015), but this work revealed it to be present in the northwestern Atlantic also. None of the sequences in this study clustered with the other eight lineages found by Laakkonen et al. (2015) (i.e. European lineages C, D, E and F, northeastern Pacific lineages G and H, Bering Sea lineage M, and northwestern Pacific lineage J) based on a 105-bp overlapping segment of the barcode region of COI. Haplotypes A and B of Laakkonen et al. (2015) are from the Southern Hemisphere (Chile and New Zealand, respectively) and were excluded from the NJ tree because they did not cluster with any Northern Hemisphere lineages.

Figure 3.

Neighbour-joining tree (COI) based on a Tamura-3-parameter substitution model with gamma distribution for Hiatella arctica s. l. NJ tree includes sequences from this study and from Laakkonen et al. (2015). Lineage names follow Laakkonen et al. (2015) and clades discovered in this study are indicated with asterisks.

Figure 3.

Neighbour-joining tree (COI) based on a Tamura-3-parameter substitution model with gamma distribution for Hiatella arctica s. l. NJ tree includes sequences from this study and from Laakkonen et al. (2015). Lineage names follow Laakkonen et al. (2015) and clades discovered in this study are indicated with asterisks.

The COI NJ tree revealed three clades in the M. balthica complex from Canada, representing M. b. balthica, M. balthica NFLD and M. petalum (Fig. 4). No Canadian sequences clustered with M. b. rubra (Fig. 4). The ITS2 NJ tree provided support for two clades; one containing M. b. balthica and M. balthica NFLD lineages and the second a clade of M. petalum. One specimen (POPMA047-14) collected from New Brunswick was assigned to M. petalum in the COI tree, but to M. b. balthica/M. balthica NFLD in the ITS2 tree (Fig. 5).

Figure 4.

Neighbour-joining tree (COI) based on a Tamura-3-parameter substitution model for Macoma balthica s. l. NJ tree includes sequences from this study and European sequences obtained from GenBank. Subspecies names suggested by Väinölä (2003) and Nikula et al. (2007) are provided and clades discovered in this work are indicated by asterisks.

Figure 4.

Neighbour-joining tree (COI) based on a Tamura-3-parameter substitution model for Macoma balthica s. l. NJ tree includes sequences from this study and European sequences obtained from GenBank. Subspecies names suggested by Väinölä (2003) and Nikula et al. (2007) are provided and clades discovered in this work are indicated by asterisks.

Figure 5.

Neighbour-joining tree (ITS2) based on a Tamura-3-parameter substitution model for Macoma balthica s. l. Subspecies names suggested by Väinölä (2003) and Nikula et al. (2007) are provided and bootstrap probabilities are shown. The specimen highlighted with an asterisk clusters with M. b. balthica/M. balthica NFLD in the ITS2 NJ tree and with M. petalum in the COI NJ tree.

Figure 5.

Neighbour-joining tree (ITS2) based on a Tamura-3-parameter substitution model for Macoma balthica s. l. Subspecies names suggested by Väinölä (2003) and Nikula et al. (2007) are provided and bootstrap probabilities are shown. The specimen highlighted with an asterisk clusters with M. b. balthica/M. balthica NFLD in the ITS2 NJ tree and with M. petalum in the COI NJ tree.

Most sequences (163 of 199) and haplotypes (66 of 83) of the Hiatella complex fell into a dominant lineage (L) that was found at all sites. All 66 haplotypes of this dominant clade showed low divergence, differing from one another by just 1 to 3 mutational steps (Fig. 6A). The remaining 17 haplotypes included representatives of three clades, two from the northeastern Pacific (Hiatella I and Hiatella N) and a third from New Brunswick (Hiatella K) (Fig. 6A, B). TCS suggested that a northwestern Atlantic haplotype was ancestral to Hiatella L (Fig. 6A). Subsequent analysis of intraspecific variation in H. arctica only examined specimens belonging to Hiatella L. The other lineages are likely to represent sibling species unrecognized in current taxonomy.

Figure 6.

A. Median-joining network of COI haplotypes for Hiatella arctica s. l. constructed using Network v. 4.6.1. All mutational steps are three or less unless indicated by a numeral. Presumptive ancestral haplotypes are marked with an asterisk. B. Proportions of four Hiatella COI lineages at each sampling location in this study (CI, Cook Inlet, Alaska; HG, Haidi Gwaii, British Columbia; VI, Vancouver Island, British Columbia; SG, Strait of Georgia, British Columbia; CH, Churchill, Manitoba; BB, Bonne Bay, Newfoundland; SA, St Andrews, New Brunswick).

Figure 6.

A. Median-joining network of COI haplotypes for Hiatella arctica s. l. constructed using Network v. 4.6.1. All mutational steps are three or less unless indicated by a numeral. Presumptive ancestral haplotypes are marked with an asterisk. B. Proportions of four Hiatella COI lineages at each sampling location in this study (CI, Cook Inlet, Alaska; HG, Haidi Gwaii, British Columbia; VI, Vancouver Island, British Columbia; SG, Strait of Georgia, British Columbia; CH, Churchill, Manitoba; BB, Bonne Bay, Newfoundland; SA, St Andrews, New Brunswick).

The three lineages in the M. balthica complex showed a maximum intraspecific divergence of 15.4% (Fig. 4). Macoma b. balthica was more broadly distributed than M. petalum and M. balthica NFLD, occurring in Alaska, British Columbia, Labrador and Manitoba (Fig. 7A, B). Macoma petalum had a narrower distribution, being found in New Brunswick, Newfoundland and Prince Edward Island, while M. balthica NFLD was only detected in Newfoundland (Fig. 7A, B). A deep divergence between M. b. balthica and M. petalum was corroborated by the high number of mutational steps that separated the two clusters (Fig. 7A). Conversely, only four mutational steps separated M. b. balthica and M. balthica NFLD in the haplotype network (Fig. 7A).

Figure 7.

A. Median-joining network of COI haplotypes for Macoma balthica s. l. constructed using Network v. 4.6.1. All mutational steps are three or less unless indicated by a numeral. Presumptive ancestral haplotypes are marked with an asterisk. B. Proportions of M. b. balthica, M. balthica NFLD and M. petalum at each sampling location, based on COI data (CI, Cook Inlet, Alaska; HG, Haidi Gwaii, British Columbia; CH, Churchill, Manitoba; MK, Makkovik, Labrador; BB, Bonne Bay, Newfoundland; PI, Pinette, Prince Edward Island; SA, St Andrews, New Brunswick).

Figure 7.

A. Median-joining network of COI haplotypes for Macoma balthica s. l. constructed using Network v. 4.6.1. All mutational steps are three or less unless indicated by a numeral. Presumptive ancestral haplotypes are marked with an asterisk. B. Proportions of M. b. balthica, M. balthica NFLD and M. petalum at each sampling location, based on COI data (CI, Cook Inlet, Alaska; HG, Haidi Gwaii, British Columbia; CH, Churchill, Manitoba; MK, Makkovik, Labrador; BB, Bonne Bay, Newfoundland; PI, Pinette, Prince Edward Island; SA, St Andrews, New Brunswick).

Patterns of genetic diversity

Table 1 provides population genetic parameters for populations of each species, excluding Hiatella lineages I, K and N. Populations in Churchill, Manitoba showed the lowest variation for Hiatella L, while populations in Labrador and Prince Edward Island had the lowest genetic diversity for M. b. balthica and M. petalum, respectively. Haplotype diversity was highest in populations of Hiatella L from Alaska and in populations of M. b. balthica and M. petalum from Newfoundland. However, nucleotide diversity was highest in Hiatella L populations from New Brunswick and in populations of M. b. balthica and M. petalum from Alaska and Newfoundland. Nucleotide and haplotype diversity did not differ significantly between species. Tajima's D values were significantly different from zero for populations of Hiatella L in New Brunswick and populations of M. b. balthica in Alaska and Labrador.

Table 1.

Genetic diversity in populations of Hiatella L, Macoma b. balthica, M. balthica NFLD and M. petalum as measured by number of sequences (n), number of haplotypes (H), haplotype diversity (h), nucleotide diversity (Π) and Tajima's D for COI data.

Species Population n H h Π D 
HiatellaAlaska 43 23 0.94 0.0072 −1.40 
British Columbia 16 10 0.92 0.0052 0.09 
Manitoba 19 0.71 0.0058 −0.55 
New Brunswick 25 18 0.93 0.0084 −1.58* 
Nova Scotia 60 25 0.88 0.0077 −1.18 
Macoma b. balthica Alaska 38 12 0.69 0.0045 −1.48* 
British Columbia 0.83 0.0035 −0.75 
Labrador 39 0.20 0.0005 −1.57* 
Manitoba 19 0.66 0.0019 −0.20 
Macoma balthica NFLD Macoma petalum Newfoundland 0.90 0.0056 
New Brunswick 51 12 0.83 0.0090 0.01 
Newfoundland 0.016 
Prince Edward Island 33 0.74 0.0053 −1.04 
Species Population n H h Π D 
HiatellaAlaska 43 23 0.94 0.0072 −1.40 
British Columbia 16 10 0.92 0.0052 0.09 
Manitoba 19 0.71 0.0058 −0.55 
New Brunswick 25 18 0.93 0.0084 −1.58* 
Nova Scotia 60 25 0.88 0.0077 −1.18 
Macoma b. balthica Alaska 38 12 0.69 0.0045 −1.48* 
British Columbia 0.83 0.0035 −0.75 
Labrador 39 0.20 0.0005 −1.57* 
Manitoba 19 0.66 0.0019 −0.20 
Macoma balthica NFLD Macoma petalum Newfoundland 0.90 0.0056 
New Brunswick 51 12 0.83 0.0090 0.01 
Newfoundland 0.016 
Prince Edward Island 33 0.74 0.0053 −1.04 

Tajima's D values significantly different from zero are marked with an asterisk.

Population structure

A comparison of genetic structure revealed more phylogeographic structure in Hiatella L (FST = 0.16) than M. b. balthica (FST = 0.09) and M. petalum (FST = 0.06). Most of the genetic variation in Hiatella L and M. b. balthica was partitioned within individual populations (Table 2). Fixation indices provided further insight into the partitioning of variation in Hiatella L and M. b. balthica. In Hiatella L, populations from Nova Scotia and New Brunswick showed little divergence (FST = 0.02), while those from British Columbia and New Brunswick were the most divergent (FST = 0.30). Populations of M. b. balthica from Alaska and British Columbia were the most similar (FST = 0.02), while those from British Columbia and Labrador were the least similar (FST = 0.34). Macoma petalum also showed a high FST (0.22) between populations in Newfoundland and Prince Edward Island, although these sites are only separated by the 40-km wide Northumberland Strait. When Slatkin's linearized FST values were plotted against geographic distance, only Hiatella L showed evidence of IBD (Fig. 8A) with a strong, positive correlation (R2 = 0.86) and a Mantel test confirmed its significance (P = 0.004). IBD results for M. petalum and M. balthica NFLD were excluded due to insufficient data.

Table 2.

Overall genetic structure measured by AMOVA for the widespread lineages Hiatella L and Macoma b. balthica using COI data.

Species Variation Df Sum of squares Variance % Variation P-value 
HiatellaAmong 59 0.41 16.3 
Within 158 335 2.12 83.7  
Macoma b. balthica Among 0.05 8.8 
Within 96 49 0.51 91.2  
Species Variation Df Sum of squares Variance % Variation P-value 
HiatellaAmong 59 0.41 16.3 
Within 158 335 2.12 83.7  
Macoma b. balthica Among 0.05 8.8 
Within 96 49 0.51 91.2  

P values <0.05 were treated as significant.

Figure 8.

FST (Slatkin's linearized) values plotted against geographic distance for populations of Hiatella L (A) and Macoma b. balthica (B) to examine the extent of isolation by distance using COI data. Results of Mantel test for significance are: A: R2 = 0.86, P = 0.004; B: R2 = 0.18, P = 0.15.

Figure 8.

FST (Slatkin's linearized) values plotted against geographic distance for populations of Hiatella L (A) and Macoma b. balthica (B) to examine the extent of isolation by distance using COI data. Results of Mantel test for significance are: A: R2 = 0.86, P = 0.004; B: R2 = 0.18, P = 0.15.

DISCUSSION

Systematics and cryptic diversity

Incorporating sequences for recognized subspecies provided insight into the taxonomic status of both bivalve complexes. Our results from both a mitochondrial (COI) and nuclear (ITS2) marker support findings that Macoma b. balthica occurs in the northeastern Pacific and the Arctic (e.g. Nikula et al., 2007), but our work has also shown that this clade occurs in the northwest Atlantic (Labrador). Moreover, a maximum intraspecific divergence of 15.4% between M. balthica NFLD and M. petalum in their zone of sympatry in eastern Canada supports their recognition as sibling species. Our results also provide evidence for a third clade in Newfoundland which is 3% divergent from M. b. balthica, although this result was only supported by COI data. Conversely, these sequences clustered into a single clade in the ITS2 NJ tree, likely due to the slower evolutionary rate of this nuclear gene. Nonetheless, these results corroborate previous work that revealed three unique clades in the M. balthica complex in Canada (Layton et al., 2014). The discovery of a separate cluster in Newfoundland (M. balthica NFLD) is also concordant with the results of Nikula et al. (2007), who found a unique subclade of M. b. balthica (C) in the northwestern Atlantic using the COIII gene. Future work should aim to determine if M. b. balthica and M. balthica NFLD show incipient or complete reproductive isolation. The discordance between COI and ITS2 for a single Macoma specimen, but without heteroplasmy, may reflect an advanced generation hybrid, so the potential for hybridization between these taxa should be investigated.

Similarly, incorporating Hiatella COI sequences from Laakkonen et al. (2015) confirmed the presence of multiple lineages in North America, including the widespread lineage L and additional lineages in the NE Pacific (I) and the NW Atlantic (K). Four unique Hiatella clades were detected by previous work in Canada (Layton et al., 2014), but the extensive molecular dataset generated by Laakkonen et al. (2015) allowed us to assign clade names to them. Sequence divergences (COI) between the four Hiatella clades ranged from 11.8 to 23.5%, suggesting that species status is warranted. Furthermore, our COI dataset provides the first evidence for Hiatella K in the northwestern Atlantic and for a novel lineage (Hiatella N) in the northeastern Pacific. These results add to earlier findings of genetic diversity in this complex (Layton et al., 2014; Laakkonen et al., 2015) and emphasize the need for a taxonomic revision of Hiatella.

Comparison of genetic diversity and structure

Although a widespread species was present in both complexes, Hiatella L showed more regional sequence variation in Canadian waters than M. b. balthica. Contrasts in genetic variation between species with similar life-history strategies have been observed in other molluscs, including the direct-developing gastropods Nucella lamellosa and N. ostrina (Marko, 2004). Hiatella L and M. b. balthica also differ in their habitat and feeding behaviour (Newell, 1965; Ali, 1970; Hines & Comtois, 1985), which may influence population structure. In fact, Marko (2004) suggested that ecological factors may often be more important than dispersal ability in determining population structure. Nevertheless, the presence of the widespread lineages Hiatella L and M. b. balthica suggests that these species have considerable gene flow among their populations. Particularly for Hiatella L, low regional divergence may reflect the use of secondary dispersal mechanisms, including transportation of juveniles in ballast water and adults in kelp holdfasts (Helmuth, Veit & Holberton, 1994). In any case, patterns of genetic structure in marine invertebrate populations of North America are largely influenced by past glacial cycles.

Patterns of population fragmentation during glacial periods are often reflected in measures of genetic diversity. Reduced genetic diversity, due to glaciation and subsequent population expansion, has been demonstrated in populations of marine fishes and in Mya arenaria (Bernatchez & Wilson, 1998; Lasota, Hummel & Wolowicz, 2004; Strasser & Barber, 2009). However, some taxa exposed to repeated glaciations exhibit high genetic diversity. For instance, the high allelic diversity of Nucella lapillus has been linked to rapid re-expansion following a severe population bottleneck (Colson & Hughes, 2004). A similar process could explain the high diversity in populations of both Hiatella L and M. petalum from New Brunswick, despite severe glacial conditions in the region (Briggs, 1970; Wares & Cunningham, 2001). Northeastern Pacific populations of the widespread lineages Hiatella L and M. b. balthica were also diverse, which may reflect their foundation through postglacial admixture (Kelly et al., 2006; Sakaguchi et al., 2011). By contrast, the low diversity of populations at Churchill (Manitoba) may be a consequence of the relatively recent (8,000 year BP) formation of Hudson Bay (Ashworth, 1996).

Glacial refugia and secondary contact

High haplotype diversity in some northern populations can be invoked as evidence for persistence in glacial refugia (Marko, 2004). Therefore, high haplotype diversity in populations from the northeastern Pacific, Newfoundland and New Brunswick suggests the possible existence of multiple refugia in these areas. In the Pacific Ocean alone, recent work has shown that populations of Atlantic Cod (Gadus morhua) survived in two refugia in the northeastern Pacific and two refugia in the northwestern Pacific (Bigg, 2014). Moreover, Nikula et al. (2007) suggested that multiple refugia existed for M. balthica, likely explaining the deep divergences in this complex.

In addition to refugial sites, genetic diversity can be elevated in zones of admixture as well (Petit et al., 2003; Kelly et al., 2006; Sakaguchi et al., 2011). The mixing of phylogenetic lineages in admixture zones combines variants derived from two or more refugia (Provan & Bennett, 2008). Secondary contact and postglacial admixture of allopatric populations previously isolated in different marine refugia may explain the pattern of haplotype distribution in Hiatella L. A section of the haplotype network for Hiatella L consists of haplotypes from only the northwestern Atlantic and Manitoba, while the remaining section of the network contains haplotypes from all localities, suggesting secondary contact and subsequent admixture of divergent lineages. The opening of the Bering Strait allowed for interchange between faunas in Canada's three oceans, with the majority of migrations occurring from the Pacific to the Atlantic (Briggs, 1970; Vermeij, 1991). This pattern of haplotype partitioning may reflect multiple trans-Arctic migrations and subsequent mixing of divergent lineages in the Atlantic, a pattern observed in both M. balthica and Mytilus edulis (Väinölä, 2003; Nikula et al., 2007). It is clear that the genetic structure of these planktotrophic bivalves has not only been influenced by their dispersal potential, but shaped by Canada's glacial history.

ACKNOWLEDGEMENTS

We thank Barry McDonald, Christina Carr, Katrin Iken, Melissa Frey, Paolo Pierossi, Rick Harbo, Sarah Hardy and Suzanne Dufour for aid in specimen collection, and staff at the CCDB for support in sequence acquisition. We also thank Elizabeth Boulding and Sean Prosser for advice on molecular techniques and feedback on data analysis. Fieldwork in Churchill, Manitoba was conducted under permits issued by Manitoba Conservation Wildlife and Ecosystem Protection to the Churchill Northern Studies Centre (CNSC) for research in the Churchill Wildlife Management Area. Collections in Alaska were conducted under a fish resource permit granted to Sarah Hardy by the State of Alaska Department of Fish and Game for scientific/educational purposes. Collections in New Brunswick and Labrador were conducted under experimental licenses from Fisheries and Oceans Canada. This research was funded, in part, by the NSERC Canadian Healthy Oceans Network and by a NSERC Discovery Grant to PDNH. Field work was aided by a Northern Scientific Training Program grant to KKSL from Indian and Northern Affairs Canada. Sequence analysis was enabled by funding from the government of Canada through Genome Canada and the Ontario Genomics Institute in support of the International Barcode of Life Project.

REFERENCES

ADDISON
J.A.
,
HART
M.W.
2005
.
Colonization, dispersal, and hybridization influence phylogeography of North Atlantic sea urchins (Strongylocentrotus droebachiensis)
.
Evolution
 ,
59
:
532
543
.
ALI
R.M.
1970
.
The influence of suspension density and temperature on the filtration rate of Hiatella arctica
.
Marine Biology
 ,
6
:
291
302
.
ALISON
R.C.
,
MARINCOVICH
L.J.
1982
.
A late Oligocene or earliest Miocene molluscan fauna from Sitkinak Island, Alaska
. In:
Jurassic (Oxfordian and late Callovian) ammonites from the western interior region of the United States
  (
Imlay
R.W.
, ed.).
Geological Survey professional paper vol. 1232, Washington
.
ARNDT
A.
,
SMITH
M.J.
1998
.
Genetic diversity and population structure in two species of sea cucumber: differing patterns according to mode of development
.
Molecular Ecology
 ,
7
:
1053
1064
.
ASHWORTH
A.
1996
.
The response of arctic Carabidae (Coleoptera) to climate change based on the fossil record of the Quaternary
.
Annales Zoologici Fennici
 ,
33
:
125
131
.
BANDELT
H-J.
,
FORSTER
P.
,
RÖHL
A.
1999
.
Median-joining networks for inferring intraspecific phylogenies
.
Molecular Biology and Evolution
 ,
16
:
37
48
.
BERNATCHEZ
L.
,
WILSON
C.C.
1998
.
Comparative phylogeography of Nearctic and Palearctic fishes
.
Molecular Ecology
 ,
7
:
431
452
.
BIGG
G.R.
2014
.
Environmental confirmation of multiple ice age refugia for Pacific cod, Gadus macrocephalus
.
Evolutionary Ecology
 ,
28
:
177
191
.
BRIGGS
J.C.
1970
.
A faunal history of the North Atlantic Ocean
.
Systematic Zoology
 ,
19
:
19
34
.
CLEMENT
M.
,
POSADA
D.
,
CRANDALL
K.
2000
.
TCS: a computer program to estimate gene genealogies
.
Molecular Ecology
 ,
9
:
1657
1660
.
COAN
E.V.
1971
.
The northwest American Tellinidae
.
Veliger
 ,
14
:
1
63
.
COLSON
I.
,
HUGHES
R.N.
2004
.
Rapid recovery of genetic diversity of dogwhelk (Nucella lapillus L.) populations after local extinction and recolonization contradicts predictions from life-history characteristics
.
Molecular Ecology
 ,
13
:
2223
2233
.
DAPPORTO
L.
2009
.
Speciation in Mediterranean refugia and post-glacial expansion of Zerynthia polyxena (Lepidoptera, Papilionidae)
.
Journal of Zoological Systematics and Evolutionary Research
 ,
48
:
229
237
.
DODSON
J.J.
,
TREMBLAY
S.
,
COLOMBANI
F.
,
CARSCADDEN
J.E.
,
LECOMTE
F.
2007
.
Trans-Arctic dispersals and the evolution of a circumpolar marine fish species complex, the capelin (Mallotus villosus)
.
Molecular Ecology
 ,
16
:
5030
5043
.
EXCOFFIER
L.
,
LISCHER
H.E.L.
2010
.
Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows
.
Molecular Ecology Resources
 ,
10
:
564
567
.
HARPER
F.M.
,
HART
M.W.
2007
.
Morphological and phylogenetic evidence for hybridization and introgression in a sea star secondary contact zone
.
Invertebrate Biology
 ,
126
:
373
384
.
HELMUTH
B.
,
VEIT
R.R.
,
HOLBERTON
R.
1994
.
Long distance dispersal of subantarctic brooding bivalve (Gaimardia trapesina) by kelp rafting
.
Marine Biology
 ,
120
:
421
426
.
HEWITT
G.
2000
.
The genetic legacy of the Quaternary ice ages
.
Nature
 ,
405
:
907
913
.
HEWITT
G.M.
1996
.
Some genetic consequences of ice ages, and their role in divergence and speciation
.
Biological Journal of the Linnaean Society
 ,
58
:
247
276
.
HINES
A.H.
,
COMTOIS
K.L.
1985
.
Vertical distribution of infauna in sediments of a subestuary of central Chesapeake Bay
.
Estuaries
 ,
8
:
296
304
.
IVANOVA
N.V.
,
FAZEKAS
A.J.
,
HEBERT
P.D.N.
2008
.
Semi-automated, membrane-based protocol for DNA isolation from plants
.
Plant Molecular Biology Reporter
 ,
26
:
186
198
.
JABLONSKI
D.
1986
.
Larval ecology and macroevolution in marine invertebrates
.
Bulletin of Marine Science
 ,
39
:
565
587
.
JI
Y.
,
ZHANG
D.
,
HE
L.
2003
.
Evolutionary conservation and versatility of a new set of primers for amplifying the ribosomal internal transcribed spacer regions in insects and other invertebrates
.
Molecular Ecology Notes
 ,
3
:
581
585
.
KEEVER
C.C.
,
SUNDAY
J.
,
PURITZ
J.B.
,
ADDISON
J.A.
,
TOONEN
R.J.
,
GROSBERG
R.K.
,
HART
M.W.
2009
.
Discordant distribution of populations and genetic variation in a sea star with high dispersal potential
.
Evolution
 ,
63
:
3214
3227
.
KELLY
D.W.
,
MUIRHEAD
J.R.
,
HEATH
D.D.
,
MACISAAC
H.J.
2006
.
Contrasting patterns in genetic diversity following multiple invasions of fresh and brackish waters
.
Molecular Ecology
 ,
15
:
3641
3653
.
KIMURA
M.
1980
.
A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences
.
Journal of Molecular Evolution
 ,
16
:
111
120
.
KYLE
C.J.
,
BOULDING
E.G.
2000
.
Comparative population genetic structure of marine gastropods (Littorina spp.) with and without pelagic larval dispersal
.
Marine Biology
 ,
137
:
835
845
.
LAAKKONEN
H.M.
,
STRELKOV
P.
,
VÄINÖLÄ
R.
2015
.
Molecular lineage diversity and inter-oceanic biogeographical history in Hiatella (Mollusca, Bivalvia)
.
Zoological Scripta
 ,
44
:
383
402
.
LASOTA
R.
,
HUMMEL
H.
,
WOLOWICZ
M.
2004
.
Genetic diversity of European populations of the invasive soft-shell clam Mya arenaria (Bivalvia)
.
Journal of the Marine Biological Association of the United Kingdom
 ,
84
:
1051
1056
.
LASSMANN
T.
,
FRINGS
O.
,
SONNHAMMER
E.L.
2009
.
Kalign2: high-performance multiple alignment of protein and nucleotide sequences allowing external features
.
Nucleic Acids Research
 ,
37
:
858
865
.
LAYTON
K.K.S.
,
MARTEL
A.L.
,
HEBERT
P.D.N.
2014
.
Patterns of DNA barcode variation in Canadian marine molluscs
.
PLoS One
 ,
9
:
e95003
.
LEE
H.J.
,
BOULDING
E.G.
2009
.
Spatial and temporal population structure of four northeastern Pacific littorinid gastropods: the effect of mode of larval development on variation at one mitochondrial and two nuclear DNA markers
.
Molecular Ecology
 ,
18
:
2165
2184
.
MAGGS
C.A.
,
CASTILHO
R.
,
FOLTZ
D.
,
HENZLER
C.
,
JOLLY
M.T.
,
KELLY
J.
,
OLSEN
J.
,
PEREZ
K.E.
,
STAM
W.
,
VÄINÖLÄ
R.
,
VIARD
F.
,
WARES
J.
2008
.
Evaluating signatures of glacial refugia for North Atlantic benthic marine taxa
.
Ecology
 ,
89
:
S108
S122
.
MANDRYK
C.A.S.
,
JOSENHANS
H.
,
FEDJE
D.W.
,
MATHEWES
W.
2001
.
Late Quaternary paleoenvironments of Northwestern North America: implications for inland versus coastal migration routes
.
Quaternary Science Review
 ,
20
:
301
314
.
MARKO
P.B.
2004
.
‘What's larvae got to do with it?’ Disparate patterns of post-glacial population structure in two benthic marine gastropods with identical dispersal potential
.
Molecular Ecology
 ,
13
:
597
611
.
MARKO
P.B.
,
HOFFMAN
J.M.
,
EMME
S.A.
,
MCGOVERN
T.M.
,
KEEVER
C.C.
,
COX
L.N.
2010
.
The ‘expansion–contraction’ model of Pleistocene biogeography: rocky shores suffer a sea change?
Molecular Ecology
 ,
19
:
146
160
.
MARSHALL
B.
,
GOFAS
S
.
2015
.
Hiatella arctica (Linnaeus, 1767). Accessed through: World Register of Marine Species at http://www.marinespecies.org/aphia.php?p=taxdetails&id=140103
.
MEEHAN
B.W.
1985
.
Genetic comparison of Macoma balthica (Bivalvia, Telinidae) from the eastern and western North Atlantic Ocean
.
Marine Ecology Progress Series
 ,
22
:
69
76
.
MEEHAN
B.W.
,
CARLTON
J.T.
,
WENNE
R.
1989
.
Genetic affinities of the bivalve Macoma balthica from the Pacific coast of North America: evidence for recent introduction and historical distribution
.
Marine Biology
 ,
102
:
235
241
.
NEI
M.
1987
.
Molecular evolutionary genetics
 .
Columbia University Press
,
New York
.
NEWELL
R.C.
1965
.
The role of detritus in the nutrition of two marine deposit feeders, the prosobranch Hydrobia ulvae and the bivalve Macoma balthica
.
Proceedings of the Royal Zoological Society of London
 ,
144
:
25
45
.
NIKULA
R.
,
STRELKOV
P.
,
VÄINÖLÄ
R.
2007
.
Diversity and trans-arctic invasion history of mitochondrial lineages in the North Atlantic Macoma balthica complex (Bivalvia: Tellinidae)
.
Evolution
 ,
61
:
928
941
.
PETIT
R.J.
,
AGUINAGALDE
I.
,
DE BEAULIEU
J-L.
,
BITTKAU
C.
,
BREWER
S.
,
CHEDDAI
R.
,
ENNOS
R.
,
FINESCHI
S.
,
GRIVET
D.
,
LASCOUX
M.
,
MOHANTY
A.
,
MÜLLER-STARCK
G.
,
DEMESURE-MUSCH
B.
,
PALMÉ
A.
,
MARTÍN
J.P.
,
RENDELL
S.
,
VENDRAMIN
G.G.
2003
.
Glacial refugia: hotspots but not melting pots of genetic diversity
.
Science
 ,
300
:
1563
1565
.
PROVAN
J.
,
BENNETT
K.D
.
2008
.
Phylogeographic insights into cryptic glacial refugia
.
Trends in Ecology and Evolution
 ,
23
:
564
571
.
RATNASINGHAM
S.
,
HEBERT
P.D.N
.
2007
.
BOLD: The Barcode of Life Data System
. .
Molecular Ecology Notes
 ,
7
:
355
364
.
REID
D.G.
1990
.
Trans-Arctic migration and speciation induced by climate change: the biogeography of Littorina (Mollusca: Gastropoda)
.
Bulletin of Marine Science
 ,
47
:
35
49
.
ROHLING
R.J.
,
FENTON
M.
,
JORISSEN
F.J.
,
BERTRAND
P.
,
GANSSEN
G.
,
CAULET
J.P.
1998
.
Magnitudes of sea-level lowstands of the past 500,000 years
.
Nature
 ,
394
:
162
165
.
SAKAGUCHI
S.
,
TAKEUCHI
Y.
,
YAMASAKI
M.
,
SAKURAI
S.
,
ISAGI
Y.
2011
.
Lineage admixture during postglacial range expansion is responsible for the increased gene diversity of Kalopanax septemlobus in a recently colonised territory
.
Heredity
 ,
107
:
338
348
.
SLATKIN
M.
1993
.
Isolation by distance in equilibrium and non-equilibrium populations
.
Evolution
 ,
47
:
264
279
.
STRASSER
C.A.
,
BARBER
P.H.
2009
.
Limited genetic variation and structure in softshell clams (Mya arenaria) across their native and introduced range
.
Conservation Genetics
 ,
10
:
803
814
.
TAJIMA
F.
1989
.
Statistical method for testing the neutral mutation hypothesis by DNA polymorphism
.
Genetics
 ,
123
:
585
595
.
TAMURA
K.
1992
.
Estimation of the number of nucleotide substitutions when there are strong transition-transversion and G+C content biases
.
Molecular Biology and Evolution
 ,
9
:
678
687
.
TAMURA
K.
,
STECHER
G.
,
PETERSON
D.
,
FILIPSKI
A.
,
KUMAR
S.
2013
.
MEGA6: Molecular Evolutionary Genetics Analysis Version 6.0
.
Molecular Biology and Evolution
 ,
30
:
2725
2729
.
TAYLOR
E.B.
,
DODSON
J.J.
1994
.
A molecular analysis of relationships and biogeography within a species complex of Holarctic fish (genus Osmerus)
.
Molecular Ecology
 ,
3
:
235
248
.
VÄINÖLÄ
R.
2003
.
Repeated trans-Arctic invasions in littoral bivalves: molecular zoogeography of the Macoma balthica complex
.
Marine Biology
 ,
143
:
935
946
.
VERMEIJ
G.
1991
.
Anatomy of an invasion: the trans-Arctic interchange
.
Paleobiology
 ,
17
:
281
307
.
WARES
J.P.
,
CUNNINGHAM
C.W.
2001
.
Phylogeography and historical ecology of the North Atlantic intertidal
.
Evolution
 ,
55
:
2455
2469
.
WARNER
B.G.
,
MATHEWES
R.W.
,
CLAGUE
J.J.
1982
.
Ice-free conditions on the Queen Charlotte Islands, British Columbia, at the height of the late Wisconsin glaciation
.
Science
 ,
218
:
675
677
.
WEIR
B.S.
,
COCKERHAM
C.C.
1984
.
Estimating F-statistics for the analysis of population structure
.
Evolution
 ,
38
:
1358
1370
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.