Geographic patterns of genetic diversity in two species complexes of Canadian marine bivalves

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 theM. 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 (Va¨ino¨la¨, 2003;Nikula, Strelkov & Va¨ino¨la¨, 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, Va¨ino¨la¨( 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 & Va¨ino¨la¨, 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.

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 220 8C 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 Va¨ino¨la¨(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.

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 8C 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 ml of ddH 2 O.After re-suspension, 2 ml of each DNA extract was placed into a well in a separate plate with 18 ml ddH 2 O 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 TAGACTTCT GGGTGCCCGAAAAACCA for H. arctica and MMacF1/LepR1: CTTTTATTAGCTGCACCTGATAT and TAAACTTCTGGA TGTCCAAAAAATCA 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 ml of diluted DNA and the following reagents to generate a 12.5 ml PCR reaction mix: 6.25 ml 10% trehalose, 2 ml ddH 2 O, 1.25 ml 10Â PCR buffer, 0.625 ml MgCl 2 (50 mM), 0.125 ml of each forward and reverse primer (10 mM), 0.0625 ml dNTP (10 mM) and 0.06 ml Platinum Taq polymerase.The thermocycling regime consisted of one cycle of 1 min at 94 8C, 40 cycles of 40 s at 94 8C, 40 s at 52 8C and 1 min at 72 8C, and finally 5 min at 72 8C.An E-GelH 96 (Invitrogen) was used to check 3 ml 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 Va¨ino¨la¨(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 & Ro¨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 (F ST ) were estimated with 1,000 permutations to determine the partitioning of variance (Weir & Cockerham, 1984).Slatkin's linearized F ST (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.

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).
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.
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).
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.
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).

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,

Population structure
A comparison of genetic structure revealed more phylogeographic structure in Hiatella L (F ST ¼ 0.16) than M. b. balthica (F ST ¼ 0.09) and M. petalum (F ST ¼ 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 (F ST ¼ 0.02), while those from British Columbia and New Brunswick were the most divergent (F ST ¼ 0.30).Populations of M. b. balthica from Alaska and British Columbia were the most similar (F ST ¼ 0.02), while those from British Columbia and Labrador were the least similar (F ST ¼ 0.34).Macoma petalum also showed a high F ST (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 F ST values were plotted against geographic distance, only Hiatella L showed evidence of IBD (Fig. 8A) with a strong, positive correlation (R 2 ¼ 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.

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 (Va¨ino¨la¨, 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.

Figure 1 .
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 3 .
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 4 .
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 Va¨ino¨la¨(2003) and Nikula et al. (2007) are provided and clades discovered in this work are indicated by asterisks.

Figure 5 .
Figure 5. Neighbour-joining tree (ITS2) based on a Tamura-3-parameter substitution model for Macoma balthica s. l.Subspecies names suggested by Va¨ino¨la¨(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 6 .
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 7 .
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 8 .
Figure 8. F ST (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: R 2 ¼ 0.86, P ¼ 0.004; B: R 2 ¼ 0.18, P ¼ 0.15.

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 (P) and Tajima's D for COI data.Tajima's D values significantly different from zero are marked with an asterisk.

Table 2 .
Overall genetic structure measured by AMOVA for the widespread lineages Hiatella L and Macoma b. balthica using COI data.
P values ,0.05 were treated as significant.