Trans-Arctic vicariance in Strongylocentrotus sea urchins

The sea urchins Strongylocentotus pallidus and S. droebachiensis first invaded the Atlantic Ocean from the Pacific following the opening of the Bering seaway in the late Miocene. While trans-Arctic dispersal during the Pleistocene is thought to have maintained species’ integrity, a recent genomic analysis identified a reproductively isolated cryptic species within S. droebachiensis. Based on previous studies, the distribution of one of these lineages (S. droebachiensis W) includes the shallow water habitats of the northwest Atlantic and Pacific, while the other (S. droebachiensis E) is found throughout the shallow habitat in the northeast but is mostly restricted to deep habitats (>65 m) in the northwest Atlantic. However, since genetic variation within S. droebachiensis has been largely unstudied in the north Pacific and Arctic oceans, the biogeography of the cryptic species is not well known, and it is difficult to identify the mechanisms driving population subdivision and speciation. Here we use population genetic analyses to characterize the distribution of each species, and to test hypotheses about the role of vicariance in the evolution of systematic and genomic divergence within the genus. We collected individuals of all three Strongylocentrotus species (n = 365) from 10 previously unsampled locations in the northeast Pacific and north Atlantic (Labrador Sea and Norway), and generated mtDNA sequence data for a 418 bp fragment of cytochrome c oxidase subunit I (COI). To assess the biogeography of all three species, we combined our alignment with five previously published data sets (total n = 789) and used statistical parsimony and maximum likelihood to identify species and characterize their distribution within and among oceans. Patterns of haplotype sharing, pairwise FST, and hierarchical analyses of molecular variance (AMOVA) identified trans-Arctic dispersal in S. pallidus and S. droebachiensis W, but other than 5 previously reported singletons we failed to detect additional mtDNA haplotypes of S. droebachiensis E in the north Pacific. Within the Atlantic, patterns of habitat segregation suggests that temperature may play a role in limiting the distribution of S. droebachiensis E, particularly throughout the warmer coastal waters along the coast of Nova Scotia. Our results are consistent with the cycles of trans-Arctic dispersal and vicariance in S. pallidus and S. droebachiensis W, but we suggest that the evolution of Atlantic populations of S. droebachiensis E has been driven by persistent trans-Arctic vicariance that may date to the initial invasion in the late Pliocene.


INTRODUCTION
The global biogeography of Strongylocentrotid sea urchins was shaped by the trans-Arctic interchange following the initial opening of the Bering seaway in the late Miocene (5.5-5.0 Mya;Marincovich & Gladenkov, 2001;Gladenkov et al., 2002), and fossil evidence reveals that Strongylocentrotus droebachiensis (along with many other Pacific taxa) reached western Europe by the late Pliocene (Durham & MacNeil, 1967). Following the initial invasion, eustatic sea level changes during the Pleistocene ice ages (2.4-0.2 Mya) periodically restricted dispersal across the Arctic Basin, causing widespread isolation and vicariance in the north Atlantic (Hewitt, 1996;Cunningham & Collins, 1998). As a result of these processes, molecular evidence from trans-Arctic taxa indicates a complex pattern of interand intra-specific divergence, with species positioned along a continuum between complete reproductive isolation and panmictic populations (see Laakkonen et al., 2021). There has been consensus that both genetic diversity and species integrity have been maintained among populations of Strongylocentrotus sea urchins from the Pacific (ancestral) and Atlantic (colonized) as a result of gene flow across the Arctic between 0.40-0.11 Mya (Palumbi & Kessing, 1991;Palumbi & Wilson, 1990;Addison & Hart, 2005;Laakkonen et al., 2021). However, the recent discovery of a cryptic species within the north Atlantic population of S. droebachiensis (Addison & Kim, 2018) suggests a more complicated role of vicariance in the evolution of the genus, demanding a re-evaluation of both the biogeography and population genetics throughout the range.
Like many marine invertebrates, Strongylocentrotus sea urchins are broadcast spawners with long-lived planktonic larvae (4 to 21 weeks; Strathmann, 1978) capable of high rates of dispersal with gene flow. Genetic studies of S. droebachiensis (Müller, 1776) in the north Atlantic detected local panmixis in the northwest (Addison & Hart, 2004;Addison & Hart, 2005), a small but significant latitudinal gradient in the northeast (Nordberg et al., 2016), and significant population substructure between the east and west coasts (Addison & Hart, 2004;Addison & Hart, 2005;Harper, Addison & Hart, 2007). Patterns of genetic variation at both microsatellites and mtDNA reveal lower levels of diversity in the eastern populations compared to the west, and pairwise F ST suggests that northwest populations are genetically more similar to the Pacific than to those from the northeast Atlantic (Addison & Hart, 2004;Addison & Hart, 2005). A more detailed study of S. droebachiensis in the northwest Atlantic indicates strong habitat segregation and reproductive isolation between distinct east and west mtDNA haplogroups, where the vast majority shallow-water coastal samples (<30 m) were a subset of haplotypes shared between the Pacific and the northwest Atlantic, and the offshore deep-water (>65 m) samples were identical to (or clustered with) populations from Norway and Iceland (Addison & Kim, 2018). These genetic patterns broadly correspond to variation in sperm morphology reported between sea urchins from the northeast Atlantic and those from the northwest Atlantic and Pacific (Manier & Palumbi, 2008;Marks et al., 2008). Together, these patterns of genetic and morphological divergence indicate that trans-Atlantic variation reflects species level differences and are not the result of limited gene flow and genetic drift between allopatric populations. Thus, there is strong evidence that the north Atlantic harbours two reproductively isolated species of Strongylocentrotus sea urchins: S. droebachiensis W whose distribution includes the shallow water habitats of the northwest Atlantic and Pacific, and S. droebachiensis E that is distributed throughout the shallow habitat in the northeast and deep habitats (>65 m) at lower latitudes in the northwest Atlantic.
Cycles of trans-Arctic dispersal and vicariance have played an important role in the evolution of new species in a variety of Pacific and Atlantic taxa, including molluscs, crustaceans, echinoderms, polychaetes, fishes, mammals, and algae (e.g., Wares, 2001;Carr et al., 2011;Layton, Corstorphine & Hebert, 2016;Neiva et al., 2018;Bringloe, Verbruggen & Saunders, 2020;Laakkonen et al., 2021). The surprising discovery of two independent lineages of S. droebachiensis suggests that allopatric speciation in this genus may have also followed the initial trans-Arctic invasion. While vicariance during the Pliocene-Early Pleistocene resulted in the speciation of the sea stars Asterias forbesi and A. rubens in the Atlantic (Wares, 2001), more recent vicariant histories during the Middle Pleistocene (1-0.2 Mya) have resulted in reciprocal monophyly and divergence among populations of Solaster endeca, Pteraster militanus, and Crosster papposus (Layton, Corstorphine & Hebert, 2016). Interoceanic divergence between allopatric Pacific and Atlantic lineages of these species ranges from 1.24% to 2.98%, and although the taxonomic status in these groups is unknown, these differences are comparable to species level differences at cytochrome c oxidase subunit I sequences (COI ) among other echinoids (Echinometra 2-3%, Palumbi et al., 1997;Leptasterias 0.4-2.2%, Hrincevich, Rocha-Olivares & Foltz, 2000; Patiriella 1.1-4.3%, Hart, Byrne & Johnson, 2003), including the cryptic lineages of S. droebachiensis (2.3%, Addison & Kim, 2018). However, the distribution and extent of the ecological segregation of both lineages of S. droebachiensis throughout the north Atlantic are not well known, particularly in the Labrador Sea, where the west flowing Greenland Current is expected to connect populations across the north Atlantic (Knutsen et al., 2007;Bringloe, Verbruggen & Saunders, 2020). Furthermore, genetic variation within S. droebacheisis has been largely unsampled in the north Pacific, making it difficult to assess the role of trans-Arctic vicariance to patterns of evolution within the genus.
Here we extend analyses of biogeography and population genetic structure in circumpolar Stronglyocentrotus sea urchins to better understand the roles that trans-Arctic and trans-Atlantic dispersal have played in the systematic and genetic divergence within the genus. We aim to establish a more complete understanding of the range limits of each species by compiling previous surveys of mtDNA sequence variation and including additional sample sites for S. pallidus and S. droebachiensis throughout the Pacific and north Atlantic. Strongylocentrotus pallidus is a circumpolar species that is abundant in shallow water (<15 m) in the north, and deeper waters of up to 1600 m at lower latitudes (Jensen, 1974;Strathmann, 1981;Gagnon & Gilkinson, 1994;Bluhm, Piepenburg & Von Juterzenka, 1998). The current known distribution of S. droebachiensis W includes the northeast Pacific and the shallow water habitat of the northwest Atlantic, and S. droebachiensis E is the only green sea urchin found in the northeast Atlantic and in the deep offshore habitat in the northwest Atlantic. However, there is less certainty about the full distribution of S. droebachiensis E, because although it appears to be circumpolar, it was only detected at low frequency in the Pacific (5/29 samples; Addison & Hart, 2005), and none of these haplotypes were shared with the Atlantic suggesting a lack of recent trans-Arctic dispersal. Although mtDNA from this lineage was extremely rare at lower latitudes in the shallow northwest Atlantic habitat (6/161 samples; Addison & Hart, 2005;Addison & Kim, 2018), a more complete understanding of the biogeography within the genus requires additional sampling at higher latitudes in both the northwest Atlantic and northeast Pacific. Late glacial and post-glacial trans-Arctic dispersal between populations of S. pallidus and S. droebachiensis W has resulted in both shallow inter-ocean divergence and widespread sharing of haplotypes (Palumbi & Kessing, 1991;Palumbi & Wilson, 1990;Addison & Hart, 2005;Laakkonen et al., 2021). Since the coastal habitat at high latitudes in the Pacific and northwest Atlantic is qualitatively similar to the northeast Atlantic where S. droebachiensis E dominates (e.g., cool water, kelp;Payne et al., 2012;Government of Canada, 2014;Bringloe & Saunders, 2019), we expect concordant biogeography and patterns of genetic diversity in S. droebachiensis E if it shared a similar history of invasion, vicariance, and secondary contact. However, an absence of haplotype sharing or discordant biogeography would suggest a lack of trans-Arctic dispersal (persistent vicariance), raising the possibility that S. droebachiensis E diverged in allopatry following the initial trans-Arctic invasion in the late Pliocene. This study will provide insight into the mechanisms driving reproductive isolation in the northern lineages of Strongylocentrotus sea urchins by defining the biogeographical distribution of allopatric and sympatric populations, and quantifying genetic subdivision both within and between oceans.

Sampling, DNA extraction, amplification, and sequencing
To examine the range distributions of all three lineages of Strongylocentrotus sea urchins we compiled COI sequence data from previous studies (GenBank accession Nos. AY504479-AY504511, Addison & Hart, 2005;EF108346-EF108365, Harper, Addison & Hart, 2007;MG098337-MG098440, Addison & Kim, 2018;MT736172-MT736220, Laakkonen et al., 2021), and collected new samples from 10 locations throughout the north Atlantic and northeast Pacific oceans ( Fig. 1; Table 1). New samples were collected from most sites at depths of 2-8 m using SCUBA, except in the Bay of Fundy (14-90 m) where samples were collected using a fixed gear dredge, and Owl's Head Nova Scotia (60 m) where collections were made using baited lobster traps as described in Filbee-Dexter & Scheibling (2014). Sea urchins were collected with permission under Section 52 permits (Department of Fisheries and Oceans Canada: 356132, NL-2619-14, TNMP-2014-16478, and S-14/15-1053, and the State of Alaska Department of Fish and Game (CF-17-004). Sea urchins from the northeast Atlantic at Kongsfjord and OsloFijord (Norway) were a subset of those analysed by Norderhaug et al. (2016) for which we generated new mtDNA sequence data. We generally observed the colour characteristics of the specimens (i.e., test, tube feet, and aboral spines) described in Jensen (1974) and Gagnon & Gilkinson (1994), but since the structures used to distinguish the species are challenging to observe under field conditions (e.g., pedicellariae, spine wedges, and pore pairs), designations were ultimately made using DNA sequence data (see Addison & Kim, 2018). We preserved gonad tissue and/or tube   Table 1  feet in 95% ethanol, and extracted total genomic DNA using DNAeasy Blood and Tissue columns (QIAGEN) following the manufacturer's recommended protocols.

Polymorphism and population genetic structure
We identified individuals as Strongylocentrotus pallidus, or one of the two reproductively isolated cryptic lineages of S. droebachiensis (see Addison & Kim, 2018) using a combination of maximum likelihood and statistical parsimony. We inferred a phylogenetic tree of unique haplotypes by maximum likelihood using PHYML 3.0 (Guindon et al., 2010), with an HKY85 substitution model, gamma distributed rate heterogeneity at sites, and an SPR tree search. Node support for the putative species clusters was estimated using nonparametric bootstrap analysis with 1,000 replicates. To visualize species assignment, we used statistical parsimony implemented in TCS v.  Rozas, 2002). Significance was assessed by 10,000 coalescent simulations. To control for the occurrence of false positives due to multiple comparisons, significance of the p-values was determined using the Bonferroni correction. To simultaneously visualize both the phylogenetic relationships and the frequency of each haplotype, we constructed separate haplotype networks for each species using statistical parsimony implemented in TCS v.1.21 (Clement, Posada & Crandall, 2000) and presented using PopART.
To evaluate the genetic subdivision among populations of each lineage within and between major oceanographic regions, we calculated global F ST and tested for pairwise genetic differences between populations. We conducted analyses of molecular variation (AMOVA) to test for hierarchical genetic structure both within and among the Pacific and Atlantic Oceans (Fig. 1). We also explored post hoc hypotheses based on patterns of pairwise F ST to further refine patterns of substructure. Indices of genetic differentiation (F ST and ) were calculated using Kimura two-parameter distances (K2P: Kimura, 1980) implemented in ARLEQUIN, and significance was assessed using 10,000 permutations of the data with Bonferroni correction for multiple tests.

RESULTS
We obtained 418 bp COI sequences (positions 6415-6832 of Jacobs et al., 1988) for 789 individual sea urchins. There were 60 variable sites and a total of 83 unique haplotypes. Based on maximum likelihood and statistical parsimony analyses (Fig. 1C) our results support the presence of three lineages reported in Addison & Kim (2018). We detected three distinct clusters of haplotypes including S. pallidus (n = 156), and both lineages of S. droebachiensis (S. droebachiensis E, n = 148; S. droebachiensis W, n = 485) ( Table 2). Mean genetic distance among all the sequences (K2P) was 2.40%, while within lineage mean genetic distance ranged from 0.33% to 0.41% (Fig. 1). Pairwise genetic distance was the greatest between S. pallidus and S. drobachiensis W, while the cryptic lineages within S. droebachiensis were 2.73% divergent. For all 3 lineages, the net between group distances (K2P) were 6-28x greater when comparing samples between oceans (i.e., trans-Arctic) than between the east and west coasts of the same ocean (trans-Pacific or trans-Atlantic; Table 3). Overall, haplotype (h) and nucleotide (π) diversity was high in all three lineages, and values ranged from 0.685-0.733 and 0.033-0.0040, respectively (Table 2). Significant negative values of Fu's F S (−6.600 to −4.400) suggest a demographic expansion (or purifying selection) in samples of S. pallidus and S. droebachiensis E from offshore sites (>60 m) in the northwest Atlantic, and S. droebachiensis W from OKH and NOM in northeast Pacific (Table 2). Sea urchin populations of all three lineages did not show an excess of alleles in shallow water habitats throughout the northwest Atlantic.

Biogeography
Strongylocentrotus pallidus was distributed in high relative abundance across all three oceanographic regions sampled. Pure populations of S. pallidus were detected above the Arctic Circle in both the northwest and northeast Atlantic (Fig. 1), and this species was relatively abundant (13-25%) in mixed aggregations at shallow sites along the coast of Labrador and western Newfoundland. However, S. pallidus was rare at all other shallow water sites in Atlantic Canada (six of 236 total samples; 2.5%), including one individual at 14 m in the Bay of Fundy where it was absent at depths >70 m. In contrast, S. pallidus was common at deeper sites, making up 28% of the samples collected at 60 m off the coast at Owl's Head, NS (19 of 69), and 69.6% of the samples collected Offshore at depths >65 m on the Scotian Shelf (39 of 56).
The green sea urchin, S. droebachiensis, was detected across all three oceanic regions sampled, but there were striking differences in the distribution of the cryptic lineages. Strongylocentrotus droebachiensis E was the only green sea urchin found in the northeast Atlantic, where pure populations were sampled throughout Iceland, Norway, and Russia ( Fig. 1). In the northwest Atlantic, S. droebachiensis E shared a distribution similar to S. pallidus where it made up ∼12% of the samples at shallow sites at higher latitudes throughout the Labrador Sea and the Gulf of St Lawrence (20/162). However, this lineage was rare in the shallow coastal samples (<60 m) at lower latitudes, comprising only ∼1% of the individuals sampled along Nova Scotia and the Bay of Fundy (two of 201), and was the only green sea urchin found offshore on the Scotian Shelf. With the exception of the     (2005), we failed to detect additional samples of S. droebachiensis E throughout the north Pacific. Strongylocentrotus droebachiensis W was distributed throughout the Pacific and coastal samples from Atlantic Canada. With the exception of the low frequency and geographically isolated haplotypes of S. droebachiensis E detected at SJI, all other green sea urchin samples collected throughout the Pacific Ocean were identified as S. droebachiensis W (Fig. 1). In the northwest Atlantic, S. droebachiensis W was detected at all coastal sites (including OWH and BOF) where it comprised an increased proportion of the samples at lower latitudes (Fig. 1)

Strongylocentrotus pallidus
Range-wide genetic structure in S. pallidus was primarily driven by differences between the Pacific and Atlantic samples, and an absence of genetic subdivision within the north Atlantic. Statistical parsimony identified a single abundant and geographically widespread genetic variant (50% of all samples) distributed across all three major oceanographic regions. In addition to this shared haplotype, there was some inter-ocean divergence as half of the haplotypes detected in both the Pacific and the shallow water northwest Atlantic were exclusive to those regions (Fig. 2). In contrast, with only two exceptions, all the S. pallidus samples from Norway (PSF) shared a single COI haplotype. Multi-locus nuclear genotypes have been scored for nine of the 20 individuals (all from ACS) harbouring the three unique northwest Atlantic S. pallidus haplotypes (Fig. 2), all of which have been identified as S. droebachiensis W (Addison & Hart, 2005;Addison & Kim, 2018;Burke, Kim & Addison, in prep.) suggesting historic hybridization and introgression (Addison & Pogson, 2009). Since these haplotypes have been segregating in S. droebachiensis W for many generations, they do not reflect the biology of S. pallidus and thus we removed them from subsequent analyses of population genetic structure. Global F ST was high (0.2816; P < 0.001), indicating strong and significant variance in the distribution of genetic variation. Pairwise F ST values were high between the Pacific (SJI) and the northwest Atlantic samples from (OWH and ACO) and northeast Atlantic (PSF; Pairwise F ST = 0.2998-0.3704, P < 0.003; Table 4), suggesting limited dispersal with gene flow across the Arctic basin. In contrast, there was not a significant difference between sample sites within the north Atlantic (Pairwise F ST = 0-0.0689, P > 0.003). Hierarchical analysis of molecular variance (AMOVA) based on the a priori grouping of populations from each oceanographic region were not significant ( CT = 0.353; P = 0.205), and we failed to detect evidence of sub-structure based on geography or depth within the Atlantic (Table 5).

Strongylocentrotus droebachiensis E
Patterns of genetic subdivision were largely driven by differences between the Pacific and Atlantic basins, and low substructure within Atlantic. We detected seven COI haplotypes S. droebachiensis E sampled from Norway and Iceland, four of which were also widespread and abundant in the northwest Atlantic (Fig. 3). Two of the three unique haplotypes in the northeast Atlantic were detected in the easternmost Arctic sites in Norway (VAR) and Russia (KOP). While S. droebachiensis E was detected throughout the northwest Atlantic, the majority of the individuals were found at deep sites (>60 m) and had rare mtDNA haplotypes. Of the 28 unique haplotypes detected for this lineage, 21 were only detected once, and of those 14 were found in the northwest Atlantic. None of the five COI haplotypes genotyped at SJI in the Pacific were found in the Atlantic populations. Global F ST was lower for this lineage compared to the others (F ST = 0.2052; P < 0.001), indicating a moderate level of genetic subdivision throughout the range. Pairwise comparisons revealed strong divergence between both SJI (Pacific), and ND2 (Norway) from most other sites (Table 6). Genetic subdivision was generally low and not significantly different from zero among most locations across the Atlantic Basin. Hierarchical AMOVA indicated strong regional grouping based on the oceanographic regions within the northeast Atlantic ( CT = 0.3287, P < 0.0001; Table 5, AMOVA), but we failed to detect significant variation across the north Atlantic or among the a priori grouping of sampling sites across oceanographic basins.

Strongylocentrotus droebachiensis W
Genetic variation within S. droebachiensis W was consistent with previous studies, with strong differences across the Arctic and genetic homogeneity among southern coastal sites in the northwest Atlantic. However, shared halpotypes and generally lower pairwise F ST values (Table 7) among new samples from the Labrador Sea and north Pacific suggest a greater influence of trans-Arctic dispersal. Unique haplotypes were found in both the Pacific and northwest Atlantic, but three high frequency genetic variants were shared throughout both oceans (Fig. 4). Based on both the distribution of haplotypes and pairwise F ST among new * 3/3 * 1/1 * 5/16 Figure 2 Sampling locations, haplotype distribution, and TCS haplotype network of COI mtDNAsequences for Strongylocentrotus pallidus (n = 156). Asterisks (*) indicate the mtDNA haplotypes removed from analyses of population genetic structure because they were recovered in individuals whose nuclear genomes (SNPs or microsatellites) were characterized as being 100% S. droebachiensis W (# tested/# individuals with the haplotype).   Addison & Hart, 2005) sites throughout the Gulf of St. Lawrence, Nova Scotia, and the Bay of Fundy were genetically homogeneous (pairwise F ST = 0 to 0.0400, P > 0.10). There was a striking correlation between latitude and genetic similarity among the northern sites in the Pacific and the Labrador Sea. Pairwise F ST was not significant between NOM (Alaska; 64.487 • N) and DUR (Nunavut; 67.038 • N), but these sites were different from the next closest sample in each region ( groups of sea urchins in the east and west Pacific, Labrador Sea, and coastal northwest Atlantic ( CT = 0.49691; P = 0.0013).

DISCUSSION
The biogeographic distribution, concordant population genetic structure, and patterns of haplotype sharing among oceanic regions suggest that cycles of vicariance and trans-Arctic gene flow has shaped diversification within circumpolar Strongylocentrotid sea urchins. While there is considerable debate about the competing contributions of both geographic isolation and divergence with gene flow to the process of speciation in the sea (e.g., Miglietta, Faucci & Santini, 2011;Faria, Johannesson & Stankowski, 2021), our results suggest that isolation across the Arctic Basin has been a driving force of genomic and systematic diversity within the genus. Consistent with earlier studies (Palumbi & Wilson, 1990;Palumbi & Kessing, 1991;Addison & Hart, 2004;Addison & Hart, 2005;Harper, Addison & Hart, 2007), we detected widespread sharing of identical haplotypes throughout the Pacific and Atlantic populations of S. pallidus and S. droebachiensis W, and patterns of population genetic subdivision among these regions suggests recent interoceanic exchange. However, there was no evidence of a similar pattern of trans-Arctic dispersal in S. droebachiensis E, as we failed to detect haplotypes from this species at additional sample sites throughout the north Pacific. Following the initial trans-Arctic invasion of the north Atlantic by Pacific ancestors, our results suggest that Strongylocentrotus droebachiensis diverged into reproductively isolated cryptic species, one of which remains connected with the pacific (S. droebachiensis W) while the other is now endemic to the Atlantic (S. droebachiensis E). Although the strong patterns of hierarchical population structure within S. droebachiensis W suggests a contribution of latitude to the distribution of genetic variation, our analysis of a putatively neutral mtDNA locus does not display a similar signature of adaptive evolution in response to temperature driven by latitudinal variation reported for populations of Atlantic cod  ( Bradbury et al., 2010) and Atlantic salmon (Jeffery et al., 2017). Our results suggest that repeated trans-Arctic gene exchange contributed to the maintenance of species integrity in two species, while vicariance may have contributed to the allopatric speciation of a third species in the north Atlantic. Like many arctic-boreal marine species, Strongylocentrotus sea urchins with a circumarctic distribution have experienced cycles of isolation and invasion throughout the Pleistocene ice ages (e.g., Laakkonen et al., 2021). Following the trans-Arctic invasion in the late Pliocene, allopatric populations of Pacific and Atlantic sea urchins diverged throughout the Pleistocene when cycles of glacial advance and retreat (ca 2.4-3.0 Mya) restricted dispersal around the Bering Strait (Einarsson, Hopkins & Doell, 1967;Herman & Hopkins, 1980;Maslin et al., 1996;Haug et al., 1999;Harris, 2005;Lisiecki & Raymo, 2005;Horikawa et al., 2015;Loeza-Quintana et al., 2019). The presence of multiple private haplotypes in Nova Scotia and New Brunswick suggests that S. droebachiensis W persisted in refugia at the southern end of their northwest Atlantic range, but moderate genetic structure and sequence similarity among oceans indicates that Pacific populations subsequently reinvaded the Atlantic during interglacial periods throughout the Pleistocene (Hewitt, 2003;Maggs et al., 2008). We detected patterns of genetic subdivision and sequence diversity within S. pallidus that are consistent with S. droebachiensis W, although the small sample size of the hierarchical analysis (three groups, six sites) suffered from low power and could not approach significance at the 5% threshold (Fitzpatrick, 2009). Genetic analysis of new samples collected at high latitudes in both oceans revealed extensive haplotype sharing and lower pairwise F ST values compared to studies conducted at more southern latitudes (Palumbi & Wilson, 1990;Palumbi & Kessing, 1991;Addison & Hart, 2004;Addison & Hart, 2005;Harper, Addison & Hart, 2007). These results indicate that trans-Arctic connectivity is likely greater than previously reported, and while qualitatively consistent with coalescent analyses (Laakkonen et al., 2021) that show the predominant migration vectors track the east flowing currents connecting the Pacific with the northwest Atlantic (Ledu et al., 2008), we cannot exclude a hypothesis of back migration from the northwest Atlantic to the Pacific (see Addison & Hart, 2005;Harper, Addison & Hart, 2007). While Pacific and northwest Atlantic populations of S. pallidus and S. droebachiensis W experienced periods of vicariance throughout the Pleistocene, late glacial and post-glacial trans-Arctic dispersal continues to maintain the integrity of these species.
Patterns of biogeography and genetic diversity suggest that S. droebachiensis E is almost exclusively limited to the Arctic and sub-Arctic in north Atlantic. Additional sampling confirmed that this lineage is the only green sea urchin found in the northeast Atlantic, and that its range in the northwest Atlantic is characterized by a clear shift to deeper habitats at lower latitudes. While S. droebachiensis E was present at some shallow sites in the Canadian Arctic, Labrador Sea, and Gulf of St. Lawrence, our sampling efforts at similar latitudes in the Pacific failed to detect additional evidence of this species beyond those reported by Addison & Hart (2005). Since the sites at Haida Gwaii (LAN, MAS), Sea of Okhotsk (OKH), and the Bering Sea (NOM) share similar coastal temperatures and macroalgal assemblages (Payne et al., 2012;Government of Canada, 2014;Bringloe & Saunders, 2019) as those supporting S. droebachiensis E in the north Atlantic, failure to detect additional representatives suggests they are not broadly distributed throughout the Pacific. In addition, the five singleton haplotypes of S. droebachiensis E reported from the San Juan Islands (Addison & Hart, 2005) were not identified in the north Atlantic, indicating very limited (or complete absence) of trans-Arctic gene flow in this species. These results suggest that the presence of S. droebachiensis E haplotypes in the Pacific could represent incomplete lineage sorting of ancestral alleles, or possibly low levels of back migration of S. droebachiensis E individuals or haplotypes (via introgression into S. pallidus or S. droebachiensis W; Addison & Hart, 2005;Harper, Addison & Hart, 2007) during interglacial periods throughout the Pleistocene. While additional analyses of both coastal and deep habitats throughout the Pacific are required before concluding that S. droebachiensis E is absent from the Pacific Ocean, our findings suggest this species may have evolved in allopatry following the initial trans-Arctic invasion of S, droebachiensis during the late Pliocene. Alternatively, the two lineages of S. droebachiensis could have initially diverged in the Pacific prior to invading the Atlantic, followed by a subsequent reduction (or possibly extirpation) S. droebachiensis E in the Pacific. At the very least, our study reveals that, unlike S. pallidus and S. droebachiensis W, Pacific and Atlantic populations of S. droebachiensis E continue to diverge in a state of persistent trans-Arctic vicariance.
The repeated trans-Arctic dispersal of S. pallidus and S. droebachiensis W following the initial period of vicariance suggests that the northwest Atlantic is a zone of secondary contact between all three species. Early studies employing microsatellites (Addison & Hart, 2005), nuclear DNA sequences (Addison & Pogson, 2009), and single nucleotide polymorphisms (SNPs;Addison & Kim, 2018) detected mitochondrial and nuclear discordance in 9 of the 305 (3.0%) individual sea urchins analyzed throughout the Pacific and Atlantic oceans. In these studies, all the hybrid individuals identified were a result of introgression of S. pallidus mtDNA into S. droebachiensis individuals. For example, Addison & Kim (2018) tested for evidence of hybridization using both COI sequences and 3,049 nuclear SNPs in a sample of 110 sea urchins collected along a depth gradient off the coast of Nova Scotia. While two S. droebachiensis W individuals from shallow sample sites harboured S. pallidus mtDNA, the lack of admixture across the nuclear genome of all samples provides evidence against widespread contemporary hybridization and suggests that reproductive isolation is complete. In Addition, patterns of endemism of the introgressed haplotypes in both oceans suggests that historic introgressive hybridization from S. pallidus into S. droebachiensis W may have occurred independently in Pacific and Atlantic populations. Previous studies have not revealed evidence of hybridization between S. droebachiensis E and the other two species. However, extensive analyses of both nuclear and mitochondrial DNA throughout the northwest Atlantic are needed to test the hypothesis that contemporary hybrids form under natural spawning conditions, particularly at sites where all 3 species co-occur (e.g., OWH and NAI).
While trans-Arctic vicariance is the dominant mechanism driving the initial divergence of S. droebachiensis E from ancestors in the Pacific, allopatry within the Atlantic has contributed to patterns of divergence in other echinoderms. Beginning in the mid-Pliocene, rapid ocean cooling and the formation of the Labrador current isolated temperate north Atlantic species where warmer mid-Atlantic and Gulf stream waters provided refuge on north American and European coasts (Berggren & Hollister, 1974;Franz, Worley & Merrill, 1981;Cronin, 1988;Wares, 2001). Genetic evidence supports this hypothesis in sea stars, where western Asterias forbesi and eastern A. rubens diverged in allopatry followed by the post-glacial recolonization and sympatry in the northwest Atlantic (Wares, 2001;Wares & Cunningham, 2001). These species now form a secondary contact zone from Nova Scotia to Cape Cod, and laboratory studies of sperm competition (Harper & Hart, 2005), morphology, and genetic surveys of natural populations (Harper & Hart, 2007) have identified hybridization and introgression. While patterns of ecological, morphological, and genetic divergence identified within S. droebachiensis are qualitatively similar to those for Asterias, our results only weakly fit the scenario of post Pliocene divergence of allopatric populations within the north Atlantic. Support for this hypothesis includes evidence of reproductive isolation between the east and west lineages (Addison & Kim, 2018), habitat segregation of the eastern lineage in the west, and a signal of range expansion in western samples of S. droebachiensis E and the co-distributed population of S. pallidus. However, we failed to detect moderate or weak population genetic structure typical of recent trans-Atlantic dispersal (Young et al., 2002;Provan, Wattier & Maggs, 2005;Chevolot et al., 2006;Jolly et al., 2006;Hoarau et al., 2007;Souche et al., 2015;Andrews et al., 2019;Neiva et al., 2020), and in contrast, we identified more private haplotypes and higher genetic diversity (h, π) in northwest Atlantic samples of both S. pallidus and S. droebachiensis E. These patterns suggest that sea urchins in the northwest Atlantic have persisted in single or multiple glacial refugia (Hewitt, 2003;Maggs et al., 2008), and were unlikely to have been extirpated during glacial maxima throughout the Pleistocene. Although repeated cycles of isolation and dispersal between the east and west coasts throughout the Pleistocene may have obscured signals of historic vicariance (Jesus et al., 2006;Maggs et al., 2008), our results suggest that lineages of S. droebachiensis have not been strictly allopatric within the north Atlantic following the initial invasion, and that vicariance within the Atlantic was not the principal driver of speciation within the genus.
Identifying the mechanisms driving speciation in the sea can be challenging because of the difficulty in identifying barriers to gene exchange, or the environmental factors driving adaptation. Addison & Kim (2018) suggest that tolerance of seasonally lower salinity may contribute to the ecological segregation of the Strongylocentrotus lineages in the southern part of their western Atlantic range (e.g., along the coast of Nova Scotia). In this study, we identified contrasting patterns of geographic distribution and habitat segregation that suggest increased water temperatures in the northwest may contribute to the near complete absence of S. droebachiensis E from shallow sites dominated by S. droebachiensis W. In the northwest Atlantic, larvae of S. droebachiensis grow rapidly at 14 • C (Hart & Scheibling, 1988), and in the Pacific and northwest Atlantic both larvae and adults can withstand temperatures up to 19.7 and 22 • C, respectively (Scheibling & Stephenson, 1984;Pearce et al., 2005). Like other species with planktonic dispersing larvae, S. droebachiensis exhibits large regional and interannual fluctuations in recruitment (e.g., Raymond & Scheibling, 1987;Scheibling & Raymond, 1990;Scheibllng, 1996), but is known to settle along the coast of Nova Scotia in July when water temperature can exceed 14 • C (Balch & Scheibling, 2000). In August and September, the nearshore water temperatures along the coast of Nova Scotia regularly reach 20 • C (Scheibling, Feehan & Lauzon-Guay, 2013). In contrast, water temperatures along the Norwegian coast are comparatively cooler (Danielssen, Svendsen & Ostrowski, 1996;Ibrahim et al., 2014), and green sea urchins experience recruitment failure in kelp beds at southern latitudes when temperatures exceed 10 • C (Fagerli, Norderhaug & Christie, 2013;Rinde et al., 2014;Nyhagen, Christie & Norderhaug, 2018). By limiting sea urchin recruitment, ocean warming is thought to be a driver of ecological change in Norway, as the southern boundary (65 • 70 N;Fagerli, Norderhaug & Christie, 2013) between kelp-dominated habitat and overgrazed urchin barren grounds continues to shift northward with corresponding increases in water temperature (Rinde et al., 2014).
Differences in thermal tolerance among lineages of S. droebachiensis may explain the habitat segregation we observed in the northwest Atlantic. The extreme rarity of S. droebachiensis E in the shallow habitat along the coast of Nova Scotia could be driven by seasonally warmer water temperatures resulting in recruitment failure, post-settlement mortality, or mortality of juveniles or adults. While summer ocean temperatures along the coast of Nova Scotia are impacted by the Gulf Stream and storm activity (Scheibling, Feehan & Lauzon-Guay, 2013), lower water temperatures in the Gulf of St Lawrence and coastal Newfoundland and Labrador are moderated by the cool south flowing Labrador Current. The increased abundance of S. droebachiensis E (and S. pallidus) at depths <15 m throughout this part of the range (i.e., north of Nova Scotia) may be explained by seasonal temperatures at or below the 10 • C threshold observed in the northeast Atlantic. The influence of temperature on the distribution of S. droebachiensis E is supported by both the decrease in recruitment success along the coast of Norway (Fagerli, Norderhaug & Christie, 2013) and the shifting population dynamics of green sea urchins in Oslofjord along the southern coast of Norway. In a response to increased sea surface temperatures (SST), Nyhagen, Christie & Norderhaug (2018) demonstrated a significant shift in population density from 10-15 m to cooler water at 20 m, and a reduction in both sea urchin size and recruitment success in 1979 and 1992 compared to 2013. Additionally, while sea urchins are present throughout southern Norway, abundant populations typically only occur at depths of 20 m or greater (e.g., site ND2 in this study; Norderhaug et al. (2016)). Although changes in coastal SST indicate a rapid warming trend in both the northwest Atlantic (∼1.0 • C per decade) and the margins of Norwegian and North Seas (between ∼0.3 and 0.7 • C per decade), particularly during the planktonic dispersal and settlement of sea urchins from late spring to autumn (Lima & Wethey, 2012), the samples analysed in our study were collected over a relatively short time scale (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) and are unlikely to have captured ongoing changes in sea urchin distribution in response to increasing SST (e.g., Hobday & Pecl, 2014). Though we suggest that temperature may be an important factor in defining the range of S. droebachiensis E, comparative analyses of the thermal tolerance of larvae, juveniles and adults of both species are required to test this hypothesis.
The evolution of gamete recognition molecules has long been viewed as an important driver of speciation in marine invertebrates (Vacquier, 1998;Palumbi, 2009;Lessios, 2011;Vacquier & Swanson, 2011). Interspecific sperm competition in the plankton is mediated by a variety of proteins and carbohydrates (sulfated polysaccharides) coating the sperm and eggs (Biermann et al., 2004), and positive selection detected at sperm Bindin (e.g., Biermann, 1998) correlates with the strength of reproductive isolation between species (Zigler et al., 2005). Palumbi & Lessios (2005) showed that, in addition to a steady accumulation of genome divergence over time, the rate of speciation in sea urchins also depends on the rate of evolution of gamete recognition proteins. In their study, Palumbi & Lessios (2005) surveyed species in eight genera and showed that the presence of sympatric species was common in genera with rapid evolution of sperm Bindin. Since studies of Bindin evolution within Strongylocentrotus only included samples of S. droebachiensis from the Pacific (Biermann, 1998;Pujolar & Pogson, 2011), it is difficult to assess patterns of positive selection and sequence divergence between the cryptic species of S. droebachiensis. However, Marks et al. (2008) detected 1.5% sequence divergence at sperm Bindin between samples S. droebachiensis from Norway, the northwest Atlantic, and northeast Pacific, and based on the conclusions of our study we suggest that this difference represents interspecific divergence. In a series of heterospecific and conspecific crosses between S. pallidus and S. droebachiensis from the Pacific and S. droebachiensis from Norway, Biermann & Marks (2000) demonstrated strong asymmetry in fertilization compatibility among allopatric populations. Our data suggests that the allopatric populations studied by Biermann & Marks (2000) represent distinct species, where eggs of S. droebachiensis E cannot be fertilized by sperm from S. pallidus, and we interpret their results as a test for reproductive isolation among the species. Consistent with Strathmann (1981), eggs of S. droebachiensis W were receptive to sperm from both S. pallidus and S. droebachiensis E, but sperm from S. droebachiensis W either failed (S. pallidus) or had very low (S. droebachiensis E) fertilization rates in heterospecific crosses. Similarly, eggs of S. droebachiensis E could not be fertilized by S. pallidus sperm, and their sperm successfully fertilized eggs of S. droebachiensis W but were not tested with eggs of S. pallidus. In addition to ecological segregation, patterns of sperm Bindin evolution and asymmetry in reproductive compatibility may contribute to the maintenance of species boundaries in sympatry in the northwest Atlantic. Detailed studies of gamete recognition molecules and sperm competition among all three species will help to further understand mechanisms driving evolution within the genera.

CONCLUSION
While previous work identified a reproductively isolated cryptic lineage of S. droebachiensis in the northwest Atlantic (Addison & Kim, 2018), this study supports the hypothesis that these species formed as a result of vicariant speciation driven by trans-Arctic isolation. Our results show widespread sharing of S. pallidus and S. droebachiensis W haplotype variants throughout the north Pacific and north Atlantic Oceans, but that S. droebachiensis E is largely restricted to the north Atlantic. We detected low genetic subdivision between S. droebachiensis W from the north Pacific and the Labrador Sea, suggesting widespread trans-Arctic gene flow in this species. There was weaker evidence of trans-Arctic dispersal in S. pallidus, which could possibly be an artefact of poor sampling of this species in the Pacific. Our analyses of biogeography and COI sequence diversity suggests that following allopatric speciation during the Pliocene or early Pleistocene, these species established a zone of secondary contact in the northwest Atlantic and the Labrador Sea. In the northwest Atlantic, we identified sites along the coast of Labrador (NAI) and Nova Scotia (OWH) where all three species of Strongylocentrotus are abundant, providing a natural laboratory for studying the ecological and molecular aspects driving the evolution of barriers to gene exchange. From a biogeographic perspective, understanding the mechanisms shaping the distribution of the S. droebachiensis species throughout the north Atlantic requires experiments to determine the physiological limits of both. We observed patterns of ecological segregation among the species that suggest temperature may play a role in habitat selection, particularly in the warmer water along the coast of Nova Scotia. In addition, while Addison & Kim (2018) provided evidence of reproductive isolation among species collected from Nova Scotia, a wider study aimed at detecting hybridization and introgression at nuclear loci is required to characterize the extent of reproductive isolation across a broader range of habitats. Viewed with a species-specific lens, both rapid sequence divergence at sperm Bindin (Marks et al., 2008) and the accumulation of interspecific gamete incompatibility (Biermann & Marks, 2000) between S. droebachiensis W and S. droebachiensis E suggests a potential role of reinforcement selection for pre-zygotic isolation (Coyne & Orr, 2004) following secondary contact. By characterizing the extent of reproductive isolation, both laboratory studies of sperm competition and interspecific